Complete left tail asymptotic for the density of branching processes in the Schr¨oder case

For the density of Galton-Watson processes in the Schr¨oder case, we derive a complete left tail asymptotic series consisting of power terms multiplied by periodic factors.


Introduction
We consider a simple Galton-Watson branching process Z t in the supercritical case with the minimum family size 1 -the so-called Schröder case.The probability of the minimum family size is 0 < p 1 < 1.Note that the case of non-zero extinction probability can usually be reduced to the supercritical case with the help of Harris-Sevastyanov transformation, see [3].We assume that the mean of offspring distribution E < +∞.Then one may define the martingal limit W = lim t→+∞ E −t Z t .It is proven in [4] that p(x), the density of W , has an asymptotic p(x) = x α V (x) + o(x α ), x → +0, with explicit α = −1 − log E p 1 and a continuous, positive, multiplicatively periodic function, V , with period E. Further references to this asymptotic are always based on the principal work [4] and do not provide any formula for V (x), see, e.g., the corresponding remark in [6] or [5] devoted to the Schröder case.Recently, some explicit expressions for V (x) in terms of Fourier coefficients of 1-periodic Karlin-McGregor function and some values of Γfunction are given in [7].The derivation is based on the complete asymptotic series for discrete relative limit densities of the number of descendants provided in [8].However, even for the first asymptotic term, the derivation is somewhat informal because the continuous martingal limit and the discrete distribution of the relative limit densities differ significantly.
In particular, only the first term in both asymptotics has a common nature, all other terms are completely different.Now, I found beautiful and independent steps allowing me to obtain the complete expansion (not only the first term) with the certain value α defined above, β = − log E p 1 > 0, and explicit multiplicatively periodic functions V 1 , V 2 , V 3 , ..., see (14) and (15).Finally, we obtain a very efficient representation of V n in terms of Fourier coefficients of Karlin-McGregor functions and some values of Γ-function, see (24).The explicit form of V 1 (x) presented in [7] significantly helped to determine this series.However, we provide an independent proof -the main result is formulated in Theorem 3.1.Summarizing above, the asymptotic consists of periodic oscillations amplified by power-law multipliers.As mentioned in, e.g., [9,10,11], such type of behavior can be important in applications in physics and biology.

Main results
Let us recall some known facts about branching processes, including some recent results obtained in [8].The Galton-Watson process is defined by where all ξ j,t are independent and identically-distributed natural number-valued random variables with the probability-generating function For simplicity, we consider the case when P is entire.In this case, the first moment is automatically finite.A polynomial P is common in practice, but the results discussed below can be extended to a wide class of non-entire P .As the Introduction mentions, we assume p 0 = 0, p 1 ̸ = 0. Another natural assumption is p 1 < 1, otherwise the case p 1 = 1 is trivial.Under these assumptions, we can define which is analytic at least for |z| < 1.The function Φ satisfies the Scröder-type functional equation The function Φ has an inverse analytic in some neighborhood of z = 0.The coefficients κ j can be determined by differentiating the corresponding Poincaré-type functional equation at z = 0.In particular, It is shown in, e.g., [12] or [7], that the density p(x), see (1), can be computed by where γ is a modification of the original contour iR discussed below and This function is entire, satisfying another Poincaré-type functional equation Note that E > 1, since p 0 = 0, 0 < p 1 < 1, and p j = 1, see ( 5).The existence of the limit in (12) follows from the standard facts that the linearization of P at 1 leads to the multiplication of the increment of the argument by E. In other words, one can use (12) and see the convergence.The same reason applies to (6).The details related to the theory of Schröder and Poincaré-type functional equations are available in, e.g., [14].Everything is ready to derive an expansion of p(x) that, in turn, gives the complete left tail asymptotic series see (11), where Using ( 7), (13), and (15), it is easy to check that all V n are multiplicatively periodic The integration contour γ in (15) should be chosen so that Π(γ) lies in the domain, where Φ is invertible, and γ must be symmetric and connects −i∞ with +i∞.One of the choices is γ = iR + ε with sufficiently large ε > 0, since Π(z) → 0 for Re z ⩾ 0 and z → ∞.For such type of contours, the second integral in (16) along γ and along γ/E is equal to the same value, due to the Cauchy integral theorem.The exact power-law decay of Π(z), when Re z ⩾ 0 and z → ∞, can be determined from the identity where K is 1-periodic function, see (18).Following [12], integral (11) exists.We assume that it exists in a strong sense meaning Π(iy) → 0 for y → ±∞ -this also means Φ(Π(iy)) → 0 for y → ±∞.Hence, K(z), see ( 18) is analytic in a neighborhood of [x + iπ 2 ln E , x + 1 + iπ 2 ln E ] for some large x ∈ R. Due to periodicity, K(z) is analytic in some strip of a positive width, a neighborhood of the line R + iπ 2 ln E .Taking into account the fact that K(z) is also analytic in the strip | Im z| < π 2 ln E , see [8], [12] and [13], we deduce that K(z) is analytic and periodic in a larger strip | Im z| < k for some k > π 2 ln E .This means that the contour log E γ lies strictly inside the domain of the definition of K, see Fig. 1, and, hence, K(log E z) is bounded and smooth for z ∈ γ.Thus, Π(z) = O(z ln p 1 ln E ), see (17), and if ln p 1 ln E < −1, integrals (15) converge absolutely.Recall that p 1 < 1 and E > 1.Hence, ln p 1 ln E < −1 is not a rare case.In fact, it can be shown that ( 15) converges anyway because ln p 1 ln E < 0 and decaying function z

