1 Introduction

The near-constancy phenomenon indicates that the tail of limit distributions of simple branching processes has a power-law asymptotic perturbed by extremely small oscillations. This phenomenon has a long history starting from the works [1, 2] and then continues to be actively studied in [3,4,5]. Note also interesting and important applications in physics and biology, see [6,7,8]. Most of the research is devoted to the martingale limit, see, e.g. [9], for branching processes, which represents a continuous function of a real argument. For this function, the power-law asymptotic multiplied by some multiplicatively periodic factor is proven. In our research, we focus on Taylor coefficients of another limit distribution, which satisfies some Schröder functional equation. The analysis of the asymptotic of the power-series coefficients of solutions of functional equations is difficult and may strongly depend on the functional equation itself, see [10,11,12,13]. For a sufficiently general type of Schröder functional equations, we provide an approximation formula for the Taylor coefficients of the solutions of these equations. The approximation formula is based on special binomial coefficients, from the asymptotic of which one may extract the leading power-law term and periodic oscillations with small amplitudes. There are two main differences between this study and other similar ones: (1) we derive a complete asymptotic; (2) all the procedures are implemented numerically.

The Galton–Watson process is defined by

$$\begin{aligned} X_{t+1}=\sum _{j=1}^{X_t}\xi _{j,t},\ \ \ X_0=1,\ \ \ t\in {{\mathbb {N}}}\cup \{0\}, \end{aligned}$$
(1)

where all \(\xi _{j,t}\) are independent and identically-distributed natural number-valued random variables with the probability generating function

$$\begin{aligned} P(z):={\mathbb {E}}z^{\xi }=p_0+p_1z+p_2z^2+p_3z^3+.... \end{aligned}$$
(2)

It is well known, see [14], that

$$\begin{aligned} {\mathbb {E}}z^{X_t}=\underbrace{P\circ ...\circ P}_{t}(z),\ \ \ t\in {{\mathbb {N}}}. \end{aligned}$$
(3)

For simplicity, let us assume that P is a polynomial of degree \(N\geqslant 2\), which is commonly the case in practice. Another assumption is \(p_0=0\) and \(p_1\ne 0\). To avoid complications due to periodicity, we assume also the existence of k such that \(p_kp_{k+1}\ne 0\), see Remark at the end of this Section. Note that the case \(p_0\ne 0\) can usually be reduced to the (supercritical) case \(p_0=0\) with the help of Harris-Sevastyanov transformation, see [15]. The so-called Böttcher case \(p_0=p_1=0\) is completely different. Both cases are discussed in the final Remark of Sect. 3.

Under the assumptions \(p_0=0\) and \(0<p_1<1\), we can define

$$\begin{aligned} \Phi (z)=\lim _{t\rightarrow \infty }p_1^{-t}\underbrace{P\circ ...\circ P}_{t}(z), \end{aligned}$$
(4)

which is analytic at least for \(|z|<1\). In fact, \(\Phi (z)\) is analytic inside an open component of the filled Julia set for P(z) containing the ball \(B_1:=\{z\in {{\mathbb {C}}}:\ |z|<1\}\), since \(|P(z)|<|z|\) for \(z\in B_1\). Recall that the points \(z\in {{\mathbb {C}}}\) with the bounded orbits \(\{\underbrace{P\circ ...\circ P}_{t}(z)\}_{t\geqslant 0}\) form the filled Julia set. The boundary of the filled Julia set is called the Julia set. Various properties and precise definitions of the objects related to holomorphic dynamics are available in a nice introductory book [16].

The function \(\Phi \) satisfies the Scröder-type functional equation

$$\begin{aligned} \Phi (P(z))=p_1\Phi (z),\ \ \ \Phi (0)=0,\ \ \ \Phi '(0)=1. \end{aligned}$$
(5)

Due to (2)–(4), the Taylor coefficients \(\varphi _n\) of

$$\begin{aligned} \Phi (z)=z+\varphi _2z^2+\varphi _3z^3+... \end{aligned}$$
(6)

describe the relative limit densities of the number of descendants

$$\begin{aligned} \varphi _n=\lim _{t\rightarrow +\infty }\frac{{\mathbb {P}}\{Z_t=n\}}{{\mathbb {P}}\{Z_t=1\}}. \end{aligned}$$
(7)

Our goal is to find simple approximations of \(\varphi _n\). The first moment

$$\begin{aligned} E:=P'(1)=\sum _{j=1}^{N}jp_j \end{aligned}$$
(8)

satisfies \(E>1\), since \(\sum _{j=1}^{N}p_j=1\) and \(0<p_1<1\). Thus, we can define

$$\begin{aligned} \Psi (z)=\lim _{t\rightarrow \infty }E^t(1-\underbrace{P^{-1}\circ ...\circ P^{-1}}_{t}(z)), \end{aligned}$$
(9)

which is analytic in some neighborhood of \(z=1\). In fact, it is not difficult to show that \(\Psi \) is analytic on (0, 1), including its neighborhood, since P is strictly monotonic on this interval and \(P(0)=0\), \(P(1)=1\). The function \(\Psi \) also satisfies the Schröder-type functional equation

$$\begin{aligned} \Psi (P(z))=E\Psi (z),\ \ \ \Psi (1)=0,\ \ \ \Psi '(1)=-1. \end{aligned}$$
(10)

Everything is ready to define

$$\begin{aligned} \Theta (z)=\Phi (z)\Psi (z)^{\frac{-\ln p_1}{\ln E}}, \end{aligned}$$
(11)

which is analytic on (0, 1). Moreover, it satisfies the functional equation

$$\begin{aligned} \Theta (P(z))=\Theta (z), \end{aligned}$$
(12)

see (5), (10), and (11). The “polynomially” periodic function \(\Theta \) can be considered as some precursor of the multiplicatively periodic Karlin–McGregor function \(K(z)=\Theta (\Psi ^{-1}(z))\), see (10) and (12). The function \(\Psi ^{-1}\) can be defined as

