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.