Charles Knessl
email: knessl@uic.edu

Wojciech Szpankowski

Enumeration of Binary Trees and Universal Types

Keywords: Binary trees, types, Lempel-Ziv'78, path length

Binary unlabeled ordered trees (further called binary trees) were studied at least since Euler, who enumerated them.

The number of such trees with n nodes is now known as the Catalan number. Over the years various interesting questions about the statistics of such trees were investigated (e.g., height and path length distributions for a randomly selected tree). Binary trees find an abundance of applications in computer science. However, recently Seroussi posed a new and interesting problem motivated by information theory considerations: how many binary trees of a given path length (sum of depths) are there? This question arose in the study of universal types of sequences. Two sequences of length p have the same universal type if they generate the same set of phrases in the incremental parsing of the Lempel-Ziv'78 scheme since one proves that such sequences converge to the same empirical distribution. It turns out that the number of distinct types of sequences of length p corresponds to the number of binary (unlabeled and ordered) trees, Tp, of given path length p (and also the number of distinct Lempel-Ziv'78 parsings of length p sequences). We first show that the number of binary trees with given path length p is asymptotically equal to Tp ∼ 2 2p/(log 2 p)(1+O(log -2/3 p)) . Then we establish various limiting distributions for the number of nodes (number of phrases in the Lempel-Ziv'78 scheme) when a tree is selected randomly among all trees of given path length p. Throughout, we use methods of analytic algorithmics such as generating functions and complex asymptotics, as well as methods of applied mathematics such as the WKB method and matched asymptotics.

Introduction

Trees are the most important nonlinear structures that arise in computer science. Applications are in abundance (cf. [START_REF] Knuth | The Art of Computer Programming[END_REF][START_REF] Knuth | Selected Papers on the Analysis of Algorithms[END_REF]); in this paper we discuss a novel application of binary unlabeled ordered trees (further called binary trees) in information theory (e.g., counting Lempel-Ziv'78 parsings and universal types). Tree structures have been the object of extensive mathematical investigations for many years, and many interesting facts have been discovered. Enumeration of binary trees, which are of principal importance to computer science, has been known already by Euler. Nowadays, the number of such trees built on n nodes is called the Catalan number.

Since Euler and Cayley, various interesting questions concerning statistics of randomly generated binary trees were investigated (cf. [START_REF] Flajolet | Analytical Combinatorics, in preparation; see also[END_REF][START_REF] Knuth | The Art of Computer Programming[END_REF][START_REF] Knuth | Selected Papers on the Analysis of Algorithms[END_REF][START_REF] Szpankowski | Average Case Analysis of Algorithms on Sequences[END_REF][START_REF] Takács | Conditional Limit Theorems for Branching Processes[END_REF][START_REF] Takács | The Asymptotic Distribution of the Total Heights of Random Rooted Trees[END_REF]). In the standard model, one selects uniformly a tree among all binary unlabeled ordered trees built on n nodes, T * n (where

|T * n | = 2n
While deep and interesting results concerning the behavior of binary trees in the standard model were uncovered, there are still many important unsolved problems of practical importance. Recently, Seroussi [START_REF] Seroussi | On Universal Types[END_REF], when studying universal types for sequences and distinct parsings of the Lempel-Ziv scheme, asked for the enumeration of binary trees with a given path length. Let T p be the set of binary trees of given path length p. Seroussi observed that the cardinality of T p corresponds to the number of possible parsings of sequences of length p in the Lempel-Ziv'78 scheme, and the number of universal types (that we discuss below). We shall first enumerate T p (cf. also Seroussi [START_REF] Seroussi | On the Number of t-ary Trees with a Given Pathlength[END_REF]), and then compute the limiting distribution of the number of nodes (phrases in the LZ'78 scheme) when a tree is selected uniformly among T p . To the best of our knowledge these problems were never addressed before, with the exception of [START_REF] Seroussi | On Universal Types[END_REF]. We show below that they are much harder than the corresponding problems in the more standard T * n model. As mentioned above, the problem of enumerating binary trees of a given path length arose in Seroussi's research on universal types. The method of types [START_REF] Csiszár | Information Theory: Coding Theorems for Discrete Memoryless Systems[END_REF] is a powerful technique in information theory, large deviations, and analysis of algorithms. It reduces calculations of the probability of rare events to a combinatorial analysis. Two sequences (over a finite alphabet) are of the same type if they have the same empirical distribution. For memoryless sources, the type is measured by the relative frequency of symbol occurrences, while for Markov sources one needs to count the number of pairs of symbols. It turns out (cf. [START_REF] Jacquet | Precise Worst Case Minimax Redundancy for Markov Sources[END_REF]) that the number of sequences of a given Markovian type can be counted by enumerating Eulerian paths in a multigraph. Recently, Seroussi [START_REF] Seroussi | On Universal Types[END_REF] introduced universal types (for individual sequences and/or for sequences generated by a stationary and ergodic source). Two sequences of the same length p are said to be of the same universal type if they generate the same set of phrases in the incremental parsing of the Lempel-Ziv'78 scheme. It is proved that such sequences have the same asymptotic empirical distribution. But, every set of phrases defines uniquely a binary tree of path length p [START_REF] Jacquet | Asymptotic Behavior of the Lempel-Ziv Parsing Scheme and Digital Search Trees[END_REF][START_REF] Seroussi | On Universal Types[END_REF] (with the number of phrases corresponding to the number of nodes in the T p model). For example, strings 10101100 and 01001011 have the same set of phrases {1, 0, 10, 11, 00} and therefore the corresponding binary trees are the same. Thus, enumeration of T p leads to counting universal types and different LZ'78 parsings of sequences of length p.

Let us now summarize our main results. It is easy to see that the generating function B(z, w) = n,p≥0 b(n, p)z n w p of the number b(n, p) of binary trees with n nodes and path length p satisfies the following functional equation [START_REF] Knuth | The Art of Computer Programming[END_REF] B(z, w) = 1 + zB 2 (zw, w).

(1.1)

Observe that this equation is asymmetric with respect to z and w. When enumerating trees in T * n , we set w = 1 to get the well known algebraic equation B(z, 1) = 1 + zB 2 (z, 1) that can be explicitly solved as B(z, 1) = 1 -√ 1 -4z /(2z) leading to the Catalan number. A randomly (uniformly) selected tree from T * n has path length L n that is asymptotically distributed as Airy's distribution [START_REF] Takács | A Bernoulli Excursion and its Various Applications[END_REF][START_REF] Takács | Conditional Limit Theorems for Branching Processes[END_REF], that is,

Pr{L n / √ 2n 3 ≤ x} → W (x)
where W (x) is the Airy distribution function defined by its moments [START_REF] Flajolet | Analytic Variations on the Airy Distribution[END_REF]. The Airy distribution arises in surprisingly many contexts, such as parking allocations, hashing tables, trees, discrete random walks, area under Brownian bridge, etc. [START_REF] Flajolet | Analytic Variations on the Airy Distribution[END_REF][START_REF] Flajolet | On the Analysis of Linear Probing Hashing[END_REF][START_REF] Knuth | Linear Probing and Graphs[END_REF][START_REF] Louchard | Kac's Formula, Lévy Local Time and Brownian Excursion[END_REF][START_REF] Louchard | The Brownian Excursion Area: A Numerical Analysis[END_REF][START_REF] Takács | A Bernoulli Excursion and its Various Applications[END_REF][START_REF] Takács | Conditional Limit Theorems for Branching Processes[END_REF][START_REF] Takács | The Asymptotic Distribution of the Total Heights of Random Rooted Trees[END_REF]. Setting z = 1 in (1.1) we arrive at

B(1, w) = 1 + B 2 (w, w)
which is not algebraically solvable. Observe that the coefficient of B(1, w) at w p of ∞ n=0 b(n, p) enumerates binary trees with given path length p. We denote the set of all such binary trees as T p and study in this paper its combinatorial and statistical properties. We shall show (and also give explicitly the error term, which involves the Airy function) that

T p := |T p | = 1 (log 2 p)
√ πp 2 2p log 2 p (1+c1 log -2/3 p+c2 log -1 p+O(log -4/3 p))

for large p, where c 1 and c 2 are explicitly computable constants. Seroussi first conjectured the form of the leading term in the exponent of the above asymptotic result, proved an upper bound of that form, and has recently obtained [START_REF] Seroussi | On the Number of t-ary Trees with a Given Pathlength[END_REF] a proof for the matching lower bound, using information-theoretic and combinatorial arguments for t-ary trees.

In this paper we further analyze the random variable N p representing the number of nodes in a randomly selected tree from the assembly T p . We show that (N p -E[N p ])/ Var [N p ] is asymptotically normal. Finally, after deriving various asymptotic expansions for B n (w) = p≥0 w p b(n, p), we analyze the number of trees b(n, p) with n nodes and path length p for various ranges of n and p. We also obtain b(n, p) in the asymptotic matching regions for various scales. In passing, we point out that T p = |T p | corresponds to the number of distinct universal types in Seroussi's sense and the number of distinct parsings of binary sequences of length p, while b(n, p) enumerates the number of Lempel-Ziv'78 parsings with n phrases. Observe that then N p represents the number of phrases in LZ'78 parsing of a sequence of length p in the T p model.

The functional equation (1.1) falls into the class of quicksort-like nonlinear functional equations (cf. [START_REF] Hennequin | Combinatorial Analysis of Quicksort Algorithm[END_REF][START_REF] Flajolet | On the Analysis of Linear Probing Hashing[END_REF][START_REF] Knuth | Selected Papers on the Analysis of Algorithms[END_REF][START_REF] Knessl | Quicksort algorithm again revisited[END_REF][START_REF] Régnier | A Limiting Distribution for Quicksort[END_REF][START_REF] Rösler | A Limit Theorem for Quicksort[END_REF]) that is still not fully analyzed (with some exceptions like the linear probing algorithm [START_REF] Flajolet | On the Analysis of Linear Probing Hashing[END_REF][START_REF] Knuth | Linear Probing and Graphs[END_REF]). Nonlinear functional equations of type (1.1) are not particularly suitable for analytic tools which work fine for linear functional equations (cf. [START_REF] Flajolet | Analytical Combinatorics, in preparation; see also[END_REF][START_REF] Szpankowski | Average Case Analysis of Algorithms on Sequences[END_REF]). Therefore, we turn to methods of applied mathematics such as matched asymptotics and the WKB method [START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers[END_REF]. These make certain assumptions about the forms of some asymptotic expansions and their asymptotic matching. When stating our main results (see Result 3 in section 2), we discuss the assumptions in more detail. The methods we use are analytic methods that are especially suitable for problems that cannot be solved exactly by transform methods (cf. [START_REF] Knessl | Quicksort algorithm again revisited[END_REF][START_REF] Knessl | Asymptotic Behavior of the Height in a Digital Search Tree and the Longest Phrase of the Lempel-Ziv Scheme[END_REF]).

The WKB method [START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers[END_REF][START_REF] Szpankowski | Average Case Analysis of Algorithms on Sequences[END_REF] was named after the physicists Wentzel, Kramers and Brillouin. It assumes that the solution, B(ξ; n), to a recurrence, functional equation or differential equation has the following asymptotic form B(ξ; n) ∼ e nφ(ξ) A(ξ) + 1 n A (1)

(ξ) + 1 n 2 A (2) (ξ) + • • • , n → ∞
where φ(ξ) and A(ξ), A (1) (ξ), . . . are unknown functions. These functions must be determined from the equation itself, often in conjunction with another tool known as the asymptotic matching principle (cf. [START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers[END_REF]). The outline for the paper is as follows: In Section 2 we present our main results and their consequences. In particular, we enumerate the number of binary trees of given path length and count the number of nodes in a randomly selected tree in the T p model. Derivations are given in Sections 3-8. Numerical studies are discussed in Section 9.

Summary of Results

We let b(n, p) denote the number of binary trees with n nodes and path length p. This function satisfies

We shall mostly analyze (2.1), and then obtain asymptotic results for b(n, p) by expanding the Cauchy integral (cf. [START_REF] Szpankowski | Average Case Analysis of Algorithms on Sequences[END_REF]) b(n, p) = 1 2πi C B n (w)w -p-1 dw.

(2.3)

Here C is any closed loop about the origin in the w-plane.

We can solve (2.2) when w = 1, noting that B(0, 1) = 1, to obtain

B(z, 1) ≡ a(z) = 1 2z 1 - √ 1 -4z and thus ∞ p=0 b(n, p) = 1 2πi C a(z) z n+1 dz = B n (1) = 1 n + 1 2n n (2.4)
is the Catalan number. This gives the total number of trees with n nodes, regardless of the total path length. By expanding (2. (2.7)

This gives the average total path length. Higher-order moments can be obtained in a similar manner. In particular, we obtain , n ≥ 0.

c n = B ′′ n (1) = ∞ p=2 p(p -
(2.8)

Asymptotically, for n → ∞, we obtain from (2.4), (2.9) and (2.10) via Stirling's formula

a n = 4 n √ πn 3/2 [1 + O(n -1 )] b n = 4 n 1 - 3 √ πn + O(n -1 ) c n = 4 n 10 3 √ π n 3/2 - 13 2 n + O( √ n) .
(2.9)

The constants b n and c n are related to the first two moments of the Airy distribution (cf. [START_REF] Flajolet | Analytic Variations on the Airy Distribution[END_REF] and [START_REF] Takács | A Bernoulli Excursion and its Various Applications[END_REF]). We can easily show that for each j B (j+1) n

(1)

B (j) n (1) = O(n 3/2 ), n → ∞. (2.10)
It is known [START_REF] Louchard | Kac's Formula, Lévy Local Time and Brownian Excursion[END_REF][START_REF] Louchard | The Brownian Excursion Area: A Numerical Analysis[END_REF][START_REF] Takács | A Bernoulli Excursion and its Various Applications[END_REF][START_REF] Takács | Conditional Limit Theorems for Branching Processes[END_REF] that the distribution of the total path length L n , that is,

Pr{L n = p} = b(n, p) ∞ p=0 b(n, p) (2.11)
follows an Airy distribution as n → ∞, and most of the mass occurs in the range p = O(n 3/2 ). In particular,

E[L n ] = √ πn 3/2 + O(n), Var [L n ] = 10 3 -π n 3 + O(n 5/2
E[A r ] = -Γ(-1/2) Γ((3r -1)/2) Ω r
where Ω r are determined by the recurrence

2Ω r = (3r -4)rΩ r-1 + r-1 j=1 r j Ω j Ω r-j , r ≥ 1
with Ω 0 = -1. Following Flajolet and Louchard [START_REF] Flajolet | Analytic Variations on the Airy Distribution[END_REF] we observe that

r≥0 Ω r w r r! = - Φ 2/3 (w) Φ -1/3 (w) , Φ ν = 2 F 0 1 2 + ν, 1 2 -ν; 3 2 w where 2 F 0 (a, b; z) = 1 + a • b 1! z + a(a + 1) • b(b + 1) 2! z 2 + • • •
is the generalized hypergeometric function (cf. [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]), and the above is a formal power series (that actually diverges).

A more difficult problem studied in this paper is to investigate the distribution of the number of nodes in trees of a fixed path length p, that is, for the ensemble T p . Let N p be the number of nodes for a tree uniformly generated from T p . It is a random variable distributed as

Pr{N p = n} = b(n, p) ∞ n=0 b(n, p)
.

(2.12)

We shall compute this distribution asymptotically, and also obtain the asymptotic structure of b(n, p) for various ranges of n and p. We note that the sums in (2.11) and (2.12) are actually finite, since b(n, p) is only non-zero in the range We now summarize our main results. Our derivations are quite complicated, and are delayed until Sec. 3-8. In passing, we should add that we use ideas of applied mathematics, such as linearization and asymptotic matching. We shall make certain assumptions about the forms of the asymptotic expansions, as well as the asymptotic matching between the various scales. In particular, we shall use the WKB method discussed in the introduction.

n J=2 ⌊log 2 J⌋ = p min (n) ≤ p ≤ p max (n) = n 2 . ( 2
We next formulate our main result concerning the cardinality of T p .

Result 1

The total number of trees of path length p is, for p → ∞

|T p | = ∞ n=0 b(n, p) = 1 (log 2 p) √ πp exp 2p log 2 log 2 p 1 - 3 2 A 0 log 2 a 1/3 Q -2/3 + M (Q)Q -1 + O(Q -4/3 ) (2.14) where Q = log p and M (Q) = (log Q)(1 + A 1 log 2) -log log 2 + (k 2 -A 1 log a) log 2, (2.15)
A 0 = 2 3 4 1/3 |r 0 | = 2.4743 . . . , A 1 = 1 log 2 - 1 3 = 1.1093 . . . , (2.16)
a = 2(log 2) 2 = .96090 . . . , r 0 = max{z : Ai(z) = 0} = -2.3381 . . . ,

Here k 2 ≈ 3.696 is obtained by numerically solving a nonlinear integral equation, and Ai(•) is the Airy function [START_REF] Andrews | Special Functions[END_REF], defined as a solution of the differential equation f ′′zf = 0 that decays as z → ∞.

It follows that the exponential growth rate of the total number of trees of path length p is

log n b(n, p) ∼ p log p 2(log 2) 2 (2.17)
with the correction terms involving the least negative root of the Airy function. This result was also recently obtained by Seroussi [START_REF] Seroussi | On the Number of t-ary Trees with a Given Pathlength[END_REF]. We will indicate how to formally obtain further terms in the asymptotic series in (2.14). Next, we discuss the random variable N p . Let

N (p) := E[N p ] = ∞ n=0 nb(n, p) ∞ n=0 b(n, p) , V(p) := Var [N p ] = ∞ n=0 (n -N (p)) 2 b(n, p) ∞ n=0 b(n, p)
.

Result 2 With the notation as above, we find that

N (p) = p Q log 2 1 - log 2 a 1/3 A 0 Q 2/3 + M (Q) -A 1 log 2 Q + O(Q -4/3 ) , (2.18)
V(p) = p Q 5/3 (log 2)A 0 6a 1/3 1 - 3A 1 A 0 a Q 1/3 + O(Q -2/3 ) , (2.19)
where A 0 , A 1 are defined in (2.16) and M (Q) is given by (2.15). Furthermore, the limiting distribution of N p is Gaussian, that is,

Pr{N p = n} = b(n, p) ∞ n=0 b(n, p) ∼ 1 2πV(p) exp - (n -N (p)) 2 2V(p) , (2.20
)

for p → ∞ and n -N (p) = O(V 1/2 (p)) = O( √ p(log p) -5/6
).

We note that while the most important scale for (2.11) is p = O(n 3/2 ), that for (2.12) is

p = n log 2 n + O[n(log n) 1/3 ],
which is close to the lower limit p min (n) (or upper limit n max (p)) of the support of b(n, p).

The above results are derived through the following main technical result. It gives detailed asymptotic results for the solution B n (w) to (2.1) as n → ∞, for various ranges of w.

This main result makes certain assumptions about the forms of various asymptotic expansions and the matching between them. Specifically, the result in item (a) below assumes that the function in (3.1) has the limiting behavior in (3.5). The result in (b) is based upon the WKB expansion in (4.9), while (2.20) assumes the asymptotic matching between the ranges in (a) and (b). To obtain the result in (c) we assumed that the function in (5.2) has an expansion of the form (5.101). The asymptotic matching assumption between ranges (b) and (c) determines a multiplicative constant in (2.21) The result in part (d) assumes the limit in (6.5), and (2.29) (respecively, (2.30)) assumes asymptotic matching between ranges (d) and (e) (respectively, (e)). The expansion in region (e) assumes the WKB form in (7.1).

Result 3 Consider binary trees with path length equal to p. Let B n (w) be its generating function satisfying (2.1). Then for n → ∞ we have the following asymptotic expansions.

(a) far right region: n → ∞, w ≫ 1

B n (w) ∼ w ( n 2 ) 2 n-1 B * (w) (2.21)
where B * (w) satisfies

B * (w) = 1 + 1 4w + 1 2w 2 + O(w -3 ), w → ∞; (2.22) B * (w) ∼ d 1 √ w -1 exp d 0 w -1 , w → 1 + , (2.23)
d 0 = log 2 0 ξ e ξ -1 dξ = .58224 . . . , d 1 = 4 √ 2π e d0/2 = 2.1350 . . . . (2.24)
The numerical calculation of B * (w) is discussed in Sections 3 and 9.

(b) right region:

w = 1 + β/n, 0 < β < ∞ B n (w) ∼ β n ĝ(β) exp[nΦ(β)], (2.25)
where

Φ(β) = log 2 + β 2 + 1 β log 2 -log(1-1 2 e -β ) ξ e ξ -1 dξ ≡ log 2 + β 2 + φ(β), ĝ(β) = 4 √ π e -β 2 /4 e -β/2 1 -e -β 2 -e -β 3/2 exp 1 2 βφ(β) + 1 2 β log 1 - 1 2 e -β .
(c) central region:

w = 1 + a/n 3/2 , -∞ < a < ∞ B n (w) = 1 n + 1 2n n + 4 n n 3/2 C(a) + 1 √ n C (1) (a) + O(n -1 ) , (2.26)
where

C(a) = (-a) D((-a) 2/3 ) = Y 3/2 D(Y ), Y = (-a) 2/3 , a < 0, D(Y ) = 1 2πi Br e sY 2 √ s + 4 2/3 Ai ′ (4 -1/3 s) Ai(4 -1/3 s) ds.
Here Br is a vertical contour on which ℜ(s) > 0, and √ s is analytic for ℜ(s) > 0 and positive for s real and positive. An alternate expression for the leading term is, for a = -Y 3/2 < 0,

B n (w) ∼ 4 n n 3/2 (-a) d dY 1 2πi Br 4 2/3 s Ai ′ (4 -1/3 s) Ai(4 -1/3 s) e sY ds = 4 n+1 n 3/2 (-a) ∞ j=0 exp(-|r j |4 1/3 Y ) (2.27)
where 0 > r 0 > r 1 > r 2 > . . . and r j are the roots of Ai(z) = 0. The correction term has the integral representation, for a < 0,

C (1) (a) = -a D1 (Y ) = Y 2
2πi Br e sY E * (s)ds, (2.28)

E * (s) = - 5 2 s + 8 h ′ (s) h(s) 2 - 4 h 2 (s) ∞ s (h ′ (v)) 3 h(v) dv = - 5 2 s + 10 h ′ (s) h(s) 2 + 4 h ′ (s) h(s) 2 log[h(s)] -s log[h(s)] - 1 h 2 (s) ∞ s h 2 (v) log[h(v)]dv, h(s) = Ai(4 -1/3 s).
For a > 0 we let a = y 3/2 with y > 0 and the leading term is

B n (w) ∼ 4 n n 3/2 a π 2 4 1/3 ∞ 0 e τ y h(ωτ )h(ω 2 τ ) dτ - 4a π ∞ 0 ℜ e πi/6 h ′ (ωτ ) h(ωτ ) e ω 2 τ y dτ (2.29)
where ω = exp(2πi/3).

(d) left region:

w = 1 -γ/n, 0 < γ < ∞ B n (w) ∼ 4 n n exp[ν 0 n 1/3 γ 2/3 + ν 1 γ log n]F 0 (γ), F 0 (γ) = 4γF 1 (γ), ν 0 = 4 1/3 r 0 = -4 1/3 |r 0 |, ν 1 = - 1 3 , (2.30)
where F 1 (•) satisfies the non-linear integral equation

e γ -1 γ F 1 (γ) = 1 0 F 1 (γx)F 1 (γ -γx)e -γH(x)/3 dx, H(x) = x log x + (1 -x) log(1 -x) (2.31)
and behaves, for γ → 0 + , as

F 1 (γ) = 1 - 2 3 γ log γ + α 1 γ + O(γ)
where

α 1 = 7 2 -γ E + log[h ′ (s 0 )] - 1 4[h ′ (s 0 )] 2 ∞ s0 h 2 (v) log[h(v)]dv = 2.9622 . . . , s 0 = 4 1/3 r 0 , γ E = Euler ′ s constant = .57721 . . . . (2.32) For γ → ∞ we have F 1 (γ) ∼ 1 √ 2π log 2 e k2γ √ γ exp 1 3 - 1 log 2 γ log γ (2.33)
where k 2 ≈ 3.696 is found numerically (cf. Sections 6 and 9).

(e) far left region: n → ∞, 0 < w ≪ 1 B n (w) ∼ e n(log 2 n) log w e n[g(w)+B * 0 (w,n)] n log 2 w (2πn) -1/2 × e g(w) w 2+ 1 log 2 e B * 0 (w,n)+B * 1 (w,n) -log 2 w -B * 1 (w, n) - 1 4 B * 2 (w, n).
(2.34)

Here B * 0 (w, n) = ∞ k=-∞ k =0 g k (w)e 2πi(log 2 n)k B * 1 (w, n) = 2πi log 2 ∞ k=-∞ kg k (w)e 2πi(log 2 n)k B * 2 (w, n) = 2πi log 2 ∞ k=-∞ 2πi log 2 k 2 -k g k (w)e 2πi(log 2 n)k (2.35)
and g(w) has the asymptotic expansion

g(w) = log 4 + 4 1/3 r 0 (1 -w) 2/3 + 1 log 2 - 1 3 (w -1) log(1 -w) -k 2 (w -1) + o(w -1), w → 1 -.
(2.36)

The numerical calculation of g(w) is discussed in Sections 7 and 9. The sum in B * 0 omits the term k = 0 and the non-constant Fourier coefficients g k (w) satisfy g k (w) = o(w -1) as w → 1 -. Numerical studies show that unless w is very small, the g k (•) are very small and we can use the approximation

B n (w) ≈ w n log 2 n e (n+1)g(w) n log 2 w n -1/2 w 2+ 1 log 2 -log 2 w 2π . (2.37)
This corresponds to neglecting B * j in (2.34) for j = 0, 1, 2. We comment that our analysis suggests that yet another scale exists, which has n → ∞ and w → 0 simultaneously, and where a different expansion for B n (w) is needed. We have not been able to analyze this scale, but it is not needed to obtain the asymptotic results for the number of trees of a given total path length. For this the important range is the asymptotic matching region between the left and far left regions, corresponding to w → 1 -, but n(1w) = γ → +∞. Since we have explicit analytic results for g(w) as w → 1 -, and g k (w) → 0 for k = 0, we can use the above results to obtain the explicit expressions in (2.14) - (2.20). To obtain the distribution of the path length in trees with n (→ ∞) nodes, the central region (c) is the most important, and the leading term corresponds to (the transform of) the Airy distribution.

We next give results for b(n, p) for n and p → ∞, and summarize the main results as items (A)-(E) below. Going from (A) to (E) corresponds to increasing n or decreasing p.

Result 4 Consider binary trees built over n nodes with path length p. Let b(n, p) denote the number of such trees. Then we have the following for p, n → ∞.

(A) n → ∞, p = n 2 -L, L = O(1), L ≥ 0 b(n, p) ∼ 2 n-1 1 2πi C w L-1 B * (w)dw (2.38)
where C is a closed loop with |w| > 1 and B * (w) is as in (a) in Result 3.

(B) p, n → ∞ with Λ = p/n 2 ∈ 0, 1 2 b(n, p) ∼ √ 2 π β * e -β * /2 1 -e -β * 2 -e -β * 3/2 1 -2Λ - 1 2e β * -1 -1/2 × 2 n+1 n 2 exp n β * (1 -2Λ) -log 1 - 1 2 e -β *
(2.39)

where β * ≡ β * (Λ) is defined implicitly by

β 2 * 1 2 -Λ + β * log 1 - 1 2 e -β * = log 2 -log(1-1 2 e -β * ) ξ e ξ -1 dξ. (C) p, n → ∞ with Ω = p/n 3/2 ∈ (0, ∞) b(n, p) ∼ - 4 n n 3 1 3Ω 1/3 ∞ j=0 56 9
4 2/3 r 2 j Ω 3 + 64 81 4 5/3 r 5 j Ω 5 Ai r 2 j 4 2/3 3 4/3 Ω 4/3 + 1 3Ω 1/3 40 3 4 1/3 r j Ω 2 + 64 27
4 4/3 r 4 j Ω 4 Ai ′ r 2 j 4 2/3 3 4/3 Ω 4/3 exp - 8|r j | 3 27Ω 2 .
(2.40)

Here r j < 0 are the roots of the Airy function.

(D) p, n → ∞ with Θ = p/n 4/3 ∈ (0, ∞) b(n, p) ∼ 4 n n 13/6 n -γ * /3 2 9 3 4 √ π |r 0 | 9/2 Θ 5 F 1 (γ * ) exp - 16n 1/3 |r 0 | 3 27Θ 2 , γ * = 32 27 |r 0 | 3 Θ 3 (2.41)
where

F 1 (•) satisfies (2.31) (cf. item (d) in Result 3). (E) p, n → ∞ with p = n log 2 n + αn, α = O(1) b(n, p) ≈ n log 2 w * 2πn w 2+ 1 log 2 * α + w 2 * g ′′ (w * ) e g(w * ) -log 2 w * × exp [ng(w * ) -nα log w * ] (2.42)
where w * = w * (α) is the solution to w * g ′ (w * ) = α.

From item (C) the Airy distribution can be recovered by dividing by the expansion of a n in (2.6). Then (2.37) becomes of the form n -3/2 times the Airy density.

In obtaining (2.42) we used (2.34) in (2.3) and neglected the non-constant terms in the Fourier series, which are numerically small. Then we evaluated the integral by the saddle point method (cf. [START_REF] Szpankowski | Average Case Analysis of Algorithms on Sequences[END_REF]). A refined approximation that uses also the non-constant terms in the Fourier series' in (2.35), can be obtained by using (2.34) in (2.3). We can also obtain an O(n -1/2 ) correction term to the Airy distribution in (2.40) by using the correction term (i.e., C (1) (a)) in (2.26) in asymptotically inverting (2.3). Our approximation(s) to b(n, p) involve the unknown functions B * (w), F 1 (γ) and g(w), whose numerical calculation we discuss in Section 9.

In view of the complexity of the results in items (A)-(E), we can get more insight into the structure of b(n, p) by giving formulas that apply in the asymptotic matching regions between the various scales. We summarize these below, with the notation (AB) denoting the asymptotic matching region between the scales in (A) and (B) , and so on. Note that the (AB) result can be obtained by either expanding (2.38) as L → ∞ (using (2.23)), or by expanding (2.39) as Λ → 1 2 -.

Result 5

The following matching asymptotics hold:

(AB) n, p → ∞; L = n 2 -p → ∞, Λ = p/n 2 → 1 2 b(n, p) ∼ 2 n πn 2 2d 0 1 -2Λ 1/2 1 √ 1 -2Λ exp n 2d 0 (1 -2Λ) - 1 2 2d 0 1 -2Λ (2.43)
where d 0 is given in (2.24).

(BC) n, p → ∞; Λ = p/n 2 → 0, Ω = p/n 3/2 → ∞ b(n, p) ∼ 4 n n 2 9 √ 3 2π Λ 2 exp - 3 4 nΛ 2 = 4 n n 3 9 √ 3 2π Ω 2 exp - 3 4 Ω 2 .
(2.44)

(CD) n, p → ∞; Ω = p/n 3/2 → 0, Θ = p/n 4/3 → ∞ b(n, p) ∼ 4 n n 13/6 |r 0 | 9/2 Θ 5 2 9 3 4 √ π exp - 16n 1/3 |r 0 | 3 27Θ 2 = 4 n n 3 |r 0 | 9/2 Ω 5 2 9 3 4 √ π exp - 16|r 0 | 5 27Ω 2 .
(2.45)

(DE) n, p → ∞; Θ = p/n 4/3 → 0, α = p/n -log 2 n → ∞ b(n, p) ∼ 1 n 13/6 |r 0 | 3 Θ 7/2 1 π √ log 2 64 9 √ 3 exp n log 4 - γ * 3 log n + 1 3 - 1 log 2 γ * log γ * + k 2 γ * - 16 27 n 1/3 |r 0 | 3 Θ 2 (2.46)
where γ * is given below (2.41).

We will show in Section 8 that the asymptotic matching region (DE) leads to the Gaussian distribution in (2.20). Note that in each of the four matching regions, our results are completely explicit functions of n and p. The result in (BC) (resp., (CD)) gives the right (resp., left) tail of the Airy distribution in (2.40). The right tail has not been characterized this precisely in previous studies [START_REF] Csőrgő | Some Asymptotic Properties of the Local Time of the Uniform Empirical Processes[END_REF] (cf. also [START_REF] Flajolet | Analytic Variations on the Airy Distribution[END_REF][START_REF] Louchard | Kac's Formula, Lévy Local Time and Brownian Excursion[END_REF][START_REF] Louchard | The Brownian Excursion Area: A Numerical Analysis[END_REF][START_REF] Takács | A Bernoulli Excursion and its Various Applications[END_REF][START_REF] Takács | Conditional Limit Theorems for Branching Processes[END_REF][START_REF] Takács | The Asymptotic Distribution of the Total Heights of Random Rooted Trees[END_REF]).

Far Right Region

We consider (2.1) for a fixed w > 1 and n → ∞. We set

B n (w) = w ( n 2 ) 2 n-1 Bn (w) (3.1)
to find that (2.1) becomes

Bn+1 (w) = 1 4 n ℓ=0 w -ℓ(n-ℓ) Bℓ (w) Bn-ℓ (w) = 1 4 [2 B0 (w) Bn (w) + 2w 1-n B1 (w) Bn-1 (w)
Bn+1 (w) = Bn (w) = 1 2 w 1-n Bn-1 (w) + 1 4 w 4-2n Bn-2 (w) + O(w -3n ) (3.4)
whose asymptotic solution is

Bn (w) = B * (w) + O(w -n ) (3.5)
for some function B * (•). From (2.1) we obtain the first few B n (w) as

B 1 (w) = 1, B 2 (w) = 2w, B 3 (w) = w 2 + 4w 3 B 4 (w) = 4w 4 + 2w 5 + 8w 6
B 5 (w) = 6w 6 + 8w 7 + 8w 8 + 4w 9 + 16w 10 (3.6) and it is easy to inductively establish that for w → ∞ and fixed n

B n (w) = 2 n-1 w ( n 2 ) + 2 n-3 w ( n 2 )-1 + 2 n-2 w ( n 2 )-2 + O w ( n 2 )-3 , (3.7)
for n ≥ 4. By comparing (3.7) to (3.1) with (3.5) we conclude that

B * (w) = 1 + 1 4w + 1 2w 2 + O(w -3 ), w → ∞. (3.8)
In section 4 we will argue, by asymptotic matching, that as w ↓ 1 we have

B * (w) ∼ d 1 √ w -1 exp d 0 w -1 , w → 1 + (3.9)
where the constants d 0 and d 1 are given in (2.24). We have not been able to determine B * (w) analytically except for its behaviors as w → ∞ and w ↓ 1. It can easily be obtained numerically by fixing w > 1 and iterating (3.2) until Bn (w) settles to a limit to some prescribed accuracy.

In Table 1, we give B * (w) for w in the range [1.04, 5]. We note that the convergence of this procedure becomes very slow (and B * (w) becomes very large) when w exceeds 1 only slightly. This is certainly consistent with the asymptotic analysis, that predicts another scale where n → ∞ with n(w -1) fixed.

Our analysis thus far assumed that w is real. However, the same arguments show that (3.1) remains valid for complex w with |w| > 1. We can also use (3.2) to calculate B * (w) for complex w with |w| > 1.

In Figure 1

(n) = n 2 . Setting L = n 2 -p = O(1)
we obtain from (2.3)

b(n, p) ∼ 2 n-1 2πi C w L-1 B * (w)dw, n → ∞, (3.10)
Tab. 1: where C is a closed contour on which |w| > 1. We can use (3.8) to compute (3.10) for L = 0, 1, 2, . . . , and in Section 4 we will obtain the asymptotic behavior of the integral as L → ∞.

w B * (w) 5

Right Region

We consider the limit w ↓ 1 and n → ∞ with

n(w -1) = β = O(1), β > 0. (4.1)
We let

B n (w) = f (β; n) = f (n(w -1); n) (4.2)
with which (2.1) becomes

f β + β n ; n + 1 = 2 1 + β n n n/2 ℓ=0 B ℓ 1 + β n f β 1 - ℓ n ; n -ℓ . (4.3)
We used the symmetry of the sum in (2.1) to truncate the upper limit on the sum in (4.3) at ℓ = n/2. To be more precise, we should distinguish the cases where n is odd or even, but it will become apparent that for n → ∞ the asymptotically dominant terms in (4.3) are those with ℓ = O(1), and the terms near ℓ = n/2 are exponentially smaller than the dominant ones. We also note that in (4.3) we re-wrote only the second factor in the convolution sum in (2.1) using the scaling in (4.1) and (4.2). The behavior of B ℓ (w) for fixed ℓ and w → 1 follows from (2.5), and we have

B ℓ 1 + β n = a ℓ + b ℓ β n + O(n -2 ) (4.4)
where

a ℓ = 1 ℓ + 1 2ℓ ℓ , b ℓ = 4 ℓ - 3ℓ + 1 ℓ + 1 2ℓ ℓ (4.5)
and these have the generating functions

∞ ℓ=0 a ℓ z ℓ = 1 2z 1 - √ 1 -4z (4.6) ∞ ℓ=0 b ℓ z ℓ = 1 z 1 - √ 1 -4z + 1 1 -4z - 3 √ 1 -4z . ( 4.7)
We also have

1 + β n n = e β 1 - β 2 2n + O(n -2 ) , n → ∞. (4.8)
We assume, for fixed β > 0 and n → ∞, that f (β; n) has an expansion of the WKB form

f (β; n) = e nΦ(β) n -1/2 g(β) + 1 n g (1) (β) + O(n -2 ) . (4.9)
The factor n -1/2 must be included in order to asymptotically match to another expansion that applies for w -1 = O(n -3/2 ). This is the "central region" analyzed in section 5. In section 9 we provide numerical justification for the ansatz (4.9). The numerical studies suggest, however, that the ratio of the correction term to the leading term in (4.9) is larger than O(n -1 ), perhaps O(n -1/2 ). Thus the series in (4.9) may actually be in powers of n -1/2 rather than n -1 . However, this will not affect our calculation of Φ(β) and g(β).

We use (4.4) and (4.9) in (4.3) and extend the limit on the sum from n/2 to ∞. We note that

exp (n + 1)Φ β + β n (n + 1) -1/2 g β + β n + 1 n + 1 g (1) β + β n + O(n -2 ) = e nΦ(β) e βΦ ′ (β)+Φ(β) n -1/2 1 + 1 n βΦ ′ (β) + 1 2 β 2 Φ ′′ (β) 1 - 1 2n × g(β) + 1 n (g (1) (β) + βg ′ (β)) + O(n -2 ) (4.10)
and

n/2 ℓ=0 B ℓ 1 + β n f β 1 - ℓ n ; n -ℓ = n/2 ℓ=0 B ℓ (1) + β n B ′ ℓ (1) + O(n -2 ) 1 √ n -ℓ exp (n -ℓ)Φ β - βℓ n × g(β) + 1 n (g (1) (β) -ℓβg ′ (β)) + O(n -2 ) = e nΦ(β) n -1/2 ∞ ℓ=0 e -βℓΦ ′ (β) a ℓ + β n b ℓ 1 + ℓ 2n e -ℓΦ(β) × 1 + β 2 ℓ 2 2n Φ ′′ (β) 1 + βℓ 2 n Φ ′ (β) g(β) + 1 n (g (1) (β) -ℓβg ′ (β)) + O(n -2 ) . (4.11)
Dividing (4.3) by 2(1 + β/n) n and letting n → ∞ we obtain, after cancelling the common factor n -1/2 exp[nΦ(β)], the limiting equation

1 2 e -β e βΦ ′ +Φ = ∞ ℓ=0 a ℓ e -ℓ(βΦ ′ +Φ) . (4.12)
The above is a non-linear ODE for the function Φ(β), which is in view of (4.9) the exponential growth rate of B n (w) on the β-scale. Using (4.6) to evaluate the sum in (4.12) leads to, after some simplification,

1 2 e -β = 1 2 1 -1 -4e -(βΦ) ′ (4.13) or (βΦ) ′ = log 4 + 2β -log(2e β -1). (4.14)
Equation (4.14) integrates to

Φ(β) = log 4 + β - 1 β β 0 log(2e v -1)dv = log 2 + β 2 - 1 β β 0 log 1 - 1 2 e -v dv = log 2 + β 2 + φ(β) (4.15)
where

φ(β) = - 1 β β 0 log 1 - 1 2 e -v dv = 1 β log 2 -log(1-1 2 e -β ) ζ e ζ -1 dζ. (4.16)
We note that Φ(β) = c/β is a homogeneous solution to (4.14), but such solutions must be discarded since asymptotic matching to the central region (w -1 = O(n -3/2 )) will force Φ(β) to be bounded as

β → 0 + (in fact Φ(0) = log 4).
We next determine g(β) in (4.9). Using (4.8), (4.10), and (4.11) in (4.3), we obtain, at order e nΦ n -1/2 • n -1 , the linear equation

1 2 e -β e (βΦ) ′ βg ′ + β 2 2 Φ ′′ + βΦ ′ + β 2 2 - 1 2 g = -βg ′ ∞ ℓ=0 ℓa ℓ e -ℓ(βΦ) ′ + βg ∞ ℓ=0 b ℓ e -ℓ(βΦ) ′ +g ∞ ℓ=0 a ℓ ℓ 2 + β 2 ℓ 2 2 Φ ′′ + ℓ 2 βΦ ′ e -ℓ(βΦ) ′ . (4.17)
Note that in view of (4.12) g (1) drops out of (4.17). The latter is a simple first order linear ODE since we have determined Φ.

To integrate (4.17) we define the sums

S j (β) = ∞ ℓ=0 ℓ j a ℓ e -ℓ(βΦ) ′ ; j = 0, 1, 2 (4.18)
T 0 (β) = ∞ ℓ=0 b ℓ e -ℓ(βΦ) ′ . (4.19)
Setting

z = e -(βΦ) ′ = 1 4 e -β (2 -e -β )
we have

∞ ℓ=0 z ℓ a ℓ = S 0 (β) = 2 2 -e -β .
(4.20)

By differentiating (4.6) we have

∞ ℓ=0 ℓz ℓ a ℓ = S 1 (β) = z √ 1 -4z 1 1 + √ 1 -4z 2 = 1 4 1 1 -e -β - 2 2 -e -β . (4.21)
Also, by differentiating S 0 with respect to β we obtain

S ′ 0 (β) = - ∞ ℓ=0 ℓ(βΦ) ′′ a ℓ e -ℓ(βΦ) ′ = -(βΦ) ′′ S 1 (β) (4.22)
and then

S ′′ 0 (β) = [(βΦ) ′′ ] 2 S 2 (β) -(βΦ) ′′′ S 1 (β). (4.23)
Using (4.20) -(4.23) in (4.17) and noting that

β 2 2 Φ ′′ + βΦ ′ = β 2 (βΦ) ′′ leads to 1 2 e -β e (βΦ) ′ βg ′ + β 2 (βΦ) ′′ + β 2 2 - 1 2 g = βgT 0 + S ′ 0 (βΦ) ′′ βg ′ - S ′ 0 2(βΦ) ′′ g + β 2 1 (βΦ) ′′ [S ′′ 0 + (βΦ) ′′′ S 1 ]g. (4.24)
But, (4.12) yields 1 2 e -β e (βΦ) ′ = S 0 and then we set

g(β) = βĝ(β) (4.25)
with which (4.24) becomes

S 0 - S ′ 0 (βΦ) ′′ ĝ′ (β) ĝ(β) = T 0 - S 0 2 [β + (βΦ) ′′ ] + S ′′ 0 + S 1 (βΦ) ′′′ 2(βΦ) ′′ . ( 4.26)
From (4.14) and (4.20) we obtain

S ′ 0 = - 2e -β (2 -e -β ) 2 = -(βΦ) ′′′ (4.27)
and

S 0 - S ′ 0 (βΦ) ′′ = 1 1 -e -β .
(4.28) Also, (4.7) yields

T 0 (β) = 4 2 -e -β + 1 (1 -e -β ) 2 - 3 1 -e -β .
ĝ′ (β) ĝ(β) = -3 + 1 1 -e -β + 4(1 -e -β ) 2 -e -β - 1 -e -β 2 -e -β β + 2 - 2 2 -e -β + 1 2 (1 -e -β )e -β 1 (1 -e -β ) 2 - 2 (2 -e -β ) 2 (4.30)
and this integrates to

ĝ(β) = (const.) e -β 2 /4 e -β/2 1 -e -β 2 -e -β 3/2 exp β 2 log 1 - 1 2 e -β + β 2 φ(β) (4.31)
where const. is an arbitrary constant and φ(β) is defined by (4.15). We next assume that (4.9) asymptotically matches to the expansion for w > 1; we write this condition symbolically as

e nΦ(β) β n ĝ(β) β→∞ ∼ w ( n 2 ) 2 n-1 B * (w) w→1 . (4.32)
The matching condition applies on some intermediate scale where

β = n(w -1) → ∞ but w → 1. For β → ∞ we have Φ(β) = β 2 + log 2 - 1 β ∞ 0 log 1 - 1 2 e -v dv + O e -β β , βφ(β) → - ∞ 0 log 1 - 1 2 e -v dv ≡ d 0 = .58224 . . . (4.33)
and thus ĝ(β) ∼ (const.) e -β 2 /4 e -β/2 e d0/2 2 -3/2 , β → ∞. (4.34)

For w → 1, we have

w ( n 2 ) = 1 + β n n(n-1) 2 = e nβ/2 e -β 2 /4 e -β/2 1 + O β 3 n . ( 4
B * (w) ∼ (const.) √ 2 √ w -1 e d0/2 exp d 0 w -1 , w → 1 + . (4.36)
We also give numerical evidence for this behavior in Section 9. We will show in Section 5 that asymptotically matching the β-scale expansion (4.9), as β → 0 + , to the central region expansion (that applies for We comment that using the following ansatz

n 3/2 (w -1) = a = O(1)),
f (β; n) ∼ e nΦ(β) n Ψ(β) g(β) (4.38) in (4.
3), which is slightly more general than (4.9), we would find that Φ(β) is as in (4.15), that

Ψ(β) = Ψ 0 is a constant,
and then

g(β) = β -Ψ0 ĝ(β)
where ĝ is as in (4.31). Then matching to the a-scale result will show that Ψ 0 = -1 2 and fix the multiplicative constant as in (4.37). Using (4.38) and matching to the expansion for w > 1 would yield

B * (w) ∼ d 1 (w -1) -Ψ0 exp d 0 w -1 , w → 1.
Finally, we use (4.9) to calculate b(n, p) from the Cauchy formula (2.3). We write

w -p-1 = 1 + β n -p-1 = exp - p n β + p 2n 2 β 2 + O pβ 3 n 3 . (4.39)
Thus if we let p, n → ∞ in such a way that p/n 2 is fixed, then e nΦ(β) and e -pβ/n grow at the same linear rate in n, for β = O(1). Then (2.3) will have saddle point(s) where

Φ ′ (β) = 1 2 + φ ′ (β) = Λ ≡ p n 2 , or 1 2 -Λ - 1 β log 1 - e -β 2 = 1 β 2 log 2 -log(1-1 2 e -β ) ζ e ζ -1 dζ. (4.40)
This transcendental equation has a unique real solution

β = β * = β * (Λ) that satisfies β * → ∞ as Λ ↑ 1 2 , β * → 0 + as Λ → 0 + . In view of p max (n) in (2.13) we need only consider Λ ∈ 0, 1 2 . The steepest descent directions at β = β * are arg(β -β * ) = ± π
2 and then (4.39) and (4.40) lead to the estimate

b(n, p) ∼ √ β * n 2 ĝ(β * )e Λβ 2 * /2 1 2πφ ′′ (β * ) exp n log 2 + β * 2 + φ(β * ) -Λβ * . (4.41)
In view of (4.40) we have

φ(β * ) = 1 2 -Λ β * -log 1 - 1 2 e -β * (4.42)
and (4.15) yields

φ ′′ (β * ) = 1 β * 1 -2Λ - 1 2e β * -1 . ( 4
(n, p) ∼ 2 n+1 n 2 √ 2 π β * e -β * /2 1 -e -β * 2 -e -β * 3/2 1 -2Λ - 1 2e β * -1 -1/2 × exp n β * (1 -2Λ) -log 1 - 1 2 e -β * , (4.44)
which establishes (2.39).

Finally we discuss the asymptotic matching between (4.44), as Λ ↑ 1 2 , and (3.10), as L → ∞. We can solve (4.40) for β * asymptotically, as β * becomes large. We have

φ(β) = d 0 β - 1 2β e -β + O R (e -2β ), β → ∞ (4.45)
and thus (4.40) becomes

Λ - 1 2 = - d 0 β 2 + 1 2β + 1 2β 2 e -β + O R (e -2β ). (4.46)
Here again O R means roughly of the order, and neglects factors algebraic in β. Inverting the relation in (4.46) we find that

β * ∼ 2d 0 1 -2Λ 1 - 1 4d 0 2d 0 1 -2Λ + 1 exp - 2d 0 1 -2Λ , Λ ↑ 1 2 . ( 4
√ 2d 0 2 n πn 2 (1 -2Λ) exp - 1 2 2d 0 1 -2Λ + n 2d 0 √ 1 -2Λ . (4.48)
We show that (4.48) agrees with (3.10) as L → ∞. To expand (3.10) for L large we argue that there is a saddle in the range w ∼ 1 and approximate B * (w) by (3.9), which leads to the integral

2 n-1 2πi d 1 C w L-1 √ w -1 exp d 0 w -1 dw. (4.49)
This has a saddle where

d dw L log w + d 0 w -1 = 0 ⇒ w = w S ≡ 1 + d 0 2L + d 0 L 1 + d 0 4L
and then the standard estimate for (4.49) is

2 n-1 d 1 √ 2π √ w S -1 2d 0 w 2 S (w S -1) 3 -L -1/2 exp L log w S + d 0 w S -1 . (4.50)
For L → ∞ we can further simplify (4.50) using

w S = 1 + d 0 L + d 0 2L + O(L -3/2 ) (4.51)
and then we note that

L = n 2 -p = n 2 2 -p - n 2 = n 2 1 2 -Λ - n 2 . (4.52)
With (4.51), (4.50) simplifies to

2 n L d 1 √ d 0 4 √ π exp 2 Ld 0 - d 0 2 = 2 n L √ d 0 √ 2π exp(2 Ld 0 ). (4.53)
But, if we use (4.52) in (4.53) and expand for Λ → 1 2 , n → ∞ with n 1 2 -Λ → ∞, we obtain precisely (4.48). This verifies the asymptotic matching between the L-scale and Λ-scale results.

Central Region

In this section we analyze (2.1) for w -1 = O(n -3/2 ) and then obtain an expansion for b(n, p) that applies for p = O(n 3/2 ). We define a by

w -1 = a n 3/2 , -∞ < a < ∞.
(5.1)

At times we will need to separately consider the cases a < 0 and a > 0.

We set

B n (w) = 4 n a n 4 -n + 1 n 3/2 Cn (a) = 1 n + 1 2n n + 4 n n 3/2 Cn (n 3/2 (w -1)) (5.2)
and note that (2.4) yields Cn (0) = 0.

(5.3)

Since a n satisfies the recurrence

a n+1 = n ℓ=0 a ℓ a n-ℓ , a 0 = 1, (5.4)
using (5.2) in (2.1) leads to

4 (n + 1) 3/2 Cn+1 1 + 1 n 3/2 a = a n+1 4 -n 1 + a n 3/2 n -1 +2 1 + a n 3/2 n n ℓ=0 a ℓ 4 -ℓ (n -ℓ) 3/2 Cn-ℓ 1 - ℓ n 3/2 a + 1 + a n 3/2 n n ℓ=0 1 ℓ(n -ℓ) 3/2 Cℓ ℓ n 3/2 a Cn-ℓ 1 - ℓ n 3/2
a .

(5.5)

We can also write B n (w) as a Taylor series about w = 1, setting

B n (w) = ∞ j=0 M j,n j! (w -1) j . (5.6)
In view of (2.5) we have M 0,n = a n , M 1,n = b n , and M 2,n = c n . By multiplying (2.1) by w -n , differentiating N times with respect to w and setting w = 1, we obtain

N i=0 N i (-1) N -i (n + N -i -1)! (n -1)! M i,n+1 = n ℓ=0 N j=0 N j M j,ℓ M N -j,n-ℓ . (5.7)
We will obtain asymptotic approximations to Cn (a), or, equivalently, M j,n . In subsection 5.1 we shall analyze (5.7), while in subsection 5.2 we analyze (5.5). Then in subsection 5.3 we will use the results to obtain the approximation to b(n, p). We can also get analogous results by analyzing the functional equation (2.2) for the double transform; we discuss this further in Appendix A.

Moment equations

We consider (5.7). For n → ∞ we write

M i,n ≡ 4 n Mi,n = 4 n n 3 2 (i-1) m i + 1 √ n mi + O(n -1 ) (5.8)
and (2.9) shows that

m 0 = 1 √ π , m0 = 0, m 1 = 1, m1 = - 3 √ π , m 2 = 10 3 √ π , m2 = - 13 2 .
(5.9)

Isolating the terms with i = N , i = N -1, and j = 0, j = N in (5.7) and using (5.8) we obtain

MN,n+1 -nN MN-1,n+1 + N -2 i=0 N i (-1) N -i (n + N -i -1)! (n -1)! Mi,n+1 (5.10)
= 1 2 n ℓ=0 M0,ℓ MN,n-ℓ + 1 4 N -1 j=1 N j n ℓ=0
Mj,ℓ MN-j,n-ℓ .

Using (5.8) in the double sum in (5.10) we obtain

N -1 j=1 N j n ℓ=0 Mj,ℓ MN-j,n-ℓ ∼ N -1 j=1 N j n ℓ=0 n 3 2 (j-1) (n -ℓ) 3 2 (N -j-1) m j m N -j (5.11) ∼ N -1 j=1 m j m N -j 1 0 x 3 2 (j-1) (1 -x) 3 2 (N -j-1) dx n 3 2 N n -2 .
Here we approximated the inner sum by an integral via the Euler-MacLaurin formula. By definition (5.8), we have M0,ℓ = a ℓ 4 -ℓ and we write

n ℓ=0 a ℓ 4 -ℓ MN,n-ℓ = n ℓ=0 a ℓ 4 -ℓ MN,n-ℓ -MN,n + MN,n n ℓ=0 a ℓ 4 -ℓ .
(5.12)

Using the generating function a(z) below (2.3) and the estimate in (2.9) we find that

n ℓ=0 a ℓ 4 -ℓ = 2 - 2 √ πn + O(n -3/2 ).
(5.13)

Again by (5.8) and the Euler-MacLaurin formula we obtain

n ℓ=0 a ℓ 4 -ℓ MN,n-ℓ -MN,n ∼ n 3 2 N n -2 m N 1 √ π 1 0 1 x 3/2 (1 -x) 3 2 (N -1) -1 dx.
(

We also have

MN,n+1 -nN MN-1,n+1 = (n + 1) 3 2 (N -1) m N + 1 √ n + 1 mN + O(n -1 ) (5.15) -nN (n + 1) 3 2 (N -2) m N -1 + 1 √ n + 1 mN-1 + O(n -1 ) = m N n 3 2 N n -3 2 -N m N -1 n 3 2 N n -2 + O n 3 2 N n -5 2 .
Using (5.11) -(5.15) in (5.10), dividing by n 3 2 N n -2 , and letting n → ∞ we obtain the limiting equation

-N m N -1 = - 1 √ π m N + m N 1 2 √ π 1 0 (1 -x) 3 2 (N -1) -1 x 3/2 dx + 1 4 N -1 j=1 N j m j m N -j 1 0 x 3 2 (j-1) (1 -x) 3 2 (N -j-1) dx, (5.16
) and replacing N by N + 2 in (5.16) we obtain

which applies for N ≥ 1. Setting u N = m N +1 (5.
-(N + 2)u N = - 3 2 √ π (N + 1)u N +1 1 0 y 3 2 N + 1 2 √ 1 -y dy + 1 4 N i=0 (N + 2)! (i + 1)!(N + 1 -i)! u i u N -i 1 0 x 3 2 i (1 -x) 3 2 (N -i) dx.
(5.18)

Here we also used

1 0 (1 -x) 3 2 (N -1) -1 x 3/2 dx = 2 -3 1 0 N -1 √ x (1 -x) 1 2 (3N -5) dx, (5.19)
which follow by integration by parts. Equation (5.18) can be recast by defining

D(y) = ∞ j=0 u j (j + 1)! y 3 2 j , y > 0.
(5.20)

Multiplying (5.18) by y 3 2 N /(N + 2)! and summing over n, we obtain

-D(y) = - 1 y √ π y 0 D ′ (v) √ y -v dv + 1 4y y 0 D(v)D(y -v)dv. (5.21)
In view of (5.6) we have the corresponding leading order approximation to B n (w) (for n → ∞ with a fixed) as

B n (w) ∼ ∞ j=0 4 n n 3 2 (j-1) m j j! (w -1) j = 4 n n 3/2 ∞ j=0 m j j! a j = 4 n n 3/2   1 √ π + a ∞ j=0 m j+1 (j + 1)! a j   = 4 n n 3/2 1 √ π + aD(a 2/3 ) , a > 0. (5.22)
We can simplify (5.18) by setting

u N = (N + 1)! Γ 3 2 N + 1 V N , V 0 = 1 (5.23)
and using the identities

1 0 x 3 2 i (1 -x) 3 2 (N -i) dx = B 3 2 i + 1, 3 2 (N -i) + 1 = Γ 3 2 i + 1 Γ 3 2 (N -i) + 1 Γ 3 2 N + 2 1 0 1 √ 1 -x x 3 2 N + 1 2 dx = B 3 2 N + 3 2 , 1 2 = √ π Γ 3 2 N + 3 2 Γ 3 2 N + 2 , ( 5.24)
where B(•, •) and Γ(•) are the Beta and Gamma functions, respectively. Using (5.23) and (5.24), (5.18) becomes

V N +1 = 3 2 N + 1 V N + 1 4 N i=0 V i V N -i , V 0 = 1. (5.25)
Since the V i are all positive, (5.25) shows that V N +1 ≥ 3 2 N + 1 V N and and thus

V N ≥ 3 2 N Γ N + 2 3 Γ 2 3
, N ≥ 0.

(5.26)

This estimate shows that the V N grow faster than exponentially, and thus (5.25) cannot be easily solved by generating functions, despite the non-linearity having the form of a convolution sum. We can also rewrite (5.18) by introducing D(Y )

D(Y ) = ∞ j=1 Y 3 2 (j-1) m j j! (-1) j = ∞ i=0 Y 3 2 i u i (i + 1)! (-1) i+1 (5.27)
and then similarly to (5.22) we find that

B n (w) ∼ 4 n n 3/2 1 √ π -a D((-a) 2/3 ) , -a = Y 3/2 > 0. (5.28)
This gives a leading order approximation to B n (w) for n → ∞ and a fixed a < 0. We also note that

∞ L=0 m L L! a L = m 0 + ∞ L=1 m L L! (-a) L-1 (-a)(-1) L = 1 √ π -a D(Y ).
(5.29) Using (5.27) in (5.18) leads to the non-linear integral equation

Y D(Y ) = - 1 √ π Y 0 D′ (v) √ Y -v dv + 1 4 Y 0 D(v) D(Y -v)dv.
(5.30)

Note that this differs from (5.21) only by the sign of the left-hand side. However, (5.30) is susceptible to solution by a Laplace transform, whereas (5.21) is not. Introducing

D * (s) = L{ D(Y )} = ∞ 0 e -sY D(Y )dY, (5.31)
where L is the Laplace transform operator, and noting that

L{Y D(Y )} = -D ′ * (s), D(0) = -m 1 = -1, L D′ (Y ) * 1 √ Y = [sD * (s) -D(0)] π s
(where * denotes convolution) we obtain

-D ′ * (s) = - √ sD * (s) + 1 4 D 2 * (s) - 1 √ s . (5.32)
This Riccatti equation is easily solved by setting

D * = 2 √ s + U (s), U (s) = 4 F ′ (s) F (s) , (5.33)
which leads to the Airy equation

F ′′ (s) = s 4 F (s).
(5.34)

The only solution with acceptable behavior as s → +∞ is given by

F (s) = (const.) Ai(4 -1/3 s) (5.35)
where Ai(•) is the Airy function. Using (5.35) in (5.33) and inverting the Laplace transform, we obtain the integral representation

D(Y ) = 1 2πi Br e sY 2 √ s + 4 2/3 Ai ′ (4 -1/3 s) Ai(4 -1/3 s) ds, Y = (-a) 2/3 (5.36)
which applies only for a < 0. Here Br is any vertical contour in which ℜ(s) > r 0 = max{z : Ai(z) = 0}. We can also write (5.36) as

D(Y ) = 1 2πi d dY Br e sY 2 √ s + 4 2/3 s Ai ′ (4 -1/3 s) Ai(4 -1/3 s) ds = d dY   2 √ πY + 4 2/3 ∞ j=0 1 r j e 4 1/3 rj Y   = 1 √ πa + 4 ∞ j=0 e 4 1/3 rj Y .
(5.37)

Here 0 > r 0 > r 1 > r 2 > . . . where r j are the roots of the Airy function Ai(•), and we evaluated the integral as a residue series, by closing the Br contour in the left half-plane. By using (5.37) in (5.28) we obtain the leading order approximation

B n (w) ∼ 4 n+1 n 3/2 (-a) ∞ j=0 e -4 1/3 |rj |Y , Y = (-a) 2/3 > 0.
(5.38)

We note that (5.38) is consistent with the fact that

B n (1) = a n ∼ 4 n /(n 3/2 √ π) as n → ∞. It is well known that as z → -∞ Ai(z) ∼ 1 √ π (-z) -1/4 sin 2 3 (-z) 3/2 + π 4 (5.39)
and thus the roots r j satisfy

|r j | ∼ 3jπ 2 2/3
, j → ∞.

(5.40) Thus, as Y → 0 + (w → 1 -), we can approximate the sum in (5.38) by the Euler-MacLaurin formula with

∞ j=0 e -4 1/3 |rj |Y ∼ 1 Y 3/2 ∞ 0 exp[-(3xπ) 2/3 ]dx = -1 4 √ πa .
We next derive a representation valid for a > 0, and show that as a → ∞ the a-scale expansion asymptotically matches to the β-scale result, as β → 0 + . We can get a rough idea of the behavior of the central region approximation as a = y 3/2 → +∞ by using the integral equation (5.21) and a WKB-type expansion. Let us assume that D(y) ∼ K(y)e Ψ(y) , y → ∞ (5.41)

where Ψ(y) ≫ log[K(y)]. We use (5.41) to approximate the various terms in (5.21). Expecting that Ψ ′ (y) > 0 and Ψ ′ (y) → ∞ as y → ∞, we have

y 0 1 √ v D ′ (y -v)dv ∼ y 0 1 √ v (K ′ e Ψ + KΨ ′ e Ψ )(y -v)dv ∼ ∞ 0 e Ψ(y) e -vΨ ′ (y) √ v 1 + 1 2 v 2 Ψ ′′ (y) (Ψ ′ (y) -vΨ ′′ (y))(K(y) -vK ′ (y)) + K ′ (y) dv ∼ Ψ ′ (y)πK(y)e Ψ(y) + K(y)e Ψ(y) ∞ 0 e -vΨ ′ √ v K ′ K -vΨ ′′ -vΨ ′ K ′ K + 1 2 v 2 Ψ ′′ Ψ ′ dv = K(y)e Ψ(y) √ π √ Ψ ′ + K ′ K 1 √ Ψ ′ - 1 2(Ψ ′ ) 3/2 Ψ ′′ + Ψ ′ K ′ K + 3 8 Ψ ′′ (Ψ ′ ) 3/2 .
(5.42)

The non-linear term in (5.21) we approximate as

y 0 D(v)D(y -v)dv = 2 y/2 0 D(v)D(y -v)dv ∼ 2 ∞ 0 D(v)e Ψ(y) e -vΨ ′ (y) K(y)dv ∼ 2D(0)K(y)e Ψ ′ (y) 1 Ψ ′ (y) . (5.43)
Recalling that D(0) = u 0 = m 1 = 1, we use (5.42) and (5.43) in (5.21). After multiplying by -y we obtain, to leading order, yK(y)e Ψ(y) ∼ Ψ ′ (y)K(y)e Ψ(y) so that

Ψ(y) = y 3 3 = a 2 3 .
(5.44)

At the next order, after cancelling the common factors Ke Ψ , we are led to

0 = - 1 2Ψ ′ + K ′ K 1 √ Ψ ′ - 1 2(Ψ ′ ) 3/2 Ψ ′′ + Ψ ′ K ′ K + 3 8 Ψ ′′ (Ψ ′ ) 3/2
(5.45) and thus

K ′ K = 1 4
Ψ ′′ Ψ ′ + 1 √ Ψ ′ = 1 2y + 1 y = 3 2y
so that K(y) = (const.) ′ y 3/2 = (const.) ′ a. We thus have formally obtained

Cn (a) ∼ (const.) ′ a 2 e a 2 /3 , a → +∞ (5.46)
for some constant (const.) ′ . We note that Cn (a) ∼ aD(a 2/3 ) dominates the first term in the right side of (5.2). In view of (5.2), (4.2) and (4.9) the asymptotic matching of the aand β-scales requires that

4 n n 3/2 Cn (a) a→∞ ∼ e nΦ(β) n -1/2 g(β) β→0 + . ( 5
Φ(β) = log 4 + β 2 3 + O(β 3 ), g(β) ∼ (const.) β 2 .
(5.48)

Since β = a/ √ n we see that the matching is indeed possible provided that const. ′ = const.

(5.49)

where the latter constant is the one that arose in Section 4. We note that our formal analysis suggests that the non-linear integral equation (5.21) may be approximated by a linear one for y → ∞. The non-linear term does not affect the exponential growth rate Ψ = y 3 /3, but it does affect the algebraic growth rate K ∝ y 3/2 . To determine the remaining constant, we continue (5.37) into the range a > 0. Since Y = (-a) 2/3 , this corresponds to arg(Y ) = ± 2π 3 . We define h(s) = Ai(4 -1/3 s).

(5.50)

By deforming the Bromwich contour in (5.37) to a piecewise linear contour that goes from s = e -2πi/3 ∞ to s = 0 and then from s = 0 to s = e 2πi/3 ∞ = ω∞, and parameterizing the two pieces, we are led to

D(Y ) - 1 √ πa = 4 2πi ∞ 0 ω h ′ (ωτ ) h(ωτ ) e ωτ Y -ω 2 h ′ (ω 2 τ ) h(ω 2 τ ) e ω 2 τ Y dτ.
(5.51)

Since ℜ(ω), ℜ(ω 2 ) < 0, these integrals converge for Y > 0. We can also write the approximation to

B n (w) for w = 1 + O(n -3/2 ) as B n (w) ∼ 4 n+1 πn 3/2 (-a) ∞ 0 ℜ e πi/6 h ′ (ωτ ) h(ωτ ) e ωτ Y dτ. (5.52)
We shall use the Wronskian identity (cf. [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF][START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers[END_REF])

Ai(ωz)ω 2 Ai ′ (ω 2 z) -ωAi ′ (ωz)Ai(ω 2 z) = i 2π
which in terms of h(•) becomes

e πi/6 h ′ (ωτ ) h(ωτ ) + e -πi/6 h ′ (ω 2 τ ) h(ω 2 τ ) = - 1 2π4 1/3 1 h(ωτ )h(ω 2 τ ) .
(5.53)

The integral

I 1 = ∞ 0 e πi/6 h ′ (ωτ ) h(ωτ ) e ωτ Y dτ, Y > 0
may be continued analytically into the range arg(Y ) ∈ -π 6 , 5π 6 , while

I 2 = ∞ 0 e -πi/6 h ′ (ω 2 τ ) h(ω 2 τ ) e ω 2 τ
I 3 = - 1 2π4 1/3 ∞ 0 e ωτ Y h(ωτ )h(ω 2 τ ) dτ (5.54) I 4 = - ∞ 0 e -πi/6 h ′ (ω 2 τ ) h(ω 2 τ ) e ωτ Y dτ. (5.55)
Since Ai(•) satisfies (cf. [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]) where a is real and positive. We have thus shown that I 1 + I 2 = I 2 + I 3 + I 4 . This continues the right side of (5.52) to a > 0 and shows that

Ai(z) ∼ 1 2 √ π z -1/4 exp - 2 3 z 3/2 , | arg z| < π (
4 π (-a) ∞ 0 ℜ e πi/6 h ′ (ωτ ) h(ωτ ) e ωτ Y dτ = a π 2 4 1/3 ∞ 0 e τ y h(ωτ )h(ω 2 τ ) dτ (5.58) - 4a π ∞ 0 ℜ e πi/6 h ′ (ωτ ) h(ωτ ) e ω 2 τ y dτ.
Next we estimate the right side as a = y 3/2 → +∞. The second integral can be expanded by Watson's lemma, and is O(a/y) = O( √ y). The first integral we evaluate by Laplace's method, first using (5.56) to approximate the integrand for τ → +∞, with

h(ωτ )h(ω 2 τ ) = |h(ωτ )| 2 ∼ (4 -1/3 τ ) -1/2 4π exp 4 3 (4 -1/3 τ ) 3/2 .
Here we also used the reflection principle, since Ai(z) is real for real z. We thus have

4 2/3 4π 2 ∞ 0 e τ y |h(ωτ )| 2 dτ ∼ 4 2/3 π ∞ 0 e τ y e -2 3 τ 3/2 2 -1/3 √ τ dτ ∼ 2 π y 2 ∞ -∞ exp y 3 3 - 1 4y (τ -y 2 ) 2 dτ = 4 √ π y 3/2 exp y 3 3 .
(5.59)

Here we used the Laplace method to find that the major contribution to the integral came from τ = y2 . Since this is asymptotically large, it is permissible to first approximate the integrand for τ large. We note that the right side of (5.59) is also the expansion of D(y) as y → +∞. Thus, this calculation verifies that obtained formally by the WKB ansatz (5.41), and also determines the constant as const. = 4 √ π .

(5.60)

The leading term on the a-scale for a > 0 is thus given by 4 n /n 3/2 times the right side of (5.58), and this is precisely the result in (2.29). We next calculate the correction term mi in (5.8). We shall show that this corresponds to

C (1) (a) if Cn (a) has the expansion Cn (a) = C (0) (a) + n -1/2 C (1) (a) + O(n -1
). We return to (5.10) and use (5.8).

For n → ∞ we have

N i=0 N i (-1) N -i (n + N -i -1)! (n -1)! Mi,n+1 = MN,n+1 -nN MN-1,n+1 + N 2 n(n + 1) MN-2,n+1 + O n 3 2 N -3 = n 3 2 (N -1) m N + 1 √ n ( mN -N m N -1 ) + 1 n mN + 3 2 (N -1)m N -N mN-1 + N 2 m N -2 + O(n -3/2 ) (5.61)
where mN is the third term in the expansion in (5.8), i.e.,

MN,n = n 3 2 (N -1) m N + 1 √ n mN + 1 n mN + O(n -3/2 ) .
We next estimate more precisely (5.12):

1 2 n ℓ=0 a ℓ [ MN,n-ℓ -MN,n ] + 1 2 MN,n n ℓ=0 a ℓ 4 -ℓ .
(5.62) From (5.13) and (5.8) we obtain

1 2 MN,n n ℓ=0 a ℓ 4 -ℓ = n
To obtain the first term in (5.62) we use an intermediate limit, writing to sum as

n ℓ=0 a ℓ 4 -ℓ [ MN,n-ℓ -MN,n ] (5.63) = n ℓ=0 a ℓ 4 -ℓ n 3 2 (N -1) 1 - ℓ n 3 2 (N -1) m N + mN √ n -ℓ -m N - mN √ n + O(n -1
) .

Now, by Euler-MacLaurin

n ℓ=0 a ℓ 4 -ℓ n 3 2 (N -1) mN √ n 1 - ℓ n 3 2 N -2 -1 ∼ 1 √ π n 3 2 N -5 2 mN (5.64) × 1 0 1 x 3/2 (1 -x)
-1 = L-1 ℓ=0 a ℓ 4 -ℓ 1 - ℓ n 3 2 (N -1)
-1 (5.65)

+ n ℓ=L a ℓ 4 -ℓ 1 - ℓ n 3 2 (N -1) -1 ≡ S 1 + S 2 .
Here L satisfies the asymptotic bounds L → ∞ but L/n → 0 as n → ∞. We use a binomial approximation in the first sum and approximate a ℓ for ℓ → ∞ in the second. Thus

S 1 = L-1 ℓ=0 a ℓ 4 -ℓ - 3 2 (N -1) ℓ n + O ℓ 2 n 2 = - 3 2 N -1 n L-1 ℓ=0 ℓa ℓ 4 -ℓ + O L 3/2 n 2 = - 3 2 N -1 n L ℓ=0 ℓa ℓ 4 -ℓ + O 1 n √ L , L 3/2 n 2 .
(5.66) Using (2.4) we have

L ℓ=0 ℓa ℓ 4 -ℓ = L ℓ=0 (ℓ + 1)4 -ℓ a ℓ - L ℓ=0 4 -ℓ a ℓ = L ℓ=0 2ℓ ℓ 4 -ℓ -2 + o L (1) = 2(L + 1) 4 L+1 2L + 2 L + 1 -2 + o L (1) = 2 √ L √ π -2 + O L 1 √ L .
(5.67)

To evaluate S 2 we use Euler-MacLaurin to find that

S 2 = n ℓ=L 1 √ π ℓ -3/2 + O(ℓ -5/2 ) 1 - ℓ n 3 2 (N +1) -1 = 1 √ π 1 √ n 1 L/n 1 x 3/2 (1 -x) 3 2 (N -1) -1 dx + O 1 n √ L = 1 √ π √ n 1 0 1 x 3/2 (1 -x) 3 2 (N -1) -1 dx + 3 √ π (N -1) √ L n + O 1 n √ L . (5.68)
Combining (5.66) -(5.68) and letting L → ∞ in such a way that L 3/2 /n → 0 we have shown that

1 2 (S 1 + S 2 ) = 1 2 √ πn 1 0 1 x 3/2 (1 -x) 3 2 (N -1) -1 dx + 3 2 N -1 n + O(n -3/2 ).
(5.69)

Next, we expand the last sum in (5.10) as we obtain, multiplying by 4,

N -1 j=1 N j n ℓ=0 Mj,ℓ MN-j,n-ℓ = N -1 j=1 N j n 3 2 N -3 m j m N -j n ℓ=0 ℓ n 3 2 (j-1) 1 - ℓ n 3 2 (N -j-1) + 2 √ n m j mN-j n ℓ=0 ℓ n 3 2 (j-1) 1 - ℓ n 3 2 (N -j)-2 + O(n -1 ) = n 3 2 N -2 N -1 j=1 m j m N -j N j 1 0 x 3 2 (j-1) (1 -x) 3 2 (N -j-1) dx +n 3 2 N -5 2 N -1 j=1 2m j mN-j N j 1 0 x 3 2 (j-1) (1 -x) 3 2 (N -j)-2 dx + O(n 3 2 N -3 ). ( 5
2N (N -1)m N -2 -4N mN-1 = - 2 √ π (3N -4) mN 1 0 (1 -x) 3 2 N -3 √ x dx (5.71) + 2 N -1 j=1 N j m j mN-j 1 0 x 3 2 (j-1) (1 -x) 3 2 (N -j)-2 dx.
This holds for N ≥ 2. Here we also rewrote the integral in (5.64) using integration by parts. Note that when N = 2 (5.71) becomes

4m 0 -8 m1 = - 8 √ π m2 + 8m 1 m1
and this is consistent with (5.9).

We solve (5.71). Note that once m N is known this is a linear equation for mN . We set

D1 (Y ) = ∞ L=1 mL L! (-1) L Y 3 2 (L-1) (5.72)
and recall the definition of D(Y ) in (5.27). Since m 0 = 1/ √ π, we have

∞ N =2 2N (N -1)m N -2 (-1) N N ! Y 3 2 (N -1) = 2 ∞ L=0 m L L! (-1) L Y 3 2 (L+1) = 2 Y 3/2 √ π + Y 3 D(Y ) . (5.73) Also, -4 ∞ N =2 N mN-1 (-1) N N ! Y 3 2 (N -1) = 4 ∞ L=1 mL L! (-1) L Y 3 2 L = 4Y 3/2 D1 (Y ), (5.74)
∞ N =2 (-1) N 2Y 3 2 (N -1) N -1 j=1 m j mN-j j!(N -j)! 1 0 x 3 2 (j-1) (1 -x) 3 2 (N -j-1) dx √ 1 -x = 2Y 3/2 1 0 D(Y x) D1 (Y -Y x) √ 1 -x dx, (5.75)
and

- 2 √ π ∞ N =2 (3N -4) (-1) N N ! Y 3 2 (N -1) mN 1 0 (1 -x) 3 2 N -3 √ x dx = - 2 √ π 1 0 1 x(1 -x) 1 1 -x ∞ N =2 3(N -1) -1 N ! (-1) N mN (Y -Y x) 3 2 (N -1) dx = - 2 √ π 1 0 1 √ x(1 -x) 3/2 [2(Y -Y x) D1 (Y -Y x) -D1 (Y -Y x) -m1 ]dx = - 4 √ π 1 0 Y x(1 -x) D′ 1 (Y -Y x)dx + 2 √ π 1 0 D1 (Y -Y x) -D1 (0) √ x(1 -x) 3/2 dx = - 4Y √ π Y 0 D′ 1 (Y -v) √ v √ Y -v dv + 2Y √ π Y 0 D1 (v) -D1 (0) √ Y -v v 3/2 dv. ( 5
2Y 1/2 √ π + 2Y 2 D(Y ) + 4Y 1/2 D1 (Y ) = 2 Y 0 D1 (v) √ v D(Y -v)dv - 4 √ π Y 0 1 √ Y -v d dv D1 (v) -D1 (0) √ v dv. (5.77)
To solve (5.77) we define

E * (s) = ∞ 0 e -sY D1 (Y ) √ Y dY (5.78)
and note that

∞ 0 e -sY Y 0 1 √ Y -v d dv D1 (v) -D1 (0) √ v dv dY = π s ∞ 0 e -sY d dY D1 (Y ) -D1 (0) √ Y dY = √ πs ∞ 0 e -sY D1 (Y ) -D1 (0) √ Y dY = √ πs E * (s) -π D1 (0). (5.79)
Since D1 (0) = -m1 = 3/ √ π, the Laplace transform of (5.77) becomes

1 s 3/2 + 2D ′′ * (s) -4E ′ * (s) = [2D * (s) -4 √ s]E * (s) + 12.
(5.80)

Recalling that h(s) = Ai(4 -1/3 s) we write (5.33) as

D * (s) = 2 √ s + 4 h ′ (s) h(s) (5.81)
whose derivatives yield

D ′ * (s) = 1 √ s + s -4 h ′ (s) h(s) 2
(5.82) and

D ′′ * (s) = - 1 2s 3/2 + 1 -2s h ′ (s) h(s) + 8 h ′ (s) h(s) 3 , (5.83)
where we used the fact that h ′′ (s) = 1 4 sh(s). Using (5.81) -(5.83) in (5.80) we obtain, after some rearrangement,

E ′ * (s) + 2 h ′ (s) h(s) E * (s) = - 5 2 -s h ′ (s) h(s) + 4 h ′ (s) h(s) 3 (5.84) or d ds [h 2 (s)E * (s)] = - 5 2 h 2 (s) -sh(s)h ′ (s) + 4 (h ′ (s)) 3 h(s) .
(5.85)

The general solution of (5.85) may be written as

E * (s) = c * h 2 (s) + 1 h 2 (s) ∞ s 5 2 h 2 (v) + vh(v)h ′ (v) -4 (h ′ (v) 3 h(v) dv.
(5.86)

But, h(s) decays exponentially as s → ∞ and so h -2 grows. Thus, if E * (s) is to be a proper Laplace transform, we must set c * = 0. Use of integration by parts and the fact that h(s) satisfies (5.34) leads to

s vh(v)h ′ (v)dv = 4 s h ′ (v)h ′′ (v)dv = 2(h ′ (s)) 2 s h 2 (v)dv = sh 2 (s) -2 s vh(v)h ′ (v)dv = sh 2 (s) -4(h ′ (s)) 2
and thus (5.86) simplifies to

E * (s) = - 5 2 s + 8 h ′ (s) h(s) 2 - 4 h 2 (s) ∞ s (h ′ (v)) 3 h(v) dv.
(5.87)

We have thus established (2.26), since for a < 0

B n (w) = ∞ j=0 M j,n (w -1) j j! = 4 n n 3/2   ∞ j=0 m j a j j! + 1 √ n ∞ j=0 mj a j j! + O(n -1 )   = 4 n n 3/2   m 0 + ∞ j=1 m j a j j! + 1 √ n ∞ j=1 mj a j j! + O(n -1 )   = 4 n n 3/2 1 √ π -a D((-a) 2/3 ) + 1 √ n (-a) D1 (Y ) + O(n -1 ) = 4 n n 3/2 1 √ π + Y 3/2 D(Y ) + 1 √ n Y 3/2 D1 (Y ) + O(n -1
) .

(5.88)

In order to asymptotically match the a-scale result to that valid for w -1 = O(n -1 ) with w < 1, we will need the behavior of D1 (Y ) as Y = (-a) 2/3 → ∞. Thus we must locate the singularity of E * (s) with the largest real part, or an appropriate saddle point. The Laplace transform of the leading term D(Y ) had simple poles at 4 1/3 r j where r j were the negative roots of the Airy function, with the pole at 4 1/3 r 0 determining the asymptotic behavior as Y → ∞ (a → -∞) (cf. (5.37)). The dominant singularity in E * (s) is also at 4 1/3 , r 0 . but (5.87) shows that the structure is more complicated, as there is a double pole combined with a logarithmic branch point at s = s 0 ≡ 4 1/3 r 0 .

We expand (5.87) near the dominant singularity at s 0 . Integrating by parts we obtain

∞ s (h ′ (v)) 3 h(v) dv = ∞ s (h ′ (v)) 2 d[log h(v)] (5.89) = -(h ′ (s)) 2 [log h(s)] - 1 2 ∞ s vh ′ (v)h(v)[log h(v)]dv and ∞ s vh(v)h ′ (v)[log h(v)]dv = ∞ s vh(v)[log h(v)]dh(v) (5.90) = -sh 2 (s)[log h(s)] - ∞ s {h 2 (v)[log h(v)] + vh(v)h ′ (v) log h(v) + vh(v)h ′ (v)}dv
which may be rewritten as

2 ∞ s vh(v)h ′ (v)[log h(v)]dv = -sh 2 (s)[log h(s)] + 2(h ′ (s)) 2 - ∞ s h 2 (v)[log h(v)]dv.
(5.91) Using (5.89) -(5.91) in (5.87) we obtain

E * (s) = - 5 2 s + 10 h ′ (s) h(s) 2 + 4 h ′ (s) h(s) 2 log h(s) -s log h(s) - 1 h 2 (s) ∞ s h 2 (v)[log h(v)]dv.
(5.92)

We let τ = ss 0 and expand (5.92) about τ = 0. We have h(s 0 ) = h ′′ (s 0 ) = 0 and

∞ s h 2 (v)[log h(v)]dv = ∞ s0 h 2 (v)[log h(v)]dv - s s0 h 2 (v)[log h(v)]dv = ∞ s0 h 2 (v)[log h(v)]dv + O(τ 3 log τ ),
and then (5.92) yields

E * (s) = 4 τ 2 log(h ′ (s 0 )τ ) + 10 τ 2 - ∞ s0 h 2 (v)[log h(v)]dv (h ′ (s 0 )) 2 τ 2 + O(log τ ), τ → 0. ( 5

.93)

Here

h ′ (s 0 ) = 4 -1/3 Ai ′ (4 -1/3 s 0 ) = 4 -1/3 Ai ′ (r 0 ). (5.94) Defining κ ≡ 10 + 4 log(h ′ (s 0 )) - 1 (h ′ (s 0 )) 2 ∞ s0 h 2 (v)[log h(v)]dv
(5.95)

we have

D1 (Y ) = √ Y 2πi Br e sY E * (s)ds = √ Y exp(-|r 0 |4 1/3 Y ) 1 2πi Br+ e τ Y 4 log τ τ 2 + κ τ 2 + O(log τ ) dτ (5.96)
where ℜ(τ ) > 0 on Br + . Evaluating explicitly the integrals in (5.96) shows that

D1 (Y ) = Y 3/2 exp(-|r 0 |4 1/3 Y ) -4 log Y + κ + 4 -4γ E + O log Y Y , Y → ∞ (5.97)
where γ E is the Euler constant. Combining (5.97) with (5.38) we have shown that as Y → ∞ the two term approximation on the a-scale behaves as

B n (w) = 4 n n 3/2 1 √ π -a D(Y ) - a √ n D1 (Y ) + O(n -1 ) (5.98) ∼ 4 n n 3/2 exp -4 1/3 |r 0 |(-a) 2/3 -4a + a 2 √ n - 8 3 log(-a) + κ + 4 -4γ E . (5.99)
This expression will be used in Section 6.

To summarize, we have used the moment equations (5.7) to get the two-term approximation (5.98) to B n (w), valid for a = n 3/2 (w -1) = O(1) and a < 0. For a > 0 we must use (5.51) and (5.58) with y = a 2/3 .

Analysis of the basic recurrence

We consider (2.1), which on the a-scale becomes (5.5). Using (2.9) we obtain

a n+1 4 -n 1 + a n 3/2 n -1 = 4 √ π a n 2 + 2 √ π a 2 n 5/2 + O(n -3 ).
(5.100)

We expand Cn (a) in (5.2) as

Cn (a) = C (0) (a) + 1 √ n C (1) (a) + O(n -1 ) (5.101)
and then estimate the various terms in (5.5). We note that C (j) (0) = 0 for all j ≥ 0, in view of (5.3).

By (5.101) and the Euler-MacLaurin formula we have

n ℓ=0 1 ℓ(n -ℓ) 3/2 Cℓ ℓ n 3/2 a Cn-ℓ 1 - ℓ n 3/2 a = 1 n 2 1 0 C (0) (x 3/2 a)C (0) ((1 -x) 3/2 a) [x(1 -x)] 3/2 dx + 2 n 5/2 1 0 C (0) (x 3/2 a)C (1) ((1 -x) 3/2 a) x 3/2 (1 -x) 2 dx + O(n -3
).

(5.102)

Note that the integrals converge since C (j) (0) = 0, and that the boundary effects near x = 0 and x = 1 do not appear until the third (O(n -3 )) term in the expansion.

Next we write the first sum in the right side of (5.5) as

n ℓ=0 a ℓ 4 -ℓ (n -ℓ) -3/2 Cn-ℓ 1 - ℓ n 3/2 a = n ℓ=0 a ℓ 4 -ℓ 1 n 3/2 1 - ℓ n -3/2 C (0) 1 - ℓ n 3/2 a -C (0) (a) + n ℓ=0 a ℓ 4 -ℓ 1 n 2 1 - ℓ n 2 C (1) 1 - ℓ n 3/2 a -C (1) (a) + C (0) (a) n 3/2 n ℓ=0 a ℓ 4 -ℓ + C (1) (a) n 2 n ℓ=0 a ℓ 4 -ℓ + C (2) (a) n 5/2 n ℓ=0 a ℓ 4 -ℓ + O(n -3 ).
(5.103)

The last three terms in (5.103) may be estimated using (5.13). By Euler-Maclaurin we have

n ℓ=0 a ℓ 4 -ℓ 1 - ℓ n -2 C (1) 1 - ℓ n 3/2 a -C (1) (a) ∼ 1 √ n 1 √ π 1 0 1 x 3/2 C (1) ((1 -x) 3/2 a) (1 -x) 2 -C (1) (a) dx (5.104)
where we approximated a ℓ for ℓ → ∞. An argument completely analogous to (5.65) shows that

n ℓ=0 a ℓ 4 -ℓ 1 - ℓ n -3/2 C (0) 1 - ℓ n 3/2 a C (0) (a) (5.105) = L ℓ=0 a ℓ 4 -ℓ 3 2 ℓ n [C (0) (a) -aC (0)′ (a)] + O ℓ 2 n 2 + n ℓ=L+1 1 √ π [ℓ -3/2 + O(ℓ -5/2 )] 1 - ℓ n -3/2 C (0) 1 - ℓ n 3/2 a -C (0) (a) = [C (0) (a) -aC (0)′ (a)] 3 2n L ℓ=0 ℓa ℓ 4 -ℓ + 1 √ πn 1 L/n x -3/2 [(1 -x) -3/2 C (0) ((1 -x) 3/2 a) -C (0) (a)]dx + O L 3/2 n 2 , 1 n √ L ∼ 1 √ πn 1 0 x -3/2 [(1 -x) -3/2 C (0) ((1 -x) 3/2 a) -C (0) (a)]dx +[C (0) (a) -aC (0)′ (a)] 3 2n L ℓ=0 ℓa ℓ 4 -ℓ -2 L π ∼ 1 √ πn 1 0 x -3/2 [(1 -x) -3/2 C (0) ((1 -x) 3/2 a) -C (0) (a)]dx - 3 n [C (0) (a) -aC (0)′ (a)].
Here we also used (5.67). Multiplying (5.5) by

w -n = 1 + a n 3/2 -n = 1 - a √ n + a 2 2n + O(n -3/2 )
and using (5.100) yields

4 (n + 1) 3/2 1 + a n 3/2 -n Cn+1 1 + 1 n 3/2 a -a n+1 4 -n 1 -1 + a n 3/2 -n = 4 n 3/2 1 - 3 2n + O(n -2 ) 1 - a √ n + a 2 2n + O(n -3/2 ) (5.106) × C (0) (a) + 1 √ n C (1) (a) + 1 n C (2) (a) + 3 2 aC (0)′ (a) - 4 √ πn 3/2 a √ n - a 2 2n + O(n -3/2 ) = 4 n 3/2 C (0) (a) + 4 n 2 -aC (0) (a) - a √ π + C (1) (a) + 4 n 5/2 -aC (1) (a) + a 2 2 √ π + a 2 2 C (0) (a) + 3 2 aC (0)′ (a) - 3 2 C (0) (a) + C (2) (a) + O(n -3 ).
Comparing (5.106) with (5.102) + 2 × (5.103), using (5.104), (5.105) and (5.13), we find at O(n -2 )

-4aC (0) (a) - 4 √ π a = - 4 √ π C (0) (a) + 2 √ π 1 0 1 x 3/2 C (0) ((1 -x) 3/2 a) (1 -x) 3/2 -C (0) (a) dx + 1 0 C (0) (x 3/2 a)C (0) ((1 -x) 3/2 a) x 3/2 (1 -x) 3/2 dx (5.107) and at O(n -5/2 ) -4aC (1) (a) + 2 √ π a 2 + 2a 2 C (0) (a) = - 4 √ π C (1) (a) + 2 √ π 1 0 1 x 3/2 C (1) ((1 -x) 3/2 a) (1 -x) 2 -C (1) (a) dx = 2 1 0 C (0) (x 3/2 a)C (1) ((1 -x) 3/2 a) x 3/2 (1 -x) 2 dx.
(5.108) Thus, the leading term in (5.101) satisfies the non-linear integral equation (5.107) and the correction term satisfies the linear equation (5.108).

We can analyze these equations for a < 0 by setting

C (0) (a) = (-a) D((-a) 2/3 ) = Y 3/2 D(Y ). (5.109)
Then (5.107) becomes

4a 2 D(Y ) - 4 √ π a = 4 √ π a D(Y ) + a 2 1 0 D(Y x) D(Y -Y x)dx - 2a √ π 1 0 D(Y -Y x) -D(Y )
x 3/2 dx.

(5.110)

Integrating by parts and using D(0) = -1 we obtain

1 0 x -3/2 [ D(Y -Y x) -D(Y )]dx = 2 + 2 D(Y ) -2Y 1 0 x -1/2 D′ (Y -Y x)dx.
(5.111) Using (5.111) in (5.110), dividing by 4a 2 /Y and noting that Y 3/2 = -a, we see that (5.110) is the same as (5.30). Similarly, we set

C (1) (a) = (-a) D1 ((a) 2/3 ) = Y 3/2 D1 (Y ) (5.112)
in (5.108). Using also (5.109) this yields

4a 2 D1 (Y ) + 2 √ π a 2 - 4 √ π a D1 (Y ) -2a 3 D(Y ) (5.113) = - 2a √ π 1 0 1 x 3/2 D1 (Y -Y x) √ 1 -x -D1 (Y ) dx + 2a 2 1 0 D(Y x) D1 (Y -Y x) √ 1 -x dx. Dividing (5.113) by a 2 / √ Y = Y 5/2 leads to 4 √ Y D1 (Y ) + 2 √ π √ Y + 4 √ πY D1 (Y ) + 2Y 2 D(Y ) (5.114) = 2 Y 0 D1 (v) √ v D(Y -v)dv + 2 √ πY Y 0 1 v 3/2 D1 (Y -v) Y Y -v -D1 (Y ) dv.
But the last term in (5.114) is

2 √ πY Z Y 0 1 v 3/2 " D1(Y -v) r Y Y -v -D1(0) r Y Y -v -D1(Y ) + D1(0) r Y Y -v -1 ! + D1(0) # dv = 4 √ πY D1(0) - 4 √ πY Z Y 0 ( D1(0) -D1(Y ) + r Y Y -v [ D1(Y -v) -D1(0)] ) d(v -1/2 ) = 4 √ πY D1(Y ) + 4 √ πY Z Y 0 1 √ v d dv » D1(Y -v) -D1(0) √ Y -v - dv (5.115)
where we integrated by parts and used

1 0 1 √ 1 -x -1 dx x 3/2 = 2.
Using (5.115) in (5.114) we regain equation (5.77). We have thus shown the equivalence of the moment method and the perturbation expansion of this subsection. Note also that (5.107) and (5.108) apply also for a ≥ 0.

Transform inversion

We now use the results of subsections 5.1 and 5.2 to obtain an approximation to b(n, p). We use (2.3) and evaluate the Cauchy integral asymptotically. We assume that

Ω ≡ pn -3/2 = O(1), 0 < Ω < ∞ (5.116)
and with (5.1) we have

w -p-1 dw = 1 + a n 3/2 -1-Ωn 3/2 da n 3/2 = e -aΩ n 3/2 [1 + O(n -3/2 )]da. ( 5
(n, p) = 4 n n 3 1 2πi Bra 1 √ π + C (0) (a) + 1 √ n C (1) (a) + O(n -1 ) e -aΩ da ∼ 4 n n 3 1 2πi Bra 1 √ π -a D(Y ) e -aΩ da = 4 n n 3 (-a) 1 (2πi) 2 Bra d dY Br 4 s h ′ (s) h(s)
e sY ds e -aΩ da.

(5.118)

Here we take ℜ(a) < 0 on the vertical contour Br a , since we used the expression (5.37) for D(Y ).

To evaluate the double integral in (5.118) we note that as -a goes from -∞i to +∞i, Y = (-a) 2/3 goes from ∞e -πi/3 to ∞e πi/3 . We let

Y = Z 2
and the contour Br a is mapped to C 1 , that goes from ∞e -πi/6 to ∞e πi/6 . Then (5.118) becomes b(n, p) ∼

4 n n 3 1 (2πi) 2 C1 e ΩZ 3 3 2 Z 4 d dZ Br 4 s h ′ (s) h(s) e sZ 2 ds dZ = 4 n n 3 -1 (2πi) 2 ∞ω ∞ω 2 e -ΩZ 3 12Z 5 Br h ′ (s) h(s)
e sZ 2 ds dZ, (5.119)

where the last integral was obtained by changing Z to -Z, and ω = e 2πi/3 . We next write d(e -ΩZ 3 ) = -3ΩZ 2 e -ΩZ 3 and do some integrating by parts, noting that

Z 4 e -ΩZ 3 = 1 9Ω 2 ∂ 2 ∂Z 2 (e -ΩZ 3 ) + 2Z 3Ω e -ΩZ 3 .
Then (5.119) becomes b(n, p) ∼ -

4 n n 3 1 (2πi) 2 ∞ω ∞ω 2 e -ΩZ 3 Br h ′ (s) h(s) (Ze sZ 2 ) ZZ 1 9Ω 2 ds (5.120) + Br h ′ (s) h(s) e sZ 2 8Z 2 Ω ds dZ = - 4 n n 3 1 (2πi) 2 ∞ω ∞ω 2 Br e -ΩZ 3 e sZ 2 h ′ (s) h(s) 40 3
Zs Ω 2 + 16 9

s 2 (1 + 2sZ 2 ) Ω 3 ds dZ.
We furthermore set

Z = s 3Ω + W, W = ζ (3Ω) 1/3
(5.121) and note that the Airy function has the integral representation

1 2πi ∞ω ∞ω 2
e ζU e -ζ 3 /3 dζ = Ai(U ).

(5.122) Differentiating (5.122) with respect to U corresponds to multiplying the integrand by ζ. Hence (5.120) becomes b(n, p) ∼ -

4 n n 3 1 2πi Br h ′ (s) h(s) exp 2s 3 27Ω 2 1 (3Ω) 1/3 40 3Ω 2 s (3Ω) 1/3 Ai ′ (U (s)) + 40 9Ω 3 + 16 9Ω 3 s 2 Ai(U (s)) + 32 9Ω 3 s 3 (3Ω) 2/3 Ai ′′ (U (s)) + 64 27Ω 4 s 4 (3Ω) 1/3 Ai ′ (U (s)) + 32 81Ω 5 s 5 Ai(U (s)) ds (5.123)
where

U (s) = s 2 (3Ω) 4/3 .
(5.124)

Using Ai ′′ (U (s)) = U (s)Ai(U (s)) to simplify (5.123) and noting that the integrand has singularities at s = 4 1/3 r j , we evaluate the integral as a residue series and the result is precisely the series in (2.40).

It is also possible to obtain and O(n -1/2 ) correction to (5.123) or (2.40), by using (5.87) and (5.96) to compute C (1) (a) in (5.118) and evaluating the integral in a manner similar to that above. Finally we discuss the expansion of the right side of (2.40) (or (5.123)) for Ω → 0 and Ω → ∞. For Ω → 0 the term in the sum with j = 0 is exponentially larger than those with j ≥ 1, and thus we get

4 n n 3 1 (3Ω) 1/3 exp - 8 27 |r 0 | 3 Ω 2 64 81 |r 0 | 5 Ω 5 4 5/3 Ai r 2 0 4 2/3 3 4/3 Ω 4/3 ∼ 4 n n 3 exp - 16 27 |r 0 | 3 Ω 2 2 9 3 4 |r 0 | 9/2 √ πΩ 5 .
(5.125)

Here we also used the asymptotic expansion (5.56) of the Airy function. We will use expression (5.125) in Section 6, when we discuss the asymptotic matching of the results for p = O(n 3/2 ) and p = O(n 4/3 ).

The limit Ω → 0 corresponds to the left tail of the Airy distribution.

The right tail corresponds to Ω → ∞ and this limit is hard to evaluate from the sum in (2.40), as there is a lot of cancelation. However, we have already shown the asymptotic matching between the expansions for w = 1 + O(n -1 ) (β-scale) and w = 1 + O(n -3/2 ) (a-scale). This matching implies also the matching in (n, p) space of the expansion for p = O(n 2 ) (Λ-scale) and for p = O(n 3/2 ), where (5.123) applies. Thus the behavior of (5.123) as Ω → ∞ is the same as that of ( 4

β * = β * (Λ) ∼ 3 2 Λ, Λ → 0. (5.127) Then 1 -2Λ - 1 2e β * -1 ∼ Λ (5.128) and β * (1 -2Λ) + log 2 -log(2 -e -β * ) ∼ log 2 - 3 4 Λ 2 .
(5.129) Using (5.127) -(5.129) in (4.44) the right side becomes

2 n+1 √ 2 n 2 π β 5/2 * √ Λ exp n log 2 - 3 4 nΛ 2 ∼ 9 √ 3 2π 4 n n 2 Λ 2 exp - 3 4 nΛ 2 = 9 √ 3 2π 4 n n 3 Ω 2 exp - 3 4 Ω 2 , Ω = p n 3/2 . (5.130)
This is the behavior of (5.123) or (2.40) as Ω → ∞. We have thus completed analysis of the scales w -1 = O(n -3/2 ) and p = O(n 3/2 ).

Left Region

We consider the scale w -1 = O(n -1 ) with w < 1, defining γ by

w = 1 - γ n , γ > 0. (6.1)
We return to (2.1), and break up the sum over [0, n] into three pieces:

[0, L] ∪ [L + 1, n -L -1] ∪ [n - L -1, n]. Furthermore, we set B n (w) = 4 n n e ν0n 1/3 γ 2/3 e ν1γ log n F (γ; n). (6.2)
Here ν 0 and ν 1 are constants that will be determined shortly. The form in (6.2) is suggested by the behavior of the a-scale result as a → -∞ (cf. (5.8) with -a = γ √ n). Using (6.2) in (2.1) leads to

4 n+1 n + 1 exp (n + 1) 1/3 ν 0 1 + 1 n γ 2/3 + ν 1 γ log(n + 1) F γ 1 + 1 n ; n + 1 = 2 1 - γ n n L ℓ=0 B ℓ 1 - γ n 4 n-ℓ n -ℓ exp ν 0 (n -ℓ) 1/3 γ 1 - ℓ n 2/3 × exp ν 1 γ 1 - ℓ n log(n -ℓ) F 1 - ℓ n γ; n -ℓ + 1 - γ n n n-L-1 ℓ=L+1 4 n ℓ(n -ℓ) exp ν 0 ℓ 1/3 γ ℓ n 2/3 + ν 0 (n -ℓ) 1/3 γ 1 - ℓ n 2/3 × exp ν 1 γ ℓ n log ℓ + 1 - ℓ n log(n -ℓ) F ℓ n γ; ℓ F 1 - ℓ n γ; n -ℓ . (6.3)
Here we used the symmetry of the sum in (2.1). Also, in the second sum in (6.3) we used (6.2) to approximate both B ℓ (w) and B n-ℓ (w), while in the first sum we only used (6.2) to approximate the latter. This careful treatment of "edge effects" is necessary to derive the limiting equation. We note that

ℓ 1/3 ℓ n 2/3 + (n -ℓ) 1/3 1 - ℓ n 2/3 = n 1/3 and ℓ n log ℓ + 1 - ℓ n log(n -ℓ) = log n + H ℓ n where H(x) = x log x + (1 -x) log(1 -x). (6.4)
Now we assume that

F (γ; n) → F 0 (γ), n → ∞. (6.5)
Then choosing L so that L → ∞ and L/n → 0 as n → ∞, the first sum in (6.3) is asymptotic to

L ℓ=0 B ℓ 1 - γ n 4 n-ℓ n -ℓ exp(ν 0 n 1/3 γ 2/3 )n ν1γ F 0 (γ) ∼ 4 n n exp(ν 0 n 1/3 γ 2/3 )n ν1γ F 0 (γ) ∞ ℓ=0 B ℓ (1)4 -ℓ = 4 n n exp(ν 0 n 1/3 γ 2/3 )n ν1γ 2F 0 (γ). (6.6)
The second sum in (6.3) we approximate by the Euler-MacLaurin formula to obtain:

4 n n exp(ν 0 n 1/3 γ 2/3 )n ν1γ 1 0 e ν1γH(x) F 0 (γx)F 0 (γ -γx) x(1 -x) dx. (6.7)
Using (6.6) and (6.7) in ( 6.3) we obtain the limiting equation

4(e γ -1)F 0 (γ) = 1 0 e ν1γH(x) x(1 -x) F 0 (γx)F 0 (γ -γx)dx. (6.8)
Here we also used

(1 -γ/n) n ∼ e -γ .
This is a non-linear integral equation that is somewhat similar to one that arises in the study of the limiting distribution of the number of comparisons in the Quicksort algorithm [START_REF] Hennequin | Combinatorial Analysis of Quicksort Algorithm[END_REF][START_REF] Knessl | Quicksort algorithm again revisited[END_REF][START_REF] Régnier | A Limiting Distribution for Quicksort[END_REF][START_REF] Rösler | A Limit Theorem for Quicksort[END_REF]. Note also that we must have F 0 (0) = 0 for the integral to converge, but this will be shown to follow by asymptotic matching to the a-scale. Setting F 0 (γ) = 4γF 1 (γ), (6.9) equation (6.8) simplifies to x) dx. (6.10)

e γ -1 γ F 1 (γ) = 1 0 F 1 (γx)F 1 (γ -γx)e ν1γH(
Setting γ = 0 in (6.10) we conclude that F 1 (0) = 1. (6.11)

We examine more closely the behavior of (6.10) as γ → 0. We expand F 1 (•) as

F 1 (γ) = 1 + α 0 γ log γ + α 1 γ + o(γ).
Using (6.11) in (6.10) and expanding for γ small we obtain at O(γ)

α 1 + 1 2 = (α 0 + ν 1 ) 1 0 H(x)dx + α 1 = 2(α 0 + ν 1 ) 1 0 x log x dx + α 1 = - 1 2 (α 0 + ν 1 ) + α 1 (6.12) so that α 0 + ν 1 = -1 (6.13)
and α 1 is arbitrary. We can compute the higher-order terms in (6.11) in terms of α 1 , but the latter cannot be determined from (6.10) alone. This shows that the solution of (6.10) is not unique.

To uniquely determine F 1 (γ) we use asymptotic matching between the a and γ scales. For γ → 0 + we obtain, from (6.2), (6.5), (6.9) and (6.11),

4 n+1 n γ[1 + α 0 γ log γ + α 1 γ + o(γ)](1 + ν 1 γ log n) exp(ν 0 n 1/3 γ 2/3 ). (6.14)
The above should agree with the a-scale approximation as a → -∞, which we obtained in (5.99). Noting that γ = -a/ √ n and comparing (5.99) to (6.14) we conclude that

ν 0 = 4 1/3 r 0 = -4 1/3 |r 0 | (6.15)
and

ν 1 γ log n + α 0 γ log γ + α 1 γ = γ - 2 3 log(γ √ n) + κ 4 + 1 -γ E . (6.16)
where γ E is the Euler constant. From (6.16) we must have We note that (6.17) is consistent with (6.13). With (6.18) the solution to (6.10) is unique and may be computed, e.g., in the form of the series (6.11).

ν 1 = - 1 3 , α 0 = - 2 3 (6.
To further analyze (6.10) (with ν 1 = -1/3) we let

F 1 (γ) = exp 1 3 γ log γ F 2 (γ) (6.19)
and obtain

(e γ -1)F 2 (γ) = γ 0 F 2 (v)F 2 (γ -v)dv. (6.20)
Introducing the Laplace transform

F (θ) = ∞ 0 e -γθ F 2 (γ)dγ (6.21)
we obtain from (6.20) the discrete logistic equation (in reverse time):

F (θ -1) = F (θ) + F 2 (θ). (6.22)
The expansion (6.11) translates to

F 2 (γ) = 1 -γ log γ + α 1 γ + o(γ), γ → 0 +
and hence

F (θ) = 1 θ + log θ θ 2 + α 1 + γ E -1 θ 2 + o(θ -2 ), θ → +∞. (6.23)
We can refine (6.23) to (using (6.22))

F (θ) = 1 θ + log θ + α * θ 2 + G(θ) θ 3 + O R (θ -4 ) (6.24)
where

α * = α 1 + γ E -1 = 2.5464 . . . and
G(θ) = (log θ + α * ) 2 -log θ -α * + 3 2 . (6.25)
We shall see in Section 9 that knowing some of the higher-order terms in expansion (6.23) will allow for more efficient numerical calculation of F (θ), and hence F 2 (γ).

While we cannot solve (6.20) or (6.22) explicitly, we can guess the behavior of F 2 (γ) as γ → +∞, which will be needed for asymptotic matching purposes. This will correspond to knowing F (θ) as θ → -∞. Let us assume that F 2 (γ) in (6.20) behaves as We conclude that

F 2 (γ) ∼ e k1γ log γ e k2γ γ k3 k 4 , γ → ∞. ( 6
k 1 = - 1 log 2 , k 3 = - 1 2 , k 4 = 1 √ 2π log 2 . (6.28)
The constant k 2 remains arbitrary, though it is uniquely determined once we know α 1 in (6.11). The numerical studies in Section 9 will show that k 2 ≈ 3.696. To summarize, we have obtained

F 1 (γ) ∼ 1 √ 2π log 2 e k2γ √ γ exp 1 3 - 1 log 2 γ log γ γ → ∞. (6.29)
Using (6.29) and (6.9) in (6.2) yields the behavior of the γ-scale result for γ → ∞, and this will be used in Section 7 to asymptotically match to an approximation valid for 0 < w < 1 and n → ∞.

We also note that for θ ≡ -τ → -∞, (6.22) is consistent with an asymptotic expansion of the form

F (θ) ∼ A 2 τ 1 + ∞ L=1 b L A -L2 τ , τ = -θ → ∞. (6.30)
Using (6.30) in (6.22) leads to the recurrence

b L-1 + L i=0 b i b L-i =      0, L odd b L/2 , L even (6.31) so that b 1 = -1/2, b 2 = -1/8
, b 3 = 0, etc. Using (6.30) to asymptotically invert the transform in (6.21) leads to

F 2 (γ) = 1 2πi Br e γθ F (θ)dθ ∼ 1 2πi Br e γθ A 2 -θ dθ = 1 2πi C 1 log 2 e -γ log 2 u e u log A du u = 1 log 2 (log A) γ/ log 2 Γ(1 + γ/ log 2) ∼ 1 √ 2πγ 1 √ log 2 exp - γ log 2 log γ log 2 exp γ log 2 (1 + log log A) , γ → ∞. (6.32)
Here the contour C goes from -∞ + i0 -to -∞ + i0 + , encircling the branch cut along the negative real axis in the u-plane. By comparing (6.26) (with (6.28)) to (6.32) we conclude that

k 2 = log log 2 log 2 + 1 + log log A log 2 . ( 6

.33)

In Section 9 we shall numerically obtain A and then (6.33) will be used to obtain k 2 . We note that (6.30) implies that

2 θ log[F (θ)] = log A - 1 2 2 -τ A -2 τ (1 + o(1)), τ = -θ → ∞ (6.34)
so that the left side should approach a constant value super-exponentially fast as θ → -∞. Finally, we invert (2.3) to get b(n, p) using the result on the γ-scale. We scale

p = Θn 4/3 , Θ = O(1), Θ > 0 (6.35)
and use (6.1) so that

w -p = exp γp n + O p n 2 ∼ exp(γn 1/3 Θ). (6.36)
With (6.2), (6.5), (6.9) and (6.36), (2.3) becomes

1 2πi C B n (w)w -p-1 dw ∼ 4 n+1 n 1 2πi Br 4γn -γ/3 F 1 (γ) exp[n 1/3 (γΘ -|r 0 |4 1/3 γ 2/3 )] dγ n . (6.37)
For n → ∞ we evaluate the last integral by the saddle point method, noting that there is a saddle point along the real axis, at

d dγ (γΘ -|r 0 |4 1/3 γ 2/3 ) = 0 ⇒ γ = γ * ≡ 32 27 |r 0 | 3 Θ 3 . (6.38)
Then the standard Laplace estimate of (6.37) yields b(n, p) ∼ 4 n n 13/6 |r 0 | 9/2 Θ 5 512 81

n -γ * /3 √ π F 1 (γ * ) exp - 16 27 n 1/3 |r 0 | 3 Θ 2 (6.39)
which applies for Θ = pn -4/3 fixed as n → ∞. Here F 1 (•) must be obtained by numerically solving (6.10) or (6.20). We have thus established (2.41).

Since F 1 (0) = 1, expanding (6.39) as Θ → ∞ corresponds simply to dropping the factors n -γ * /3 and F 1 (γ * ). Thus (6.39) asymptotically matches to the result valid for p = O(n 3/2 ), in view of (5.125) and the fact that Ω = Θn -1/6 . We can also infer the behavior of the right side of (6.39) as Θ → 0, which corresponds to γ * → ∞ and p = o(n 3/2 ). Using (6.29) we obtain

4 n n 13/6 512 81π 1 √ 2 log 2 |r 0 | 9/2 Θ 5 Θ 3/2 |r 0 | 3/2 3 √ 3 4 √ 2 exp - 16 27 n 1/3 |r 0 | 3 Θ 2 × exp - 1 3 γ * log n + 1 3 - 1 log 2 γ * log γ * + k 2 γ * = 1 n 13/6 |r 0 | 3 Θ 7/2 64 9 √ 3 1 π √ log 2 exp n log 4 -4 1/3 9 16 |r 0 |n 1/3 γ 2/3 * - 1 3 γ * log n + 1 3 - 1 log 2 γ * log γ * + k 2 γ * . (6.40)
We will use (6.40) for asymptotic matching purposes in Sections 7 and 8.

Far Left Region

We consider (2.1) for 0 < w < 1 and n → ∞, and

(2.3) for n → ∞ with p -n log 2 n = O(n). Note that p -n log 2 n = O(n) contains p min (n) in (2.
13) (see also (7.29)). We assume first that B n (w) has an expansion of the form

B n (w) = e -n log nf (w) e ng(w) n h(w) q(w)[1 + o(1)], (7.1)
for 0 < w < 1 and n → ∞. If f (w) > 0 this means that B n (w) will decay faster than exponentially as

n → ∞. Since (n + 1) log(n + 1) = n log n + log n + 1 + O(n -1 )
we see that with (7.1) the left side of (2.1) becomes e -n log nf (w) e ng(w) n -f (w)+h(w) q(w)e g(w)-f (w

) [1 + o(1)]. (7.2)
To evaluate the right side we treat the sum as an implicit Laplace-type integral. Using (7.1) the sum becomes asymptotic to

e n log w n ℓ=0 ℓ h(w) (n -ℓ) h(w) q 2 (w)e ng(w) exp{-[ℓ log ℓ + (n -ℓ) log(n -ℓ)]f (w)}. (7.3)
We approximate the sum via Euler-MacLaurin to get e n log w n 2h(w)+1 q 2 (w)e ng(w) e -n log nf (w)

1 0 x h(w) (1 -x) h(w) e -nf (w)H(x) dx (7.4)
where H(x) = x log x+(1-x) log(1-x) is as in (6.4). Assuming for now that f (w) > 0 we can evaluate (7.4) as a Laplace integral for n → ∞, with the major contribution coming from

x = 1 2 + O(n -1/2 ). Using H 1 2 = -log 2, H ′ 1 2 = 0, H ′′ 1 2 = 4, (7.5)
(7.4) becomes asymptotic to

e ng(w) e -n log nf (w) q 2 (w)n 2h(w)+1 π 2f (w)n 1 2 2h(w) exp[n(f (w) log 2 + log w)]. (7.6) Comparing (7.1) to (7.6)
we conclude that

f (w) = - log w log 2 = -log 2 w > 0, (7.7)
h(w) = log 2 w - 1 2 , (7.8)
and

q(w) = w 2+ 1 log 2 -log 2 w 1 √ 2π e g(w) . (7.9)
However, the function g(w) remains undetermined. To summarize, assuming the ansatz (7.1) we have obtained

B n (w) = w n log 2 n e (n+1)g(w) n log 2 w -log 2 w 2πn w 2+ 1 log 2 [1 + o(1)]. (7.10)
The numerical studies in Section 9 suggest that (7.10) is approximately correct for w > 0, but there are oscillations present that become numerically significant when w becomes small. We thus re-examine (2.1) with the more general ansatz

B n (w) = e -n log nf (w) e nB(w,n) n h(w) Q(w, n)[1 + o(1)]. (7.11)
We allow B(w, n) and Q(w, n) to depend weakly upon n, in such a way that B(w, n + 1) ∼ B(w, n) and Q(w, n + 1) ∼ Q(w, n) as n → ∞. Repeating the calculation we did with (7.1) we now find that for n → ∞ the left side of (2.1) is e -n log nf (w) e (n+1)B(w,n+1) n h(w)-f (w) Q(w, n + 1)e -f (w) . (7.12)

The sum in the right side of (2.1) becomes

e n log w n ℓ=0 ℓ h(w) (n -ℓ) h(w) Q(w, ℓ)Q(w, n -ℓ)e ℓB(w,ℓ)+(n-ℓ)B(w,n-ℓ) × exp{-[ℓ log ℓ + (n -ℓ) log(n -ℓ)]f (w)}. (7.13)
Again evaluating this by the Laplace method, we use

ℓB(w, ℓ) + (n -ℓ)B(w, n -ℓ) = nB w, n 2 + 1 n ℓ - n 2 2 n 2 2 B nn w, n 2 + 2nB n w, n 2 + . . . , (7.14)
which follows upon expanding in Taylor series about ℓ = n/2. Then (7.13) becomes, for n → ∞, e n log w e -n log nf (w) e n(log 2)f (w) n 1+2h(w) 1 2

2h(w) × e nB(w, n 2 ) Q 2 w, n 2
π 2n f (w) -nB n w, n 2 - n 2 4 B nn w, n 2 -1/2
. Comparing (7.12) to (7.14) we again obtain (7.7). Then using

(n + 1)B(w, n + 1) = nB(w, n) + B(w, n) + nB n (w, n) + • • • we find that B(w, n) satisfies B(w, n) = B w, n 2 . (7.15)
The most general solution of (7.15) is

B(w, n) = g(w) + ∞ k=-∞ k =0 g k (w)e 2πi(log 2 n)k , (7.16)
where the sum represents an arbitrary, zero-mean periodic function of log 2 n, with period one. Note that (7.15) and (7.16) are indeed consistent with our assumption of slow variation. Unfortunately, we cannot determine explicitly the Fourier coefficients g k (w) using purely the recurrence (2.1). With (7.16) we write

nB n w, n 2 = nB n (w, n) = B * 1 (w, n), n 2 B nn w, n 2 = n 2 B nn (w, n) = B * 2 (w, n)
where

B * 1 (w, n) = 2πi log 2 ∞ k=-∞ kg k (w)e 2πi(log 2 n)k B * 2 (w, n) = 2πi log 2 ∞ k=-∞ 2πi log 2 k 2 -k g k (w)e 2πi(log 2 n)k . (7.17)
Returning to (7.12) and (7.14) we regain (7.8) and instead of (7.9) we obtain the equation

Q 2 w, n 2 = Q(w, n)w 2+ 1 log 2 √ 2πe g(w) × e B * 0 (w,n)+B * 1 (w,n) -log 2 w -B * 1 (w, n) - 1 4 B * 2 (w, n). (7.18)
Here we wrote (7.16) as B(w, n) = g(w) + B * 0 (w, n). A particular solution to (7.18) is

Q p (w, n) = w 2+ 1 log 2 1 √ 2π e g+B * 0 +B * 1 -log 2 w -B * 1 - 1 4 B * 2 (7.19)
and setting

Q(w, n) = Q p (w, n) Q(w, n) we find that Q satisfies Q2 w, n 2 = Q(w, n). (7.20)
Setting

G(w, n) = log[ Q(w, n)]
we find that

2 G w, n 2 = G(w, n) (7.21)
whose most general solution is

G(w, n) = n × [periodic function of log 2 n]. (7.22)
However, then Q(w, n) can be incorporated into the factor exp[nB(w, n)] in (7.11). Thus in (7.11) we can use (7.7), (7.8), (7.16) and (7.19). We have thus established (2.34) and (2.35).

We next examine the asymptotic matching between the results for 0 < w < 1 and w = 1 -O(n -1 ). For γ → ∞ we expand (6.2) to get

4 n n+1 exp[4 1/3 r 0 n 1/3 γ 2/3 - 1 3 γ log n] e k2γ √ 2π log 2 √ γ exp 1 3 - 1 log 2 γ log γ (7.23)
where we also used (6.5), (6.9) and (6.29). By asymptotic matching, (7.23) should agree with the expansion of (7.11) as w → 1. Since (7.23) contains no oscillatory terms, we conclude that nB * 0 (w, n) → 0 as w → 1 and thus g k (w) = o(1w) as w → 1, for each k. Then (7.23) must match to the right side of (7.10) (which neglected the oscillations), as the latter is expanded for w → 1. Noting that

w n log 2 n = 1 - γ n n log 2 n ∼ exp - γ log 2 log n and n log 2 w ∼ 1, (7.10) becomes 1 √ 2π log 2 1 -w n exp - γ log 2
log n e (n+1)g(1-γ/n) .

Thus, the matching is possible provided that

g(w) = log 4 + 4 1/3 r 0 (1 -w) 2/3 + 1 log 2 - 1 3 (w -1) log(1 -w) + k 2 (1 -w) + o(1 -w), w ↑ 1. (7.24)
We have thus used asymptotic matching to infer the behavior of g(w) as w → 1, and this will play a major role in Section 8. Finally we discuss briefly the limit w → 0 with n fixed. From the discussion in Section 3, each B n (w) is a polynomial of the form

B n (w) = C n w pmin(n) + • • • + 2 n-1 w ( n 2 ) (7.25)
where

p min (n) = n J=2 ⌊log 2 J⌋ = (n + 1)⌊log 2 n⌋ -2(2 ⌊log 2 n⌋ -1).
By using (7.25) in (2.1) we find that the dominant terms for w → 0 come from near the midpoint of the sum ℓ ≈ n/2. More precisely, for n even we let n = 2N and obtain

C 2N +1 = j∈S(N ) C N +j C N -j (7.26)
where

S(N ) = {ℓ : p min (N + ℓ) + p min (N -ℓ) = 2p min (N )}.
For n odd we set n = 2N + 1 and get

C 2N +2 = j∈T (N ) C N +j C N +1-j (7.27)
where

T (N ) = {ℓ : p min (N + ℓ) + p min (N + 1 -ℓ) = p min (N ) + p min (N + 1)}.
Note that 0 ∈ S(N ) for all N , and 0, 1 ∈ T (N ). If ℓ ∈ S(N ) then -ℓ ∈ S(N ), and if ℓ ∈ T (N ) then 1ℓ ∈ T (N ). For example, we have S(1) = {0}, S(2) = {-1, 0, 1}, S(3) = {0}, T (0) = {0, 1}, T (1) = {0, 1}, T (2) = {0, 1}. Thus for n fixed we have

B n (w) ∼ C n w pmin(n) , w → 0. (7.28)
We next show that (7.28) cannot be a special case of (7.14), as the latter is expanded for w → 0. Using the Fourier series'

{x} = 1 2 - ∞ k=-∞ k =0 e 2πikx 2πik = 1 2 - ∞ k=1 sin(2πkx) πk , 2 -{x} = 1 2 ∞ k=-∞ 1 log 2 + 2kπi e 2πikx ,
where the first sum omits the term k = 0, we can represent p min (n) as

p min (n) = n log 2 n + n - 1 2 - 1 log 2 = n ∞ k=-∞ k =0 log 2 2πik(log 2 + 2πik) e 2πik log 2 n = log 2 n + 3 2 + ∞ k=-∞ k =0 e 2πik log 2 n 2πik . (7.29)
This form yields the asymptotic behavior as n → ∞. Now, if (7.14) were to contain (7.28) as a special case (at least for n → ∞), it would need to behave as C asy n exp(p min (n) log w) for w → 0, where C asy n is understood as the asymptotic behavior of C n as n → ∞. By comparing (7.29) to (2.34) we see that the largest factors, i.e., exp(n log 2 n log w), match automatically, and the factors that are exp(O(n)) agree if

g(w) ∼ - 1 2 - 1 log 2 log w, w → 0, g k (w) ∼ log 2 2πik(log 2 + 2πik) log w, w → 0. (7.30)
With (7.30) we then have

e g(w) w 2+ 1 log 2 ∼ w 3/2 , w → 0,
and the term log 2 n in p min (n) corresponds to the factor n log 2 w = w log 2 n in (2.34). Furthermore

exp[B * 0 (w, n) + B * 1 (w, n)] = exp      ∞ k=-∞ k =0 1 + 2πik log 2 g k (k)e 2πik log 2 n      ∼ exp      log w    ∞ k=-∞ k =0 e 2πik log 2 n 2πik         , w → 0
which corresponds to the last sum in (7.29). However, (2.34) still has the factor

1 √ 2πn -log 2 w -B * 1 (w, n) - 1 4 B * 2 (w, n)
which for w → 0 would look likelog 2 w times a function of n, and such a factor can have no analog in (7.28), regardless of the form of C n for n → ∞. This shows that yet another scale is needed in order to completely understand the asymptotic behavior of B n (w). This scale would have w → 0 and n → ∞ simultaneously. We have not been able to identify this scale precisely nor to analyze it. It is, however, not important for understanding the asymptotic for n b(n, p), as this requires w → 1, as shown in Section 8. Finally, we use the form (7.10) to obtain an approximation to b(n, p). We set For n → ∞ we can use the saddle point method. There is a saddle in the range w = O(1) if we can solve the equation g ′ (w) = α/w. The numerical studies of g(w) in Section 9 suggest this is indeed possible, and then the standard saddle point approximation leads to (2.42).

p = n log 2 n + αn, α = O(1)

The Matching Region Between the Left and Far Left Scales

In Section 7 we showed that the expansions for 0 < w <

1 (p = n log 2 n + O(n)) and w = 1 -O(n -1 ) (p = O(n 4/3
)) can match in some intermediate limit. Here we examine in more detail the solution in this matching region, as this is the key to understanding the distribution of the number of nodes in trees of a given (large) total path length p (cf. (2.12)). For a fixed p, as we increase n we move from right to left. The result in (6.39) shows that for a fixed Θ = pn -4/3 , b(n, p) still grows with n, due to the dominant exponential factor of 4 n . However, for p = n log 2 n + O(n), (2.42) shows that for fixed p/nlog 2 n, b(n, p) grows (decays) with n according as g(w * )α log w * > 0 (< 0). However, numerical studies (cf. Section 9) show that this quantity is negative unless w * → 1. Thus to find the limiting distribution of the number of nodes we need the maximum of b(n, p) over n, and this occurs precisely in the matching region, where this function will achieve a Gaussian peak.

As (pn log 2 n)/n = α → ∞ we have w * → 1 and the non-constant Fourier coefficients in (2.34) (i.e., g k (w), k = 0) vanish as w → 1. Thus we use (2.37) and (2.42). For α → ∞ we can solve the equation

wg ′ (w) = α (8.1)
for the saddle w = w * (α) asymptotically, using the relations

g(w) = log 4 -4 1/3 |r 0 |(1 -w) 2/3 + 1 log 2 - 1 3 (w -1) log(1 -w) -k 2 (w -1) + o(w -1), (8.2)
g ′ (w) = 4 1/3 2 3 |r 0 | 1 1 -w 1/3 1 log 2 - 1 3 log(1 -w) + - 1 3 + 1 log 2 -k 2 + o(1), (8.3)
and

g ′′ (w) ∼ 2 9 4 1/3 |r 0 |(1 -w) -4/3 , w → 1 -. (8.4)
To facilitate the calculation, we set

A 0 = 2 3 4 1/3 |r 0 |, A 1 = 1 log 2 - 1 3 , A 2 = - 1 3 + log 2 -k 2 = A 1 -k 2 (8.5)
and

w * = 1 -∆. (8.6)
We rewrite (8.1) as

α = [A 0 ∆ -1/3 + A 1 log ∆ + A 2 + o(1)](1 -∆), ∆ → 0 (8.7)
which can be rearranged to give

∆ A 3 0 = [α -A 1 log ∆ -A 2 + α∆ + o(1)] -3 = α -3 1 - A 1 log ∆ + A 2 α + ∆ + o 1 α -3 = α -3 1 + 3 α (A 1 log ∆ + A 2 ) + o(α -1 ) . (8.8)
It follows that

∆ ∼ A 3 0 α 3 , α → ∞ (8.9)
and this can be refined to the expansion

∆ = A 0 α 3 1 + δ 1 α + δ 2 α 2 + • • • , α → ∞ (8.10)
where

δ 1 = 3A 2 -9A 1 log α A 0 . (8.11)
We next evaluate g(w * )α log w * as α → ∞ and w * → 1. We note that this is the exponential growth rate, in n, of (2.42). Using (8.6), (8.10) and (8.11) we obtain

g(w * ) -α log w * = g(1 -∆) -α log(1 -∆) = log 4 + 1 2 A 1 ∆ log ∆ + 3 2 A 2 + k 2 ∆ - 1 2 α∆ + O(α∆ 2 ) = log 4 - A 3 0 2 1 α 2 + A 3 0 α 3 3A 1 log α A 0 + k 2 + O R (α -4 ). (8.12)
Here we also used (8.7) in the form α∆ We will see that Ψ 0 (p), ñ(p) and V 0 (p) will provide asymptotic approximations to the exponential growth rate of the total number of trees of path length p, and the mean and variance of the Gaussian distribution of the number of nodes in such trees. We use (8.2) and (8.10) in (8.13), which as w * → 1 yields log 4 -log n + 1 log 2

∼ A 0 ∆ 2/3 + A 1 ∆ log ∆ + A 2 ∆.
A 3 0 α 3 + 3A 3 0 α 4 A 2 -3A 1 log α A 0 - 3 2 A 0 A 2 0 α 2 1 + 2 α A 2 -3A 1 log α A 0 + k 2 A 3 0 α 3 -A 1 A 3 0 α 3 3 log A 0 α + o(α -3 ) = 0.
To leading order, log 4 must be balanced by the term proportional to α -3 log n, and we have

α ∼ A 0 (log 2) 1/3 (log n) 1/3 (log 4) 1/3 = A 0 (log n) 1/3 2 1/3 (log 2) 2/3 , n = ñ(p). (8.19)
This shows that for a fixed p, the maximum of b(n, p) occurs in the range p = n log 2 n + O[n(log n) where we approximated the sum by the Laplace method. Now,

∆ ∼ 2(log 2) 2 log n ∼ 1 -w, n = ñ(p) → ∞. ( 8
n log 2 w = exp[log 2 n log(1 -∆)] ∼ exp[-2 log 2] = 1 4 , e g( w) ∼ exp[log 4] = 4 and α + w2 g ′′ ( w) ∼ α + 1 3 A 0 ( ∆) -4/3 ∼ 1 3
A 0 2 -4/3 (log 2) -8/3 (log ñ) we obtain

0 = log 4 - S log 2 ∆ - S 2 log 2 ∆ 2 - 3 2 A 0 ∆ 2/3 -A 1 ∆ log ∆ + k 2 ∆ - 1 log 2 ∆ + O R (S -2 ) (8.27)
where we also used (8.2). From (8.27) it follows that as S → ∞, ∆ has the expansion

∆ = a S + b S 4/3 + c S 5/3 + d S 2 + O R (S -7/3 ) (8.28)
where a = 2(log 2) 2 , b = 0, (8.29)

c = - 3 2 (log 2)A 0 a 2/3 = - 3 2 2 2/3 A 0 (log 2) 7/3 , (8.30) d = log 2 k 2 - 1 log 2 a - a 2 2 log 2 + A 1 a log S a (8.31)
and we note that d depends weakly on S (as log S).

Next we use (8.3) and rewrite (8.1) as

p ñ -log 2 ñ = (1 -∆)[A 0 ∆ -1/3 + A 1 log ∆ + A 2 + o(1)]. (8.32)
We set ñ(p) = pn * (p), Q = log p (8.33)

with which (8.32) becomes

1 n * - S log 2 = A 0 S a 1/3 -A 0 c 3a 4/3 1 S 1/3 -A 0 d 3a 4/3 1 S 2/3 + A 1 -log S a + c a 1 S 2/3 + A 2 -A 0 a 2/3 S 2/3 + O R (S -1 ) (8.34)
where we expanded ∆ in (8.32) using (8.28). Noting that

S = Q + log n * ∼ Q and S 1/3 ∼ Q 1/3 + 1 3 Q -2/3 log n *
we rewrite (8.34) as

1 n * = Q log 2 + log n * log 2 + A 0 a 1/3 Q 1/3 + A 2 + A 1 log a -A 1 log Q + O R (Q -1/3 ). (8.35)
It follows that n * has the expansion

n * = log 2 Q + ν Q 5/3 + ν ′ Q 2 + O R (Q -7/3 ) (8.36) where ν = - A 0 a 1/3 (log 2) 2 and ν ′ = log 2[log Q -log log 2] + A 1 log Q(log 2) 2 -(A 2 + A 1 log a)(log 2) 2 .
Hence, as p → ∞ (with Q = log p) we have

ñ(p) = p Q log 2 1 - A 0 a 1/3 log 2 Q 2/3 + Z Q + O R (Q -4/3 ) , Z = (A 1 log 2 + 1) log Q -log log 2 -(A 2 + A 1 log a) log 2. (8.37)
Recalling that n = ñ(p) corresponds to the maximum of b(n, p) for a fixed large p, (8.37) is also the expansion of N (p), so we have established (2.18). We also note that

n max (p) -ñ(p) ∼ A 0 (log 2) 2 a 1/3 p (log p) 5/3 = o(p), p → ∞.
Here n max (p) is the inverse of p min (n) in (2.13).

Next we refine (8.25) and also obtain the variance in (8.18). From (8.17), (8.24) and (8.28) we have

Ψ 0 (p) = p + ñ log 2 a S + c S 5/3 + d + a 2 2 1 S 2 + O R (S -7/3 ) = p + p log p + O p Q 5/3 a Q + c Q 5/3 + d + a 2 2 + a(log Q -log log 2) 1 Q 2 + O R (Q -7/3 ) = p Q a + c Q 2/3 + a + a 2 2 + d -a log log 2 + a log Q 1 Q + O R (Q 4/3
= Q log 2 + O(Q 1/3 ) 2 = Q 2 (log 2) 2 + O(Q 4/3 ).
From (8.19) and (8.37) we obtain

α ∼ A 0 2 1/3 (log 2) 2/3 [log p + log log 2 -log Q + O(Q -2/3 )] 1/3 ∼ A 0 2 1/3 (log 2) 2/3 Q 1/3 .
Using the above in (8.39) we obtain

F nn (ñ(p); p) ∼ - 1 ñ(p) Q 2 (log 2) 2 3 A 0 2 4/3 (log 2) 8/3 Q -4/3 ∼ - Q 5/3 p 6a 1/3 A 0 log 2 . ( 8

.40)

We have thus obtained the leading term in the variance in (2.19). We will obtain the correction term using a somewhat different method (cf. (8.46) and (8.47)).

Using (8.40) and (8.37) in (8.23) leads to (2.14). To summarize, in the limit α → ∞ and pn -4/3 → 0 we have obtained the approximation b(n, p) ∼ 3 π

2 1/6 √ log 2|r 0 | 1/2 1 n A 0 α 7/2 exp - A 3 0 α 3 log n log 2 × 4 n exp - A 3 0 n 2α 2 + A 3 0 n α 3 (k 2 + 3A 1 log α -3A 1 log A 0 ) . ( 8
H(α) = - A 3 0 2α 2 + A 3 0 α 3 (k 2 + 3A 1 log α -3A 1 log A 0 ) Φ(n, p) ≡ n log 4 + nH(α). (8.42)
Then ∂Φ/∂n = 0 is equivalent to solving

log 4 + H(α) - 1 log 2 H ′ (α) - p n H ′ (α) = 0. (8.43)
We denote the solution to (8.43) by

n = n(p), α(p) = p n(p) -log 2 [n(p)]
where we are viewing the problem for a fixed p. We also have, since α = p/nlog 2 n,

∂ 2 Φ ∂n 2 = - 1 n log 2 H ′ (α) + 1 n p n + 1 log 2 2 H ′′ (α)
and we get Φ = Φ(n(p), p) and

V(p) = - 1 Φ nn (n(p), p) .
Solving for n asymptotically as p → ∞ regains the expansion of the mean in (8.37), though using the more implicit form (8.43) may be advantageous numerically (cf. Section 9). Thus we can use ( Φ, n, V) as the approximations to the growth rate, mean and variance, and these are asymptotically equivalent to the results in (2.14) -(2.19). For example, solving (8.43) for p → ∞ we have

α = A 0 Q a 1/3 + A 2 + A 1 log a -A 1 log Q + O R (Q -1/3 ) (8.44)
and we note that ). Note that while the leading terms in the expansion of the growth rate (log n b(n, p)) and the mean do not involve the root r 0 of the Airy function, the leading term for the variance does involve this root. Also, the correction terms for the growth rate and mean are smaller than the leading terms by factors of Q -2/3 = (log p) -2/3 , but the correction term to the variance is smaller by only an O(Q -1/3 ) = O((log p) -1/3 ) factor. This completes the proof of Result 2.

H ′′ (α) = - 3A 3 0 α4 + A 3 0 α5 36A 1 log α A 0 -12A 2 -9A 1 . ( 8
H ′′ (α) = - 3 A 0 a Q 4/3 - 9A 1 A 2 0 a Q 5/3 + O R (Q -2 ). ( 8

Numerical Studies

We give numerical results that determine some unknown constants/functions, and also provide back-up to the various assumptions we made.

For w > 1 we gave in Table 1, numerical values of

B * (w) (∼ B n (w)2 1-n w -( n 2 )
, n → ∞) for various w, and sketched this function in Figure 1 (cf. Section 2). We test the accuracy of the asymptotic relation in (2.23). Let Bn (w) be the exact solution to (3.2). Then our analysis predicts that Bn (w) → B * (w) as n → ∞ for each w > 1, and then as w ↓ 1. In Table 2, we compute (w -1) log[ B∞ (w)] and also (w -1) log[ B∞ (w)] -1 2 log(w -1) for various w > 1. Here, for each fixed w, B∞ (w) is computed by iterating (3.2) for n sufficiently large until Bn (w) settles to a constant to 3 digits. Both of the tabulated functions should approach d 0 as w ↓ 1, with the latter function approaching more rapidly, since it includes information from the algebraic factor √ w -1 in (2.23). We see that both functions are indeed approaching d 0 as w → 1, with the second function having a faster rate of approach. Next we attempt to verify the constant d 1 in the asymptotic relation (2.23), by computing We next consider the β-scale, where w = 1 + β/n = 1 + O(n -1 ) with β > 0. Here the asymptotic result is given by (2.25). In Tables 4 and5, we compare

D 1 (w) ≡ exp - d 0 w -1 B∞ (w) √ w -1 in
1 n log B n 1 + β n (9.2) to Φ(β) = log 2 + β 2 + φ(β) (cf. (2. 25
)) for β = .25, .5, 1, 2 and 4, respectively, and for various n. Our WKB expansion predicts that (9.2) should approach Φ(β) as n → ∞. The tables clearly demonstrate this convergence. The data are also consistent with an O(n -1 log n) error term, which is predicted by (2.25). The convergence is slowest when β = .25 (cf. Table 4), which is to be expected since our analysis predicts that once β becomes O(n -1/2 ), the expansion becomes invalid and we must use the a-scale result.

In Table 6, we compare

√ ne -nΦ(β) B n 1 + β n (9.3)
to √ βĝ(β), which should be the limit of (9.3) as n → ∞. Table 6 has β = .5, 1, and 2, respectively. Again the tables show the convergence predicted by the WKB ansatz (4.9). However, the convergence is much slower than one would expect with the O(n -1 ) error term in (4.9). The numerical data suggest that the error term may be O(n -1/2 ) (i.e., a term n -1/2 g (1/2) (β) should be included in the series in (4.9)). To resolve this issue more conclusively, we would have to examine higher-order asymptotic matching between the a-scale and β-scale. This would necessitate, among other things, continuing (2.28) into the range a > 0 and the evaluating the result for a → +∞.

We next consider the a-scale result. We recall that obtaining the right tail of the Airy distribution involved showing that

D(y) = D(a 2/3 ) ∼ 4 √ π ae a 2 /3 , a → +∞. (9.4)
In Table 7 we give D(a 2/3 )a -1 e -a 2 /3 for various values of a ≥ 1. This clearly shows that the function is approaching the constant 4/ √ π = 2.2567 . . . . To compute the exact values of D(y) we used the recurrence (5.25), (5.23) and (5.20). Now consider the expansion on the γ-scale, where w = 1 -O(n -1 ). Here we must solve (6.20) or, equivalently, (6.22). Given the asymptotic behavior of F (θ) in (6.24) as θ → ∞, we used the following numerical scheme: (i) Fix a large N , (ii) let F (J) satisfy the recurrence

F (J -1) = F (J) + F 2 (J), J = N, N -1, . . . , -M + 1, (9.5)
subject to the "terminal" condition

F (N ) = 1 N + log N + α * N 2 + (log N + α * ) 2 -log N -α * + 3 2 N 3 , (iii) use the approximation F (-M ) ≈ F (-M ).
Tab. 4: For M large and positive, (6.30) shows that

β = .25 n 1 n log(B n ) Φ(β) 10
2 -M log[F (-M )] = log A - 1 2 2 -M A -2 M + • • • , M → ∞. (9.6)
Thus the left side of (9.5) should approach the constant log A very rapidly (as a double exponential) as M → ∞, if A > 1. We have found that taking M = 5 is sufficient to get log A to several decimal places.

In Table 8 we take M = 5, solve (9.5) for various N and give 2 -5 F (-5), which is an approximation to log A that should become exact as N → ∞. This shows that log A ≈ 6.880 (9.7) and then (6.33) yields k 2 ≈ 3.696. (9.8)

The latter constant corresponds to the exponential growth rate of F 1 (γ) in (6.29), and thus that of the approximation on the γ-scale in (6.2) or (2.30) -(2.33).

Next we consider n → ∞ with 0 < w < 1, where the approximation (2.37) applies (or its refinement (2.34), that allows for oscillations). If we were to ignore the oscillations and define g (n) num (w) ≡ num (w) → g(w) as n → ∞ for each fixed 0 < w < 1. Of course we should also have n -1 log[B n (w)w -n log 2 n ] → g(w), but including the other factors in the right side of (9.9) should help with the convergence, as shown by the asymptotic analysis of Section 7. In Table 9, we plot g (n) num (w) large n, for several w < 1. These data show that for n sufficiently large, g (n) num (w) is constant to several decimal places, except when w ≤ .001, where this function oscillates (with n) over a certain range. This range is indicated in Table 9 for w = .001, .0001 and .00001. Our study suggests that while the oscillations corresponding to g k (w) in (2.35) are indeed present, they are numerically very small, unless w itself is very small. We also recall that the analysis in Section 7 predicted that g k (w) has a logarithmic singularity, for each k, as w → 0 (cf. (7.30)). For practical purposes we can use (2.37) and (2.42) as the approximations to B n (w) and b(n, p), though g(w) must still be determined numerically, using, e.g., Table 9. We also note from the table that g(w) is approaching, as w → 1, the theoretical value g(1) = log 4 = 1.3862 . . . .

In Table 10 we examine the matching region between the scales w = 1 -O(n -1 ) and 0 < w < 1. By asymptotic matching we showed that g(w) behaves as (8. For w not small, we can choose n sufficiently large to make (9.10) independent of n to several decimal places. Then as w → 1 our analysis shows that (9.9) should approach the constant k 2 . The convergence becomes very slow with n when w is only slightly less than one, but the data shows that CONST is decreasing to some value < 3.792, and this is consistent with our theoretical value in (9.8). This provides numerical evidence both for the form (8.2), as well as the numerical value of k 2 . Finally we do some numerical studies of our main results, which consist of the growth rate of the total number of trees of a given path length, and the distribution of these trees as characterized by their numbers of nodes. In Figure 2, we plot the exact b(n, p) for p = 190 and p = 250 with n ∈ [n min (p), n max (p)] = [START_REF] Régnier | A Limiting Distribution for Quicksort[END_REF]49] and [START_REF] Seroussi | On the Number of t-ary Trees with a Given Pathlength[END_REF]61]. While this graph possibly resembles a Gaussian near the peak, the value of p is too small for us to be certain. In Table 11 we compare the exact values of log[ n b(n, p)] to the one-term, two-term and three-term asymptotic approximations to the exponent in (2.14). The one-term approximation is ap/Q = 2(log 2) 2 p/ log p, the two-term includes the O(Q -2/3 ) correction, and the three-term includes also the O R (Q -1 ) correction. Table 11 shows that the one-term approximation underestimates the true value by about 20%, the two-term approximation is very poor, while the three-term approximation is quite accurate. In Table 12, we give the exact mean total path length N (p) (cf. the one-two-and three-term asymptotic approximations in (2.18). We again see that the two-term approximation is actually worse than the one-term, and the three-term result is better than the one-term. In Table 13, we give the exact variance V(p) along with the one-and two-term approximations in (2.19). We see that the one-term approximation is fairly accurate (10% errors), while the two-term underestimates significantly the true value.

To explain these numerical trends we note that the ratios of coefficients of the second term to the first term are respectively given by: growth rate: R 0 ≡ All of these numbers exceed unity. To compensate for the ratio being large we need to have Q -2/3 (growth rate, mean) or Q -1/3 (variance) sufficiently small. However, when p = 250, Q 1/3 = (log p) 1/3 = 1.767 . . . and Q 2/3 = 3.123 . . .; increasing p to 500 gives Q 1/3 = 1.838 . . . and Q 2/3 = 3.380 . . . . Thus to make the second term in the growth rate numerically smaller than the first, we need p ≈ 67.3. Indeed, Table 11 shows that the two-term approximation is actually negative for p = 20 -50. To make the second term one half as large as the leading term, we would need Q 2/3 = 2R 0 or p = exp[(2R 0 ) 3/2 ] ≈ 148, 000. This is well outside the range of Table 11, where the maximum value of p is 250. Thus, the fact that the two-term approximations tend to actually be worse than the one-term ones is not surprising, in view of the relatively large values of R j and the fact that p must be extremely large before (log p) -1/3 or (log p) -2/3 is sufficiently small. The numerical results are not particularly bad, considering that the asymptotic series involves powers of Q -1/3 . We have also tested the accuracy of the approximations ( Φ, n, V), defined by (8.43) and below. However, we found these to be worse than the (more explicit) results in (2.15), (2.18) and (2.19). We would guess that having an explicit expression for g(w) (rather than only the local behavior in (8.2)) would lead to more accurate approximations. However, refining (8.2) (i.e., computing explicitly the O(1w) error term) would seem to require a lot of work, as we would need to compute higher-order terms in the expansions on the aand γ-scales, and use higher-order asymptotic matching. The γ-scale would involve numerically solving certain integral equations.

Another idea to approximating log ( n b(n, p)), N and V would be to define ( ŵa , na ) as the solution to the system .

Of course, asymptotically these approximations have the same order of accuracy as the explicit results in (2.15), (2.18) and (2.19), as both are based on the approximation (9.13) for g(w). However, numerically ( Φa , na , V a ) lead to somewhat different results. In Table 14, we give these approximations for various p in the range [START_REF] Régnier | A Limiting Distribution for Quicksort[END_REF]250]. These should be compared to the exact values in Tables 111213. In doing the numerics we first replaced n by n = exp g a (w) log 2 log w -1 (9.13)

where E * (s) is given by (5.87). Noting that Y 2 = -a √ Y , we have thus regained the correction term C (1) We have thus shown that a perturbation analysis of (2.2) yields results essentially equivalent to those obtained by the methods of Section 5. However, here we had to interpret certain integrals as distributions; this issue never arose when we analyzed (2.1) or the moment equations (5.7).

1 .

 1 the recurrence relation b(n, p) = k+ℓ=n-1 r+s+n-1=p b(k, r)b(ℓ, s), n ≥ 1 with the boundary conditions b(0, 0) = 1; b(0, p) = 0, p ≥The generating functionB n (w) = ∞ p=0 b(n, p)w p becomes B n+1 (w) = w n n ℓ=0 B ℓ (w)B n-ℓ (w), n ≥ 0 (2.1)with B 0 (w) = 1. Furthermore, the double generating functionB(z, w) = ∞ n=0 ∞ p=0 b(n, p)w p z n = ∞ n=0 z n B n (w)satisfies the functional equation B(z, w) = 1 + zB 2 (zw, w).

2 )

 2 about w = 1, with B n (w) = a n + b n (w -1) + 1 2 c n (w -1) 2 + O((w -1) 3 ) (2.5) B(z, w) = a(z) + b(z)(w -1) + 1 2 c(z)(w -1) 2 + O((w -1) 3 ) (2.6)we are led tob(z) = B w (z, 1) = 2z 2 B(z, 1)B z (z, 1) 1 -2zB(z, 1) = 2z 2 a(z)a ′ (z) 1 -2za(z)

+2w 4

 4 -2n B2 (w) Bn-2 (w) + • • • ]. (3.2) But B 0 (w) = B 1 (w) = 1 and B 2 (w) = 2w so that (3.1) yields B0 (w) = 2, B1 (w) = 1, B2 (w)

Fig. 1 :

 1 Fig. 1: Plots of ℜ[B * (w)] and ℑ[B * (w)] for 1.5 < |w| < 4.

  we plot ℜ[B * (w)] and ℑ[B * (w)] for 1.5 < |w| < 4. These surfaces are viewed from the +u perspective, where w = u + iv. With (3.1) and (3.5) we can infer the behavior of b(n, p) for p close to p max

  27) -(4.29) in (4.26) and multiplying the equation by 1e -β leads to

  .47) Using (4.47) in (4.44) the right side becomes

17

 17

  )

2 .

 2 Y dτ, Y > 0 may be similarly continued into the range arg(Y ) ∈ -5π 6 , π 6 in this range arg(ω 2 Y ) ∈ -3π 2 , -π Using (5.53) we write I 1 = I 3 + I 4 where

, 5π 6 or

 6 5.56) we see that both h(ωτ ) and h(ω 2 τ ) grow faster than exponentially as τ → +∞, and thus I 3 defines an entire function of Y . The integral I 4 may be viewed as defining an analytic function in the range arg(Y ) ∈ -π 6 arg(Y ) ∈ -13π 6 , -7π 6 . Now we let arg(Y ) = -2π 3 and set Y = ω 2 y = ω 2 a 2/3(5.57)

3 2 N 3 2 N - 2 ) 3 2 N - 5 2

 23235 -2 -1 dx where we used a ℓ 4 -ℓ ∼ 1/[ √ πℓ 3/2 ]. The remaining terms in (5.63) must be estimated more precisely, as we need both the O(n leading term (as in (5.14)) and the O(n ) correction. Breaking up the sum into two pieces we obtain

3 2 N - 5 2

 35 .70) Combining (5.61), (5.62) -(5.64), (5.69) and (5.70) we see that at O n O n

  .76) Combining (5.72) -(5.76) in (5.71) the equation transforms to (after dividing by Y )

  .117) Thus the scale w -1 = O(n -3/2 ) in the transform space translates to (5.116) in the (n, p) space, and (5.117) and (5.98) lead to b

2 +

 2 .44) as Λ = p/n 2 = Ω/ √ n → 0. The latter is easily obtained from (4.40) and (4.44). As Λ → 0, we have β * → 0 and from (4.16) we obtain φ(β) = log 2 -O(β 3 ), β → 0. (5.126) It follows from (4.40) since φ ′ (β) + 1 2 ∼ 2 3 β that

17

 17

) and α 1 = κ 4 + 1 -

 41 γ E = 2.9692 . . . .(6.18)

. 26 )γ γ e k1γ log γ e k2γ γ k3 k 4 ∼ 1 0γ∼ γ 2k3 1 4 k3 k 2 4 ee k2γ e k1γ log γ e -k1γ log 2 π

 26142 We use(6.26) in(6.20) and evaluate the integral by Laplace's method, to yield e 2k3 x k3 (1x) k3 k 2 4 e k2γ e k1[γ log γ+γH(x)] dx k2γ e k1γ log γ e k1γH(

so that w n log 2

 2 n w -p-1 = w -αn-1 and we have, for α fixed, b(n, p) ≈ 1 2πi C 1 w e g(w) n log 2 wlog 2 w 2πn w 2+ 1 log 2 exp[n(g(w)α log w)]dw. (7.31)

. 20 ) 2 F

 202 Here we write w * = w * (n, p) via (8.1) and letw = w(p) = w * (ñ(p), p), ∆ = ∆(p) = 1w * (ñ(p), p), α = pñ(p) log 2 [ñ(pg ′′ ( w) n log 2 we g( w)log 2 w × ∞ n=-∞exp ng( w)nα log w + 1 nn (ñ(p); p)(nñ(p)) 2(8.22)

(w - 1 )

 1 log[B * (w)] → d 0 =

2 ,

 2 n (w)w -n log 2 n n log 2 w √log ww 2+ 1 log

  2) for w ↑ 1. Let us define CONST = CONST(w; n) 1) log(1w) . (9.10)

3 2 A 0 log 2 a 1 / 3 3 = 1 .

 21331 = 2 1/3 (log 2) 1/3 |r 0 | = 2.6070 . . . mean: R 1 ≡ 2 3 2 1/3 (log 2) 1/3 |r 0 | = 1.7380 . . . 3272 . . . .

(a) p = 190 with n in the range [ 20 Fig. 2 :

 202 Fig. 2: Plots of b(n, p).

  wg ′ a (w) = α = p n log 2 n, g a (w) + log n + 1 log 2 log w = 0 (9.11)whereg a (w) ≡ log 4 + 4 1/3 r 0 (1w) 1) log(1w)k 2 (w -1) (9.12)is the truncation of the asymptotic series (8.2), up to the point that we explicitly computed it. For a given (large) p, we solve (9.11) (with (9.12)) for (w, n) numerically. Then our approximations to the growth rate, mean and variance are Φa = Φa (p) = na [g a ( ŵa )αa log( ŵa )],na and V a (p) = = -

  (a) in (2.26) and (2.28). Here we must again interpret1 2πi Br ζe -ζ dζ = Y 2 2πi Br se Y s ds = Y 2 δ ′ (Y ) = 0as a distribution.

  ). in distribution and A is the random variable possessing the Airy distribution. It is characterized by moments[START_REF] Flajolet | Analytic Variations on the Airy Distribution[END_REF]

	More precisely,	Pr{L n / √	2n 3 ≤ x}

d →Pr{A ≤ x} where d → denotes convergence

  Asymptotically, for n → ∞, [p min , p max ] ∼ n log 2 n, n 2 2 and, for p → ∞, [n min , n max ] ∼

	2p,	p log 2 p	.

.13) Here p min and p max are the minimal and maximal total path lengths possible in a tree with n nodes. If we view the problem as having p fixed and varying n, then b(n, p) is non-zero in the range n ∈ [n min (p), n max (p)] where n min (p) = min{n : n 2 ≥ p} and n max (p) = max n : n J=2 ⌊log 2 J⌋ ≤ p .

  as a → +∞, determines the constant in (4.31) and (4.36) as

	With (4.36) and (4.37) we obtain (2.25).
	const. =	4 √ π	.	(4.37)

  To find the maximum of ng(w * )nα log w * we need to asymptotically solve Here we note that w * depends on n through α, but g ′ (w * ) = α/w * by (8.1). Also, nα = pn log 2 n so that ∂(nα)/∂n =log 2 n -1/ log 2. Let us write the solution of (8.13) as n = ñ = ñ(p). Setting

	∂ ∂n	[ng(w ) +	log n + 1 log 2	log w * = 0.	(8.13)
			F (n; p) = ng(w * ) -nα log w *	(8.14)
	we have
		F n (n; p) = g(w * ) +	log n + 1 log 2	log w *	(8.15)
	and
	F nn (n; p) = g ′ (w * )	∂w * ∂n	+	log w * n log 2	+	log n + 1 w * log 2	∂w * ∂n
		=	α +	log n + 1 log 2	1 w *	∂w * ∂n	+	1 n log 2	log w * .	(8.16)
	We also define
			Ψ 0 (p) = F (ñ(p); p)	(8.17)
	and
			V 0 (p) = -1/F nn (ñ(p); p).	(8.18)

* )nα log w * ] = g(w *

  1/3 ]. We compare this to p min (n) = n log 2 n + O(n) (cf. (2.13)). From (8.9) and(8.19) we conclude that

  4/3 .To refine(8.25) and make (8.23) more explicit, we consider (8.1) and (8.13) as a simultaneous system to determine w and ñ (or, equivalently, α) as functions of p. By expanding (8.13) about w

	Thus (8.22) becomes
	n	b(n, p) ∼	1 2πñ	1 log ñ	7/6	3 A 0	2 7/6 (log 2) 11/6
				×	2π -F nn (ñ(p); p)	exp[F (ñ(p); p)].	(8.23)
	In view of (8.13) and (8.14) we have
	F (ñ(p); p)) = -p +	ñ log 2	log w = p +	ñ log 2	[ ∆ + O( ∆2 )].	(8.24)
	From (8.23), (8.24), (8.19) and (8.20) we obtain the growth rate
		log	n	b(n, p) ∼	p log p	(2 log 2 2), p → ∞.	(8.25)
								= 1 and
	setting
					S = S(p) = log ñ	(8.26)

*

  .41) This follows by expanding (2.42) for w * → 1, and also matches to (6.40), as can be seen by replacing α in (8.41) by p/nlog 2 n ∼ p/n. We note that in this limit Expanding (8.41) about n = ñ(p) (or α = α) leads to the Gaussian form in(2.20).For numerical purposes it sometimes proves useful to base our approximations on (8.41). To this end we define

	γ * =	32 27	|r 0 | 3 Θ 3 =	32 27	|r 0 | 3 n 4 p 3 ∼	A 3 0 α 3 n.

Table 3 ,

 3 for various w > 1. Our analysis predicts that D 1 (w) → d 1 = 2.1350 . . . as w ↓ 1. It became very difficult to compute B∞ (w) for w ≤ 1.04, since for w 1 the convergence of Bn to B∞ is very slow. We see from Table3that D 1 (w) still changes appreciably as w goes from 1.08 to 1.06 to 1.04. The inconsistent with the theoretical value, but w would have to be made much closer to one for us to reach a more definite conclusion.

		Tab. 2:
	w	(w -1) log( B∞ ) (w -1) log( B∞ ) -1 2 log(w -1)
	2	.358	.358
	1.8	.372	.461
	1.6	.390	.544
	1.4	.417	.601
	1.2	.462	.623
	1.18	.468	.622
	1.16	.475	.621
	1.14	.482	.620
	1.12	.490	.618
	1.10	.499	.614
	1.08	.510	.611
	1.06	.521	.606
	1.04	.535	.600
	1.02	.552	.592
	1	.582	.582

(N -1) m N + 1 √ n mN + 1 n mN -1 √ πn m N + O(n -3/2 ) .

N -2 we regain(5.16). At

Acknowledgements

We thank Dr. Gadiel Seroussi (HPL, Palo Alto) for motivating the problem, and sharing with us his early results concerning the number of trees of a given path length.

† The work was supported by NSF Grant DMS 05-03745 and NSA Grant MDA 904-03-1-0036. ‡ The work was supported by NSF Grants CCR-0208709 and DMS 05-03745, NIH Grant R01 GM068959-01, and NSA Grant MDA 904-03-1-0036. 1365-8050 c 2005 Discrete Mathematics and Theoretical Computer Science (DMTCS)

A APPENDIX

Here we briefly re-examine the central region by using the functional equation (2.2). We introduce the scaling w = 1 + a n 3/2 , z =

we obtain from (A.3)

we obtain from (A.4) at the first two orders (O(n -1 ) and O(n -3/2 )) the equations

and

To solve (A.5) we set

to obtain the Airy equation 4a 2 H ζζ + ζH = 0 and thus

We use (A.8) in (A.7) and note that

We have

where again the first integral in (A.11) must be interpreted using distributions. We have thus regained the leading term for the scale w = 1 + O(n -3/2 ) and w < 1. However, the present derivation is not as clean as those in Section 5, where all the functions that arose may be interpreted in the classical sense.

We next obtain the correction term B (1) (ζ, a). Differentiating (A.5) yields

We rewrite the linear equation (A.6) as

ζ + 2

By multiplying (A.13) by H 2 and using (A.7) and (A.12), we obtain

which follows from (A.8) and the exponential decay rate of the Airy function Ai(•). It follows that

and the right side of (A.14) will be asymptotic to

To avoid this growth we set

we see that the approximation to B n (w) becomes