$$\begin{aligned} \Pi (z):=\Psi ^{-1}(z)=\lim _{t\rightarrow +\infty }\underbrace{P\circ ...\circ P}_{t}(1-E^{-t}z), \end{aligned}$$
(13)

see (9). The function \(\Pi (z)\) is an entire function satisfying the Poincaré-type functional equation

$$\begin{aligned} P(\Pi (z))=\Pi (Ez),\ \ \ \Pi (0)=1,\ \ \ \Pi '(0)=-1. \end{aligned}$$
(14)

The next step is to introduce

$$\begin{aligned} K^*(z):=K(E^z)=\Theta (\Pi (E^z))=\Phi (\Pi (E^z))p_1^{-z}, \end{aligned}$$
(15)

which is an 1-periodic function defined on some strip \(\{z:\ |\mathop {\textrm{Im}}\nolimits z|<\vartheta ^*\}\) with \(\vartheta ^*>0\). Thus, we have the Fourier series

$$\begin{aligned} K^*(z)=\sum _{m=-\infty }^{+\infty }\theta _me^{2\pi {\textbf{i}} mz},\ \ \ |\mathop {\textrm{Im}}\nolimits z|<\vartheta ^* \end{aligned}$$
(16)

with the exponentially rapidly decreasing coefficients \(\theta _m=\overline{\theta _{-m}}\). We denote \({\textbf{i}}=\sqrt{-1}\). Using the definitions \(\Psi ^{-1}=\Pi \) and (15), along with (16) we obtain

$$\begin{aligned} \Theta (z)=\sum _{m=-\infty }^{+\infty }\theta _m\Psi (z)^{\frac{2\pi {\textbf{i}} m}{\ln E}}, \end{aligned}$$
(17)

which with (11) leads to

$$\begin{aligned} \Phi (z)=\sum _{m=-\infty }^{+\infty }\theta _m\Psi (z)^{\frac{2\pi {\textbf{i}} m+\ln p_1}{\ln E}}. \end{aligned}$$
(18)

There are two main ideas: (1) the amplitudes of oscillations \(\theta _m\) are usually extremely small; (2) \(z=1\) is a unique singular point for \(\Phi (z)\) on the unit circle, see Remark at the end of this Section. Thus, we can take few terms to obtain a good approximation

$$\begin{aligned} \Phi (z)\approx \sum _{m=-M}^{M}\theta _m\Psi (z)^{\frac{2\pi {\textbf{i}} m+\ln p_1}{\ln E}}. \end{aligned}$$
(19)

After that, we substitute the expansion of \(\Psi (z)\) near \(z=1\), see (9) and (10),

$$\begin{aligned} \Psi (z)= & {} 1-z+O((1-z)^2)\ \ \textrm{or}\ \ \Psi (z)=1-z\nonumber \\{} & {} +\frac{P''(1)}{2(E^2-E)}(1-z)^2+...+O((1-z)^R) \end{aligned}$$
(20)

into (19) to obtain

$$\begin{aligned} \Phi (z)\approx \sum _{m=-M}^{M}\theta _m(1-z+...)^{\frac{2\pi {\textbf{i}} m+\ln p_1}{\ln E}}. \end{aligned}$$
(21)

In comparison to the computations of \(\varphi _n\), the Taylor coefficients of RHS in (21) can be computed easily. For example, if we take the first approximation in (20) only, then (21) becomes

$$\begin{aligned} \Phi (z)\approx \sum _{m=-M}^{M}\theta _m(1-z)^{\frac{2\pi {\textbf{i}} m+\ln p_1}{\ln E}}=\sum _{n=0}^{+\infty }z^n\sum _{m=-M}^M(-1)^n\left( {\begin{array}{c}\frac{2\pi {\textbf{i}} m+\ln p_1}{\ln E}\\ n\end{array}}\right) \theta _m, \end{aligned}$$
(22)

which gives an approximation

$$\begin{aligned} \varphi _n\approx {{\tilde{\varphi }}}_n^M:=(-1)^n\sum _{m=-M}^M\left( {\begin{array}{c}\frac{2\pi {\textbf{i}} m+\ln p_1}{\ln E}\\ n\end{array}}\right) \theta _m, \end{aligned}$$
(23)

where the generalized binomial coefficients

$$\begin{aligned} \left( {\begin{array}{c}k\\ n\end{array}}\right) =\frac{k(k-1)(k-2)...(k-n+1)}{n!},\ \ \ n\geqslant 1 \end{aligned}$$
(24)

are Taylor coefficients of \((1+z)^k\) for arbitrary k. Using the asymptotic

$$\begin{aligned} \left( {\begin{array}{c}k\\ n\end{array}}\right) \approx \frac{(-1)^n}{\Gamma (-k)n^{k+1}},\ \ \ n\rightarrow +\infty , \end{aligned}$$
(25)

and the Stirling approximation

$$\begin{aligned} \Gamma (-k)\approx C_{m}e^{-\frac{\pi ^2|m|}{\ln E}}|m|^{-\frac{\ln p_1}{\ln E}-\frac{1}{2}} \end{aligned}$$
(26)

for some bounded constants \(C_{m}\), and \(k=\frac{2\pi {\textbf{i}} m+\ln p_1}{\ln E}\), see (23), with \(m\rightarrow \infty \), we obtain the approximation

$$\begin{aligned} (-1)^n\left( {\begin{array}{c}\frac{2\pi {\textbf{i}} m+\ln p_1}{\ln E}\\ n\end{array}}\right) \theta _m\approx C_m^{-1}\theta _me^{\frac{\pi ^2|m|}{\ln E}}|m|^{\frac{\ln p_1}{\ln E}+\frac{1}{2}}n^{-k-1}, \end{aligned}$$
(27)

for large n and m. While this approximation is very rough and we do not use it directly, it allows us to estimate the possible order of growths of terms in (23). Below, we discus some other details related to (27).

