Complete monotonicity of a sequence related to tetration
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.
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.