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$.
