A STIRLING-TYPE FORMULA FOR THE DISTRIBUTION OF THE LENGTH OF LONGEST INCREASING SUBSEQUENCES, APPLIED TO FINITE SIZE CORRECTIONS TO THE RANDOM MATRIX LIMIT

. The discrete distribution of the length of longest increasing subsequences in random permutations of order n is deeply related to random matrix theory. In a seminal work, Baik, Deift and Johansson provided an asymptotics in terms of the distribution of the largest level of the large matrix limit of GUE. As a numerical approximation, however, this asymptotics is inaccurate for small lengths and has a slow convergence rate, conjectured to be just of order n − 1 / 3 . Here, we suggest a diﬀerent type of approximation, based on Hayman’s generalization of Stirling’s formula. Such a formula gives already a couple of correct digits of the length distribution for n as small as 20 but allows numerical evaluations, with a uniform error of apparent order n − 3 / 4 , for n as large as 10 12 ; thus closing the gap between a table of exact values (that has recently been compiled for up to n = 1000 ) and the random matrix limit. Being much more eﬃcient and accurate than Monte-Carlo simulations for larger n , the Stirling-type formula allows for a precise numerical understanding of the ﬁrst few ﬁnite size correction terms to the random matrix limit, a study that has recently been initiated by Forrester and Mays, who visualized the form of the ﬁrst such term. We display also the second one, of order n − 2 / 3 , and derive (heuristically) expansions of expected value and variance of the length, exhibiting several more terms than previously known.