Remark on exponential decaying \(\theta _m\). If \(z=x+{\textbf{i}}y\) for \(x\rightarrow -\infty \) and \(y\in {{\mathbb {R}}}\), then \(\Pi (E^z)=1-E^xe^{{\textbf{i}}y\ln E}+o(E^x)\), since \(\Pi (0)=1\) and \(\Pi '(0)=-1\). Thus,

$$\begin{aligned} K^*(z)=\Phi (1-E^xe^{{\textbf{i}}y\ln E}+o(E^x))p_1^{-x-{\textbf{i}}y}. \end{aligned}$$
(28)

As is mentioned above, the function \(\Phi \) is defined in the open filled Julia set for the polynomial P which contains the unit ball. Hence, the maximal y for which \(K^*(z)\) is defined is at least \(\frac{\pi }{2\ln E}\) for \(x\rightarrow -\infty \), see (28), since the unit circle intersects real axis at the angle \(\frac{\pi }{2}\). Thus, the maximal width of the strip where the 1-periodic function \(K^*\) is defined satisfies \(2\vartheta ^*\geqslant \frac{\pi }{\ln E}\). Moreover, at the edges of the smaller strip with the width of \(\frac{\pi }{\ln E}\), the function \(K^*\) is bounded by some value \(C_P\), since the open filled Julia set for the polynomial P contains a segment of the line \(\{1+{\textbf{i}}y,\ y\in {{\mathbb {R}}}\}\) located near \(z=1\), except this unique point \(z=1\). For instance, this can be verified by applying the polynomial P iteratively to the points of this segment. For small y, each next iteration decreases the real part of the points

$$\begin{aligned} P(1+{\textbf{i}}y)=1+E{\textbf{i}}y-\frac{P''(1)}{2}y^2+O(y^3), \end{aligned}$$
(29)

and at some step the points are inside the unit ball. The modified Karlin-McGregor function

$$\begin{aligned} K^*\biggl (\frac{\ln z}{2\pi {\textbf{i}}}\biggr )=\sum _{m=-\infty }^{+\infty }\theta _mz^m \end{aligned}$$
(30)

is analytic in the ring \(\{z:\ e^{-2\pi \vartheta ^*}<|z|<e^{2\pi \vartheta ^*}\}\). Using \(\vartheta ^*\geqslant \frac{\pi }{2\ln E}\) and applying Cauchy estimates for this smaller ring we obtain

$$\begin{aligned} |\theta _{m}|\leqslant C_Pe^{-\frac{|m|\pi ^2}{\ln E}},\ \ \ m\in {{\mathbb {Z}}}, \end{aligned}$$
(31)

which shows rapid exponential decaying of \(\theta _m\). Unfortunately, the binomial factor in (23) may compensate this exponential decaying, see (27). However, the factor \(|m|^{\frac{\ln p_1}{\ln E}+\frac{1}{2}}\) in (27) and, especially, the exact value of \(\vartheta ^*\), which is often greater than \(\frac{\pi }{2\ln E}\), allow us to hope for a very rapid decrease of terms in (23). (Note that (23) converges anyway, see (38) below.) As an example, we draw the Julia set for the polynomial \(P(z)=0.1z+0.5z^2+0.4z^3\), see Fig. 1, where the so-called “critical angle” is strictly less than \(\pi \), which leads to \(\vartheta ^*>\frac{\pi }{2\ln E}\), see the explanations above (29). The “critical angle” is the minimal angle in the neighborhood of 1, the exterior of which asymptotically belongs to the filled Julia set when the radius tends to 0.

For practical implementations, it is useful to compute Fourier coefficients of \( K^{*} \left( {x - \frac{{{\textbf{i}}\pi }}{{2\ln E}}} \right) \) instead of \(K^*(x)\), because they contain already the exponential factor \(e^{\frac{m\pi ^2}{\ln E}}\). After that, \(\frac{m\pi ^2}{\ln E}\) should be subtracted from \( \ln \left( {\frac{{2\pi {\textbf{i}}m + \ln p_{1} }}{{\mathop {\ln E}\limits _{n} }}} \right) \) before taking the exponent. This procedure leads to stable and accurate computations of the terms in (23) for positive m. It is not necessary to compute the terms with negative m, since they are conjugations of “positive” ones. We actively use this trick for the computations in the examples Section.

Fig. 1
figure 1

Julia set for the polynomial \(P(z)=0.1z+0.5z^2+0.4z^3\). It is the boundary of the open filled Julia set, which is the domain of definition for the analytic function \(\Phi \), see (4) and (5). The unit ball belongs to the filled Julia set. The critical angle \({{\tilde{\vartheta }}}^*=2(\pi -\vartheta ^*\ln E)\) is computed approximately

Full size image

Remark on computation of \(\theta _m\) Approximations (19)-(23) are based on the computation of \(\theta _m\), \(m\in {{\mathbb {Z}}}\), or also \(e^{\frac{\pi ^2m}{\ln E}}\theta _m\), \(m\in {{\mathbb {N}}}\). Thus we should compute \(K^*(z)\), since it is periodic and

$$\begin{aligned} \theta _m=\int _0^1K^*(x)e^{-2\pi \textbf{i} m x}\ \ \textrm{or}\ \ e^{\frac{\pi ^2m}{\ln E}}\theta _m=\int _0^1K^*\left( x-\frac{{\textbf{i}}\pi }{2\ln E}\right) e^{-2\pi \textbf{i}mx}dx. \end{aligned}$$
(32)

In the end, the calculation of

$$\begin{aligned} K^*(z)=\Phi (\Pi (E^z)){p_1}^{-z} \end{aligned}$$
(33)

is reduced to the calculation of \(\Phi \) and \(\Pi \), see (15). Recall that

$$\begin{aligned} P(z)=p_1z+p_2z^2+...+p_Nz^N,\ \ \ 0<p_1<1, \end{aligned}$$
(34)

see below (2). By induction, let us define

$$\begin{aligned} \Phi _0(z)=z,\ \ \ \Phi _{t+1}(z)= & {} p_1^{-t-1}\underbrace{P\circ ...\circ P}_{t+1}(z) = p_1^{-t}p_1^{-1}P(\underbrace{P\circ ...\circ P}_{t}(z))\nonumber \\= & {} \Phi _t(z)+\sum _{j=2}^Np_jp_1^{(j-1)t-1}\Phi _t(z)^j, \end{aligned}$$
(35)

see (34). The recurrence formula (35) provides an exponentially fast convergence \(\Phi _t(z)\rightarrow \Phi (z)\), since \(p_1^{(j-1)t-1}\rightarrow 0\) for \(t\rightarrow +\infty \) exponentially fast when \(j\geqslant 2\).

To compute \(\Pi (z)\), we rewrite (34) as

$$\begin{aligned} P(z)=1-E(1-z)+q_2(1-z)^2+...+q_N(1-z)^N, \end{aligned}$$
(36)

see (8). By induction, let us define

$$\begin{aligned} \Pi _0(z)=z,\ \ \ \Pi _{t+1}(z)= & {} \underbrace{P\circ ...\circ P}_{t+1}(1-E^{-t-1}z) =\underbrace{P\circ ...\circ P}_{t}(P(1-E^{-t-1}z))\nonumber \\= & {} \underbrace{P\circ ...\circ P}_{t}(1-E^{-t}(z-\sum _{j=2}^Nq_jE^{-(j-1)t-j}z^j))\nonumber \\= & {} \Pi _t(z-\sum _{j=2}^Nq_jE^{-(j-1)t-j}z^j). \end{aligned}$$
(37)

The type of the recurrence formula (37) differs from that of (35), but (37) also provides the exponentially fast convergence \(\Pi _t(z)\rightarrow \Pi (z)\) for \(t\rightarrow +\infty \), since \(E^{-(j-1)t-j}\rightarrow 0\) for \(t\rightarrow +\infty \) exponentially fast when \(j\geqslant 2\). In the examples below, we use numerical procedures based on (35) and (37).

Remark on the approximations for \(M\rightarrow +\infty \). For fixed n, the series (23) and (42) converge when \(M\rightarrow +\infty \), since \(\theta _m\) decay exponentially fast, while the binomial coefficients are polynomials of m, see (24). Using (16) and (23), one can express the best among \(\{{{\tilde{\varphi }}}_n^M\}_{M\geqslant 0}\) approximation \({{\tilde{\varphi }}}_n^{\infty }\) through the derivatives of \(K^*\), namely

$$\begin{aligned} {{\tilde{\varphi }}}_n^{\infty }=(-1)^n\left( {\begin{array}{c}\frac{\frac{\partial }{\partial z}+\ln p_1}{\ln E}\\ n\end{array}}\right) K^*(z)|_{z=0}, \end{aligned}$$
(38)

since the derivative of \(K^*\) leads to the multiplication of each coefficient \(\theta _m\) by \(2\pi {\textbf{i}}m\).

Remark on the periodicity It is necessary to take into account the well known fact that, in some cases, we should be careful using the approximations. For example, if all \(p_{2n}=0\), then all \(\varphi _{2n}=0\) as well. In such cases, the approximation procedure (19)-(23) should be modified. The modification is not very difficult. One of the options is to take into account all the singularities of \(\Phi (z)\) on the unit circle, not only \(z=1\). These singularities coincide with the intersections of the Julia set with the unit circle. Suppose that \(p_kp_{k+1}\ne 0\) for some k. If \(|z|=1\) and \(z\ne 1\), then

$$\begin{aligned} |P(z)|\leqslant & {} |p_1z+...+p_{k-1}z^{k-1}|+|p_kz^k+p_{k+1}z^{k+1}|+|p_{k+2}z^{k+2}+...+p_{N}z^{N}|\\\leqslant & {} p_1+...+p_{k-1}+|p_k+zp_{k+1}|+p_{k+2}+...+p_N\\= & {} 1-p_k-p_{k+1}+|p_k+zp_{k+1}|\\< & {} 1, \end{aligned}$$

and, hence, there are no points of the Julia set belonging to the unit circle except one unique \(z=1\). Thus, this unique point is the only singularity of \(\Phi (z)\) in the closed unit ball. The same arguments lead to the inequality \(|P(z)|<|z|\) for \(|z|\leqslant 1\) and \(z\ne 1\), which means that \(z=0\) is the unique attractor for P in the unit ball.

We now move on to examples in Sect. 2, the complete asymptotic expansion will be considered in Sect. 3, see identities (58)–(61), the conclusion will be given in Sect. 4.

2 Examples

If \(\frac{\ln p_1}{\ln E}>-2\), then \(M=0\) or \(M=1\) in (23) already give a good approximation of \(\varphi _n\). The more interesting case is when \(p_1\) is small. Let us assume that

$$\begin{aligned} -3<\frac{\ln p_1}{\ln E}\leqslant -2. \end{aligned}$$
(39)

In this case, we take two terms in the Taylor expansion for \(\Psi (z)\), since \((1-z)^{\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}}\) and \((1-z)^{1+\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}}\) give power-series with non-attenuating coefficients, while the coefficients of \((1-z)^{2+\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}}\) already tend to 0. Thus, we have

