Complete monotonicity of a sequence related to tetration Let $\Delta$ denote the forward difference operator on a sequence:
$$\Delta s_n = s_{n+1} - s_n,$$
and $\Delta^m$ denote the forward difference of the order $m$:
$$\Delta^0 s_n = s_n, \quad \Delta^{m+1} s_n = \Delta\left(\Delta^m s_n\right).$$
We say that a sequence $s_n$ is completely monotone iff
$$(-1)^m \Delta^{m} s_n > 0,$$
i.e. the sequence $s_n$ itself is positive and decreasing, its first differences $\Delta s_n$ are negative and increasing (that is, decreasing in absolute values), its second differences $\Delta^2 s_n$ are positive and decreasing  (like the sequence $s_n$ itself), and so on. So, the sequence of second differences of a completely monotone sequence is also completely monotone (and so is the sequence of its differences of any even order).
Let $a$ be a real number in the interval $1<a<e^{1/e}$. Let ${^n a}$ denote tetration:
$${^0 a}=1, \quad {^{(n+1)} a} = a^{\left({^n a}\right)}.$$
The sequence ${^n a}$ is strictly increasing and converges to the limit
$${^\infty a} = \lim_{n\to\infty}\left({^n a}\right) = \frac{W(-\ln a)}{-\ln a},$$
where $W(x)$ is the Lambert $W$-function. Let $s_n = {^\infty a} - {^n a}$. It appears that the sequence $s_n$ is completely monotone. How can we prove it?
 A: So from my comments above, $\lambda$ is the multiplier at the fixed point.  $a^{L+x}=L+\lambda x + O x^2\;\;\ln_a(L+x) = L +\frac{x}{\lambda}  + O x^2 $
