In the same spirit as DeTemple–Wang for a series expansion of harmonic numbers, I tried to approach the problem as $$H_n\sim\frac 12 \log(n^2+n+a)+\gamma-\frac 1{b(n^2+n+a)+\Delta}\tag 1$$ hoping to minimize $\Delta$ in order to have a good approximation of the inverse of harmonic number $$n^2+n+a=\frac{2}{b\, W\left(\frac{2 }{b}e^{2 (\gamma - H_n)}\right)}\tag 2$$
Cancelling as many coefficient as possible, one obtains very quickly $$b=\frac{6}{3 a-1}\qquad \text{and} \qquad \Delta=-\frac{3 \left(15 a^2-10 a+2\right)}{5 (3 a-1)^2}$$
Then the approximation $$\color{blue}{n(a) \sim -\frac 1 2+\sqrt{\frac 14+\frac{\left(a-\frac{1}{3}\right)}{ W\Big[\left(a-\frac{1}{3}\right) e^{2 (\gamma - H_n)}\Big]}-a}} \tag 3$$
$$\text{rhs - lhs}= \frac{1575 a^4-2100 a^3+1260 a^2-480 a+76}{12600 \,(3 a-1)\, n^6}+O\left(\frac{1}{n^7}\right)\tag 4$$
At this point, we could select $\color{red}{a=\frac 1{10}}$ (remember that in Ramanujan's expansion $a=0$) which corresponds to $\Delta_{\text{min}}=\frac 9{7}$ and gives $\text{rhs - lhs}=-\frac{2209}{504000 n^6}$ or try to cancel the expression given in $(4)$; but here starts the problem because of the quartic.
If we cancel this term, that is to say using $$a=0.64646984692343004209959681070738106406762333935250\cdots$$ $$\text{rhs - lhs}= \frac{-1125 a^6+6000 a^5-8200 a^4+5600 a^3-2340 a^2+560 a-56}{2000 \,(3 a-1)^2\,n^8}+O\left(\frac{1}{n^9}\right)$$ $$\text{rhs - lhs}\sim \frac 1 {420\,n^8}+O\left(\frac{1}{n^9}\right)$$ which would be a significant improvement.
Making the root of the quartic rational, we could use any of the following value $$ \left\{\frac{2}{3},\frac{9}{14},\frac{11}{17},\frac{64}{99},\frac{1099}{1700},\cdots\right\}$$
Using the middle one $a=\frac{11}{17}$, $b=\frac{51}{8}$ and $\Delta=-\frac{1569}{1280}$ which seems to be sufficiently small to justify the propsed inverse. But, alas, $$\text{rhs - lhs} =\frac{10771}{990460800\, n^6}+O\left(\frac{1}{n^7}\right)$$
Using the simplest $a=\frac{1}{10}$, the results seem to be already decent. Using $n=2^p$, the results are $$\left( \begin{array}{cc} p & z\left(\frac{1}{10}\right) & \log_{10} \left(z'\left(\frac{1}{10}\right)\right)\\ 0 & 0.994119772820588 & -1.374 \\ 1 & 1.99876796576233 & -2.099 \\ 2 & 3.99978963553326 & -2.884 \\ 3 & 7.99996879285659 & -3.719 \\ 4 & 15.9999957334094 & -4.585 \\ 5 & 31.9999994416678 & -5.469 \\ 6 & 63.9999999285723 & -6.362 \\ 7 & 127.999999990967 & -7.260 \\ 8 & 255.999999998864 & -8.160 \\ 9 & 511.999999999858 & -9.062 \\ 10 & 1023.99999999998 & -9.965 \\ 11 & 2048.00000000000 & -10.87 \end{array} \right)$$
My qestions
- Did I introduced any bias from the very beginning ?
- Could it be possible to minimize at the same time $\Delta$ and $(\text{rhs - lhs})$ ?
- Would another approach be possible to have at the same time, with a very limited expansion, a good approximation of $H_n$ and its analytical inverse (without resorting to cubic equations).