A curious series related to the asymptotic behavior of the tetration

Question

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

Answer 1

The coefficients of the series (6) can be found symbolically in $\lambda$ and rational numbers.

Consider some base $b$ with $1 \le b \le e^{1/e}$ , the fixpoint $t$ for $f(x)=b^x $ and $f(t)=t$. For bases $b$ in this interval we have a real fixpoint $t$ with $1 \le t \le e$ and $\lambda = \log(t) $ with $0 \le \lambda \le 1$ with $b^t = t$ (which is also attracting when $\lambda,t,b$ are not on the bounds of the resp. intervals).

$\ $

We assume now a conjugate function $d(z)=z/t-1$ and $d°^{-1}(z)=(z+1)\cdot t$ which I write here for shortness in sup-suffix notation $z´ = d(z) $ and $z`=d°^{-1}(z)$ .

The Schröder-mechanism, using this conjugacy, allows an invertible "schlicht" function $\sigma(x)$ such that $$ f°^h(x) = \left(\sigma°^{-1}( \lambda^h \cdot \sigma(x´))\right)`$$

Now the given problem assumes $x=1$ such that $x´ = 1/t-1 = \exp(-\lambda)-1$ and we have explicitely $$ f°^h(1) = \left(\sigma°^{-1}( \lambda^h \cdot \sigma(\exp(-\lambda)-1))\right) \cdot \exp(\lambda) + t$$ and because $\sigma(x) = x + O(x^2) $ (meaning to be called "schlicht" ) also its formal power series inverse is of this form and we get $$ f°^h(1) = t + \lambda^h c_\lambda + O(\lambda^{2h})$$ where $c_\lambda = \sigma(1/t-1) $
Interestingly, leaving $\lambda$ indeterminate and $t=\exp(\lambda)$ unevaluated as formal power series Pari/GP gives $c_\lambda$ as formal power series in the form as the OP in def. (6) provides its leading coefficients.

---

Derivation/implementation in Pari/GP

Here I show, illustrated by sample calculations with Pari/GP, that there is no numerical approximation needed to get the coefficients at $\lambda^k$ in definition (6). Instead there is a derivation on formal power series which provides exact rational coefficients. The method below reproduces perfectly your finding.

---

Preliminaries: a Carleman-matrix B associated to a function $f(x)$ which has a power series with more than zero radius of convergence (and generalized: can be Euler-summed in a certain interval of its argument) is initially only a collection of column-vectors $B[,c]$ ($c$ indexing the c-olumn), which contain the coefficients of the formal power series of $f(x)$ and more precisely of $f(x)^c$ for $c=0 \ldots \infty$ in the column of index $c$ , and is thus in both direction of infinite size. It is an important special case when the power series of some function $f(x)$ has no constant term: then the shape of $B$ is lower triangular and the analytic handling of such Carlemanmatrices has been substantially established. We'll need this in the derivation below.
_{It is not trivial to show, that such a collection of columns-vectors can indeed be used as Matrix, applicable to concepts of matrix-algebra and -analysis, for instance raising to powers and especially diagonalization - but there is reliable literature about this so that we can omit such a step here. One pioneer is surely Eri Jabotinsky, see for instance this 1963-article "Analytic iteration" (introducing the term "representation matrix")}

After the Carlemanmatrix B contains the coefficients of formal power series we need also a construct $V(x)$ as row vector $V(x)=[1,x,x^2,x^3,\ldots]$ which I call of "Vandermonde"-type here. Using $V(x)$ as diagonal-matrix I pre-superscribe this like $ \;^dV(x)$

---

For instance the Carlemanmatrix $B_e$ for the function $f(x) = \exp(x)$ looks like

   1      1       1        1         1       ...
   0      1       2        3         4       ...
   0     1/2!   2^2/2!   3^2/2!    4^2/2!    ...
   0     1/3!   2^3/3!   3^3/3!    4^3/3!    ...
   0     1/4!   2^4/4!   3^4/4!    4^4/4!    ...
  ...    ...     ...      ...       ....     ...  