$$\begin{aligned} \Psi (z)^{\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}}= & {} \biggl (1-z+\frac{P''(1)}{2(E^2-E)}(1-z)^2+O((1-z)^3)\biggr )^{\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}}\nonumber \\= & {} (1-z)^{\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}} \biggl (1+\frac{P''(1)}{2(E^2-E)}(1-z)+O((1-z)^2)\biggr )^{\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}}\nonumber \\= & {} (1-z)^{\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}}+ \frac{P''(1)(2\pi {\textbf{i}}m+\ln p_1)}{2(E^2-E)\ln E}(1-z)^{1+\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}}\nonumber \\{} & {} +O \left( {\left( {1 - z} \right) ^{{2 + \frac{{2\pi {\textbf{i}}m + \ln p_{1} }}{{\ln E}}}} } \right) , \end{aligned}$$
(40)

which leads to the approximation

$$\begin{aligned} \Phi (z)\approx \sum _{m=-M}^M\biggl ((1-z)^{\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}} +\frac{P''(1)(2\pi {\textbf{i}}m+\ln p_1)}{2(E^2-E)\ln E}(1-z)^{1+\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}}\biggr )\theta _m, \end{aligned}$$
(41)

see (19), which, in turn, leads to

$$\begin{aligned} \varphi _n\approx \varphi _n^M:= & {} (-1)^n\sum _{m=-M}^M\biggl (\left( {\begin{array}{c}\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}\\ n\end{array}}\right) \nonumber \\{} & {} +\frac{P''(1)(2\pi {\textbf{i}}m+\ln p_1)}{2(E^2-E)\ln E}\left( {\begin{array}{c}1+\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}\\ n\end{array}}\right) \biggr )\theta _m. \end{aligned}$$
(42)

