A curious series related to the asymptotic behavior of the 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$$
where power towers are evaluated from top to bottom.
Let $a$ be a real number in the interval $1<a<e^{1/e}$. It is convenient to introduce a notation $$\lambda = -W(-\ln a),\tag3$$
where $W(z)$ denotes the principal branch of the Lambert W-function, satisfying $W(z) \, e^{W(z)}=z$. Note that the restrictions on $a$ imply $0<\lambda<1$. We may also observe that $\ln\ln a = - \lambda + \ln \lambda$.
The sequence $\{^n a\}$ converges to a limit
$$\lim_{n\to\infty} {^n a} = e^\lambda = \frac\lambda{\ln a}.\tag4$$
Its asymptotic behavior can be represented as
$${^n a} = e^\lambda - c_{\lambda} \cdot \lambda^n + O(\lambda^{2n}),\tag5$$
where $c_{\lambda}$ is some positive coefficient depending on $\lambda$ (or, equivalently, on $a$). It is known that the coefficients in $O(\lambda^{2n})$ and all higher-order terms can be expressed in a closed form via $c_{\lambda}$. But, apparently, the dependency between $\lambda$ and $c_{\lambda}$ has no known simple representation. Here is how its graph looks:

It appears to have a maximum near (or exactly at?) $\lambda = \ln 2$ that corresponds to $a=\sqrt2$. Numerical evidence suggests that $c_\lambda$ has a Taylor–Maclaurin series expansion with rational coefficients
$$c_\lambda = \frac{\lambda }{1!}+\frac{\lambda ^2}{2!}-\frac{2 \lambda ^3}{3!}-\frac{11 \lambda ^4}{4!}-\frac{44 \lambda ^5}{5!}-\frac{89 \lambda ^6}{6!}-\frac{636 \lambda
   ^7}{7!}-\frac{615 \lambda ^8}{8!}-...\tag6$$
(more numerators can be found here)
I could not find a formula for the coefficients (apparently, they are not yet in the OEIS), and I do not even have a proof that the coefficients given above are exact, so I am asking for your help with it. Is it possible to sum this series in terms of known special functions? Is it actually unimodal as is suggested by its graph? What is the exact location of its maximum?
Related questions: [1][2][3][4][5][6].
 A: Let $$c_n(\lambda)=\frac{e^\lambda-{^n a}}{ \lambda^n},\quad d_n(\lambda)=e^{-\lambda}c_n(\lambda)$$ ($d_n$ seems to be simpler). In particular initial functions $c_0(\lambda)=e^\lambda-1$ and $d_0(\lambda)=1-e^{-\lambda}\sim\lambda$ (as $\lambda\to 0$) have good power series of the variable  $\lambda$. We are intersted in the limits $$c_\lambda=\lim\limits_{n\to\infty}c_n(\lambda),\quad d_\lambda=\lim\limits_{n\to\infty}d_n(\lambda).$$  The recurrence
\begin{gather*}
d_{n+1}(\lambda)=\frac{1-e^{-\lambda^{n+1}d_n(\lambda)}}{\lambda^{n+1} }
\end{gather*}
follows  from the definition. By induction $d_n(\lambda)\sim\lambda$. Expanding the recurrence above we get
\begin{gather*}
d_{n+1}(\lambda)=d_n(\lambda)-\frac{\lambda^{n+1}d_n^2(\lambda)}{2}+\ldots\equiv d_n(\lambda)\pmod{\lambda^{n+3} }.
\end{gather*}
It means that the limit $d_\lambda$ (as a power series) does exist.
EDT. Numerators $a(k)$ of $c_\lambda$-coefficients are not very random. At least first $41$ coefficients from OP satisfy the congruences
\begin{gather*}
a(k)\equiv\begin{cases}
0,&\text{if } k\equiv 1\pmod{2 },\\
1,&\text{if }k\equiv 0\pmod{4 },\\
3,&\text{if }k\equiv 2\pmod{4 },\\
\end{cases}\pmod{4 }\quad(k\ge 4); \\
a(k)\equiv\begin{cases}
0,&\text{if }k\equiv 1,2\pmod{6 },\\
1,&\text{if }k\equiv0,3,4,5\pmod{6 },\\
\end{cases}\pmod{3 }\quad(k\ge 3).
\end{gather*}
A: This is not an answer, only an extended comment at V. Reshetnikov's request in his comment at my answer 
@VladimirReshetnikov - not so far now. I think this problem of a closed form
is also intimately related to the so far unsolved question of the "mag"-numbers as discussed in 
"Is there an “elegant” non-recursive formula for these coefficients?" 
which is focusing the series expansion for the Schröder-function $\sigma()$ as I've used it here.
Also your approach using q-binomials for a series "An explicit series representation for the analytic tetration with complex height" was/is much interesting for me - but although I had similar analyses involving q-binomials I didn't arrive at any conciser formulae than the formal algebraic descriptions following from the matrix-based concepts as discussed, for instance, here:  "Exponential polynomial interpolation" and also here  "Relation between binomial- and diagonalization method" but at least showing at its end the series (5) in your OP assuming $\lambda=\log 2$ with evaluated coefficients in real numbers and which gives also a short resume of the -as I christened it- "exponential polynomial interpolation" method which involves the q-binomials.
In general: the -in my opinion- most sophisticated explication of the coefficients of the function $t^x-1$ (and by conjugacy as well of $b^x$ ) in terms of the logarithm $u$ of $t$ and of $u^h$ comes via the diagonalization-approach as shown in "Eigendecomposition of triangular matrix-operators (here: Ut)" from where the coefficients of your formulae (5) and (6) can be derived, and even symbolically
But as well as with the q-binomial approach I did not yet find any simplifications/closed forms besides that (recursive)
basic generating schemes for that formulae/coefficients, I am really sorry. 
Of course, in case I'll find something worth sometime I'll surely link to that discussion here.