The clue with such Carlemanmatrices is, that the left-multiplication of some $B$ with a vector $V(x)$ as in $$V(x) \cdot B = V(f(x)) \tag 1$$ gives a resulting vector of the same structure $V(f(x)) = [ 1 , f(x), f(x)^2, f(x)^3, ... ]$ .
Of course, using the only "relevant" column in $B$ alone , then this gives the scalar expression $f(x)$: $$V(x) \cdot B [,1] = f(x) \tag {1.1}$$

This means: a composition of a function with itself $f(f(x))$ (in terms of its power series) is mapped to a repeated multiplication with its Carleman-matrix or in other words: is mapped to multiplication with powers of its Carlemanmatrix.

The further clue is now, that for appropriate matrices a fractional power can be determined and for that cases we can thus easily define power series for the fractional iterates of the associated function.

---

Technically this means that the concept of diagonalization can be transferred to that Carleman-ansatz (of course for suitable Carlemanmatrices $B$). A case, for which this is widely accepted is the case where the Carleman-matrices for some function $f(x)$ is lower triangular.

We have the function $f(x) = b^x $ and by the definition in the OP $1 \lt b \lt e^{1/e}$ . With this restriction, $f(x)$ has a fixpoint $t$ such that $f(t)=t$ (or: $b^t=t$ ) and also this fixpoint is attracting (which is important here). We can, as I wrote in the first part of text above, apply conjugacy to arrive at the function $g(x)=f(x`)´ = t^x-1 = ux + (ux)^2/2! + O(x^3) $ (for notational convenience I write $u$ instead of $\lambda$ in the following).

Given the lower triangular matrix of Stirlingnumbers $2$'nd kind $S2$ and its factorially similarity-scaled transform which I denote here as $fS2F$ we have the lower triangular Carlemanmatrix $U_t$ for $g(x)$ by $$ U_t = \; ^dV(u) \cdot fS2F \tag {2.1}$$

---

The top-left edge of $U_t$ looks like

  1       .          .             .             .             .      
  0    u/1!          .             .             .             .      
  0  u^2/2!     u^2*2!/2!          .             .             .      
  0  u^3/3!   3*u^3*2!/3!     u^3*3!/3!          .             .      
  0  u^4/4!   7*u^4*2!/4!   6*u^4*3!/4!     u^4*4!/4!          .      
  0  u^5/5!  15*u^5*2!/5!  25*u^5*3!/5!  10*u^5*4!/5!     u^5*5!/5!   

or with cancelled entries

  1          .         .        .      .    .
  0          u         .        .      .    .
  0    1/2*u^2       u^2        .      .    .
  0    1/6*u^3       u^3      u^3      .    .
  0   1/24*u^4  7/12*u^4  3/2*u^4    u^4    .
  0  1/120*u^5   1/4*u^5  5/4*u^5  2*u^5  u^5

and taking the coefficients from the second column $$ g(x) = \sum_{k=1}^\infty u^k \cdot \frac{x^k}{ k!} \tag {2.2}$$

---

The formal power series for $h$'th fractional iterates $g°^h(x)$ are now generated by fractional powers of $U_t$ instead and this can be done by diagonalization:

$$ U_t = M \cdot D \cdot W \qquad \qquad \text{where I denote $W$ for } M^{-1} \tag {3.1}$$ where also the diagonal entries in $D$ are the consecutive powers of $u$ $$ D = \; ^d V(u) \text{ and } D^h = \; ^dV(u^h) \tag {3.2}$$ such that $$ U_t^h = M \cdot \; ^dV(u^h) \cdot W \tag {3.3} $$

Unfortunately the current Pari/GP-version's diagonalization procedure cannot exploit the advantage of $U_t$ being triangular and also cannot keep the parameter $u$ symbolically. But it is not difficult to have a small recursive routine to generate $M$ , $D$ and $W (=M^{-1})$ for the triangular Carleman-case; the additional norming of $M$ to have a unit-diagonal is also included and with this the matrix $M$ is also of the Carleman-type (for details and an implementation in Pari/GP see the explanation in my linked article). The matrix M is now the Carlemanmatrix for the above introduced Schröder-function $\sigma(x)$ which takes its coefficients from the second column in M.

The top-left of the matrix $M$ looks like

  1                                    .             .      .
  0                                    1             .      .
  0                           u/(-2*u+2)             1      .
  0      (2*u^3+u^2)/(6*u^3-6*u^2-6*u+6)      u/(-u+1)      1

The entries in the (only relevant) $2$'nd column have a partially discernable pattern; I've described this in earlier discussions, see for instance here, but which is too much to show at this place.

---

Having this we can write for the $h$'th fractional iterate $g(x)$ of some $x$:

$$ \begin{array}{rll} V(x) \cdot U_t^h &= V( g°^h(x)) \\ &= V(x) \cdot (M \cdot \;^dV(u^h) \cdot W) \qquad \qquad \text{ and by associativity}\\ &= (V(x) \cdot M) \cdot \;^dV(u^h) \cdot W \end{array} \tag {4.1}$$

Because $M$ is of Carlemantype we have also $$ V(x) \cdot M = V(\sigma(x)) \tag {4.2}$$ or -using also conjugacy- $$ V(x´) \cdot M = V(\sigma(x´)) \tag {4.3}$$ and it is relevant here, that $\sigma(x´)$ is identical to the (normed) Schröder-function for the iteration of $b^x$.

In the problem as given in the OP we assume $x=1$ or $x´=1/t-1 = \exp(-u)-1$ and to have the Schröder-function symbolically we have to look at $$ V(1´) \cdot M = V(\exp(-u)-1) \cdot M = V(\sigma(\exp(-u)-1)) \tag {4.4} $$ and this can still be done symbolically in $u$; having M in its symbolical representation, Pari/GP is able to give the leading part of $\sigma(\exp(-u)-1)$ by formal expansion of $\exp(-u)$ as power series taken as argument of $\sigma(\cdot)$ with still rational coefficients.

Having the left part of $(4.2)$ evaluated we need now $$ V(\sigma(\exp(-u)-1)) \cdot \; ^dV(u^h) \cdot W[,1] = y =\sigma°^{-1}(u^h \cdot \sigma(1´))\tag {4.5}$$ by the relevant $2$'nd column of $W$ only.

Note that this gives a power series in $u^h$ with the (constant and rational) coefficients in $W[,1]$ and powers of the expression $\sigma(\exp(-u)-1)$ . To help Pari/GP to keep $u^h$ as a single symbol, we need to substitute it by the name of a unevaluated variable, say $u_h$ or $v$ and re-substitute it later after series-expansions.

The final result $f°^h(1) $ occurs by the inverse conjugacy on $y$ $$f°^h(1) = y` = (y+1)\cdot t = y\cdot t + t =t + y \cdot \exp(u) \tag {4.6} $$ which can still be expressed as formal power series in the indeterminate $u$ and $u^h$.

The OP asks now for the first coefficient in $y \cdot \exp(u) $ and namely for its series expression when $u$ and $v (=u^h) $ are both kept as indeterminate parameter. Pari/GP is still able to evaluate the last expression symbolically and arrives at the same set of coefficients as given in eq (6) in the OP.

---

Additional computation: a short recalculation gives also the leading terms of the powerseries of the following coefficients in (5) of the OP, which I index here as $c_{\lambda,k}$ :

$$\small \begin{array} {r} c_{\lambda,1}=& - 1 u^1/1!& - 1 u^2/2!& + 2 u^3/3!& + 11 u^4/4!& + 44 u^5/5!& + 89 u^6/6!& + 636 u^7/7!& + O(u^8) \\ c_{\lambda,2}=&&& - 3 u^3/3!& - 12 u^4/4!& - 5 u^5/5!& + 150 u^6/6!& + 1463 u^7/7!& + O(u^8) \\ c_{\lambda,3}=&&&&& - 40 u^5/5!& - 240 u^6/6!& - 840 u^7/7!& + O(u^8) \\ c_{\lambda,4}=&&&&&&& - 1260 u^7/7!& + O(u^8) \\ \end{array} $$ and recalling that $\; ^hb $ is my notation for your $ \; ^na $ we get $$ \; ^na = \; ^hb = t +c_{\lambda,1} u^h + c_{\lambda,2} u^{2h} + c_{\lambda,3} u^{3h} + \cdots \tag 5$$ Remark: I don't think that a closing term $O(u^{4h})$ makes really sense here because we assume $h>-2$ to be of any real value, especially fractional, and I'm not firm with using that notation in such cases.

Answer 2

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*}

Answer 3

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.