It is useful to compare (23) with (42), where the new correction terms appear. These new terms are not leading, but they improve the approximation significantly.

It is possible to extract the conceptual behavior - the power factor and the oscillations - of \(\varphi _n^M\) from (42). It is enough to take two terms in the asymptotic of binomial coefficients

$$\begin{aligned} \left( {\begin{array}{c}k\\ n\end{array}}\right) \simeq \frac{(-1)^n}{\Gamma (-k)n^{k+1}}+\frac{(-1)^nk(k+1)}{2\Gamma (-k)n^{k+2}}+..., \end{aligned}$$
(43)

since the next terms already tend to zero for \(n\rightarrow \infty \) and \(\Re k>-3\). Using (43) in (42) along with \(\Gamma (k+1)=k\Gamma (k)\), we obtain

$$\begin{aligned} \varphi _n\approx \varphi _n^{\Gamma ,M}:= & {} n^{\frac{-\ln p_1}{\ln E}-1}\sum _{m=-M}^M\frac{\theta _me^{2\pi {\textbf{i}}m\frac{-\ln n}{\ln E}}}{\Gamma (-\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E})}\nonumber \\{} & {} +\frac{n^{\frac{-\ln p_1}{\ln E}-2}(P''(1)+E-E^2)}{2(E^2-E)\ln E}\sum _{m=-M}^M\frac{(2\pi {\textbf{i}}m+\ln p_1)\theta _me^{2\pi {\textbf{i}}m\frac{-\ln n}{\ln E}}}{\Gamma (-\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}-1)}.\qquad \end{aligned}$$
(44)

Indeed, (44) consists of the terms \(n^{\alpha }\textrm{Periodic}(\log _E n)\) that is very useful for understanding the conceptual behavior of \(\varphi _n\).

Let us discuss how to compute exact values \(\varphi _n\) for the cubic polynomial \(P(z)=p_1z+p_2z^2+p_3z^3\). Substituting (6) into (5), we obtain

$$\begin{aligned} \sum _{n=1}^{+\infty }\varphi _n(p_1z+p_2z^2+p_3z^3)^n=p_1\sum _{n=1}^{+\infty }\varphi _nz^n,\ \ \ \varphi _1=1, \end{aligned}$$
(45)

which can be rewritten in the form

$$\begin{aligned} \sum _{n=1}^{+\infty }\varphi _n\sum _{i+j+k=n}\frac{n!p_1^ip_2^jp_3^kz^{i+2j+3k}}{i!j!k!}=p_1\sum _{n=1}^{+\infty }\varphi _nz^n,\ \ \ \varphi _1=1, \end{aligned}$$
(46)

which, in turn, leads to

$$\begin{aligned} \sum _{n=1}^{+\infty }\varphi _n\sum _{r=n}^{3n}c_{n,r}z^r= & {} p_1\sum _{n=1}^{+\infty }\varphi _nz^n,\ \ \ \varphi _1=1,\ \ \ \nonumber \\ c_{n,r}:= & {} \sum _{k=0}^{\frac{r-n}{2}}\frac{n!p_1^{2n+k-r}p_2^{r-n-2k}p_3^k}{(2n+k-r)!(r-n-2k)!k!}, \end{aligned}$$
(47)

which, finally, gives the recurrence formula

$$\begin{aligned} \varphi _n=\frac{1}{p_1-p_1^n}\sum _{j=1}^{\frac{2n}{3}}c_{n-j,n}\varphi _{n-j},\ \ \ \varphi _1=1 \end{aligned}$$
(48)

that allows us to compute all \(\varphi _n\) recurrently. Yes, even for cubic polynomials, the computation of exact values \(\varphi _n\) is much harder than the computation of their approximations, as is seen from the comparison of (47) and (48) with (42) and (44).

2.1 Example \(p_1=0.2\), \(p_2=0.6\), \(p_3=0.2\)

We plot exact values (48) and their approximations (42) and (23) for the cubic polynomial \(P(z)=0.2z+0.6z^2+0.2z^3\), see Figs. 2 and 3. The simple approximations \(\varphi _n^0\) and \({{\tilde{\varphi }}}_n^0\) give already good results, see Fig. 2a. They show the dominant power-law growth \(\frac{\theta _0n^{-1-\ln p_1/\ln E}}{\Gamma (-\ln p_1/\ln E)}\), see (23) and (25). However, \(\varphi _n^1\) and, especially, \(\varphi _n^2\) are much better than \(\varphi _n^0\) and, of course, than \({{\tilde{\varphi }}}_n^0\). Namely, comparing Fig. 3.(a) and (b), it is seen that the additional terms in (42) play a very important role. The accurate \(\Gamma \)-approximations are even better than others, see Fig. 4. Let us provide some important characteristics for this example

$$\begin{aligned}{} & {} \frac{-\ln p_1}{\ln E}\approx 2.32,\ \ \ \theta _0\approx 1.94,\ \ e^{\frac{\pi ^2}{\ln E}}\theta _1\approx 0.16-1.3{\textbf{i}}, \\{} & {} e^{\frac{2\pi ^2}{\ln E}}\theta _2\approx -0.0036+0.04{\textbf{i}},\ \ e^{\frac{\pi ^2}{\ln E}}\approx 1.5\cdot 10^6. \end{aligned}$$
Fig. 2
figure 2