Introduction
As witnessed by a number of outstanding surveys and monographs (see, e.g., [1,5,44,48]), a surprisingly rich topic in combinatorics and probability theory, deeply related to representation theory and to random matrix theory, is the study of the lengths L n (σ) of longest increasing subsequences of permutations σ on the set [n] = {1, 2, . . ., n} and of the behavior of their distribution in the limit n → ∞.Here, L n (σ) is defined as the maximum of all k for which there are Writing permutations in the form σ = (σ 1 σ 2 • • • σ n ) we get, e.g., L 9 (σ) = 5 for σ = (4 1 2 7 6 5 8 9 3), where one of the longest increasing subsequences has been highlighted.Enumeration of the permutations with a given L n can be encoded probabilistically: by equipping the symmetric group on [n] with the uniform distribution, L n becomes a discrete random variable with cumulative probability distribution (CDF) P(L n l) and probability distribution (PDF) P(L n = l).
Constructive combinatorics.Using the Robinson-Schensted correspondence [45], one gets the distribution of L n in the following form (see, e.g., [49, § §3.3-3.7]): (1) Here λ n denotes an integer partition λ 1 λ 2 • • • λ l λ > 0 of n = l λ j=1 λ j and d λ is the number of standard Young tableaux of shape λ, as given by the hook length formula.By generating all partitions λ n, in 1968 Baer and Brock [3] computed tables of P(L n = l) up to n = 36; in 2000 Odlyzko and Rains [40] for n = 15, 30, 60, 90, 120 (the tables are online, see [38]), reporting a computing time for n = 120 of about 12 hours (here p n , the number of partitions, is of size 1.8 × 10 9 ).This quickly becomes infeasible, 1 as p n is already as large as 2.3 • 10 14 for n = 250.Another use of the combinatorial methods is the approximation of the .Discrete distribution of Ln for n = 10 5 near its mode vs. the random matrix limit given by the leading order terms in ( 3) and (42a) (solid red line); here and in the figures below, discrete distributions are shown as blue bars centered at the integers.Left: CDF P(Ln l); right: PDF P(Ln = l).The discrete distributions were computed using the Stirling-type formula (5), with additive errors estimated to be smaller than 10 −5 , cf.Fig. 3, which is well below plotting accuracy.
distribution of L n by Monte Carlo simulations [3,40]: one samples random permutations σ and calculates L n (σ) by the Robinson-Schensted correspondence. 2nalytic combinatorics and the random matrix limit.For analytic methods of enumeration the starting points is a more or less explicit representation of a suitable generating function; here, the suitable one turns out to be the exponential generating function of the CDF P(L n l), when considered as a sequence of n with the length l fixed: (We note that f l is an entire function of exponential type.)In fact, Gessel [30, p. 280] obtained in 1990 the explicit representation (2) f l (z 2 ) = D l (z), D l (z) := l det j,k=1 in terms of a Toeplitz determinant of the modified Bessel functions I m , m ∈ Z, which are entire functions of exponential type themselves.By relating, first, the Toeplitz determinant to the machinery of Riemann-Hilbert problems to study a double-scaling limit of the generating function and by using, next, a Tauberian theorem3 to induce from that limit an asymptotics of the coefficients, Baik, Deift and Johansson [4] succeeded 1999 in establishing 4(3) uniformly in l ∈ N; it will be called the random matrix limit of the length distribution throughout this paper since F 2 (t) denotes the Tracy-Widom distribution for β = 2 (that is, the probability that in the soft-edge scaling limit of the Gaussian unitary ensemble (GUE) t = F2(t).
However, since the limit distribution F2 is continuous, by a standard Tauberian follow-up [54, Lemma 2.1] of the Portmanteau theorem in probability theory, this convergence is known to hold, in fact, uniformly in t.
Table 1.The values of P(Ln l) and various of its approximations for n = 20.The Monte Carlo simulation was run with T = 5 • 10 6 samples.Clearly observable is a multiplicative (i.e., relative) error of already less than 1% for the Stirling-type formula (5) as well as an additive (i.e., absolute) error of order 4 • 10 −4 ≈ T −1/2 for Monte Carlo and an additive error of order 10 −1 ≈ n −1/3 for the limit (3) from random matrix theory.the scaled largest eigenvalue is bounded from above by t).This distribution can be evaluated numerically based on its representation either in terms of the Airy kernel determinant [23] or in terms of the Painlevé-II transcendent [51]; see Remark 3.1 and [9] for details.
As impressive as the use of the limit (3) might look as a numerical approximation to the distribution of L n near its mode for larger n, cf.Fig. 1, there are two notable deficiencies: first, since the error term in (3) is additive (i.e., w.r.t.absolute scale), the approximation is rather poor for l 2 √ n; second, the convergence rate is rather slow, in fact conjectured to be just of the order O(n −1/3 ); see [27] and the discussion below.Both deficiencies are well illustrated in Table 1 for n = 20 and in Fig. 2 for n = 1000.
A Stirling-type formula.In this paper we suggest a different type of numerical approximation to the distribution of L n that enjoys the following advantages: (a) it has a small multiplicative (i.e., relative) error, (b) it has faster convergence rates, apparently even faster than (3) with its first finite size correction term added, and (c) it is much faster to compute than Monte Carlo simulations.In fact, the distribution of L n for n = 10 5 , as shown in Fig. 1, exhibits an estimated maximum additive error of less than 10 −5 and took just about five seconds to compute; whereas Forrester and Mays [27] have recently reported a computing time of about 14 hours to generate T = 5 • 10 6 Monte Carlo trials for this n; the error of such a simulation is expected to be of the order 1/ √ T ≈ 4 • 10 −4 .Specifically, we use Hayman's generalization [33] of Stirling's formula for H-admissible functions; for expositions see [22,39,55].For simplicity, as is the case here for f (z) = f l (z), assume that is an entire function with positive coefficients a n and consider the real auxiliary functions If f is H-admissible, then for each n ∈ N the equation a(r n ) = n has a unique solution r n > 0 such that b(r n ) > 0 and the following generalization of Stirling's formula 5 holds true: (4) As in Table 1, the error is already below 1% for n as small as n = 20.).The discrete distribution of Ln is shown near its mode vs. the random matrix limit given by the leading order terms in ( 9) and (42a) (solid red line).Left: CDF P(Ln l); right: PDF P(Ln = l).The exact values of the distribution of Ln and their approximation by the Stirling-type formula (5) differ just by additive errors of the order 10 −4 (see Fig. 3), which is well below plotting accuracy.
We observe that the error is multiplicative here.In Thm.2.2 we will prove, using some theory of entire functions, the H-admissibility of the generating functions f l .Hence the Stirling-type formula (4) applies without further ado to their coefficients P(L n l)/n!.Since the error is multiplicative, nothing changes if we multiply the approximation by n! and we get ( 5) where we have labeled all quantities when applied to f = f l by an additional index l.An approximation to P(L n = l) is then obtained simply by taking differences.The power of these approximations, if used as a numerical tool even for n as small as n = 20, is illustrated in Table 1 and Figs.1-3, as well as in Table 2 below.
Numerical evaluation of the generating function.For the Stirling-type formula (5) to be easily accessible in practice, we require an expression for f l (r) that can be numerically evaluated, for r > 0, in a stable, accurate, and efficient fashion.Since the direct evaluation of the Toeplitz determinant (2) is numerically highly unstable, and has a rather unfavorable complexity of O(l 3 ) for larger l, we look for alternative representations.One option-used in [11] to numerically extract P(L n l) from f l (z) by Cauchy integrals over circles in the complex plane that are centered at the origin with the same radius r l,n as in (5)-is the machinery, cf., e.g., [15], to transform Toeplitz determinants into Fredholm determinants which are then amenable for the numerical method developed in [10].However, since we need the values of the generating function f l (r) for real r > 0 only, there is a much more efficient option, which comes from yet another connection to random matrix theory.
To establish this connection we first note that an exponentially generating function f (r) of a sequence of probability distributions has a probabilistic meaning if multiplied by e −r : a process called Poissonization.Namely, if the draws from the different permutation groups are independent and if we take N r ∈ N 0 := {0, 1, 2, 3, . ..} to be a further independent random variable with Poisson distribution of intensity r 0, we see that ( 6) is, for fixed r 0, the cumulative probability distribution of the composite discrete random variable L Nr .On the other hand, for fixed l ∈ N, also e −r f l (r) turns out to be a probability to the points in display.Red +: error of leading order terms (random matrix limit) in ( 9), (42a); α = ( 1 3 , 2 3 , 1). Green •: error of expansions ( 9), (42a) truncated after the first finite size correction term; α = ( 23 , 1, 4 3 ), where F2,1 has been approximated as in Fig. 4. Blue •: error of the Stirling-type formula (5); α = ( (0; [0, 4r], l) = e −r f l (r).
Here, E (0; [0, t], l) denotes the probability that, in the hard-edge scaling limit of the Laguerre unitary ensemble (LUE) with parameter l, the smallest eigenvalue is bounded from below by t 0. Now, the point here is that this distribution can be evaluated numerically, stable and accurate with a complexity that is largely independent of l, based on two alternative representations: either in terms of the Bessel kernel determinant [23] or in terms of the Jimbo-Miwa-Okamoto σ-form of the Painlevé-III transcendent [52]; see [9] for details.We will show in Sect. 3 that the auxiliary functions a l (r) and b l (r) fit into both frameworks, too.
Finite size corrections to the random matrix limit.In a double logarithmic scaling, a plot of the additive errors (taking the maximum w.r.t.l ∈ {1, 2, . . ., n}) in approximating the distribution P(L n l) by either the random matrix limit (3) or by the Stirling-type formula (5) exhibits nearly straight lines; see Fig. 3 for n between 10 and 1000.Fitting the data in display to a model of the form c 1 n −α 1 + c 2 n −α 2 + c 3 n −α 3 with simple triples α = (α 1 , α 2 , α 3 ) of rationals strongly suggests that, uniformly in l ∈ {1, 2, . . ., n} as n → ∞, The approximation order (and the size of the implied constant) in (8b) is much better than the one in (8a) so that the Stirling-type formula can be used to reveal the structure of the O(n −1/3 ) term in the random matrix limit.In fact, as the error plot in Fig. 3 suggests and we will more carefully argue in Sect.4.1, this can even be iterated yet another step and we are led to the specific conjecture6 (9) Rescaled differences between the distributions of Ln and their expansions truncated after the leading order term (i.e., the random matrix limit)-see (9) for the CDF resp.(42a) for the PDF; data points (to avoid clutter just every 5 th is displayed) have been calculated using the Stirling-type formula (5) for n = 10 6 (red +), n = 10 8 (green •), n = 10 10 (blue •).Left: CDF errors rescaled by [27,Fig. 7] for a similar figure with data from Monte Carlo simulations for n = 2 • 10 4 and n = 10 5 .The solid line is a polynomial F2,1(t) of degree 64 fitted to the 836 data points for n = 10 10 with −8 t 10; it approximates F2,1(t) in that interval.Right: PDF errors rescaled by . The solid line displays the function F 2,1 (t) + F 2 (t)/24 as an approximation of F 2,1 (t) + F 2 (t)/24, with the polynomial F2,1(t) taken from the left panel.The dotted line shows the term F 2,1 (t) only.
as n → ∞, uniformly in l ∈ N. Compelling evidence for the existence of the functions F 2,1 and F 2,2 is given in the left panels of Figs. 4 and 6. 7  We note that a corresponding expansion 8 for the Poissonization (6) of the length distribution was studied by Baik and Jenkins [6, Eq. ( 25)] (using the machinery of Riemann-Hilbert problems up to an error of order O(r −1/2 )) and by Forrester and Mays [27, Eqs.(1.18), (2.29)] (using Fredholm determinants), who obtained the expansion, as r → ∞ for bounded t * l : (10a) , with the explicit functional form (identified by means of Painlevé representations) Though (10a) adds to the plausibility of the expansion ( 9), the de-Poissonization lemma of Johansson [35,Lemma 2.5] and its commonly used variants (see [4,5,44]) would not even allow us to deduce from (10) the existence of the term F 2,1 (t), let alone to obtain its functional form.
Remark (added in proof).On the other hand, by inserting the Poissonized expansion (10) (and the induced expansions of the quantities b(r) and r l,n ) into the Stirling-type formula (8b) with its conjectured error of 7 Based on Monte-Carlo simulations, Forrester and Mays [27, Eq. (1.10) and Fig. 7] were recently also led to conjecture an expansion of the form However, the substantially larger errors of Monte Carlo simulations as compared to the Stirling-type formula (5) would inhibit them from getting, in reasonable time, much evidence about the next finite size correction term.

8
Expansions of probability distributions are sometimes called Edgeworth expansions in reference to the classical one for the central limit theorem.In random matrix theory a variety of such expansions have been studied: e.g., for the soft-edge scaling limits of GUE/LUE [16] and GOE/GSE [17], for the hard-edge scaling limit of LUE [12] and LβE [28], for the bulk scaling limit of CUE/COE/CSE [14].For the the hard-to-soft edge transition limit of LUE see the expansion (10) and its discussion.
order O(n −2/3 ), we are led to the conjecture This functional form is in perfect agreement with the data displayed in Fig. 4; see Footnote 29 and Remark 4.4 for further numerical evidence.Details will be given in a forthcoming paper of the author [13], where the expression (11) for F2,1(t) (as well as one for F2,2(t)) is also obtained by a complex-analytic modification (related to H-admissibility) of the de-Poissonization process.
We will argue in Sect.4.3 that the expansion (9) of the length distribution allows us to derive an expansion of the expected value of L n , specifically Similarly, we will derive in Sect.4.4 an expansion of the variance of L n of the form The values of µ 0 and ν 0 are the known values of mean and variance of the Tracy-Widom distribution F 2 , cf. [9, Table 10].(The leading parts of ( 12) up to µ 0 n 1/6 and of ( 13) up to ν 0 n 1/3 had been established previously by Baik, Deift and Johansson [4, Thm.1.2].) 2. H-Admissibility of the Generating Function and its Implications 2.1.H-admissible functions.For simplicity we restrict ourselves to entire functions.We refrain from displaying the rather lengthy technical definition of H-admissibility, 9 which is difficult to be verified in practice and therefore seldomly directly used.Instead, we start by collecting some usefuls facts and criteria from Hayman's original paper [33]: 10   Theorem 2.1 (Hayman 1956).Let f (z) = ∞ n=0 a n z n and g(z) be entire functions and let p(z) denote a polynomial with real coefficients.
a. f (r) > 0 for all sufficiently large r > 0, so that in particular the auxiliary functions 11 b. for r > 0 as in I.a there is log f (r) strictly convex in log r, a(r) strictly monotonically increasing, and b(r) > 0 such that a(r), b(r) → ∞ as r → ∞; in particular, for large integers n there is a unique c. if the coefficients a n of f are all positive, then I.b holds for all r > 0; 9 Since we consider entire functions only, H-admissibility is here understood to hold in all of C.

10
Interestingly, the powerful criterion in part III (which is [33, Thm.XI]) is missing from the otherwise excellent expositions [22,39,55] of H-admissibility. 11In terms of differential operators we have r d dr = d d log r .
II.If f and g are H-admissible, then: e. f (z)g(z), e f (z) and f (z) + p(z) are H-admissible; f. if the leading coefficient of p is positive, f (z)p(z) and p(f (z)) are H-admissible; g. if the Taylor coefficients of e p(z) are eventually positive, e p(z) is H-admissible.
III.If f has genus zero 12 with, for some δ > 0, at most finitely many zeros in the sector | arg z| π 2 + δ and satisfies I.a such that b(r) → ∞ as r → ∞, then f is H-admissible.Obviously, the Stirling-type formula (4) is obtained from the approximation result (15) by just inserting the particular choice r = r n .
Remark 2.1.We observe that, if a n 0, Eq. ( 15) has an interesting probabilistic content: 13  as a distribution in the discrete variable n ∈ N 0 , the Boltzmann 14 probabilities a n r n /f (r) associated with an H-admissible entire function f are, for large intensities r > 0, approximately normal with mean a(r) and variance b(r); see the right panel of Fig. 5 for an illustrative example using the generating function f 5 (z).The additional freedom that is provided in the normal approximation (15) by the uniformity w.r.t.n will be put to good use in Sect.2.3.
12 By definition, an entire function f has genus zero if it is a polynomial or if it can be represented as a convergent infinite product of the form where c ∈ C is a constant, m ∈ N0 is the order of the zero at z = 0, and z1, z2, . . . is the sequence of the non-zero zeros, where each one is listed as often as multiplicity requires. 13Note that the particular case f (z) = e z (which is H-admissible by Thm.2.1.II.g) specifies to the well-known normal approximation of the Poisson distribution for large intensities-which is a simple consequence of the central limit theorem if we observe that the sum of k independent Poisson random variates of intensity ρ is one of intensity r = ρk. 14We follow the terminology in the theory of Boltzmann samplers [20], a framework for the random generation of combinatorial structures.Note that mean and variance of the Boltzmann propabilities anr n /f (r) are exactly the auxiliary functions a(r) and b(r) as defined in (14), cf.[20, Prop.2.1].
The classification of entire functions (by quantities such as genus, order, type, etc.) and their distribution of zeros is deeply related to the analysis of their essential singularity at z = ∞.For the purposes of this paper, the following simple criterion is actually all we need.The proof uses some theory of entire functions, which can be found, e.g., in [37].
Lemma 2.1.Let f (z) be an entire function of exponential type with positive Maclaurin coefficients.If there are constants c, τ, ν > 0 such that there holds, for the principal branch of the power function and for each 0 < δ π 2 , the asymptotic expansion15 (16) For r → ∞ the associated auxiliary functions a(r) and b(r) satisfy and the solution r n of a(r n ) = n satisfies Proof.The expansion ( 16) is equivalent to which readily implies: • since δ > 0 is arbitrary, f has the Phragmén-Lindelöf indicator so that f has order 1 2 and type τ , hence genus zero; Since the Maclaurin coefficients of f are positive we have f (r) > 0 for r > 0 and the auxiliary functions a(r), b(r) in ( 14) are well-defined for r > 0. In fact, both functions can be analytically continued into the domain of uniformity of the expansion (19) and by differentiating this expansion (which is, because of analyticity, legitimate by a theorem of Ritt, cf.[41]) we obtain (17); this implies, in particular, b(r) → ∞ as r → ∞.Thus, all the assumptions of Thm.2.1.III are satisfied and f is shown to be H-admissible.

2.2.
Singularity analysis of the generating function at z = ∞.Establishing an expansion of the form (16) for D l (z) = f l (z 2 ) as given by (2), that is to say, for the Toeplitz determinant suggests to start with the expansions (valid for all 0 < δ π 2 , see [41, p. 251]) This does not yield (16) at once, as there could be, however unlikely it would be, eventually a catastrophic cancellation of all of the expansion terms when being inserted into the determinant expression defining D l (z).For the specific cases l = 1, 2, . . ., 8 a computer algebra system shows that exactly the first l − 1 terms of the expansion (21) mutually cancel each other in forming the determinant, and we get by this approach 16 the expansions All of them, inherited from (21), are valid as z → ∞ while |arg z| π 2 − δ with the uniformity content implied by the symbol O(z −1 ).From these instances, in view of ( 21) and the multilinearity of the determinant, we guess that and observe Though this is very likely to hold for all l ∈ N-a fact that would at once yield the Hadmissiblity of all the generating function f l by Lemma 2.1-a proof seems to be elusive along these lines, but see Remark 2.3 for a remedy.
Inspired by the fact that the one-dimensional Laplace's method easily gives the leading order term in (21) when applied to the Fourier representation (22) I m (2z) = 1 2π we represent the Toeplitz determinant D l (z) in terms of a multidimensional integral and study the limit z → ∞ by the multidimensional Laplace method discussed in the Appendix.In fact, (22) shows that the symbol of the Toeplitz determinant D l (z) is exp(2z cos θ) and a classical formula of Szegő's [50, p. 493] from 1915, thus gives, without further calculation, the integral representation 18 denotes the Vandermonde determinant of the complex numbers w 1 , . . ., w l . 16Odlyzko [39, Ex. 10.9] reports that he and Wilf had used this approach, before 1995, for small l in the framework of the method of "subtraction of singularities" in asymptotic enumeration.He states the expansions for D4(z) and D5(z), cf.[39, Eqs.(10.30)/(10.39)],with a misrepresented constant factor in D5(z), though.
No attempt, however, was made back then to guess the general form. 17https://oeis.org/A000178 18This induces, see (7), an integral representation of the distribution E (hard) 2 (0; [0, s], l) which has been derived in 1994 by Forrester [24] using generalized hypergeometric functions defined in terms of Jack polynomials.
Remark 2.2.By Weyl's integration formula on the unitary group U (l), cf.[46, Eq. 1.5.89], the integral (23) can be written as where the expectation E is taken with respect to the Haar measure.Without any reference to (2), this form was derived in 1998 by Rains [42,Cor. 4.1] directly from the identity which he had obtained most elegantly from the representation theory of the symmetric group.
We are now able to prove our main theorem.
Theorem 2.2.For each 0 < δ π/2 and l ∈ N there holds the asymptotic expansion Thus, by Lemma 2.1, the generating functions f l (z) are H-admissible and their auxiliary functions satisfy, as r → ∞, ). Proof.We write (23) in the form The phase function of this multidimensional integrand, that is to say ∆(e iθ 1 , . . ., e iθ l ) where the degree of the homogeneous polynomial ∆(θ 1 , . . ., θ l ) 2 is l(l − 1).Therefore, by the multidimensional Laplace method as given in Corollary A.1 (see also formula (54) following it) we obtain immediately where the evaluation of this multiple integral is well-known in random matrix theory, e.g., as a consequence of Selberg's integral formula, cf.[2, Eq. (2.5.11)]. 19This factor is zero at θ * = 0, so that Hsu's variant (48) of Laplace's method, which is the one predominantly found in the literature, does not yield the leading term of D l (z).That we have to expand the non-exponential factor up to degree l(l−1) for the first non-zero contribution to show up, corresponds to the mutual cancellation of the leading terms of ( 21) when being inserted into the Toeplitz determinant that defines D l (z).
Remark 2.3.By (7), Thm.2.2 implies an asymptotics first rigorously proven, using Riemann-Hilbert problem machinery, by Deift, Krasovsky and Vasilevska [19] in 2010.Besides that our proof is much simpler, their result, which is for real s > 0 only, would by itself not suffice to establish the H-admissibility of the generating function f l (z); one would have to complement it with the arguments given above for expanding the Toeplitz determinant (20) based on the expansions (21) of the modified Bessel functions.However, their result is more general in another respect: it covers parameters α ∈ C of the LUE with α > −1 instead of just l ∈ N; the superfactorial factor ! is then to be replaced by G(1 + α), where G(z) is the Barnes G-function. 20  We complement the large r expansion (24) of the auxiliary functions with their expansions as r → 0 + , which are simple consequences of elementary combinatorics.
Lemma 2.2.The auxiliary functions of the generating function f l satisfy, as r → 0 + , Proof.Because of L n n and since there is just one permutation σ with L n = n, we get This implies, by truncating the power series of f l at order l + 1, Logarithmic differentiation of the power series thus yields, as z → ∞, ), and the results follow from specializing to z = r > 0.
The left panel of Fig. 5 visualizes the expansions ( 24) and (26) for a 5 (r).Apparently, as a function of log r, the auxiliary log a l (r) interpolates monotonically, concavely, and from below between the following two extremal regimes: • log r as r → 0 + , which reflects, by the proof of Lemma 2.2, the regime l n, and • 1 2 log r + log l as r → ∞, which reflects the regime l n 1/4 (see Sect. 2.3).
It is this seamless interpolation between the two regimes l n and l n 1/4 that helps to understand the observed uniformity of the Stirling-type formula w.r.t.l, cf.(8b).

