Improving an approximation for the inverse of the Riemann–Siegel θ-function Recall the Riemann–Siegel θ-function:
$$\theta(z) = \arg\Gamma\left(\frac{1}{4}+\frac{i\,z}{2}\right) - \frac{z\,\log \pi}{2},$$
that describes the complex phase of the Riemann $\zeta$-function on the critical line.
There is a known approximation for its inverse:
$$\theta^{\small(-1)}(x)=\frac{\pi+8{\tiny\text{ }}x}{4\,W\!\left(\frac{\pi+8{\tiny\text{ }}x}{8{\tiny\text{ }}\pi{\tiny\text{ }}e}\right)}+o(1),$$
where $W(x)$ is the Lambert W-function, which becomes more precise as $x$ grows.
I wonder if it is possible to improve this approximation by including higher-order terms, so that the remaining error term decays as $o(x^{-1})$, $o(x^{-2})$, etc. Can those higher-order terms be expressed using only elementary functions and $W(x)$?
 A: We start with the asymptotics
$$
\theta (t) = \frac{t}{2}\log \frac{t}{{2\pi }} - \frac{t}{2} - \frac{\pi }{8} + \frac{1}{{48t}} + \mathcal{O}\!\left( {\frac{1}{{t^3 }}} \right),
$$
i.e.,
$$
\frac{{\theta (t)}}{\pi } + \frac{1}{8} = \frac{t}{{2\pi }}\log \frac{t}{{2\pi }} - \frac{t}{{2\pi }} + \frac{1}{{48\pi t}} + \mathcal{O}\!\left( {\frac{1}{{t^3 }}} \right).
$$
This may be re-written in the form
$$
\frac{{\theta (t)}}{\pi } + \frac{1}{8} = \left( {\frac{t}{{2\pi }} + g(t)} \right)\log \left( {\frac{t}{{2\pi }} + g(t)} \right) - \left( {\frac{t}{{2\pi }} + g(t)} \right),
$$
where
$$
g(t) = \frac{1}{{48\pi t\log \frac{t}{{2\pi }}}} + \mathcal{O}\!\left( {\frac{1}{{t^3 \log t}}} \right).
$$
Thus,
$$
\frac{1}{e}\left( {\frac{{\theta (t)}}{\pi } + \frac{1}{8}} \right) = \frac{{\frac{t}{{2\pi }} + g(t)}}{e}\log \frac{{\frac{t}{{2\pi }} + g(t)}}{e},
$$
i.e.,
$$
\frac{{\frac{{\theta (t)}}{\pi } + \frac{1}{8}}}{{W\!\left( {\frac{1}{e}\left( {\frac{{\theta (t)}}{\pi } + \frac{1}{8}} \right)} \right)}} = \frac{t}{{2\pi }} +g(t)= \frac{t}{{2\pi }} + \frac{1}{{48\pi t\log \frac{t}{{2\pi }}}} + \mathcal{O}\!\left( {\frac{1}{{t^3 \log t}}} \right).
$$
Iterating this once yields
$$
\frac{{\frac{{\theta (t)}}{\pi } + \frac{1}{8}}}{{W\!\left( {\frac{1}{e}\left( {\frac{{\theta (t)}}{\pi } + \frac{1}{8}} \right)} \right)}} = \frac{t}{{2\pi }} + \frac{1}{{96\pi ^2 \left[ {\frac{{\frac{{\theta (t)}}{\pi } + \frac{1}{8}}}{{W\left( {\frac{1}{e}\left( {\frac{{\theta (t)}}{\pi } + \frac{1}{8}} \right)} \right)}}} \right]\log \left[ {\frac{{\frac{{\theta (t)}}{\pi } + \frac{1}{8}}}{{W \left( {\frac{1}{e}\left( {\frac{{\theta (t)}}{\pi } + \frac{1}{8}} \right)} \right)}}} \right]}} \\ + \mathcal{O}\!\left( {\frac{{\log ^2 \theta (t)}}{{\theta ^3 (t)}}} \right).
$$
By solving for $t$, simplifying and introducing the inverse function, we find
$$
\theta ^{ - 1} (t) = \frac{{8t + \pi }}{{4W\!\left( {\frac{{8t + \pi }}{{8\pi e}}} \right)}} - \frac{{W\!\left( {\frac{{8t + \pi }}{{8\pi e}}} \right)}}{{6 (8t + \pi )\left( {\log \left( {\frac{{8t + \pi }}{{8\pi }}} \right) - \log W\!\left( {\frac{{8t + \pi }}{{8\pi e}}} \right)} \right)}} + \mathcal{O}\!\left( {\frac{{\log ^2 t}}{{t^3 }}} \right).
$$
For $t=100$ this, without the error term, gives $108.5639773824\ldots$ whereas the exact value is $108.5639773815\ldots$. It is possible to obtain higher terms by using more terms from the asymptotics of $\theta(t)$, obtaining more terms for $g(t)$ and so on. But this leads to elaborate computations once one starts iterating.
A: (this is not an answer but too long for a comment)
(+1) Interesting discussion and answers! Three years earlier I searched the best constant $C$ in following approximate value of the imaginary part of the $n$-th non-trivial zero (from your initial expression of course) : $$\;t_n\approx 2\pi\,\exp(W((n-7/8-C)/e)+1)=2\pi\dfrac{n-7/8-C}{W((n-7/8-C)/e)}$$
and conjectured that $C$ had to be exactly $\dfrac 12$ (computing different moving averages and so on). Further the actual error doesn't exceed $\pm 1$ for the first $2$ million zeros as illustrated  :

Notice the vertical symmetry around $0$ and the slow decrease of the variance of the error with $n$ (a correction term depending of $n$ appears less interesting than in your question, if needed at all, since the mean error remains near $0$ for values as large as $10^{22}$ using Andrew Odlyzko's tables ).
Anyway I found this a neat illustration of the gentle statistical distribution of the zeros.
We seem further able to find the position of the $n$-th zero for $n$ as large as we want with an error of less than one (the error for the $10^4$ zeros following $10^{22}$ is less than $0.21$).
For $\,n=10^{22}+1\,$ for example the formula gives us
$t_n\approx 1370919909931995308226.770224\ $ while the actual zero is at :
$t_n= 1370919909931995308226.680161\cdots$
