How to compute this constant with high precision $\sum_{n=1}^\infty \left(\frac{1}{a_n}-\frac{1}{(n+1) \ln (n+1)} \right)$

$\DeclareMathOperator{\li}{li}$ I can confirm the convergence of the series, but it is the slowest convergence I have ever encoutered in a series not constructed for this purpose. Here the approximation of $a_n$ using the logarithmic integral function is used. The tail of the series can be evaluated using the Euler-MacLaurin formula and the value $0.02225375...$ is found using the first 100 terms of the series and evaluation of the tail. All 8 decimals are correct. The approximation of $a_n$ is again crucial as well as the fact that the integrals occuring can be expressed using elementary functions.

For a proof of the convergence, we only need that the sequence $a_n$ defined by $a_1=2$ and $a_{n+1}=a_n+\log(a_n)$ satisfies $a_n=\li^{-1}(n)+O(\log(n)^2)=(n+1)(\log(n+1)+O(\log(\log(n))))$ which follows form this answer. It implies that $$\frac1{a_n}-\frac1{(n+1)\log(n+1)}= O\left(\frac{\log(\log(n))}{n\log(n)^2}\right).$$ Now, as shown in Yurij S' answer, $\sum_{n=3}^\infty\frac{\log(\log(n))}{n\log(n)^2} $ converges because the corresponding integral $\int_{3}^\infty\frac{\log(\log(t))}{t\log(t)^2}\,dt=\int_{\log(3)}^\infty\frac{\log(s)}{s^2}\,ds$ converges. This proves the convergence of the series. It also permits to show the slow convergence: If we sum the tail of the above series from $n=N$ to infinity then the integral is of the order $\log(\log(N))/\log(N)$. By summing $10^8$ terms, the tail of the given series then still is (a multiple of) $0.158...$. It is hopeless to try to evaluate the series by summing its terms even with a very fast computer.

To evaluate the given series, we use in a first step the approximation of this answer. There is a constant $C$ such that $$g(a_n)=n+C+O\left(\frac1{n\log(n)}\right),\mbox{ where } g(x)=\li(x)+\frac12\log(\log(x)).$$ Here $\li$ denotes the logarithmic integral function (see this page.) As $$g'(x)=\frac1{\log(x)}+\frac1{2x\log(x)},$$ we have in particular $(g^{-1})'(y)=\frac1{g'(g^{-1}(y))}\sim \log(g^{-1}(y))\sim \log(y)$ and hence $a_n=g^{-1}(n+C)+O\left(\frac1{n}\right)$. This yields $$\frac1{a_n}-\frac1{g^{-1}(n+C)}=O\left(\frac1{n^3\log(n)^2}\right).$$ Hence for $M>N$ sufficiently large $$\sum_{n=N}^M\left(\frac1{a_n}-\frac1{g^{-1}(n+C)}\right)=O\left(\frac1{N^2\log(N)^2}\right).$$

$\newcommand{\ds}{\displaystyle}$ In a secong step, we apply the Euler-MacLaurin formula with $h(t)=1/g^{-1}(t+C)$: $$\begin{array}{rcl}\ds\sum_{n=N}^M h(n)&=&\ds\int_N^M h(t)\,dt+\frac12h(N)+\frac12h(M)+\frac1{12}h'(M)-\frac1{12}h'(N) +O\left(\int_N^M|h^{(2)}(t)|\,dt\right)\\&=&\ds\int_N^M h(t)\,dt+\frac12h(N)+\frac12h(M)+O\left(\frac1{N^2\log(N)}\right).\end{array}$$ The proof of the estimate for the error terms is omitted. We now evaluate each of the terms in the above formula. Using the substitution $t=g(s)-C$, we find $$\begin{array}{ccl}\ds\int_N^M h(t)\,dt&=&\ds\int_N^M \frac1{g^{-1}(t+C)}\,dt= \ds\int_{g^{-1}(N+C)}^{g^{-1}(M+C)} \frac{1}sg'(s)\,ds=\\ &=&\ds\int_{g^{-1}(N+C)}^{g^{-1}(M+C)} \left( \frac1{s\log(s)}+\frac1{2s^2\log(s)} \right)\,ds\\ &=&\ds \left( \log(\log(s))+\frac12\li(1/s) \middle)\right|_{s=g^{-1}(N+C)}^{s=g^{-1}(M+C)}=\ds \left( \log(\log(s))+\frac12\li(1/s) \middle)\right|_{a_N}^{a_M}+O\left(\frac1{N^2\log(N)}\right). \end{array}$$ In the last equality, $a_n=g^{-1}(n+C)+O\left(\frac1{n}\right)$ has been used again. This gives altogether $$\ds\sum_{n=N}^M h(n)=\ds \left( \log(\log(s))+\frac12\li(1/s) \middle)\right|_{a_N}^{a_M}+\frac1{2a_N}+\frac1{2a_M}+O\left(\frac1{N^2\log(N)}\right)$$

Similarly, but simpler since we are dealing with an elementary function, we find $$\ds\sum_{n=N}^M\frac1{(n+1)\log(n+1)}=\log(\log(s))|_{N+1}^{M+1}+\frac12\frac1{(N+1)\log(N+1)}+ \frac12\frac1{(M+1)\log(M+1)}+O\left(\frac1{N^2\log(N)}\right).$$

Now we combine the above estimates and let $M$ tend to $\infty$. Using that $ \log(\log(a_M)) - \log(\log(M+1))\to0$ since $a_M\sim M\log(M)$, we obtain the wanted approximations of the tails $$\begin{array}{rcl}\ds\sum_{n=N}^\infty\left(\frac1{a_n}-\frac1{(n+1)\log(n+1)}\right)&=& \ds- \log(\log(a_N))-\frac12\li(1/a_N)+\frac1{2a_N}+\\&&+\ds \log(\log(N+1))-\frac12\frac1{(N+1)\log(N+1)} +O\left(\frac1{N^2\log(N)}\right).\end{array}$$

Finally we can evaluate the given series by summing its first $N-1$ terms and using the above approximation of the tail. Here $N\in\mathbb N$ is a parameter. We obtain $$\begin{array}{rcl}b=\ds\sum_{n=1}^\infty\left(\frac1{a_n}-\frac1{(n+1)\log(n+1)}\right)&=& \ds\sum_{n=1}^{N-1}\left(\frac1{a_n}-\frac1{(n+1)\log(n+1)}\right)\\&& \ds- \log(\log(a_N))-\frac12\li(1/a_N)+\frac1{2a_N}\\&&+\ds \log(\log(N+1))-\frac12\frac1{(N+1)\log(N+1)}\\&&\ds +O\left(\frac1{N^2\log(N)}\right).\end{array}$$

Using this formula for $N=100$ gives the value $0.02225375...$ for the series; comparison with the results for larger $N$ show that all 8 decimals after the dot are correct. Using $N=10^7$ as in the question and a refined (more complicated) formula having an error term $O\left(\frac1{N^3\log(N)}\right)$ provides a value $0.0222\ 5375\ 6202\ 8220\ 6538\ 1016...$ with $24$ correct decimals.