2.3.
A new proof of Regev's asymptotic formula for fixed l.A simple closed form expression in terms of n and l is obtained by studying the asymptotics, as n → ∞ for fixed l, of the Stirling-type formula (5) itself-or even easier yet, because of its added flexibility, of Hayman's normal approximation (15) for a suitable choice of r.As tempting as it might appear, however, this stacking of asymptotics leads, first, to a considerable loss of approximation power for small n, cf.Table 2, and, second, to a lack of uniformity w.r.t.l since the result is effectively conditioned to the constraint l n 1/4 . 20For real α > −1, Tracy and Widom [52] had conjectured this asymptotics in 1994 based on a guess of the connection formula (33) for the Painlevé transcendent (32) and a numerical exploration of the constant factor.
In the same year Forrester [24] confirmed this to be true for α ∈ N by sketching an argument that, basically, uses the idea underlying the multidimensional Laplace method in the proof of Thm.2.2.So, Corollary A.1 can be used to spell out the details there, and for the generalization to β-ensembles sketched in [25, p. 608].
Table 2.The Stirling-type formula (5) vs. Regev's formula (29) for l = 5.The exact values of #{σ : Ln(σ) 5} = n! • P(Ln 5) were computed using the holonomic 4-term recursion w.r.t.n (see [7, p. 468]).Note that the Stirling-type formula gives 4 correct digits already for n = 160; Regev's formula starts to deliver 4 correct digits only for n as large as n ≈ 2 • 10 5 , where one is just about to enter the regime l n 1/4 .Including the Gaussian factor (28) helps to improve the accuracy a little, relaxing the bound on l to about l 2n To avoid notational clutter, we suppress the index l from the generating function f l , its auxiliaries a l , b l and from the radius r l,n .Solving a(r n ) = n yields, by (18), the expansion ( 27) which suggests to plug its leading order term r * n := n 2 /l 2 into (15).Thm.2.2 gives, as n → ∞, ).The Gaussian term in (15) has thus the expansion (28) exp which indicates that we can expect r * n to deliver a quality of approximation comparable to r n (which corresponds to using the Stirling-type formula) only if l n 1/4 ; see Table 2 for an illustrative example.Altogether Hayman's normal approximation (15) gives, choosing r = r * n , as n → ∞.Wrapping up by using Stirling's formula in the form we have thus given a new proof of Regev's formula [43, Eq. (4.5. 2)]: Remark 2.4.The fixed l asymptotics (29) was first proved by Regev [43] in 1981, cf.[48,Thm. 7].His delicate and rather long 21 proof proceeds, first, by identifying the leading 21 Though Regev studies, with a real parameter β > 0, the more general combinatorial sums this generality adds only marginally to the complexity of his proof.In its final step he refers to the same instance of Selberg's integral, cf.[2, Eq. (2.5.11)], that we have used to obtain (25) in the specific case β = 2.
contributions to the finite sum (1) using Stirling's formula, and then, after trading exponentially decaying tails (the basic idea of Laplace's method), by approximating the sum by a multidimensional integral which, finally, leads to the evaluation of Selberg's integral (25).