Representation of V n (x) through the Fourier coefficients of Karlin-McGregor function
Recall that the Karlin-McGregor function, see [1] and [2], is 1-periodic function given by where the corresponding Fourier coefficients θ m decay exponentially fast.The decay rate is determined by the width of the strip, where the Karlin-McGregor function is defined, see details in, e.g., [8].We recall only that the width is at least log E π.Let us rewrite (15) in terms of K: where θ * n m are Fourier coefficients of 1-periodic function They are convolution powers of original Fourier coefficients θ m .We already know the form of V 1 (x): see [7].Comparing ( 19) with ( 21) and taking into account κ 1 = 1, see (10), we guess the main identity A direct independent proof of ( 22) can be based on well-known Hankel integral representations of the Γ-function.It is easy to check that, after changes of variables, (22) becomes equivalent to the formula presented in Example 12.2.6 on page 254 of [15].This is a standard, but interesting, exercise in complex analysis and special functions.Substituting p n 1 instead of p 1 into (22), we obtain also Combining ( 19) with (23), we obtain the final result The formula (24) is more convenient for computations than (15) because it is similar to (21) for which efficient numerical schemes developed in [8] were already applied in [7].
To justify changing the order of summation and integration in (19), it is enough to take into account ln p 1 ln E < −1 and remember that K is smooth and bounded in the symmetric strip, parallel to R, of the width log E π, see details at the end of Section 2. Hence, K(z) is well approximated by its Fourier series for z ∈ log E γ, with the necessary rate of convergence.Gathering all these remarks, let us formulate the main result.Theorem 3.1.If P is entire, integral (11) exists in the strong sense Π(iy) → 0 for y → ±∞, and ln p 1 ln E < −1 then (2) along with (24) hold true.As it is mentioned at the end of Section 2, the assumption ln p 1 ln E < −1 can be usually omitted.The convergence of ( 24) is quite fast because θ * n m = O(e −2π|m|k ) for some k > π 2 ln E , which more than compensates small values of Γ in the denominator, see, e.g., the corresponding analysis based on Stirling approximations of Gamma function explained in [8].Thus, Fourier coefficients in (24) decay exponentially fast.The condition Π(iy) → 0 for y → ±∞ can be replaced with a more simple condition, which follows directly from (13) and from the definition of the Julia set.Proof.Indeed, recall that the filled Julia set related to P (z) contains the open unit disk because all the coefficients of the polynomial P (z) are positive and the maximal value at the boundary of the unit disk is P (1) = 1.Thus, if Π(iy) lies in the interior of the unit disk then Π(iy) lies in the interior of the corresponding component of the filled Julia set as well.Thus, P -iterations performed by the first formula in (13) converges to 0 by the standard properties of the Julia sets, since 0 is the unique attracting point inside the unit disc for P (z).
Using fast exponentially convergent algorithms for the computation of Π(z) developed in, e.g., [8], we can check the condition formulated in Proposition 3.2 numerically.On the other hand, if the interior of the corresponding component of the filled Julia set contains a little bit more than the open unit disk, namely a small sector near z = 1 of the angle ⩾ π, then the condition formulated in Proposition 3.2 is automatically satisfied, since the point w, defined by lies inside the filled Julia set and after a few, say m, P -iterations, see (13), Π(E m iy) = P • ... • P m (w) will be small enough by the properties of Julia sets.All the filled Julia sets we tested in various examples contain such sectors, see, e.g., the right panel in Fig. 2 computed with the help of this site 1 .

Examples
As examples, we consider the same cases as in [7], see Fig. 3.The difference is that we focus on the second asymptotic term.We use the same numerical procedures based on fast algorithms developed in [8] plus FFT to compute p(x).At the moment, FFT is a bottleneck because it is applied to the infinite integral, see the first formula in (11), with the use of some tricks.The tricks are necessary because FFT is developed for the finite integralsfrom −π to π.This leads to insufficient accuracy for small x.Now, I am working on the improvement of the situation.However, this fact shows how important the asymptotic is, because the computation of V 1 and V 2 is sufficiently accurate and does not involve infinite integrals.

ln p 1 Figure 1 :
Figure 1: Gray area is the domain of definition of K(z).Paths γ and log E γ are plotted with the blue color.