For the polynomial \(P(z)=0.2z+0.6z^2+0.2z^3\), Tailor coefficients \(\varphi _n\) (black curve) and their approximations \(\varphi _n^M\) (red, blue, and green curves) for \(M=0,1,2\) are plotted in (a), the differences \(\varphi _n^M-\varphi _n\) (red, blue, and green curves) are plotted in (b) (Color figure online)

Full size image
Fig. 3
figure 3

For the polynomial \(P(z)=0.2z+0.6z^2+0.2z^3\), the differences \({{\tilde{\varphi }}}_n^M-\varphi _n\) and \(\varphi _n^M-\varphi _n\) (blue and green curves) are plotted in (a) and (b) respectively (Color figure online)

Full size image
Fig. 4
figure 4

For the polynomial \(P(z)=0.2z+0.6z^2+0.2z^3\), the differences \({{\tilde{\varphi }}}_n^{\Gamma ,M}-\varphi _n\) and \(\varphi _n^{\Gamma ,M}-\varphi _n\) (blue and green curves) are plotted in (a) and (b) respectively. The \(\Gamma \)-approximations \({{\tilde{\varphi }}}_n^{\Gamma ,M}\) and \(\varphi _n^{\Gamma ,M}\) are given by (23) and (43) respectively, where, in the first case, the binomial coefficients are replaced by their approximations (24) (Color figure online)

Full size image

2.2 Example \(p_1=0.1\), \(p_2=0.5\), \(p_3=0.4\)

As another example, we compute \(\varphi _n\) and their approximations for the polynomial \(P(z)=0.1z+0.5z^2+0.4z^3\), see Figs. 5 and 6. Some important characteristics for this example are

$$\begin{aligned}{} & {} \frac{-\ln p_1}{\ln E}\approx 2.76,\ \ \ \theta _0\approx 2.887,\ \ e^{\frac{\pi ^2}{\ln E}}\theta _1\approx 9.94-0.62{\textbf{i}},\ \\{} & {} e^{\frac{2\pi ^2}{\ln E}}\theta _2\approx 1.21+1.66{\textbf{i}},\ \ e^{\frac{\pi ^2}{\ln E}}\approx 1.4\cdot 10^6. \end{aligned}$$
Fig. 5
figure 5

For the polynomial \(P(z)=0.1z+0.5z^2+0.4z^3\), Tailor coefficients \(\varphi _n\) (black curve) and their approximations \(\varphi _n^M\) (red, blue, and green curves) for \(M=0,1,2\) are plotted in (a), the differences \(\varphi _n^M-\varphi _n\) (red, blue, and green curves) are plotted in (b) (Color figure online)

Full size image
Fig. 6
figure 6

For the polynomial \(P(z)=0.1z+0.5z^2+0.4z^3\), the differences \({{\tilde{\varphi }}}_n^M-\varphi _n\) and \(\varphi _n^M-\varphi _n\) (blue and green curves) are plotted in (a) and (b) respectively (Color figure online)

Full size image

2.3 Example \(p_1=0.25\), \(p_2=0.5\), \(p_3=0.25\)

In fact, I made a lot of tests, not only for cubic polynomials. All of them confirm the main observations obtained in the previous examples. To finish this Section, I made one example, which complements some results from [17] (see Fig. 3 there). Namely, the approximations \(\varphi _n^1\) and \(\varphi _n^2\) improve \(\varphi _n^0\) significantly, see Fig. 7. This example is my main personal motivation for this work. Some important characteristics for this example are

$$\begin{aligned}{} & {} \frac{-\ln p_1}{\ln E}=2,\ \ \ \theta _0\approx 1.46,\ \ e^{\frac{\pi ^2}{\ln E}}\theta _1\approx 0.12+0.01{\textbf{i}}, \\{} & {} e^{\frac{2\pi ^2}{\ln E}}\theta _2\approx 0.002+0.0001{\textbf{i}},\ \ e^{\frac{\pi ^2}{\ln E}}\approx 1.5\cdot 10^6. \end{aligned}$$
Fig. 7
figure 7

For the polynomial \(P(z)=0.25z+0.5z^2+0.25z^3\), Tailor coefficients \(\varphi _n\) (black curve) and their approximations \(\varphi _n^M\) (red, blue, and green curves) for \(M=0,1,2\) are plotted in (a), the differences \(\varphi _n^M-\varphi _n\) (red, blue, and green curves) are plotted in (b) (Color figure online)

Full size image

3 General case

In this section, we consider some methods, allowing us to determine the complete asymptotic of \(\varphi _n\) for \(n\rightarrow +\infty \), up to terms that tend to 0. Following the ideas presented in the first part of this article, one may write

$$\begin{aligned} \Psi (z)^{\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}}=\sum _{j=0}^{+\infty }R_j(2\pi {\textbf{i}}m)(1-z)^{\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}+j}, \end{aligned}$$
(49)

where the polynomials \(R_j(x)\) have real coefficients determined from the functional equation for \(\Psi \), see (10) and (40),

$$\begin{aligned} R_0(x)=1,\ \ \ R_1(x)=\frac{P''(1)(x+\ln p_1)}{2(E^2-E)\ln E},...\ \end{aligned}$$
(50)

The next step is the use of expansions

$$\begin{aligned} (1-z)^{\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}+j}=\sum _{n=0}^{+\infty }(-1)^n\left( {\begin{array}{c}\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}+j\\ n\end{array}}\right) z^n \end{aligned}$$
(51)

in (49), which leads to

$$\begin{aligned} \Psi (z)^{\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}}=\sum _{n=0}^{+\infty }z^n\sum _{j=0}^{+\infty }R_j(2\pi {\textbf{i}}m)(-1)^n\left( {\begin{array}{c}\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}+j\\ n\end{array}}\right) . \end{aligned}$$
(52)