Numerical Evaluation of the Generating Function and its Auxiliaries
The numerical evaluation of the Stirling-type formula (5) requires the evaluation of the generating function f l (r) and its auxiliaries a l (r), b l (r) for real r > 0. This will be based on the representation (7).That is to say, by writing for the functions from random matrix theory, we obtain As is common in the discussion of the LUE, we generalize this by replacing l ∈ N with a real parameter α > −1.Dropping the index altogether we write, briefly, just g(s), v(s), and u(s).
subject to the following initial condition, which is consistent with (26) for α = l ∈ N: A numerical integration of the initial value problem gives approximations to v(s) and v (s), thus also to u(s) = sv (s).As explained in [9] a direct numerical integration has stability issues as the solution v is a separatrix solution of the σ-Painlevé-III equation.It is therefore advisable to solve the differential equation numerically as an asymptotic boundary problem by supplementing the initial condition by its connection formula, that is the corresponding expansion for x → ∞: Note the consistency with (24) for α = l ∈ N; for general α > −1 this connection formula was conjectured by Tracy and Widom [52, Eq. (3.1)], a rigorous proof is given in [19], cf.Remark 2.3.

Compiling a table of exact rational values. As observed recently by Forrester
and Mays [27, Sec.4.2], substituting a truncated power series expansion of v(x) into the σ-Painlevé-III equation ( 31) is a comparatively cost-efficient way 22 to compile a table of the exact rational values of the distribution P(L n l); they report to have done so up to n = 700.We note that instead of dealing directly with (31) in this fashion, it is of advantage to use an equivalent third-order differential equation belonging to the Chazy-I class, 23 namely (34)  2. For l > 5 the polynomial coefficients quickly become unwieldy, though. 23In fact, this equation is obtained as the particular choice c1 which is obtained from differentiating (31) w.r.t.x and dividing the result by 2x 2 v 2 l .Note that the Chazy-I equation ( 34) is linear in the highest order derivative of v l and quadratic in the lower orders, whereas the σ-Painlevé-III equation ( 31) is quadratic in the highest order derivative and cubic in the lower orders.Therefore, substituting the expansion into the Chazy-I equation ( 34) yields a much simpler recursive formula for the a n , n = l+1, . . .: uniquely determining the coefficients a n from the initial value (32), that is, from It is now a simple matter of truncated power series calculations in a modern computer algebra system to expand the generating function itself, Avoiding the overhead of reducing fractions and computing common denominators in exact rational arithmetic, we have used significance arithmetic with 2.5 log 10 (1000!) = 6420 digits and subsequent rational reconstructions to compile a table 24 of P(L n = l), 1 l n, up to n = 1000 (in just about 1.5 hours CPU time using one core of a 3GHz Xeon server).This table is used in Figs. 3 and 6 as well as Sect.4.3 (note that Table 1 could have been compiled with the values for up to n = 36 that were tabulated in the work of Baer and Brock [3]).

