Finding a summation form of this sum - and its value

Let $s_n=\sum\limits_{k=1}^{n-1}a_k$, then $s_1=0$ and $s_{n+1}=s_n+n^{-s_n}$. We're looking for $s=\lim\limits_{n\to\infty}s_n$. $\newcommand{\ntail}[1]{\color{\LightBlue}{#1\cdots}}$

A closed form for $s$ looks elusive. But the asymptotics of $d_n=s-s_n$ attracts my curiosity. Moreover, it can be used to accelerate the convergence of $s_n$, hence to compute a more precise approximation of $s$. Namely, if we have a good enough $d_n\asymp d(n,s)$ as $n\to\infty$, and take a large $N$ such that we're still able to compute $s_N$, then $s_N^{(1)}=s_N+d(N,s_N)$ is a closer approximation of $s$ than $s_N$ itself, and the solution $s=s_N^{(2)}$ of $s=s_N+d(N,s)$ is even a much better one.

Let's illustrate this with the simplest possible $d(n,s)$. Since $n^s(d_n-d_{n+1})=n^{d_n}$ tends to $1$ as $n\to\infty$, it looks plausible to suppose that $d_n\asymp n^{1-s}/(s-1)$ (in fact this is not hard to prove). Now the computation suggested above, with $d(n,s)=n^{1-s}/(s-1)$, gives:

A more elaborate asymptotics of $d_n$ has the form $$d_n\asymp\sum_{p,q,r\geqslant 0}c_{p,q,r}n^{-p-q(s-2)}\log^r n$$ and can be obtained by successive refinement. Starting with $d_n\asymp n^{1-s}/(s-1)$ as above, we repeatedly compute the asymptotic expansion of $n^s(d_n-d_{n+1})-e^{d_n\log n}$, and replace $d_n$ with $d_n+t_n$, where $t_n$ is chosen so that the main term of the computed expansion vanishes.

I have implemented this for PARI/GP, representing $d_n$ by $f(n^{-1},n^{2-s},\log n)$ where $f(x,y,z)$ is a power series with $x\ll y\ll z$; then, for instance, $d_{n+1}=g(n^{-1},n^{2-s},\log n)$ where $$g(x,y,z)=f\big(x(1+x)^{-1},y(1+x)^{s-2},z+\log(1+x)\big).$$

The result of this computation, with $a=s-1$ denoted for an easier look, is \begin{align*} d_n\asymp\frac{n^{-a}}{a}&+\frac{n^{-a-1}}{2}+\frac{n^{-2a}}{4a^3}(2a\log n+1)+\frac{a+1}{12}n^{-a-2}\\&+\frac{n^{-2a-1}}{2a(2a+1)^2}\big((2a+1)(3a+1)\log n+a\big)\\&+\frac{n^{-3a}}{108a^5}(36a^2\log^2 n+33a\log n+11)+\ldots \end{align*} With this approximation used as $d(n,s)$, the computation of $s_N^{(2)}$ gives more digits: