A limit related to asymptotic growth of tetration
The tetration is denoted $^n a$, where $a$ is called the base and $n$ is called the height, and is defined for $n\in\mathbb N\cup\{-1,\,0\}$ by the recurrence $$ {^{-1} a} = 0, \quad {^{n+1} a} = a^{\left({^n a}\right)},\tag1$$ so that $${^0 a}=1, \quad {^1 a} = a, \quad {^2 a} = a^a, \quad {^3 a} = a^{a^a}, \, \dots \quad {^n a} = \underbrace{a^{a^{{.^{.^{.^a}}}}}}_{n\,\text{levels}}.\tag2$$ Let $a$ be a real number in the interval $e^{-e} < a < e^{1/e}$. It is known that the following limit exists $$L(a) = \lim_{n\to\infty} {^n a},\tag3$$ where $L(a)$ satisfies $a^{L(a)}=L(a)$. For example, $L\left(\!\sqrt2\right)=2$.
It is also known that $$\lim_{n\to\infty} \, \frac{L(a) - {^{n+1} a}}{L(a) - {^n a}} = \ln L(a).\tag4$$ Finally, it is known that the following limit exists $$C(a) = \lim_{n\to\infty} \, \frac{L(a) - {^n a}}{\left(\ln L(a)\right)^n}.\tag5$$ Apparently, no closed form for the function $C(a)$ is known. But the numerical evidence suggests the following conjecture (basically, this is the coefficient of the linear term in the Taylor series expansion of $C(a)$ near $a=1$):
$$C'(1) = \lim_{a\to1} \, \frac{C(a)}{a-1} \stackrel?= 1.\tag{$\diamond$}$$
How can we prove it? Can we find values of some higher-order derivatives? Are they all integers? Is there a general formula, recurrence or an efficient algorithm to compute them?
Related questions: [1][2][3].
I am aware of this answer as well as of this question in MO and its answer. There are, of course, some overlaps, but I proceed differently and answer the questions of the OP.
To fix notation, we consider $a$ close to 1, $b=b(a)=\log(a)$ and $L=L(a)$ the solution close to 1 of $L=a^{L}(=\exp(b\,L))$. Keep in mind that $L$ is close to 1 and $b$ is close to $0$. The sequence ${^n\!a}$, $n\in\mathbb N$ is defined by $^0 a=1$ and ${^{n+1}\!a}=\exp(b\ {^n\! a})$. We want to study the function defined by $$C=C(a)=\lim_{n\to\infty}\frac{L-{^n\! a}}{(bL)^n}\mbox{ for }a\neq1.$$ We show
Theorem: The function $C(a)$ has a power series around $a=1$. Its derivatives $C^{(n)}(1)$ at $a=1$ are integers for all $n$. These numbers can be computed using computer algebra.
Algorithms to compute the numbers $C^{(n)}(1)$ are also contained in other answers, but I have not seen a proof that they are all integers.
Part 1. Since $b$ is close to 0, it is more convenient to consider it as the basic variable and have $a=\exp(b)$ etc. For convenience, we do not indicate the dependence upon $b$ when unnecessary. As also used in other answers, the idea is to linearise the iteration ${^{n+1}\!a}=\exp(b\ {^n\! a})$ around its fixed point $L$. So we introduce the sequence $z_n={^n\!a}-L$ and the function $f(z)=f(b,z)=L\left(\exp(bz)-1\right)$. Then the iteration can be written $z_0=1-L$, $z_{n+1}=f(z_n)$ and we have $f(0)=0$, $f'(0)=bL$. The OP asks about properties of $C=-\lim_{n\to\infty}(bL)^{-n}z_n$.
Lemma 1: There is a uniquely determined series $T(b,z)=z+\sum_{k=2}^{\infty}d_k(b)z^k$ having a convergence radius larger than 1 for all sufficiently small $b$ such that $$T(b,f(b,z))=bL T(b,z)\mbox{ for all small }b\mbox{ and }|z|<1.$$
As a consequence, we have $C=-T(1-L)$: Indeed, on the one hand $(bL)^{-n}T(z_n)=(bL)^{-n}z_n(1+o(1))$ tends to $-C$, on the other hand, by induction using the functional equation, we have $(bL)^{-n}T(z_n)=T(z_0)=T(1-L)$. This shows that $C$ is a holomorphic function of small $b$ and hence in a small disk around $a=1$.
For the proof, we rewrite the functional equation for $T$. We first put $T(z)=z+V(b,z)$. The functional equation for $V$ is then \begin{equation}\label{eqv}\tag{1}V(b,z)=\frac1b(\exp(bz)-1-bz)+\frac1{bL}V(b,L(\exp(bz)-1)). \end{equation} Now we put $V(b,z)=z^2U(b,z)$ and obtain the functional equation \begin{equation}\label{equ}\tag{2}\begin{array}{rcl} U(b,z)=(\varphi(U))(b,z)&:=& b\,\frac{\exp(bz)-1-bz}{b^2z^2}\\&&+bL\left(\frac{\exp(bz)-1}{bz}\right)^2U(b,L(\exp(bz)-1)). \end{array}\end{equation} Now consider some $R>1$, $\delta>0$ and the Banach space $\mathcal B$ of all bounded holomorphic functions $U$ defined for $|z|<R$ and $|b|<\delta$, equipped with the maximum norm. If $\delta$ is sufficiently small then $\varphi:{\mathcal B}\to{\mathcal B}$ is a contraction and hence has a unique fixed-point by the Banach fixed-point Theorem.
The sequence of series $U_0=0$, $U_{n+1}=\varphi(U_n)$ also converges with respect to the $b$-valuation: Since we have $U_{n+1}(b,z)-U_n(b,z)= bL\left(\frac{\exp(bz)-1}{bz}\right)^2(U_n-U_{n-1})(b,L(\exp(bz)-1))$, the $b$-valuation of $U_{n+1}-U_n$ is at least $n+1$ and hence $U_n$ form a Cauchy sequence. Therefore it converges to a limit $U$, a series which we identify with the function $U$. As a consequence, the sequence defined by $V_0=0,$ $V_{n+1}(b,z)=\frac1b(\exp(bz)-1-bz)+\frac1{bL}V_n(b,L(\exp(bz)-1))$ converges in the $b$-valuation to the solution of (\ref{eqv}).
Part 2. As a consequence of the convergence of $U_n$ with respect to the $b$-valuation, the coefficient series $d_k(b)$ of $U$ have rational coefficients. Hence $C=-T(1-L)=-1+L-(1-L)^2U(1-L)$ also can be written as a series in powers of $b$ with rational coefficients. It is a bit tricky to prove that $C(a)=\sum_{k=1}^\infty \frac{c_k}{k!}(a-1)^k$ with integers $c_k$ which implies that the derivatives $C^{(m)}(1)=c_m$ are integers. This will be done in this second part of the answer. Observe that the convergence of the series for small $b$ or $z$ is not important here.
$\renewcommand{\AA}{{\mathcal A}}$ Let us introduce the set $\AA$ of all formal series $s(b)=\sum_{k=0}^\infty \frac{s_k}{k!}b^k$ where the $s_k$ are integers. Equivalently, this is the set of all formal series $s$ such that $s^{(m)}(0)=s_m$ is an integer for all integers $m\geq0$. It is easily seen that $\AA$ is closed under sums and products and that it is complete for the topology generated by the $b$-valuation. Hence the units are the elements $s$ of $\AA$ with $s(0)=\pm1$: we can write $\pm\frac1s=1+\sum_{k\geq1}(1\mp s)^k$. Important elements of $\AA$ are $e^b$ and $L$. For the latter observe that it can be written $L=-W(-b)/b=\sum_{m=0}^{\infty} \frac{(m+1)^{m-1}}{m!}b^m$ with Lambert's function $W$. We will need later that $1/L\in\AA$ and $\frac1b(e^b-L)=\sum_{m=2}^\infty \frac{1-(m+1)^{m-1}}m\frac{b^{m-1}}{(m-1)!}\in\AA$.
We will also use the set $\AA[[z]]$ of all formal power series with coefficients in $\AA$. A series $S=S(b,z)\in\AA[[z]]$ if and only if for all integers $k\geq0$ we have $\frac{\partial^k S}{\partial b^k}(0,z)\in{\mathbb Z}[[z]]$. $\AA[[z]]$ is closed under sums, products and partial derivatives and complete with respect to both the $b$- and the $z$-valuation. If $S(b,z),R(b,z)\in\AA[[z]]$, $R(b,0)=0$, then also $S(b,R(b,z))\in\AA[[z]]$ because it can be written as $\sum_{k\geq0} s_k(b)R(b,z)^k$ with the coefficients $s_k(b)\in\AA$ of $S(z)$.
Claim 1 If $S(b,z)\in\AA[[z]]$ with $S(0,z)=0$ then $P_k(b,z)=\frac1{k!}S(b,z)^k\in\AA[[z]]$ for all $k\geq1$.
For the proof, we proceed by induction. For $k=1$, there is nothing to prove. Now suppose that $P_{k-1}(b,z)\in\AA[[z]]$ and hence $\frac{\partial^m P_{k-1}}{\partial b^m}(0,z)\in{\mathbb Z}[[z]]$ for all $m\geq0$. Then $P_k(0,z)=0\in{\mathbb Z}[[z]]$ and since $\frac{\partial P_{k}}{\partial b}(b,z)=P_{k-1}(b,z)\frac{\partial S}{\partial b}(b,z) \in\AA[[z]]$, we have $\frac{\partial^m P_{k}}{\partial b^m}(0,z)\in{\mathbb Z}[[z]]$ for $m\geq1$.
As a consequence of Claim 1, we note
Claim 2 If $S(b,z)\in\AA[[z]]$ with $S(0,z)=0$ then $\exp(bS(b,z))-1\in b\AA[[z]]$.
For a proof, it is sufficient to write $\exp(bS(b,z))-1=\sum_{k\geq1}b^k P_k(b,z)$ with the $P_k$ of Claim 1.
Now we are in a position to prove the most important statement of part 2.
Lemma 2 The series $V$ solution of (\ref{eqv}) satisfies $bV(b,z)\in\AA[[z]]$.
Proof: We have seen that the sequence $V_n$ defined by $V_0=0,$ $V_{n+1}(b,z)=\frac1b(\exp(bz)-1-bz)+\frac1{bL}V_n(b,L(\exp(bz)-1))$ converges in the $b$-valuation to the solution $V$ of (\ref{eqv}). Therefore it is sufficient to prove that $bV_n\in\AA[[z]]$ for all $n$. For $n=0,1$ this is evident. Surprisingly, it is now better to combine two steps of the recursion, that is express $V_{n+2}$ in terms of $V_{n}$. We use again the series $f(b,z)=L(\exp(bz)-1)\in \AA[[z]]$ satisfying $f(0,z)=0$. We find \begin{equation}\nonumber\begin{array}{rcl} bV_{n+2}(b,z) & = & \exp(bz)-1-bz + \frac1{L^2}g(b,z)-\frac1L\,f(b,z)+\\ & &+ \frac1{b^2L^2} bV_n(b,b\,g(b,z))\mbox{ where } g(b,z)=\frac Lb(\exp(b\,f(b,z))-1). \end{array}\end{equation} Here, by Claim 2, $g(b,z)\in\AA[[z]]$ because $f(0,z)=0$. In order to prove that $bV_n\in\AA[[z]]$ implies $bV_{n+2}\in\AA[[z]]$, it remains to show that $bV_n(b,b\,g(b,z))\in b^2\AA[[z]]$ if $bV_n(b,z)\in\AA[[z]]$. If we write $bV_n(b,z)=\sum_{k\geq2}v_k^{(n)}(b)z^n$ with $v_k^{(n)}\in\AA$, then $$bV_n(b,b\,g(b,z))=\sum_{k=2}^\infty v_k^{(n)}(b) b^k g(b,z)^k$$ and it becomes evident that $bV_n(b,b\,g(b,z))\in b^2\AA[[z]]$. This completes the proof of Lemma 2.
Finally, we study $C=-T(b,1-L)$ which is regarded a series in powers of $b$ for the moment.
Claim 3 $C\in\AA$.
We have by (\ref{eqv}) $$\begin{array}{rcl}T(b,1-L)&=&1-L+V(1-L)\\&=& 1-L+\frac1b(\exp(b(1-L))-1)+\frac1{bL}V(b,L(\exp(b(1-L))-1)). \end{array}$$ Now we calculate $\exp(b(1-L))=e^b\exp(-bL)=e^b/L$ and hence $L(\exp(b(1-L))-1)=e^b-L=b\,q(b)$ with some element $q(b)\in\AA$ as we have seen at the beginning of part 2. Now we have, if $bV(b,z)=\sum_{k\geq2}v_k(b)z^k$ with $v_k(b)\in\AA$ according to Lemma 2, $$\begin{array}{rcl}T(b,1-L)&=&1-L+\frac1L q(b)+\frac1{b^2L}bV(b,bq(b))\\&=& 1-L+\frac1L q(b)+\frac1L\displaystyle\sum_{k=2}^\infty v_k(b)b^{k-2}q(b)^k\in\AA. \end{array}$$ This proves that $C=-T(1-L)\in\AA$, say $C=\sum_{k=1}^\infty\tilde c_k\frac1k!b^k$ with integers $\tilde c_k$.
In a last step, we show that $C=C(a)=\sum_{k=1}^\infty\tilde c_k\frac1k!(\log(a))^k$ can be written as $C=\sum_{k=1}^\infty c_k\frac1k!(a-1)^k$ with integers $c_k$. For this it is sufficient to show that all derivatives $u_k^{(m)}(1)$, $m\geq0$, are integers for the functions $u_k(a)=\frac1{k!}\log(a)^k$, $k\geq1$. This follows by induction in the same way as Claim 2.
This proves the second statement of the Theorem that all derivatives $C^{(m)}(1)$ are integers.
Part3. I indicate one way to calculate the above coefficients $\tilde c_k$ and hence also the $c_k$ of $C=\sum_{k=1}^\infty c_k\frac1k!(a-1)^k$, but this had also been done in other answers.
I use pari/gp where one can conveniently manipulate truncated series with remainder: $\sum_{k=0}^{K-1}d_k b^k+O(b^K)$ with rational coefficients $d_k$. As there seems to be no recursion for the coefficients $c_k$, I work with truncated series of two variables which pari/gp writes $\sum_{\ell=0}^{L-1} D_\ell z^\ell+O(z^L)$ with truncated series $D_\ell$ of one variable as above.
Then doubly truncated series for $V$ satisfying (\ref{eqv}), which we also denote $V$, can be found by the iteration $V_0=0$, $V_{n+1}(b,z)=\frac1b(\exp(bz)-1-bz)+\frac1{bL}V_n(b,L(\exp(bz)-1))$. Since the full series converge in the $b$-valuation to $V$, the truncated series will become stationary after $K$ steps. We use $L=\sum_{k=0}^{K-1}\frac{(k+1)^{k-1}}{k!}b^k+O(b^K)$ and also truncate the exponential. Then $C=-1+L-V(b,1-L)$ and the series in $a-1$ is obtained by substituting $b=\log(1+d)$, where $d$ stands for $a-1$. This method is probably not the fastest, but it works for small $K$.
We find as values of $c_k=C^{(k)}(1)$, $k=1,\dots,25$: $$\begin{array}{l} 1, 2, 6, 26, 120, 474, - 3500, - 169744, - 4739628, - 122528220, \\- 3244006128, - 89971866744, - 2643601630488, - 82449886989120,\\ - 2730313541889120, - 95853665484598656 , - 3561107748108889344,\\ - 139703010646898138688, - 5774800668716738596896,\\ - 250987866830927324395200, - 11446384958969565652481952,\\ - 546691342237064262980473248, - 27294921436196504583147875328,\\ - 1422140544716132719283821967232, - 77201569725471588510395303978880. \end{array}$$ It is not surprising that they are growing fast in modulus. It seems that the first six numbers are positive and all the following are negative, but I have no idea how to prove this. Maybe a study of the singularities of $C(a)$ on its circle of convergence might help.
Update: The natural variable in this context is not $b$, but $\lambda=bL=\log(L)$ since the recursion $z_{n+1}=L(\exp(b\,z_n)-1)$ converges for $z_0$ in a neighborhood of the fixed-point $L$ if and only $|\lambda|<1$. In the sequel we regard $L=\exp(\lambda)$, $b=\lambda\exp(-\lambda)$ and $a=\exp(b)$ as functions of $\lambda$ defined for all complex $\lambda$.
Conjecture: $C=C(\lambda)$ has a power series expansion having a radius of convergence $1$.
Remark: If this is true, then it is easy to show that the radius of convergence of $C=C(a)$ around $a=1$ is $e^{1/e}-1$ and $e^{1/e}$ is the only singularity on the circle of convergence. This in turn is an important step in showing that the coefficients $c_k$ have the same sign with finitely many exceptions.
The evidence supporting the conjecture is strong.
a) The functional equation for $U$ (from the proof of Lemma 1) as a function of $\lambda=bL$ and $z$
permits to construct it as a holomorphic bounded function of complex $\lambda$ , $|\lambda|<\rho<1$ and complex $z$, $|z|<\delta$ small. So $U(\lambda,z)$ and therefore also $T$ can be continued analytically
to $|\lambda|<1$ and small $z$.
b) For 4000 points $\lambda$ with $|\lambda|=0.99$, there is a path (actually invariant by $z\mapsto L(e^{bz}-1)$) connecting $^1\!a=a$ and $^{400}a$ such that the latter is at a distance $<0.04$ from the fixed point $L$. Therefore $T(b,^{400}\!a)$ can be defined using the series $T(b,z)$ and hence also $C(\lambda)=\lambda^{-400}T(b,^{400}a)$ is well defined for these points. This suggest that $C(\lambda)$ can be defined in this way for all $\lambda$, $|\lambda|=0.99$ and then, by the maximum principle, also for all $|\lambda|\leq0.99$.
The fact that it works for the radius $0.99$ suggest that it works for all radii $<1$.
This is just an attempt to come nearer to the problem, perhaps the last step can be seen as a proof for the right equality in your last equation. Maybe this relatively easy-to-see point of view has escaped your attention and you might do something more with this ansatz. (For the cursory reader the link to an answer of mine in MO(jun 2017) might be interesting, where a much related problem has also been discussed)
Let's first rewrite $a = 1+x$ and let's tend to look at our formula with some $x \to 0$ instead of $a \to 1$.
Then for the iterates $ ^n (1+x)$ we find power series with a remarkable pattern when increasing $n$:
Table 1: taylor-coefficients for (1+x)^^n (rowwise)
n|x^0 x x^2 x^3 x^4 x^5 x^6 x^7 x^8 x^9 x^10
---+-------------------------------------------------------------------------- ---------
1| 1 1 0 0 0 0 0 0 0 0 0
2| 1 1 1 1/2 1/3 1/12 3/40 -1/120 59/2520 -71/5040 131/10080
3| 1 1 1 3/2 4/3 3/2 53/40 233/180 5627/5040 2501/2520 8399/10080
4| 1 1 1 3/2 7/3 3 163/40 1861/360 33641/5040 8363/1008 22391/2160
5| 1 1 1 3/2 7/3 4 243/40 3421/360 71861/5040 54371/2520 69281/2160
6| 1 1 1 3/2 7/3 4 283/40 4321/360 102941/5040 85871/2520 61333/1080
7| 1 1 1 3/2 7/3 4 283/40 4681/360 118061/5040 106661/2520 81583/1080
8| 1 1 1 3/2 7/3 4 283/40 4681/360 123101/5040 115481/2520 93013/1080
9| 1 1 1 3/2 7/3 4 283/40 4681/360 123101/5040 118001/2520 97333/1080
10| 1 1 1 3/2 7/3 4 283/40 4681/360 123101/5040 118001/2520 98413/1080
This pattern has already been discussed in the tetration forum and I think also here in MSE.
The first interesting aspect is, that for increasing $n \to \infty$ the leading coefficients stabilize and we might formulate the powerseries for the fixpoint $L(a)$ or better $L(1+x)$ as limit of this sequence of powerseries, like
$$ L(1+x)= 1 + x + x^2 + 3/2 x^3 + 7/3 x^4 + 4 x^5 + 283/40 x^6 + 4681/360 x^7 \\ \qquad + 123101/5040 x^8 + 118001/2520 x^9 + 98413/1080 x^{10} \\ + 13580899/75600 x^{11} + O(x^{12}) \qquad \qquad \tag 1 $$
In your eq(5) you have in the numerator the difference between the fixpoint ($(1+x)$ to the infinite height) and $(1+x)$ to the $n$'th height.
Thus we might look in the differences of the above coefficients-table and the fixpoint's coefficients:
Table 2: differences (1+x)^^oo - (1+x)^^n (coefficients of taylor series)
n | x^0 x x^2 x^3 x^4 x^5 x^6 x^7 x^8 x^9 x^10
---+-----------------------------------------------------------------------------------
1| 0 0 1 3/2 7/3 4 283/40 4681/360 123101/5040 118001/2520 98413/1080
2| 0 0 0 1 2 47/12 7 1171/90 17569/720 78691/1680 2755171/30240
3| 0 0 0 0 1 5/2 23/4 281/24 2797/120 275/6 2730367/30240
4| 0 0 0 0 0 1 3 47/6 71/4 9247/240 11629/144
5| 0 0 0 0 0 0 1 7/2 61/6 101/4 8503/144
6| 0 0 0 0 0 0 0 1 4 51/4 103/3
7| 0 0 0 0 0 0 0 0 1 9/2 187/12
8| 0 0 0 0 0 0 0 0 0 1 5
9| 0 0 0 0 0 0 0 0 0 0 1
10| 0 0 0 0 0 0 0 0 0 0 0
Next, for the expression $L(1+x)$ we can find a power series for its logarithm which, by Pari/GP, is
$$ \ln(L(1+x))= x + 1/2 x^2 + 5/6 x^3 + 13/12 x^4 + 28/15 x^5 + 187/60 x^6+ 1781/315 x^7 \qquad \\ \qquad + 52247/5040 x^8 + 16487/840 x^9 + 569879/15120 x^{10} \\ \qquad + 17486729/237600 x^{11} + O(x^{12}) \qquad \qquad \tag 2 $$
and powers of it can easily be determined. For instance we find
$$ \ln(L(1+x))^8 = x^8 + 4 x^9 + 41/3 x^{10} + 39 x^{11} + O(x^{12}) \qquad \qquad \tag 3 $$
Then, if we put this in the numerator and denominator of eq 5, we get -for instance for iteration-height $n=8$
$$ \text{eq $5$ with $n=8$ and $x=a-1$ inserted as taylor series}\\ C_8(1+x)= x + x^2 + x^3 + 13/12 x^4 + x^5 + 79/120 x^6 - 25/36 x^7 + O(x^8) \tag 4$$ Checking with $L(a)=1.1$ and $a \approx 1.09051015407$ (computed from $L(a)$ using the LambertW-function) gives
(L(a)- a^^8)/ln(L(a))^8
%213 = 0.0995227877076
and the Pari/GP method for substituting a value into a truncated powerseries
subst(Pol(Ser(difflist[8,])/log(La)^8),x,a-1)
%216 = 0.0995228135777
Going one step further, the rhs in your last equation claims that for $n \to \infty$
$$ \lim_{a \to 1} C'(a) = {C(a) \over a-1} = {C_n(1+x)\over x} = 1 \tag 5$$
and we get, for instance with $n=8$ ,
$$ \begin{eqnarray}{ C_8(1+x)\over x} &= 1 &+ x + x^2 + 13/12 x^3 + x^4 + 79/120 x^5 - 25/36 x^6 + O(x^7) \\ \lim_{ x \to 0} {C_8(1+x)\over x} &= 1 \end{eqnarray}\tag 6 $$ which is -in principle- your conjecture.
Of course we should consider this example, where I use the taylor series based on $n=8$, with the limiting series for $n \to \infty$ as well. However, it occurs, that the leading part of this series up to the $n$'th coefficient does not change when $n$ is increased, so we get the same truncation of the taylor series if $n \to \infty$ were used
That series gives also the coefficients for the higher derivatives as asked (and as hinted to in the comments of the OP). We get by multiplication by factorials the sequence of coefficients
$$ \{1, 2, 6, 26, 120, 474, -3500, -169744, -4739628, -122528220,... \}
\tag 7$$
which is also known to OEIS (contributed 2017 by the OP's author). A simpler construction-pattern seems to be not known so far.
I'm a bit surprised about one aspect here. From the coefficients in the difference-table I'd expect, that, for instance, we divide the $6$'th row by the $8$'th power of the logarithm because the first two leading coefficients are then equal, and only one exponent of $x$ must be shifted. Which would give $$ { L(a) - ^{n-2} a \over \ln(L(a))^n} \cdot (a-1) = 1 - O( x^2) $$ instead. But, well, I don't see then the relation to your derivative-expression...