3.3.
Evaluation in terms of Bessel kernel determinants and traces.In [10] the author has shown that Nyström's method for integral equations can be generalized to the numerical evaluation of Fredholm determinants.Thus, as advocated in [9], there is a stable and efficient numerical method to directly address the representation derived by Forrester [23] in 1993, in terms of the Bessel kernel This numerical method was extended in [14,Appendix] to the evaluation of general terms involving determinants, traces, and resolvents of integral operators.The evaluation of the auxiliary functions v(s), u(s), as defined in (30), is thus facilitated by the following theorem.
This way we can keep the space L 2 (0, 1) fixed while the kernels become dependent on the parameter s > 0; in particular, then, there is no need to distinguish in notation between kernels and their induced integral operators.Now, if we denote differentiation w.r.t. to the parameter s by a dot, we get Ks (x, y) = K (sx, sy) = s −1 K s (x, y), and thus, by [14, Lemma 1], the logarithmic derivative which proves the asserted formula for v(s) = −sg (s)/g(s).Since, cf.[52, Eq. (2.4)], d ds we get by the linearity of the trace which finishes the proof after a back-transformation to L 2 (0, s).
For the Bessel kernel (36b) at hand we get the derived kernels We observe that both, K and K , induce integral operators of finite rank, namely where we have put Hence the results of Theorem 3.1 simplify considerably: first, we obtain 25 next, by observing we get, because of symmetry and linearity, (36d) Both formulae for the auxiliary functions u and v can now be easily implemented in the author's Matlab toolbox for Fredholm determinants (which provides also commands to evaluate traces and inner products of general operator terms including resolvents; cf.[9,10] and [14,Appendix]).Since all the numerical evaluations come with an estimate of the 25 Correcting an obvious typo, (36c) is precisely [12, Eq. ( 6)].As it was noted there, (36c) can also be found, though not explicitly, in [52].On the other hand, formula (36d) seems to be new.
(absolute) error there, the implied approximation errors in computing the generating function f l (r) and its auxiliaries a l (r) and b l (r) can straightforwardly be assessed.
Remark 3.1.A result similar to Theorem 3.1 holds for smooth integral kernels K(x, y), with sufficient decay at ∞, which induce integral operators K on L 2 (s, ∞) for all s ∈ R.Here we define the derived kernel as and get, if g(s) = det(I − K)| L 2 (s,∞) > 0 for all s ∈ R, the logarithmic derivative The proof goes by considering K s (x, y) = K(s + x, s + y) and transforming (s, ∞) to (0, ∞) by a shift.As an example, the Tracy-Widom distribution F 2 (s) used in ( 3) is known to be given in terms of the Airy kernel determinant [23], x − y .
Here we have K (x, y) = − Ai(x) Ai(y), i.e., K = − Ai ⊗ Ai, and thus The last formula was used for the calculations shown in Table 3.In the same manner we get as well as similar formulae for higher order derivatives of F 2 . 26  3.4.Implementation details.First, by uniqueness, solving a l (r) = n for r = r l,n can easily be accomplished by an iterative solver.In view of the left panel in Fig. 5 and the expansion (27) we take as initial guess r 0 := max(n, (n/l) 2 + n/2).
Second, the numerical evaluation of the Stirling-type formula (5) for larger values of n requires to avoid severe overflow of intermediate terms.Based on the representations in (30), and by rearranging terms, we can write (5) equivalently as follows: 27 where g l , v l and u l are evaluated at s = 4r l,n and there is In IEEE hardware arithmetic we take the definition of τ n until n = 100 and only switch to the shown Stirling expansion for larger n-thus seamlessly providing full accuracy. 28This allows us to approximate the PDF P(L n = l) near its mode for up to n = 10 12 and larger.For accurate tails, such as in Table 2, we have to resort to higher precision arithmetic, though. 26Though the inner products in (37b) and (37c) appear in the work of Tracy and Widom [52, Eq. (1. 3)], the formula for F 2 (s) is not given there. 27We have to stabilize the numerical evaluation of the expression h − log(1 + h) for small h := v l /n ≈ 0. This is done, first, by using h-log1p(h) and, second, by switching to a suitable Taylor expansion for very small h. 28In fact, nothing of substance would change if we just replaced τn by 1 since the thus committed error would be in the same ballpark as the one of the Stirling-type formula (5) itself.We did not bother to do so, though.Left: CDF errors rescaled by n 2/3 , horizontal axis is t = (l − 2 √ n)/n 1/6 .The solid line is a polynomial F2,2(t) of degree 48 fitted to all the 152 data points with −7.5 t 9.5; it approximates F2,2(t) in that interval.Right: PDF errors rescaled by n 5/6 , horizontal axis is t   [27] have recently initiated the study of finite size corrections to the random matrix limit (3), which is (39) P(L n l) = F 2 (t l ) + o(1), t l := l − 2 √ n n 1/6 , as n → ∞, uniformly in l ∈ N. We will refine their study by using the much more accurate and efficient Stirling-type formula (5) instead.Looking at the error δ 0 (n) := max l∈{1,...,n} for n up to 1000 (see the red crosses in the left panel of Fig. 3 in a double logarithmic scaling) suggest that δ 0 (n) ≈ c 1 n −1/3 + c 2 n −2/3 and yields the conjecture ( 40) ) for some function F 2,1 (t).Numerically, the conjecture has been convincingly checked against the data obtained by the Stirling-type formula (5) for n = 10 6 , n = 10 8 , and n = 10 10 ; see the left panel of Fig. 4. We have fitted a polynomial F2,1 (t) of degree 64 to the 836 data points obtained for n = 10 10 in the interval −8 t 10, thus approximating the putative function F 2,1 (t) there.
Remark 4.2.To validate conjecture (41) against numerical data, obtained by replacing F 2,1 (t) by the approximation F2,1 (t), we have to be careful with an effective choice of n, though.
On the one hand, as a perturbation of F 2,2 (t) the error of about 10 −4 in F2,1 , as estimated in Remark 4.1, would get amplified by n 1/3 .On the other hand, an extrapolation of the errors displayed in the left panel of Fig. 3 shows that the Stirling-type formula would induce an additional perturbation of size ≈ (0.031n −2/3 + 0.058n −1 )n 2/3 .The sweet spot of both perturbations combined is at n ≈ 1.4 • 10 4 with a minimum error of about 3.6 • 10 −2 .Thus, we better stay with the tabulated exact values of the distribution of L n up to n = 1000, which restricts the size of the perturbation to just less than the order of n 1/3 10 −4 | n=1000 = 10 −3 .
Thus, staying with the tabulated values of the distribution of L n for n = 250, n = 500, and n = 1000 we get a convincing picture; see the left panel of Fig. 6.We have fitted a polynomial F2,2 (t) of degree 48 to all of the 152 data points in the interval −7.5 t 9.5, approximating the putative function F 2,2 (t) there. .
Note the shift by 1/2 in the numerator of tl as compared to t l (defined in (39)).There is compelling numerical evidence for the expansion (42a); see the right panels of Figs. 4 and 6.where the integral over R n was evaluated by spherical symmetry and the resulting gamma integral by Eq. ( 51); ω n−1 denotes the surface area of the sphere S n−1 .
Lemma A.3.There is η → 0 for → 0 such that, for 0 < 0 with 0 sufficiently small and z ∈ C with |arg z| if we define the term given by the large bracket to be η .