Then, using (52) along with (18) and (6), we deduce that

$$\begin{aligned} \varphi _n=\sum _{j=0}^{+\infty }\sum _{m=-\infty }^{+\infty }\theta _mR_j(2\pi {\textbf{i}}m)(-1)^n\left( {\begin{array}{c}\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}+j\\ n\end{array}}\right) . \end{aligned}$$
(53)

The next step is the approximation of the binomial coefficients

$$\begin{aligned} (-1)^n\left( {\begin{array}{c}k\\ n\end{array}}\right) \simeq \frac{n^{-k-1}}{\Gamma (-k)}\sum _{r=0}^{+\infty }\frac{S_{2r}(k)}{n^r}, \end{aligned}$$
(54)

with polynomials \(S_{2r}\) of the degree 2r. These polynomials can be determined explicitly from the identity

$$\begin{aligned} \biggl (1-\frac{k+1}{n+1}\biggr )(-1)^n\left( {\begin{array}{c}k\\ n\end{array}}\right) =(-1)^{n+1}\left( {\begin{array}{c}k\\ n+1\end{array}}\right) , \end{aligned}$$
(55)

which leads to the formal identities

$$\begin{aligned}{} & {} \biggl (1-\frac{k+1}{n}+\frac{k+1}{n^2}-...\biggr )\sum _{r=0}^{+\infty }\frac{S_{2r}(k)}{\Gamma (-k)n^{r+k+1}}\nonumber \\{} & {} \quad =\sum _{r=0}^{+\infty }\frac{S_{2r}(k)}{\Gamma (-k)(n+1)^{r+k+1}}\nonumber \\{} & {} \quad =\sum _{r=0}^{+\infty }\frac{S_{2r}(k)}{\Gamma (-k)n^{r+k+1}}\sum _{p=0}^{+\infty }\frac{1}{n^p}\left( {\begin{array}{c}-r-k-1\\ p\end{array}}\right) . \end{aligned}$$
(56)

Equating the terms with the same power of n and using \(S_0\equiv 1\), we can determine \(S_{2r}(k)\):

$$\begin{aligned} S_2(k)=\frac{k(k+1)}{2},\ \ \ S_{4}(k)=\frac{k(k+1)(k+2)(3k+1)}{24},...\ \ \ \end{aligned}$$
(57)

These polynomials resemble those from the ratio of two Gamma functions, see [18]. Using (54) in (53) leads to

$$\begin{aligned} \varphi _n\simeq \sum _{j=0}^{+\infty }\sum _{m=-\infty }^{+\infty }\theta _mR_j(2\pi {\textbf{i}}m) \sum _{r=0}^{+\infty }\frac{S_{2r}(\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}+j)}{\Gamma (-\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}-j)n^{r+\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}+j+1}}. \end{aligned}$$
(58)

To obtain the main part of the asymptotic for \(n\rightarrow +\infty \), we should take the terms with non-negative powers of n:

$$\begin{aligned}{} & {} \varphi _n\approx \varphi ^A_n=\sum _{0\leqslant r+j\leqslant \frac{-\ln p_1}{\ln E}-1}n^{\frac{-\ln p_1}{\ln E}-1-r-j}\nonumber \\{} & {} \sum _{m=-\infty }^{m=+\infty }\frac{S_{2r}(\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}+j)R_j(2\pi {\textbf{i}}m)\theta _me^{2\pi {\textbf{i}}m\frac{-\ln n}{\ln E}}}{\Gamma (-\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}-j)}. \end{aligned}$$
(59)

Introducing periodic functions

$$\begin{aligned} K_j(z)=\sum _{m=-\infty }^{+\infty }\frac{\theta _me^{2\pi {\textbf{i}}mz}}{\Gamma (-\frac{2\pi {\textbf{i}}m+\ln p_1}{\ln E}-j)}, \end{aligned}$$
(60)

we rewrite (59) as

$$\begin{aligned} \varphi _n\approx \varphi ^A_n=\sum _{0\leqslant r+j\leqslant \frac{-\ln p_1}{\ln E}-1}n^{\frac{-\ln p_1}{\ln E}-1-r-j}S_{2r}\biggl (\frac{\frac{\partial }{\partial x}+\ln p_1}{\ln E}+j\biggr )R_j\left( \frac{\partial }{\partial x}\right) K_j(x)|_{x=\frac{-\ln n}{\ln E}},\nonumber \\ \end{aligned}$$
(61)

since \(\partial /\partial x\) leads to the multiplication of the m-th Fourier coefficient by \(2\pi {\textbf{i}}m\). The asymptotic (61) contains explicitly the power-terms and periodic factors. As an interesting exercise, one may rewrite (60) as a convolution of two functions. In this exercise, it is useful to take into account the symmetry \(\theta _m=\overline{\theta _{-m}}\), and the exponential decaying of \(\theta _m\), when \(m\rightarrow \pm \infty \). As is already mentioned above, (61) is very useful in the conceptual context, since it consists of a finite number of explicit terms of the form \(n^{\alpha }\textrm{Periodic}(\log _E n)\).

Remark on the cases \(p_0\ne 0\) and \(p_0=p_1=0\). If \(p_0\ne 0\), then the role of \(\Phi (z)\), see (4), plays

$$\begin{aligned} \Phi (z):=\lim _{t\rightarrow \infty }P'(q)^{-t}(\underbrace{P\circ ...\circ P}_{t}(z)-q), \end{aligned}$$
(62)

where the extinction probability q is the minimal root of \(P(z)-z\). If \(E>1\), see (8), then the ideas presented in the main part of this work are applicable to the case \(p_0\ne 0\) as well. If \(E\leqslant 1\), then the situation is more complex since \(z=1\) is no longer a repelling point. However, if E is strictly less than 1, then there are still two marked points: attracting \(z=q=1\) and repelling \(z=\widetilde{q}>1\), where \(\widetilde{q}=P(\widetilde{q})\) is the second positive root of \(P(z)-z\). Thus, the presented ideas are still applicable if we use