$$L= \frac{W(-\ln a)}{-\ln a}\;\;\;\lambda=\ln(L)$$
For the bases in question $1<L<e,\;\;0<\lambda<1$.  Then an approximate solution for the superfunction for large values of z is $$S_1(z)=L-\lambda^z=L-e^{z \ln(\ln(L))}$$
$S_1(x)$ is a completely monotonic analytic function at the real axis, where the odd derivatives are all >0, and the even derivatives are all <0, which we will use as the definition of a completely monotonic analytic function.   This is because $\ln(\ln(L))<0$ for all bases in question; and we substitute that into the Taylor series for exp(z).  $\ln_a(x)$ is also completely monotonic from the Taylor series of $\ln(x)$, for the Op's bases, in the domain $0<x<L$, which is the domain we will use in the iterations below.  From that initial $S_1(z)$ estimate, we generate a sequence of better and better estimates, which converges to the Schröder function solution by lemma1 which is described below.
$$S_n(z)=\log_a(S_{n-1}(z+1))$$
We can show that $S_n(x)$ is completely monotonic if $Re(x+1)>S^{-1}_{n-1}(0)>$; by lemma2, which is that the composition of two completely monotonic analytic functions is completely monotonic in the appropriate range. In the limit, $S(x)$ is completely monotonic at the real axis if $Re(x+1)>S^{-1}(0)$
$$S(z) = \lim_{n\to\infty} S_n(z)\;\;\;S(z+1)=a^{S(z)} \;\;\; \lim_{x\to\infty}S(x)=L$$
lemma 1, the sequence of $S_n(z)$ functions; $S(z)$ converges to the   Superfunction generated from the Schröder equation.  I would prove this by generating error estimates from the Superfunction generated from the formal inverse Schröder equation which is periodic; $p=\frac{2\pi i}{\ln{\lambda}}$ so $S(z)$ can be expressed as an analytic Fourier series $S(z)=L+\sum a_n(-1)^n\lambda^{n z}$ where $a_1=1$ and $a_n$ are the formal inverse Schröder function Taylor series coefficients and $S(z)$ converges if $\Re(z+1)>S^{-1}(0)\;\;$ end lemma 1
Then we shift by a constant to get the desired sexp function.
$$\text{sexp}_a(z)=S(z+k)\;\;\;k=S^{-1}(0)+1\;\;\;S(k)=1$$
Finally, $\text{sexp}_a(z)$ is completely monotonic at the real axis if z>-2.  Then the sequence $^n a=\text{sexp}_a(n)$ is also a completely monotonic sequence by lemma3 that if an analytic function is completely monotonic over a given range, then a sequence of equally spaced samples of that function in that range is also completely monotonic.
So, if we prove lemma(1,2,3), then we have the Ops result.  Lemma1 would come from a good book on complex dynamics, Lennart Carlson's book, would be the place to start, but I have seen it mentioned before.  The other two Lemmas should be well known in the theory of monotonic functions.  Real valued Schröder superfunction solutions will often be fully monotonic; so presumably this is known.
There is a possible alternative approach to lemma1, that $S(x)$ is completely monotonic.  As noted above, the formal superfunction solution can be expressed as a Fourier series: $S(z)=L+\sum a_n(-1)^n\lambda^{n z}$ Then $a_1=1;$ and by observation $a_n(-1)^n<0 \;\forall n$.  This sequence leads immediately to a completely monotonic solution for $S(z)$ at the real axis although I haven't tried proving the observation.
A: This answer is an attempt to explain lemma1.  Why is $S_1(z)$ the correct initial approximation in the previous answer, and why does the sequence converge?
$$S_1(z)=L-\lambda^z=L-e^{z \ln(\ln(L))}\;\;\;S_n(z)=\log_a(S_{n-1}(z+1))$$
Starting with ideas from the previous answer, along with the definition of the inverse Schröder function, which we will call $h(z)$, where the desired $S(z)$ function follows directly from $h(z)$ as follows: $S(z)=\sum a_n(-1)^n\lambda^{n z}\;\;$ where $a_n$ are the Taylor series coefficients of $h(z)$.  We start with the fixed point L from the previous answer.  $\ln_a(L+z) = L +\frac{z}{\lambda} + O z^2\; $ and the multiplier at the fixed point is $\lambda$.
$$L= \frac{W(-\ln a)}{-\ln a}\;\;\;\lambda=\ln(L) = \ln(a) \cdot L; \;\:1<L<e;\;\;0<\lambda<1$$
$$h(z)=\sum_{n=0}^{\infty}a_n z^n\;\;\;a_0=L;\;a_1=1\;\;\;h(\lambda z) = a^{h(z)}\;\;\;\ln_a(h(\lambda z))= h(z)$$
$$S(z) = h(\lambda^z)\;\;\;S(z+1)=a^{S(z)}\;\;\;\lim_{\Re(z) \to \infty}S(z)=L$$ 
lemma4: if all of the odd Taylor series coefficients of $h(z)$ at z=0 are positive, and all of the even Taylor series coefficients are negative, then $S(z)$ is fully monotonic, since $0<\lambda<1$, and $\lambda^{nz}=-\exp(n\ln(\lambda) z)$ is fully monotonic for all n, since $\ln(\lambda)$ is a negative number, and $a_n(-1)^n$ is also a negative number, so $S(z)$ is fully monotonic under that case.
There is a formal solution for $h(z)$, as the inverse of the Schröder equation.  But alternatively, we start with $h_1(z)=L+z$ and iterate as follows.
$$h_{m+1}(z)=\ln_a(h_m(\lambda z))\;\;\; h(z)=\lim_{m\to\infty}h_m(z)\;\;\;h(z)= \ln_a(h(\lambda z))$$
$$h_2(z)=\ln_a(h_1(z)) =\ln_a(L+\lambda z) $$
$$\ln_a(L+z) =\ln_a(L(1+z/L)) = L + \sum_{n=1}^{\infty}\frac{(-1)^{n+1}z^n}{n \cdot \ln(a) \cdot L^n}$$
$$h_2(z)= \ln_a(L+\lambda z) = L + \sum_{n=1}^{\infty}\frac{(-1)^{n+1} \cdot \lambda \cdot z^n}{n \cdot \ln(a) \cdot L^n} $$
$$h_2(z)= \ln_a(L+\lambda z) =  L + \sum_{n=1}^{\infty}\frac{(-1)^{n+1} \cdot z^n}{n \cdot L^{n-1}} \;\;\; a_0=L; a_1=1$$
Notice that $\ln_a(z)$ has all odd Taylor series coefficients positive, and all even Taylor series coefficients negative, just as desired.  Notice that $h_2(z)$ has all odd Taylor series coefficients positive, and all even Taylor series coefficients negative, just as desired. Similar to lemma2, the composition $h_3(z)=\ln_a \circ h_2(\lambda z)$ has the same property, so $h_3(z)$ has all odd Taylor series coefficients positive, and all even Taylor series coefficients negative, just as desired.  And by induction, the function $h_m(z)$ has all odd Taylor series coefficients positive, and all even Taylor series coefficients negative, just as desired, so this is also true of $h(z)$.  lemma5: the $h_m(z)$ sequence converges, which has been shown in the literature of complex dynamics.  Also, notice that iterating from $h_{m+1}$ from $h_m$ that $a_0=L$, and $a_1=1$, for the Taylor series of each version of $h_m$.  This approach is identical but a little bit more rigorous than the iterations in the previous answer, since it directly shows the stability in iterating $h_n(z)$, which makes proving convergence of lemma5 easier.  Let us say the golden values of $a_n$ for the formal inverse Schröder function are $g_n$, and we are doing m iterations to generate $h_m(z)$. Then by observation, $\forall n>2, \;a_n=g_n+O\lambda^{m+n}\;$ after m iterations. Formalizing this error estimate would be a way of showing that the $h_m(z)$ sequence converges to $h(z)$, and that lemma5 is valid.
So if lemma2, lemma4, and lemma5 hold, the conclusion is that $S(z)$ is fully monotonic at the real axis in its range of analyticity, where $\Re(z+1)>S^{-1}(0)$
Finally, as in the previous answer we generate the desired sexp function.
$$\text{sexp}_a(z)=S(z+k)\;\;\;k=S^{-1}(0)+1\;\;\;\text{sexp}_a(z)=a^{\text{sexp}(z-1)}$$
$\text{sexp}_a(z)$ is completely monotonic at the real axis if z>-2.  Then the sequence $^n a=\text{sexp}_a(n)$ is also a completely monotonic sequence by lemma3 that if an analytic function is completely monotonic over a given range, then a sequence of equally spaced samples of that function in that range is also completely monotonic.