Figure 1
Figure1.Discrete distribution of Ln for n = 10 5 near its mode vs. the random matrix limit given by the leading order terms in (3) and (42a) (solid red line); here and in the figures below, discrete distributions are shown as blue bars centered at the integers.Left: CDF P(Ln l); right: PDF P(Ln = l).The discrete distributions were computed using the Stirling-type formula(5), with additive errors estimated to be smaller than 10 −5 , cf.Fig.3, which is well below plotting accuracy.

Figure 2 .
Figure2.Display of the notable inaccuracy of the random matrix limit (3) for n = 1000 (see the contrast with Fig.1for n = 10 5 ).The discrete distribution of Ln is shown near its mode vs. the random matrix limit given by the leading order terms in (9) and (42a) (solid red line).Left: CDF P(Ln l); right: PDF P(Ln = l).The exact values of the distribution of Ln and their approximation by the Stirling-type formula (5) differ just by additive errors of the order 10 −4 (see Fig.3), which is well below plotting accuracy.

3. 1 .
Evaluation in terms of σ-Painlevé-III.The work of Tracy and Widom [52] shows g(s) = exp − where v(s) satisfies a Jimbo-Miwa-Okamoto σ-form of the Painlevé-III equation (related to the Hamiltonian formulation PIII in the work of Okamoto; cf.[25, §8.2]), i.e., the nonlinear second order differential equation

Figure 6 .
Figure 6.Rescaled differences between the distributions of Ln and their expansions truncated after the first finite size correction term-see(9) for the CDF resp.(42a) for the PDF; the data points have been calculated with the polynomial F2,1(t) from Fig.4for the exact values of the distribution for n = 250 (red +), n = 500 (green •), n = 1000 (blue •).Left: CDF errors rescaled by n 2/3 , horizontal axis is t = (l − 2 √ n)/n 1/6 .The solid line is a polynomial F2,2(t) of degree 48 fitted to all the 152 data points with −7.5 t 9.5; it approximates F2,2(t) in that interval.Right: PDF errors rescaled by n 5/6 , horizontal axis is t = (l − 1 2 − 2

4 .
First and Second Finite Size Corrections to the Random Matrix Limit 4.1.The CDF of the distribution of L n .Based on data from Monte-Carlo simulations, Forrester and Mays