$$\begin{aligned} \Psi (z):=\lim _{t\rightarrow \infty }P'(\widetilde{q})^{t}(\widetilde{q}-\underbrace{P^{-1}\circ ...\circ P^{-1}}_{t}(z)) \end{aligned}$$
(63)

instead of (9). Expansions (20), (49), etc. are assumed near \(z=\widetilde{q}\) instead of \(z=1\). Recall that we should always check the condition: is the repelling point a single singularity for \(\Phi \) in the corresponding ball? Another general note is that the choice of the extinction probability q as the center of the ball for the power series of \(\Phi \), see (62), leads to an analysis very similar to the one presented in the main part of the article. However, in this case, (7) is not fulfilled for \(\varphi _n\) given by

$$\begin{aligned} \Phi (z)=(z-q)+\varphi _2(z-q)^2+\varphi _3(z-q)^3+.... \end{aligned}$$

The case \(p_0=p_1=0\) is very different. There are no more direct analogs of (7), but we can try to analyze power-series expansions of another type of functions. Let \(L(\geqslant 2)\) be the minimal number for which \(p_L\ne 0\). Then we can define

$$\begin{aligned} \Phi (z):=(\underbrace{P\circ ...\circ P}_{t}(z))^{L^{-t}}, \end{aligned}$$
(64)

which satisfies the functional equation

$$\begin{aligned} \Phi (P(z))=\Phi (z)^L,\ \ \ \Phi (0)=0. \end{aligned}$$
(65)

The first serious issue is that the domain of definition for \(\Phi \) can be smaller than the filled Julia set for P, and, hence, \(z=1\) can be not a minimal singular point for \(\Phi (z)\). In principle, we can define an analog of the periodic function \(K^*\), namely,

$$\begin{aligned} K^*(z):=\Phi (\Pi (E^z))^{L^{-z}}, \end{aligned}$$
(66)

see (15). Then, taking another periodic function

$$\begin{aligned} K_1^*(z):=\ln K^*(z)=\sum _{m=-\infty }^{m=+\infty }\theta _{1m}e^{2\pi {\textbf{i}}mz}, \end{aligned}$$
(67)

we can rewrite (66) as

$$\begin{aligned} \ln \Phi (z)=\sum _{m=-\infty }^{m=+\infty }\theta _{1m}\Psi (z)^{\frac{2\pi {\textbf{i}}m+\ln L}{\ln E}}, \end{aligned}$$
(68)

which is similar to (18). However, the difference is that \(\ln \Phi (z)\) has a singular part \(\ln z\) when \(z\rightarrow 0\). Possibly more promising will be a direct analysis of another expansion, which follows from (64):

$$\begin{aligned} \ln \Phi (z)=\ln z+\frac{\ln p_L}{L-1}+\frac{\ln Q(z)}{L}+\frac{\ln Q(P(z))}{L^2}+\frac{\ln Q(P(P(z)))}{L^3}+..., \end{aligned}$$
(69)

where the polynomial Q satisfies \(P(z)=p_Lz^LQ(z)\) and \(Q(0)=1\). Note that

$$\begin{aligned} \underbrace{P\circ ...\circ P}_{t}(z)\sim p_L^{\frac{L^t-1}{L-1}}z^{L^t},\ \ \ z\rightarrow 0 \end{aligned}$$
(70)

and series (69) converge very fast. Let

$$\begin{aligned} \ln \Phi (z)=\ln z+\frac{\ln p_L}{L-1}+\sum _{n=1}^{+\infty }\varphi _{1n}z^n \end{aligned}$$
(71)

be a power-series expansion of LHS in (69). Denoting all the zeros of the corresponding polynomials in RHS of (69) as

$$\begin{aligned} \{z_{tj}\}_{j=1}^{N^t(N-L)}:=\underbrace{P^{-1}\circ ...\circ P^{-1}}_{t}\circ Q^{-1}(0), \end{aligned}$$
(72)

and using (70), we write an explicit expression for the power-series coefficients

$$\begin{aligned} \varphi _{1n}=\frac{-1}{n}\sum _{0\leqslant t\leqslant \frac{\ln n}{\ln L}}\frac{1}{L^{t+1}}\sum _{j=1}^{N^t(N-L)}z_{tj}^{-n}. \end{aligned}$$
(73)

Further consideration of (73) consists of different cases. Let us assume that the root of Q with the minimal norm is unique, say it is \(z_{01}\). We assume also that its norm \(|z_{01}|<1\). Taking into account the inequality \(|P(z)|<|z|\) for \(|z|<1\), see Remark at the end of Sect. 1, we can write \(|z_{01}|<r<|z_{tj}|\) for some \(r>0\) and all t, j except \(t=0\) and \(j=1\). Thus,

$$\begin{aligned} \biggl |\varphi _{1n}-\frac{-z_{01}^{-n}}{nL}\biggr |\leqslant N\frac{r^{-n}}{nL}\cdot n^{\frac{\ln N-\ln L}{\ln L}}, \end{aligned}$$
(74)

see (73), that means \(\varphi _{1n}\sim \frac{-z_{01}^{-n}}{nL}\). A similar analysis is possible when Q has several roots closest to 0 with the same norm less than 1. More asymptotic terms can also be obtained. The difficulty appears when this norm is equal to or greater than 1. The asymptotic analysis of harmonic measures related to the corresponding Julia set may be required in this case.

4 Conclusion

For the solutions of Schröder-type functional equations, we provide simple approximations of the corresponding power-series coefficients. A good agreement between the exact values of Taylor coefficients and their approximations is demonstrated by several numerical examples. The complete asymptotic of Taylor coefficients consisting of a finite number of terms \(n^{\alpha }\textrm{Periodic}(\log _E n)\) is also discussed. However, at the moment, I do not have rigorous estimates showing how good the asymptotic is. The analysis of the Bötcher case \(p_0=p_1=0\), presented in the current work, is incomplete.