Connection between Bernoulli numbers and Riemann-Siegel theta function? I have come across a strange approximation for the Riemann-Siegel theta function involving the Bernoulli numbers - namely that
$$\frac{1}{2} \log \left| B_{2 n}\right|\approx \vartheta (2n)\ ,\quad n \in \mathbb{N}$$
where $B_n$ is the $n$th Bernoulli number, and $\vartheta,$ the Riemann-Siegel theta function. More accurately, it seems that
$$\frac{1}{2} \left(\log \left(\frac{\left| B_{2 n}\right| }{\sqrt{2 n}}\right)-\left(\frac{2}{3}\right)^{2/3} \pi \right)$$
is closer to the actual value of $\vartheta (2n),$ and it may be true to say that 
$$\lim_{n\rightarrow\infty}\frac{1}{2} \left(\log \left(\frac{\left| \sqrt{2}n \zeta(1-2n)\right| }{\sqrt{ n}}\right)-\left(\frac{2}{3}\right)^{2/3} \pi \right)=\vartheta (2n),$$
and it seems likely that large values of $|B_{2n}|$ can be estimated with
$$\sqrt{2 n} \exp \left(2 \vartheta (2 n)+\pi  \left(\frac{2}{3}\right)^{2/3}\right).$$
I haven't seen any references anywhere alluding to any connection between the absolute values of the Bernoulli numbers at even $n,$ and the Riemann-Siegel theta function previously. Is the above statement correct, and if it is, what is the connection?
 A: 
Summary: Your approximation is very good, but not quite perfect. The absolute error does not go to $0$ as $n\to\infty$. But it is astonishingly close: the error is approximately $1.3\times 10^{-8}+o(1)$.

Detailed answer:
we can easily check your proposed approximation for $\vartheta(2n)$, since both $|B_{2n}|$ and $\vartheta(t)$ have well-known asymptotic expansions. In the case of $|B_{2n}|$, we have
$$
B_{2n} = (-1)^{n-1} \zeta(2n) \frac{2(2n)!}{(2\pi)^{2n}} = (1+O(2^{-2n})) (-1)^{n-1} \frac{2(2n)!}{(2\pi)^{2n}},
$$
so, using this together with Stirling's approximation in the form
$$
\log(m!) = m\log m-m+\frac12 \log(2\pi m)+O(m^{-1})
$$
we have that
\begin{align}
\frac12\log \frac{|B_{2n}|}{\sqrt{2n}} &= 
\frac12 \Big(\log 2+
2n\log(2n)-2n+\frac12\log(4\pi n)+O(n^{-1}) - (2n)\log(2\pi) \\ & \qquad - \frac12\log(2n)
+O(2^{-2n}) \Big)
\\ &= n \log\frac{n}{\pi}-n + \frac14 \log(8\pi) + O(n^{-1}).
\end{align}
On the other hand, for $\vartheta(t)$ we have the asymptotic formula
$$
\vartheta(t) = \frac{t}{2}\log\frac{t}{2\pi} - \frac{t}{2} - \frac{\pi}{8} + O(t^{-1}).
$$
which for $t=2n$ gives
$$
\vartheta(2n) = n\log\frac{n}{\pi} - n - \frac{\pi}{8} + O(n^{-1}).
$$
It follows that the difference $\frac12\log \frac{|B_{2n}|}{\sqrt{2n}}-\vartheta(2n)$ satisfies
$$
\lim_{n\to\infty}\left( \frac12\log \frac{|B_{2n}|}{\sqrt{2n}}-\vartheta(2n)\right)
= \frac{\pi}{8} + \frac14 \log(8\pi) \neq \frac12\left(\frac23\right)^{2/3}\pi.
$$
Now note that, amazingly,
$\frac{\pi}{8} + \frac14 \log(8\pi)\approx 1.19874193858$ and $\frac12\left(\frac23\right)^{2/3}\pi \approx 1.19874195162$. The difference between these two constants is $\approx 1.304\times 10^{-8}$. 
To summarize, there is certainly something interesting going on here, and I am very curious to know how you obtained your approximation. Was it just a lucky guess? If not, there seems to be some mysterious mechanism at play causing the constant
$$
\frac12\left(\frac23\right)^{2/3}\pi - \frac{\pi}{8} - \frac14 \log(8\pi)
$$
to assume a surprisingly small numerical value.
Note however that my approach of using the standard asymptotic expansions for $|B_{2n}|$ and $\vartheta(t)$ says nothing about a possible connection between the Bernoulli numbers and the Riemann-Siegel theta function, other than maybe the trivial fact that the asymptotics of the Bernoulli numbers is governed by factorials (aka, the gamma function), and the Riemann-Siegel $\vartheta(t)$ is also defined in terms of the gamma function.
