A Stirling-Type Formula for the Distribution of the Length of Longest Increasing Subsequences

Abstract

The discrete distribution of the length of longest increasing subsequences in random permutations of n integers 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 scaled largest level of the large matrix limit of GUE. As a numerical approximation, however, this asymptotics is inaccurate for small n and has a slow convergence rate, conjectured to be just of order $n^{-1/3}$. Here, we suggest a different 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^{-2/3}$, for n as large as $10^{12}$, thus closing the gap between a table of exact values (compiled for up to $n=1000$) and the random matrix limit. Being much more efficient and accurate than Monte Carlo simulations, the Stirling-type formula allows for a precise numerical understanding of the first few finite size correction terms to the random matrix limit. From this we derive expansions of the expected value and variance of the length, exhibiting several more terms than previously put forward.

1 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(\sigma )$ of longest increasing subsequences of permutations $\sigma$ on the set $[n] = \{1,2,\ldots ,n\}$ and of the behavior of their distribution in the limit $n\rightarrow \infty$. Here, $L_n(\sigma )$ is defined as the maximum of all k for which there are $1\leqslant i_1< i_2< \cdots < i_k \leqslant n$ such that $\sigma _{i_1}< \sigma _{i_2}< \cdots < \sigma _{i_k}$. Writing permutations in the form $\sigma = (\sigma _1\, \sigma _2\,\cdots \, \sigma _n)$ we get, e.g., $L_9(\sigma )= 5$ for $\sigma = (4\,\textbf{1}\, \textbf{2}\, 7\, 6\, \textbf{5}\, \textbf{8}\,\textbf{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) ${\mathbb {P}}( L_n \leqslant l)$ and probability distribution (PDF) ${\mathbb {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]):

$$\begin{aligned} {\mathbb {P}}( L_n \leqslant l) = \frac{1}{n!} \sum _{\lambda \,\vdash n\,:\, l_\lambda \leqslant l} d_\lambda ^2. \end{aligned}$$

(1)

Here $\lambda \,\vdash n$ denotes an integer partition $\lambda _1\geqslant \lambda _2 \geqslant \cdots \geqslant \lambda _{l_\lambda }> 0$ of $n=\sum _{j=1}^{l_\lambda } \lambda _j$ and $d_\lambda$ is the number of standard Young tableaux of shape $\lambda$, as given by the hook length formula. By generating all partitions $\lambda \,\vdash n$, in 1968 Baer and Brock [3] computed tables of ${\mathbb {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\times 10^9$). This quickly becomes infeasible,¹ as $p_n$ is already as large as $2.3\times 10^{14}$ for $n=250$. Another use of the combinatorial methods is the approximation of the distribution of $L_n$ by Monte Carlo simulations [3, 40]: one samples random permutations $\sigma$ and calculates $L_n(\sigma )$ by the Robinson–Schensted correspondence.²

Analytic 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 ${\mathbb {P}}( L_n \leqslant l)$, when considered as a sequence of n with the length l fixed:

$$\begin{aligned} f_l(z) := \sum _{n=0}^\infty \frac{{\mathbb {P}}( L_n \leqslant l)}{n!}z^n \qquad (z\in {\mathbb {C}},\, l \in {\mathbb {N}}). \end{aligned}$$

(We note that $f_l$ is an entire function of exponential type.) In fact, Gessel [30, p. 280] obtained in 1990 the explicit representation

$$\begin{aligned} f_l(z^2) = D_l(z),\qquad D_l(z) :=\det _{j,k=1}^l I_{j-k}(2z), \end{aligned}$$

(2)

in terms of a Toeplitz determinant of the modified Bessel functions $I_m$, $m\in {\mathbb {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 theorem³ to induce from that limit an asymptotics of the coefficients, Baik, Deift and Johansson [4] succeeded 1999 in establishing⁴

$$\begin{aligned} {\mathbb {P}}( L_n \leqslant l) = F_2\left( \frac{l-2\sqrt{n}}{n^{1/6}}\right) + o(1) \qquad (n\rightarrow \infty ), \end{aligned}$$

(3)

uniformly in $l\in {\mathbb {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 $\beta =2$ (that is, the probability that in the soft-edge scaling limit of the Gaussian unitary ensemble (GUE) 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.

Fig. 1

Discrete distribution of $L_n$ 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 ${\mathbb {P}}( L_n \leqslant l)$; right: PDF ${\mathbb {P}}( L_n = 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 (Color figure online)

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 \ll 2\sqrt{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$.

Table 1 The values of ${\mathbb {P}}( L_n \leqslant l)$ and various of its approximations for $n=20$

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\times 10^6$ Monte Carlo trials for this n; the error of such a simulation is expected to be of the order $1/\sqrt{T} \approx 4\cdot 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

$$\begin{aligned} f(z) = \sum _{n=0}^\infty a_n z^n \qquad (z \in {\mathbb {C}}) \end{aligned}$$

is an entire function with positive coefficients $a_n$ and consider the real auxiliary functions

$$\begin{aligned} a(r) = r \frac{\textrm{d}}{\textrm{d}r} \log f(r),\qquad b(r) = r \frac{\textrm{d}}{\textrm{d}r} a(r)\qquad (r>0). \end{aligned}$$

If f is H-admissible, then for each $n\in {\mathbb {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⁵ holds true:

$$\begin{aligned} a_n = \frac{f(r_n)}{r_n^n \sqrt{2\pi \, b(r_n)}}(1+ o(1)) \qquad (n\rightarrow \infty ). \end{aligned}$$

(4)

We observe that the error is multiplicative here. In Theorem 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 ${\mathbb {P}}(L_n \leqslant l)/n!$. Since the error is multiplicative, nothing changes if we multiply the approximation by n! and we get

$$\begin{aligned} {\mathbb {P}}(L_n \leqslant l) = \frac{n! \cdot f_l(r_{l,n})}{r_{l,n}^n \sqrt{2\pi \, b_l(r_{l,n})}}(1+ o(1)) \qquad (n\rightarrow \infty ), \end{aligned}$$

(5)

where we have labeled all quantities when applied to $f=f_l$ by an additional index l. An approximation to ${\mathbb {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, 2 and 3, as well as in Table 2.

Fig. 2

Display of the notable inaccuracy of the random matrix limit (3) for $n=1000$ (see the contrast with Fig. 1 for $n=10^5$). The discrete distribution of $L_n$ 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 ${\mathbb {P}}( L_n \leqslant l)$; right: PDF ${\mathbb {P}}( L_n = l)$. The exact values of the distribution of $L_n$ 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 (Color figure online)

Fig. 3

Maximum absolute (i.e., additive) errors of various approximations to (left panel) the CDF ${\mathbb {P}}(L_n\leqslant l)$ and (right panel) the PDF ${\mathbb {P}}(L_n= l)$ in a double logarithmic scaling, based on tabulated exact values up to $n=1000$, cf. Sect. 3.2; solid lines are fits of the form $c_1 n^{-\alpha _1} + c_2 n^{-\alpha _2} + c_3 n^{-\alpha _3}$ (CDF) and $n^{-1/6}\cdot (c_1 n^{-\alpha _1} + c_2 n^{-\alpha _2} + c_3 n^{-\alpha _3})$ (PDF) to the points in display. Red $+$: error of leading order terms (random matrix limit) in (9), (42a); $\alpha = (\tfrac{1}{3},\tfrac{2}{3},1)$. Green $\circ$: error of expansions (9), (42a) truncated after the first finite size correction term; $\alpha =(\tfrac{2}{3}, 1, \tfrac{4}{3})$, where $F_{2,1}$ has been approximated as in Fig. 4. Blue $\bullet$: error of the Stirling-type formula (5); $\alpha = (\tfrac{2}{3}, 1, \tfrac{4}{3})$ (Color figure online)

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 ${\mathbb {P}}(L_n\leqslant 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 \in {\mathbb {N}}_0 :=\{0,1,2,3,\ldots \}$ to be a further independent random variable with Poisson distribution of intensity $r\geqslant 0$, we see that

$$\begin{aligned} {\mathbb {P}}(L_{N_r} \leqslant l) = \sum _{n=0}^\infty \frac{e^{-r} r^n}{n!} {\mathbb {P}}(L_n \leqslant l) = e^{-r} f_l(r) \end{aligned}$$

(6)

is, for fixed $r\geqslant 0$, the cumulative probability distribution of the composite discrete random variable $L_{N_r}$. On the other hand, for fixed $l\in {\mathbb {N}}$, also $e^{-r} f_l(r)$ turns out to be a probability distribution w.r.t. the continuous variable $r\geqslant 0$: specifically, in terms of precisely the Toeplitz determinant (2), Forrester and Hughes [26, Eq. (3.33)] arrived in 1994 at the representation

$$\begin{aligned} E^{\text {(hard)}}_2(0; [0,4r],l) = e^{-r} f_l(r). \end{aligned}$$

(7)

Here, $E^{\text {(hard)}}_2(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\geqslant 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 $\sigma$-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.

Table 2 The Stirling-type formula (5) versus Regev’s formula (29) for $l=5$

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\in \{1,2,\ldots ,n\}$) in approximating the distribution ${\mathbb {P}}(L_n\leqslant 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^{-\alpha _1} + c_2 n^{-\alpha _2} + c_3 n^{-\alpha _3}$ with simple triples $\alpha =(\alpha _1, \alpha _2, \alpha _3)$ of rationals strongly suggests that, uniformly in $l \in \{1,2,\ldots ,n\}$ as $n\rightarrow \infty$,

$$\begin{aligned} {\mathbb {P}}( L_n \leqslant l)&= F_2\left( \frac{l-2\sqrt{n}}{n^{1/6}}\right) + O(n^{-1/3}), \end{aligned}$$

(8a)

$$\begin{aligned} {\mathbb {P}}( L_n \leqslant l)&= \frac{n! \cdot f_l(r_{l,n})}{r_{l,n}^n \sqrt{2\pi \, b_l(r_{l,n})}} + O(n^{-2/3}). \end{aligned}$$

(8b)

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 conjecture⁶

$$\begin{aligned} {\mathbb {P}}(L_n\leqslant l) = F_2(t_l) + n^{-1/3} F_{2,1}(t_l) + n^{-2/3} F_{2,2}(t_l) + O(n^{-1}), \quad t_l := \frac{l-2\sqrt{n}}{n^{1/6}}, \end{aligned}$$

(9)

as $n\rightarrow \infty$, uniformly in $l\in {\mathbb {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.⁷

Fig. 4

Rescaled differences between the distributions of $L_n$ 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^\text {th}$ is displayed) have been calculated using the Stirling-type formula (5) for $n=10^6$ (red $+$), $n=10^8$ (green $\circ$), $n=10^{10}$ (blue $\bullet$). Left: CDF errors rescaled by $n^{1/3}$, horizontal axis is $t=(l-2\sqrt{n})/n^{1/6}$, cf. [27, Fig. 7] for a similar figure with data from Monte Carlo simulations for $n=2\times 10^4$ and $n=10^5$. The solid line is a polynomial $\tilde{F}_{2,1}(t)$ of degree 64 fitted to the 836 data points for $n=10^{10}$ with $-8\leqslant t \leqslant 10$; it approximates $F_{2,1}(t)$ in that interval. Right: PDF errors rescaled by $n^{1/2}$, horizontal axis is $t=(l-\frac{1}{2}-2\sqrt{n})/n^{1/6}$. The solid line displays the function ${{\tilde{F}}}_{2,1}'(t)+F_2'''(t)/24$ as an approximation of $F_{2,1}'(t)+F_2'''(t)/24$, with the polynomial $\tilde{F}_{2,1}(t)$ taken from the left panel. The dotted line shows the term ${{\tilde{F}}}_{2,1}'(t)$ only (Color figure online)

We note that a corresponding expansion⁸ 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\rightarrow \infty$ for bounded $t_l^*$:

$$\begin{aligned} {\mathbb {P}}(L_{N_r} \leqslant l) = F_2(t_l^*) + r^{-1/3} F_{2,1}^*(t_l^*) + O(r^{-2/3}), \quad t_l^* := \frac{l-2\sqrt{r}}{r^{1/6}}, \end{aligned}$$

(10a)

with the explicit functional form (identified by means of Painlevé representations)

$$\begin{aligned} F_{2,1}^*(t) = -\frac{1}{10} \left( F_2''(t) + \frac{t^2}{6} F_2'(t)\right) . \end{aligned}$$

(10b)

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 order $O(n^{-2/3})$, we are led to the conjecture

$$\begin{aligned} F_{2,1}(t) = F_{2,1}^*(t) - \frac{1}{2} F_2''(t) = -\frac{1}{10} \left( 6 F_2''(t) + \frac{t^2}{6} F_2'(t)\right) . \end{aligned}$$

(11)

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 $F_{2,1}(t)$ (as well as one for $F_{2,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

$$\begin{aligned}{} & {} {\mathbb {E}}(L_n) = 2\sqrt{n} + \mu _0 n^{1/6} + \frac{1}{2} + \mu _1 n^{-1/6} + \mu _2 n^{-1/2} + O(n^{-5/6}),\nonumber \\{} & {} \mu _0 = \int _{-\infty }^\infty t \, F_2'(t)\,\textrm{d}t = -1.77108\,68074\cdots ,\nonumber \\{} & {} \mu _1 = \int _{-\infty }^\infty t \, F_{2,1}'(t)\,\textrm{d}t = 0.06583\,238\cdots , \nonumber \\{} & {} \mu _2 = \int _{-\infty }^\infty t \, F_{2,2}'(t)\,\textrm{d}t = 0.26122\,27\cdots . \end{aligned}$$

(12)

Similarly, we will derive in Sect. 4.4 an expansion of the variance of $L_n$ of the form

$$\begin{aligned} \textrm{Var}(L_n)&= \nu _0 n^{1/3} + \nu _1 + \nu _2 n^{-1/3} + O(n^{-2/3}),\nonumber \\ \nu _0&= \int _{-\infty }^\infty t^2 F_2'(t)\,\textrm{d}t - \mu _0^2 = 0.81319\,47928\cdots ,\nonumber \\ \nu _1&= \int _{-\infty }^\infty t^2 F_{2,1}'(t)\,\textrm{d}t +\frac{1}{12} -2\mu _0\mu _1 = -1.20720\,507\cdots ,\nonumber \\ \nu _2&= \int _{-\infty }^\infty t^2 F_{2,2}'(t)\,\textrm{d}t -\mu _1^2 - 2\mu _0\mu _2 = 0.56715\,6\cdots . \end{aligned}$$

(13)

The values of $\mu _0$ and $\nu _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 $\mu _0n^{1/6}$ and of (13) up to $\nu _0n^{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,⁹ which is difficult to be verified in practice and therefore seldom directly used. Instead, we start by collecting some useful facts and criteria from Hayman’s original paper [33]:¹⁰

Theorem 2.1

(Hayman 1956) Let $f(z) = \sum _{n=0}^\infty a_nz^n$ and g(z) be entire functions and let p(z) denote a polynomial with real coefficients.

I. If f is H-admissible, then:

• a. $f(r) > 0$ for all sufficiently large $r>0$, so that in particular the auxiliary functions¹¹

  $$\begin{aligned} a(r) = r \frac{\textrm{d}}{\textrm{d}r} \log f(r),\qquad b(r) = r \frac{\textrm{d}}{\textrm{d}r} a(r) \end{aligned}$$

  (14)

  are well defined there;

• 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)\rightarrow \infty$ as $r\rightarrow \infty$; in particular, for large integers n there is a unique $r_n>0$ that solves $a(r_n)=n$, it is $r_n\rightarrow \infty$ as $n\rightarrow \infty$;

• c. if the coefficients $a_n$ of f are all positive, then I.b holds for all $r>0$;

• d. as $r\rightarrow \infty$, uniformly in $n \in {\mathbb {N}}_0$,

  $$\begin{aligned} \frac{a_n r^n}{f(r)} = \frac{1}{\sqrt{2\pi b(r)}}\left( \exp \left( -\frac{(n-a(r))^2}{2b(r)}\right) + o(1)\right) . \end{aligned}$$

  (15)

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¹² with, for some $\delta >0$, at most finitely many zeros in the sector $|\arg z|\leqslant \frac{\pi }{2}+\delta$ and satisfies I.a such that $b(r)\rightarrow \infty$ as $r\rightarrow \infty$, 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\geqslant 0$, Eq. (15) has an interesting probabilistic content:¹³ as a distribution in the discrete variable $n\in {\mathbb {N}}_0$, the Boltzmann¹⁴ 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.

Fig. 5

Left: the auxiliary function $a_5(r)$ (blue solid line) associated with the generating function $f_5(z)$, together with the asymptotics $a_5(r)= r + O(r^6)$ as $r\rightarrow 0$ (green dotted line) and $a_5(r)= 5 r^{1/2} + O(1)$ as $r\rightarrow \infty$ (red dashed line). Right: Illustration of the approximation (15) of the Boltzmann probabilities (blue bars) associated with the generating function $f_5(z)$ for intensity $r_{5,15}\approx 18.23$, cf. the notation in (5). The normal distribution (red solid line) has mean $a_5(r_{5,15}) = 15$ (cf. the left panel) and variance $b_5(r_{5,15})\approx 10.80$ (Color figure online)

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=\infty$. 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, \tau , \nu > 0$ such that there holds, for the principal branch of the power function and for each $0 < \delta \leqslant \frac{\pi }{2}$, the asymptotic expansion¹⁵

$$\begin{aligned} D(z):=f(z^2) = cz^{-\nu } e^{\tau z}(1 + O(z^{-1}))\qquad (z\rightarrow \infty ,\; |\!\arg z|\leqslant \tfrac{\pi }{2}-\delta ), \end{aligned}$$

(16)

then f is H-admissible. For $r\rightarrow \infty$ the associated auxiliary functions a(r) and b(r) satisfy

$$\begin{aligned} a(r) = \frac{\tau }{2}r^{1/2} - \frac{\nu }{2} + O(r^{-1/2}), \qquad b(r) = \frac{\tau }{4}r^{1/2} + O(r^{-1/2}), \end{aligned}$$

(17)

and the solution $r_n$ of $a(r_n)=n$ satisfies

$$\begin{aligned} r_n =\frac{4n}{\tau ^{2}} (n+\nu ) + O(1) \qquad (n\rightarrow \infty ). \end{aligned}$$

(18)

Proof

The expansion (16) is equivalent to

$$\begin{aligned} f(z) = cz^{-\nu /2} e^{\tau z^{1/2}}(1 + O(z^{-1/2}))\qquad (z\rightarrow \infty ,\; |\!\arg z|\leqslant \pi -2\delta ), \end{aligned}$$

(19)

which readily implies:

• since $\delta >0$ is arbitrary, f has the Phragmén–Lindelöf indicator

  $$\begin{aligned} {{\,\mathrm{\text {lim sup}}\,}}_{r\rightarrow \infty } \frac{\log |f(r e^{i\theta })|}{r^{1/2}} = \tau \cos (\theta /2) \qquad (-\pi< \theta < \pi ), \end{aligned}$$

  so that f has order $\tfrac{1}{2}$ and type $\tau$, hence genus zero;

• for sufficiently large $R = R_{\delta }>0$, there are no zeros z of f with

  $$\begin{aligned} |z|\geqslant R,\qquad |\!\arg z| \leqslant \pi -2\delta . \end{aligned}$$

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) \rightarrow \infty$ as $r\rightarrow \infty$. Thus, all the assumptions of Theorem 2.1.III are satisfied and f is shown to be H-admissible. $\square$

2.2 Singularity Analysis of the Generating Function at $z=\infty$

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

$$\begin{aligned} D_l(z) =\det _{j,k=1}^l I_{j-k}(2z), \end{aligned}$$

(20)

suggests to start with the expansions (valid for all $0<\delta \leqslant \frac{\pi }{2}$, see [41, p. 251])

$$\begin{aligned} I_m(z)&\sim \frac{e^{z}}{(2\pi z)^{1/2}} \sum _{n=0}^\infty (-1)^n \frac{A_n(m)}{z^n}\qquad \qquad (z\rightarrow \infty ,\; |\!\arg z|\leqslant \tfrac{\pi }{2}-\delta ), \end{aligned}$$

(21a)

$$\begin{aligned} A_n(m)&= \frac{(4m^2-1^2)(4m^2-3^2)\cdots (4m^2-(2n-1)^2)}{n! \,8^n}. \end{aligned}$$

(21b)

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,\ldots ,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¹⁶ the expansions

$$\begin{aligned} D_1(z)&= \frac{e^{2 z}}{2 \pi ^{1/2} z^{1/2}}(1+O(z^{-1})),&D_2(z)&= \frac{e^{4 z}}{8 \pi z^2}(1+O(z^{-1})),\\ D_3(z)&= \frac{e^{6z}}{32 \pi ^{3/2} z^{9/2}}(1+O(z^{-1})),&D_4(z)&= \frac{3 e^{8 z}}{256 \pi ^2 z^8}(1+O(z^{-1})),\\ D_5(z)&= \frac{9 e^{10 z}}{1024 \pi ^{5/2}z^{25/2}}(1+O(z^{-1})),&D_6(z)&= \frac{135 e^{12 z}}{8192 \pi ^3z^{18}}(1+O(z^{-1})),\\ D_7(z)&= \frac{6075 e^{14 z}}{65536 \pi ^{7/2}z^{49/2}}(1+O(z^{-1})),&D_8(z)&= \frac{1913625 e^{16 z}}{1048576 \pi ^4z^{32}}(1+O(z^{-1})). \end{aligned}$$

All of them, inherited from (21), are valid as $z\rightarrow \infty$ while $|\!\arg z|\leqslant \tfrac{\pi }{2}-\delta$ 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

$$\begin{aligned} D_l(z) = c_l \frac{e^{2lz}}{(4\pi z)^{l/2} (2z)^{l(l-1)/2}} (1+ O(z^{-1})) \qquad \qquad (z\rightarrow \infty ,\; |\!\arg z|\leqslant \tfrac{\pi }{2}-\delta ) \end{aligned}$$

and observe

$$\begin{aligned} c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8, \ldots = 1, 1, 2, 12, 288, 34560, 24883200, 125411328000,\ldots \;. \end{aligned}$$

Consulting the OEIS¹⁷ (On-Line Encyclopedia of Integer Sequences) suggests the coefficients to be generally of the form

$$\begin{aligned} c_l = 0!\cdot 1!\cdot 2! \, \cdots \, (l-1)!\,. \end{aligned}$$

Though this is very likely to hold for all $l\in {\mathbb {N}}$—a fact that would at once yield the H-admissibility 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

$$\begin{aligned} I_m(2z) = \frac{1}{2\pi } \int _{-\pi }^{\pi }e^{2z\cos \theta } e^{-im \theta }\,\textrm{d}\theta \qquad (z\in {\mathbb {C}}, m \in {\mathbb {Z}}), \end{aligned}$$

(22)

we represent the Toeplitz determinant $D_l(z)$ in terms of a multidimensional integral and study the limit $z\rightarrow \infty$ 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 \theta )$ and a classical formula of Szegő’s [50, p. 493] from 1915, thus gives, without further calculation, the integral representation¹⁸

$$\begin{aligned} D_l(z) = \frac{1}{(2\pi )^l \,l!} \int _{-\pi }^{\pi }\cdots \int _{-\pi }^{\pi } e^{2z \sum _{j=1}^l \cos \theta _j}\cdot \big |\Delta (e^{i\theta _1},\ldots ,e^{i\theta _l})\big |^2 \,\textrm{d}\theta _1 \cdots \textrm{d}\theta _l, \end{aligned}$$

(23)

where

$$\begin{aligned} \Delta (w_1,\ldots ,w_l) := \prod _{j > k} (w_j- w_k) \end{aligned}$$

denotes the Vandermonde determinant of the complex numbers $w_1,\ldots ,w_l$.

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

$$\begin{aligned} D_n(z) = {\mathbb {E}}_{U\in U(l)} e^{z {{\,\mathrm{\text {tr }}\,}}(U+U^*)}, \end{aligned}$$

where the expectation ${\mathbb {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

$$\begin{aligned} {\mathbb {P}}(L_n\leqslant l) = {\mathbb {E}}_{U\in U(l)} (|{{\,\mathrm{\text {tr }}\,}}U|^{2n}), \end{aligned}$$

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< \delta \leqslant \pi /2$ and $l\in {\mathbb {N}}$, there holds the asymptotic expansion

$$\begin{aligned}{} & {} D_l(z) = f_l(z^2) =\frac{0!\cdot 1!\cdot 2! \, \cdots \, (l-1)!\cdot e^{2l z}}{(2\pi )^{l/2}(2z)^{l^2/2}}(1+ O(z^{-1}))\\ {}{} & {} \qquad (z\rightarrow \infty ,\; |\!\arg z|\leqslant \tfrac{\pi }{2}-\delta ). \end{aligned}$$

Thus, by Lemma 2.1, the generating functions $f_l(z)$ are H-admissible and their auxiliary functions satisfy, as $r\rightarrow \infty$,

$$\begin{aligned} a_l(r) = lr^{1/2} - \tfrac{1}{4} l^2 + O(r^{-1/2}),\qquad b_l(r) = \tfrac{1}{2} l r^{1/2} + O(r^{-1/2}). \end{aligned}$$

(24)

Proof

We write (23) in the form

$$\begin{aligned} e^{-lz} D_l(z/2) = \frac{1}{(2\pi )^l \,l!} \int _{-\pi }^{\pi }\cdots \int _{-\pi }^{\pi } e^{-z \sum _{j=1}^l (1-\cos \theta _j)}\cdot \big |\Delta (e^{i\theta _1},\ldots ,e^{i\theta _l})\big |^2 \,\textrm{d}\theta _1 \cdots \textrm{d}\theta _l. \end{aligned}$$

The phase function of this multidimensional integrand, that is to say

$$\begin{aligned} S(\theta _1,\ldots ,\theta _l) := \sum _{j=1}^l (1-\cos \theta _j), \end{aligned}$$

takes it minimum $S(\theta _*)=0$ at $\theta _* = 0$ with the expansion $S(\theta _1,\ldots ,\theta _l) = \frac{1}{2}\theta ^T\theta + O(|\theta |^4)$ as $\theta \rightarrow 0$, where $|\cdot |$ denotes Euclidean length. Likewise we get for the non-exponential factor¹⁹

$$\begin{aligned} \big |\Delta (e^{i\theta _1},\ldots ,e^{i\theta _l})\big |^2&= \prod _{j>k} |e^{i\theta _j}-e^{i\theta _j}|^2 = \prod _{j>k} \left| i\theta _j-i\theta _j + O(|\theta |^2)\right| ^2 \\&= \Delta (\theta _1,\ldots ,\theta _l)^2 + O(|\theta |^{l(l-1)+1})\qquad (\theta \rightarrow 0), \end{aligned}$$

where the degree of the homogeneous polynomial $\Delta (\theta _1,\ldots ,\theta _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

$$\begin{aligned} e^{-lz} D_l(z/2) = c_l \frac{(2\pi )^{l/2} z^{-\frac{l+l(l-1)}{2}}}{(2\pi )^l} (1+O(z^{-1})) \qquad (z\rightarrow \infty ,\; |\!\arg z|\leqslant \tfrac{\pi }{2}-\delta ) \end{aligned}$$

with

$$\begin{aligned} c_l := \frac{1}{l!}\frac{1}{(2\pi )^{l/2}} \int _{{\mathbb {R}}^l} e^{- \theta ^T\theta /2}\cdot \big |\Delta (\theta _1,\ldots ,\theta _l)\big |^2 \,\textrm{d}\theta = 0!\cdot 1!\cdot 2! \, \cdots \, (l-1)!, \end{aligned}$$

(25)

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)]. $\square$

Remark 2.3

By (7), Theorem 2.2 implies

$$\begin{aligned} E^{\text {(hard)}}_2(0; [0,s],l) = \frac{0!\cdot 1!\cdot 2! \, \cdots \, (l-1)!}{(2\pi )^{l/2}} \cdot s^{-l^2/4} e^{-s/4 + l s^{1/2}}(1+O(s^{-1/2}))\qquad (s\rightarrow \infty ), \end{aligned}$$

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 $\alpha \in {\mathbb {C}}$ of the LUE with $\Re \alpha > -1$ instead of just $l \in {\mathbb {N}}$; the superfactorial factor $0!\cdot 1!\cdot 2! \, \cdots \, (l-1)!$ is then to be replaced by $G(1+\alpha )$, where G(z) is the Barnes G-function.²⁰

We complement the large r expansion (24) of the auxiliary functions with their expansions as $r\rightarrow 0^+$, which are simple consequences of elementary combinatorics.

Lemma 2.2

The auxiliary functions of the generating function $f_l$ satisfy, as $r\rightarrow 0^+$,

$$\begin{aligned} a_l(r) = r - \frac{r^{l+1}}{l! \cdot (l+1)!} + O(r^{l+2}),\qquad b_l(r) = r-\frac{r^{l+1}}{l!\cdot l!} + O(r^{l+2}). \end{aligned}$$

(26)

Proof

Because of $L_n \leqslant n$ and since there is just one permutation $\sigma$ with $L_n=n$, we get

$$\begin{aligned} {\mathbb {P}}(L_n \leqslant l) = {\left\{ \begin{array}{ll} 1 &{} \quad n \leqslant l,\\ 1- \frac{1}{(l+1)!} &{}\quad n = l+1. \end{array}\right. } \end{aligned}$$

This implies, by truncating the power series of $f_l$ at order $l+1$,

$$\begin{aligned} f_l(z)= & {} 1 + z + \frac{z^2}{2!} + \cdots + \frac{z^l}{l!} + \frac{1-\frac{1}{(l+1)!}}{l!}z^{l+1} + O(z^{l+2}) \\= & {} e^z - \frac{z^{l+1}}{((l+1)!)^2} + O(z^{l+2}). \end{aligned}$$

Logarithmic differentiation of the power series thus yields, as $z\rightarrow \infty$,

$$\begin{aligned} z \frac{f'(z)}{f(z)} = z - \frac{z^{l+1}}{l! (l+1)!} + O(z^{l+2}), \quad z \frac{\textrm{d}}{\textrm{d}z} \left( z \frac{f'(z)}{f(z)}\right) = z - \frac{z^{l+1}}{(l!)^2} + O(z^{l+2}), \end{aligned}$$

and the results follow from specializing to $z = r>0$. $\square$

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\rightarrow 0^+$, which reflects, by the proof of Lemma 2.2, the regime $l \geqslant n$, and

• $\frac{1}{2}\log r + \log l$ as $r\rightarrow \infty$, which reflects the regime $l \ll n^{1/4}$ (see Sect. 2.3).

It is this seamless interpolation between the two regimes $l \geqslant n$ and $l \ll 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 \rightarrow \infty$ 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\ll n^{1/4}$.

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

$$\begin{aligned} r_n = \frac{n^2}{l^2} + \frac{n}{2} + O(1) \qquad (n\rightarrow \infty ); \end{aligned}$$

(27)

which suggests to plug its leading order term $r^*_n := n^2/l^2$ into (15). Theorem 2.2 gives, as $n\rightarrow \infty$,

$$\begin{aligned}{} & {} f(r_n^*) = \frac{0!\cdot 1!\cdot 2! \, \cdots \, (l-1)!}{(2\pi )^{l/2} 2^{l^2/2}} \cdot e^{2n} \left( \frac{l}{n}\right) ^{l^2/2} (1+ O(n^{-1})),\\{} & {} a(r_n^*) = n - \tfrac{1}{4}l^2 + O(n^{-1}),\qquad b(r_n^*) = \tfrac{1}{2} n + O(n^{-1}). \end{aligned}$$

The Gaussian term in (15) has thus the expansion

$$\begin{aligned} \exp \left( -\frac{(n-a(r_n^*))^2}{2b(r_n^*)}\right) = e^{-l^4/16n}\big (1+ O(n^{-2})\big ) = 1 - \frac{l^4}{16n} + O(n^{-2}), \end{aligned}$$

(28)

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\ll n^{1/4}$; see Table 2 for an illustrative example. Altogether Hayman’s normal approximation (15) gives, choosing $r=r_n^*$,

$$\begin{aligned} \#\{\sigma : L_n(\sigma ) \leqslant l\} = n!\cdot {\mathbb {P}}(L_n\leqslant l) = \frac{0!\cdot 1!\cdot 2! \, \cdots \, (l-1)!}{(2\pi )^{l/2} 2^{l^2/2}} \cdot \frac{(n!)^2 \left( \frac{e}{n}\right) ^{2n} l^{2n + l^2/2}}{\sqrt{\pi }\, n^{(l^2+1)/2}}(1+ o(1)) \end{aligned}$$

as $n\rightarrow \infty$. Wrapping up by using Stirling’s formula in the form

$$\begin{aligned} (n!)^2 \left( \frac{e}{n}\right) ^{2n} = 2\pi n \,(1+O(n^{-1})) \end{aligned}$$

we have thus given a new proof of Regev’s formula [43, Eq. (4.5.2)]:

$$\begin{aligned} \#\{\sigma : L_n(\sigma ) \leqslant l\} = \frac{0!\cdot 1!\cdot 2! \, \cdots \, (l-1)! \cdot l^{2n + l^2/2}}{(2\pi )^{(l-1)/2}(2n)^{(l^2-1)/2}}(1+o(1))\qquad (n\rightarrow \infty ). \end{aligned}$$

(29)

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²¹ proof proceeds, first, by identifying the leading 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).

3 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

$$\begin{aligned} g_l(s) :=E^{\text {(hard)}}_2(0; [0,s],l), \qquad v_l(s) := - s \frac{\textrm{d}}{\textrm{d}s} \log g_l(s), \qquad u_l(s):= s v_l'(s), \end{aligned}$$

(30a)

for the functions from random matrix theory, we obtain

$$\begin{aligned} f_l(r) = e^r g_l(4r),\qquad a_l(r) = r - v_l(4r),\qquad b_l(r) = r - u_l(4r). \end{aligned}$$

(30b)

As is common in the discussion of the LUE, we generalize this by replacing $l\in {\mathbb {N}}$ with a real parameter $\alpha >-1$. Dropping the index altogether we write, briefly, just g(s), v(s), and u(s).

3.1 Evaluation in Terms of $\sigma$-Painlevé-III

The work of Tracy and Widom [52] shows

$$\begin{aligned} g(s) = \exp \left( -\int _0^s v(x)\,\frac{\textrm{d}x}{x}\right) \qquad (s\geqslant 0), \end{aligned}$$

where v(s) satisfies a Jimbo–Miwa–Okamoto $\sigma$-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

$$\begin{aligned} (x v'')^2 - (\alpha v')^2 + v'(v-x v')(4v'-1) = 0\qquad (x>0), \end{aligned}$$

(31)

subject to the following initial condition, which is consistent with (26) for $\alpha =l\in {\mathbb {N}}$:

$$\begin{aligned} v(x) = \frac{1}{\Gamma (\alpha +1)\Gamma (\alpha +2)} \left( \frac{x}{4}\right) ^{\alpha +1} (1+ O(x)) \qquad (x\rightarrow 0^+). \end{aligned}$$

(32)

A numerical integration of the initial value problem gives approximations to v(s) and $v'(s)$, thus also to $u(s)= s v'(s)$. As explained in [9], a direct numerical integration has stability issues as the solution v is a separatrix solution of the $\sigma$-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\rightarrow \infty$:

$$\begin{aligned} v(x) = \frac{x}{4} - \frac{\alpha }{2}x^{1/2} + \frac{\alpha ^2}{4} + \frac{\alpha }{16}x^{-1/2} + \frac{\alpha ^2}{16}x^{-1} + O(x^{-3/2})\qquad (x\rightarrow \infty ). \end{aligned}$$

(33)

Note the consistency with (24) for $\alpha = l\in {\mathbb {N}}$; for general $\alpha >-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.

3.2 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 $\sigma$-Painlevé-III equation (31) is a comparatively cost-efficient way²² to compile a table of the exact rational values of the distribution ${\mathbb {P}}(L_n \leqslant 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,²³ namely

$$\begin{aligned} v_l''' +\frac{1}{x}v_l'' -\frac{6}{x}v_l'^2 +\frac{4}{x^2}v_lv_l' +\frac{x-l^2}{x^2} v_l'-\frac{1}{2 x^2}v_l = 0, \end{aligned}$$

(34)

which is obtained from differentiating (31) w.r.t. x and dividing the result by $2x^2v_l''^2$. 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 $\sigma$-Painlevé-III equation (31) is quadratic in the highest order derivative and cubic in the lower orders. Therefore, substituting the expansion

$$\begin{aligned} v_l(x) = \sum _{n=l+1}^\infty a_n x^n \end{aligned}$$

(35a)

into the Chazy-I equation (34) yields a much simpler recursive formula for the $a_n$, $n=l+1,\ldots$ :

$$\begin{aligned} (n+1)(n^2-l^2)a_{n+1} + (n-\tfrac{1}{2})a_n - 2\sum _{m=l+1}^{n-l} m a_m \cdot (3(n-m)+1) a_{n+1-m} = 0, \end{aligned}$$

(35b)

uniquely determining the coefficients $a_n$ from the initial value (32), that is, from

$$\begin{aligned} a_{l+1} = \frac{1}{4^{l+1}l!(l+1)!}. \end{aligned}$$

(35c)

It is now a simple matter of truncated power series calculations in a modern computer algebra system to expand the generating function itself,

$$\begin{aligned} f_l(r) = \exp \left( r - \int _0^{4r} v_l(x)\frac{\textrm{d}x}{x}\right) = \exp \left( r - \sum _{n=l+1}^K \frac{4^n a_n}{n}r^n + O(r^{K+1})\right) . \end{aligned}$$

Avoiding the overhead of reducing fractions and computing common denominators in exact rational arithmetic, we have used significance arithmetic with $\lceil 2.5\log _{10} (1000!) \rceil = 6420$ digits and subsequent rational reconstructions to compile a table²⁴ of ${\mathbb {P}}(L_n = l), 1\leqslant l \leqslant 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

$$\begin{aligned} g(s) = \det (I- K)|_{L^2(0,s)}, \end{aligned}$$

(36a)

derived by Forrester [23] in 1993, in terms of the Bessel kernel

$$\begin{aligned} K(x,y) := \frac{\smash [b]{J_\alpha (\sqrt{x}) \sqrt{y} J_{\alpha -1}(\sqrt{y}) - \sqrt{x} J_{\alpha -1}(\sqrt{x}) J_\alpha (\sqrt{y})} }{2(x-y)}. \end{aligned}$$

(36b)

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.

Theorem 3.1

Let K(x, y) be a smooth kernel that induces an integral operator K on $L^2(0,s)$ for all $s>0$ and define the derived kernel as

$$\begin{aligned} K'(x,y) := K(x,y) + xK_x(x,y) + yK_y(x,y). \end{aligned}$$

Then, if we assume $g(s) = \det (I- K)|_{L^2(0,s)}>0$ for all $s>0$, there holds

$$\begin{aligned} v(s)&= - s \frac{\textrm{d}}{\textrm{d}s} \log g(s) = {{\,\mathrm{\text {tr }}\,}}\big ((I-K)^{-1} K' \big )|_{L^2(0,s)},\\ u(s)&= s v'(s) = {{\,\mathrm{\text {tr }}\,}}\big ( (I-K)^{-1} K'' + ((I-K)^{-1} K')^2\big )|_{L^2(0,s)}. \end{aligned}$$

Proof

Rescaling integrals w.r.t. the measure $d\mu (y)=K(x,y)\,dy$ from being taken over the interval (0, s) to (0, 1) induces a transformation of the kernel K(x, y) according to

$$\begin{aligned} K_s(x,y) = s K(sx,sy). \end{aligned}$$

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

$$\begin{aligned} {\dot{K}}_s(x,y) = K'(sx,sy) = s^{-1} K_s'(x,y), \end{aligned}$$

and thus, by [14, Lemma 1], the logarithmic derivative

$$\begin{aligned} g'(s)/g(s) = - {{\,\mathrm{\text {tr }}\,}}\big ((I-K_s)^{-1} \dot{K}_s \big ) |_{L^2(0,1)} = - s^{-1} {{\,\mathrm{\text {tr }}\,}}\big ((I-K)^{-1} K' \big )|_{L^2(0,s)}, \end{aligned}$$

which proves the asserted formula for $v(s) = -s g'(s)/g(s)$. Since, cf. [52, Eq. (2.4)],

$$\begin{aligned} \frac{\textrm{d}}{\textrm{d}s} (I-K_s)^{-1} = (I-K_s)^{-1} \dot{K}_s (I-K_s)^{-1}, \end{aligned}$$

we get by the linearity of the trace

$$\begin{aligned} u(s) = s v'(s) = {{\,\mathrm{\text {tr }}\,}}\big ( (I-K_s)^{-1} K_s'' + ((I-K_s)^{-1} K_s')^2\big )|_{L^2(0,1)}, \end{aligned}$$

which finishes the proof after a back-transformation to $L^2(0,s)$. $\square$

For the Bessel kernel (36b) at hand, we get the derived kernels

$$\begin{aligned} K'(x,y)&= \frac{1}{4}J_\alpha (\sqrt{x})J_\alpha (\sqrt{y}),\\ K''(x,y)&= (1-\alpha ) K'(x,y) + \frac{1}{8}\Big (J_\alpha (\sqrt{x}) \sqrt{y} J_{\alpha -1}(\sqrt{y}) + \sqrt{x} J_{\alpha -1}(\sqrt{x}) J_\alpha (\sqrt{y})\Big ). \end{aligned}$$

We observe that both, $K'$ and $K''$, induce integral operators of finite rank, namely

$$\begin{aligned} K' = \phi \otimes \phi ,\quad K'' = (1-\alpha )\, \phi \otimes \phi + \frac{1}{2} (\phi \otimes \psi + \psi \otimes \phi ), \end{aligned}$$

where we have put

$$\begin{aligned} \phi (x):= \frac{1}{2} J_\alpha (\sqrt{x}), \quad \psi (x):= \frac{1}{2} \sqrt{x} J_{\alpha -1}(\sqrt{x}). \end{aligned}$$

Hence, the results of Theorem 3.1 simplify considerably: first, we obtain²⁵

$$\begin{aligned} v(s) = {{\,\mathrm{\text {tr }}\,}}\big ((I-K)^{-1} \phi \otimes \phi \big )|_{L^2(0,s)} = \langle (I-K)^{-1}\phi , \phi \rangle _{L^2(0,s)}; \end{aligned}$$

(36c)

next, by observing

$$\begin{aligned} K' (I-K)^{-1} K' = \langle (I-K)^{-1}\phi , \phi \rangle _{L^2(0,s)} \cdot \phi \otimes \phi = v(s)\cdot K', \end{aligned}$$

we get, because of symmetry and linearity,

$$\begin{aligned} u(s) = (1-\alpha )v(s) + \langle (I-K)^{-1}\phi , \psi \rangle _{L^2(0,s)} + v(s)^2. \end{aligned}$$

(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 (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 $\infty$, which induce integral operators K on $L^2(s,\infty )$ for all $s \in {\mathbb {R}}$. Here we define the derived kernel as

$$\begin{aligned} K'(x,y):= K_x(x,y) + K_y(x,y) \end{aligned}$$

and get, if $g(s) =\det (I-K)|_{L^2(s,\infty )}>0$ for all $s\in {\mathbb {R}}$, the logarithmic derivative

$$\begin{aligned} \frac{\textrm{d}}{\textrm{d}s} \log g(s) = -{{\,\mathrm{\text {tr }}\,}}\big ((I-K)^{-1} K'\big ) |_{L^2(s,\infty )}. \end{aligned}$$

The proof goes by considering $K_s(x,y)=K(s+x,s+y)$ and transforming $(s,\infty )$ to $(0,\infty )$ 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],

$$\begin{aligned} F_2(s) = \det (I-K)|_{L^2(s,\infty )}, \qquad K(x,y) = \frac{{{\,\mathrm{\text {Ai}}\,}}(x){{\,\mathrm{\text {Ai}}\,}}'(y)-{{\,\mathrm{\text {Ai}}\,}}'(x){{\,\mathrm{\text {Ai}}\,}}(y)}{x-y}. \end{aligned}$$

(37a)

Here we have $K'(x,y) = -{{\,\mathrm{\text {Ai}}\,}}(x){{\,\mathrm{\text {Ai}}\,}}(y)$, i.e., $K'=-{{\,\mathrm{\text {Ai}}\,}}\otimes {{\,\mathrm{\text {Ai}}\,}}$, and thus

$$\begin{aligned} F_2'(s) = -F_2(s)\cdot {{\,\mathrm{\text {tr }}\,}}\big ((I-K)^{-1} K' \big )|_{L^2(s,\infty )} = F_2(s)\cdot \langle (I-K)^{-1}{{\,\mathrm{\text {Ai}}\,}},{{\,\mathrm{\text {Ai}}\,}}\,\rangle _{L^2(s,\infty )}. \end{aligned}$$

(37b)

The last formula was used for the calculations shown in Table 3. In the same manner we get

$$\begin{aligned} F_2''(s)&= 2 F_2(s) \cdot \langle (I-K)^{-1} {{\,\mathrm{\text {Ai}}\,}}, {{\,\mathrm{\text {Ai}}\,}}'\,\rangle _{L^2(s,\infty )}, \end{aligned}$$

(37c)

$$\begin{aligned} F_2'''(s)&= 2 F_2(s) \cdot \left( \langle (I-K)^{-1} {{\,\mathrm{\text {Ai}}\,}}', {{\,\mathrm{\text {Ai}}\,}}'\,\rangle _{L^2(s,\infty )} + \langle (I-K)^{-1} {{\,\mathrm{\text {Ai}}\,}}, {{\,\mathrm{\text {Ai}}\,}}''\,\rangle _{L^2(s,\infty )}\right) , \end{aligned}$$

(37d)

as well as similar formulae for higher order derivatives of $F_2$.²⁶

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

$$\begin{aligned} r_0 := \max (n, (n/l)^2 + n/2). \end{aligned}$$

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:²⁷

$$\begin{aligned} {\mathbb {P}}(L_n \leqslant l) = \tau _n \cdot g_l \cdot \frac{ \exp \left( n \left( \dfrac{v_l}{n}-\log \left( 1+ \dfrac{v_l}{n}\right) \right) \right) }{\sqrt{1 + \dfrac{v_l-u_l}{n}}}(1+ o(1)) \qquad (n\rightarrow \infty ), \end{aligned}$$

(38)

where $g_l$, $v_l$ and $u_l$ are evaluated at $s=4r_{l,n}$ and there is

$$\begin{aligned} \tau _n{} & {} := \frac{n!}{\sqrt{2\pi n}} \left( \frac{e}{n}\right) ^n = 1+\frac{n^{-1}}{12}+\frac{n^{-2}}{288}-\frac{139n^{-3}}{51840}\\ {}{} & {} \quad -\frac{571n^{-4}}{2488320}+\frac{163879n^{-5}}{209018880}+O(n^{-6}). \end{aligned}$$

In IEEE hardware arithmetic we take the definition of $\tau _n$ until $n=100$ and only switch to the shown Stirling expansion for larger n—thus seamlessly providing full accuracy.²⁸ This allows us to approximate the PDF ${\mathbb {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.

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 [27] have recently initiated the study of finite size corrections to the random matrix limit (3), which is

$$\begin{aligned} {\mathbb {P}}( L_n \leqslant l) = F_2(t_l) + o(1),\qquad t_l:=\frac{l-2\sqrt{n}}{n^{1/6}}, \end{aligned}$$

(39)

as $n\rightarrow \infty$, uniformly in $l\in {\mathbb {N}}$. We will refine their study by using the much more accurate and efficient Stirling-type formula (5) instead. Looking at the error

$$\begin{aligned} \delta _0(n) := \max _{l\in \{1,\ldots ,n\}}\big |{\mathbb {P}}( L_n \leqslant l) - F_2(t_l)\big | \end{aligned}$$

for n up to 1000 (see the red crosses in the left panel of Fig. 3 in a double logarithmic scaling) suggest that $\delta _0(n) \approx c_1 n^{-1/3} + c_2 n^{-2/3}$ and yields the conjecture

$$\begin{aligned} {\mathbb {P}}( L_n \leqslant l) = F_2(t_l) + n^{-1/3} F_{2,1}(t_l) + O(n^{-2/3}) \end{aligned}$$

(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 ${{\tilde{F}}}_{2,1}(t)$ of degree 64 to the 836 data points obtained for $n=10^{10}$ in the interval $-8 \leqslant t \leqslant 10$, thus approximating the putative function $F_{2,1}(t)$ there.

Remark 4.1

The error of approximating $F_{2,1}(t)$ by this procedure can be estimated as follows. Extrapolating the errors displayed in the left panel of Fig. 3 shows that the Stirling-type formula induces a perturbation of size $\approx (0.031 n^{-2/3} + 0.058 n^{-1}) n^{1/3}|_{n=10^{10}} = 1.4\times 10^{-5}$. On the other hand, the finite size effect of the next order term $n^{-2/3} F_{2,2}(t)$, displayed in the left panel of Fig. 6, induces a perturbation of size $\approx 0.25 n^{-1/3}|_{n=10^{10}} = 1.2\cdot 10^{-4}$. Thus, altogether ${{\tilde{F}}}_{2,1}$ approximates $F_{2,1}$ up to an error²⁹ of the order $10^{-4}$.

Fig. 6

Rescaled differences between the distributions of $L_n$ 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 ${{\tilde{F}}}_{2,1}(t)$ from Fig. 4 for the exact values of the distribution for $n=250$ (red $+$), $n=500$ (green $\circ$), $n=1000$ (blue $\bullet$). Left: CDF errors rescaled by $n^{2/3}$, horizontal axis is $t=(l-2\sqrt{n})/n^{1/6}$. The solid line is a polynomial ${{\tilde{F}}}_{2,2}(t)$ of degree 48 fitted to all the 152 data points with $-7.5 \leqslant t \leqslant 9.5$; it approximates $F_{2,2}(t)$ in that interval. Right: PDF errors rescaled by $n^{5/6}$, horizontal axis is $t=(l-\frac{1}{2}-2\sqrt{n})/n^{1/6}$. The solid line displays ${{\tilde{F}}}_{2,2}'(t)+{{\tilde{F}}}_{2,1}'''(t)/24 + F_2^{(5)}(t)/1920$ as an approximation of $F_{2,2}'(t)+ F_{2,1}'''(t)/24 + F_2^{(5)}(t)/1920$, with the polynomial ${{\tilde{F}}}_{2,2}(t)$ taken from the left panel and ${{\tilde{F}}}_{2,1}(t)$ as in Fig. 4. The dotted line displays the term $\tilde{F}_{2,2}'(t)$ only (Color figure online)

If we iterate this approach yet another step, by looking at the error

$$\begin{aligned} \delta _1(n) := \max _{l\in \{1,\ldots ,n\}}\big |{\mathbb {P}}( L_n \leqslant l) - F_2(t_l) - n^{-1/3} {{\tilde{F}}}_{2,1}(t_l)\big | \end{aligned}$$

for n up to 1000, then a double logarithmic plot (the green circles in the left panel of Fig. 3) suggests that $\delta _1(n) \approx c_1 n^{-2/3} + c_2 n^{-1}$. As stated in the introduction, this yields the refinement (9) of conjecture (40), namely that there is further a function $F_{2,2}(t)$ such that

$$\begin{aligned} {\mathbb {P}}(L_n\leqslant l) = F_2(t_l) + n^{-1/3} F_{2,1}(t_l) + n^{-2/3} F_{2,2}(t_l) + O(n^{-1}) \qquad (n\rightarrow \infty ), \end{aligned}$$

(41)

uniformly in $l\in {\mathbb {N}}$.

Remark 4.2

To validate conjecture (41) against numerical data, obtained by replacing $F_{2,1}(t)$ by the approximation ${{\tilde{F}}}_{2,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 ${{\tilde{F}}}_{2,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 $\approx (0.031 n^{-2/3} + 0.058 n^{-1}) n^{2/3}$. The sweet spot of both perturbations combined is at $n\approx 1.4\times 10^4$ with a minimum error of about $3.6\times 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 ${{\tilde{F}}}_{2,2}(t)$ of degree 48 to all of the 152 data points in the interval $-7.5 \leqslant t \leqslant 9.5$, approximating the putative function $F_{2,2}(t)$ there.

4.2 The PDF of the Distribution of $L_n$

If we apply the central differencing formula, for smooth functions F(x) and increments $h\rightarrow 0$, that is to say

$$\begin{aligned} F(x+h)-F(x)= & {} h F'(x+h/2) +\frac{h^3}{24} F'''(x+h/2)\\ {}{} & {} + \frac{h^5}{1920} F^{(5)} (x+ h/2) + O(h^7), \end{aligned}$$

with the increment $h=n^{-1/6}$ to the conjectured expansion (41), now written in the form

$$\begin{aligned} {\mathbb {P}}(L_n=l)= & {} {\mathbb {P}}(L_n\leqslant l) - {\mathbb {P}}(L_n\leqslant l-1)\\= & {} \big (F_2(t_l)-F_2(t_{l-1}) \big ) + n^{-1/3} \big (F_{2,1}(t_l)-F_{2,1}(t_{l-1}) \big ) \\{} & {} + n^{-2/3} \big (F_{2,2}(t_l) - F_{2,2}(t_{l-1}) \big ) + \cdots , \end{aligned}$$

we get, assuming some uniformity, an induced expansion of the PDF:

$$\begin{aligned} {\mathbb {P}}(L_n= l)= & {} n^{-1/6} F_2'({{\hat{t}}}_l) + n^{-1/2} \Big (F_{2,1}'({{\hat{t}}}_l)+\frac{1}{24}F_2'''({{\hat{t}}}_l)\Big )\nonumber \\{} & {} + n^{-5/6} \Big (F_{2,2}'({{\hat{t}}}_l) + \frac{1}{24}F_{2,1}'''({{\hat{t}}}_l)+\frac{1}{1920}F_2^{(5)}({{\hat{t}}}_l)\Big ) + O(n^{-7/6})\qquad (n\rightarrow \infty )\nonumber \\ \end{aligned}$$

(42a)

uniformly in $l\in {\mathbb {N}}$, where we have briefly written

$$\begin{aligned} {{\hat{t}}}_l := \frac{l-\frac{1}{2}-2\sqrt{n}}{n^{1/6}}. \end{aligned}$$

(42b)

Note the shift by 1/2 in the numerator of ${{\hat{t}}}_l$ 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.

4.3 The Expected Value of $L_n$

By a shift and rescale, the expected value of the discrete random variable $L_n$ can be written in the form

$$\begin{aligned} {\mathbb {E}}(L_n) = \sum _{l=1}^n l \cdot {\mathbb {P}}(L_n=l) = 2\sqrt{n} + \frac{1}{2} + n^{1/6} \sum _{l=1}^n {{\hat{t}}}_l \cdot {\mathbb {P}}(L_n=l). \end{aligned}$$

If we write (42a), with obvious definitions of the functions $G_j(t)$, in the form

$$\begin{aligned} {\mathbb {P}}(L_n= l) = n^{-1/6} G_0({{\hat{t}}}_l) + n^{-1/2} G_1({{\hat{t}}}_l) + n^{-5/6} G_2({{\hat{t}}}_l) + O(n^{-7/6}) \end{aligned}$$

(43)

we get the induced expansion

$$\begin{aligned} {\mathbb {E}}(L_n)&= 2\sqrt{n} + n^{1/6} \mu _0^{(n)} + \frac{1}{2} + n^{-1/6} \mu _1^{(n)} + n^{-1/2} \mu _2^{(n)} + \cdots \\ \mu _j^{(n)}&:= n^{-1/6} \sum _{l=1}^n {{\hat{t}}}_l \, G_j({{\hat{t}}}_l). \end{aligned}$$

Now, if we assume (a) that the decay $G_j(t) \rightarrow 0$ (and likewise for all the derivatives) is exponentially fast as $t\rightarrow \pm \infty$ (see Figs. 4 and 6) and (b) that the $G_j$ can be extended analytically to a strip containing the real axis, we obtain

$$\begin{aligned} \mu _j^{(n)} \doteq n^{-1/6} \sum _{l=-\infty }^\infty {{\hat{t}}}_l \, G_j({{\hat{t}}}_l) \doteq \int _{-\infty }^\infty t\, G_j(t)\,\textrm{d}t =:\mu _j, \end{aligned}$$

(44)

where “$\doteq$” denotes equality up to terms that are exponentially small for large n. Here, in the first step the series was obtained by adding, under assumption (a), the exponentially small tail, and in the next step we have identified the series as the trapezoidal rule with step-size $h = n^{-1/6}$—a quadrature rule known to converge, under assumption (b), exponentially fast to the integral, cf. [53].

Table 3 Exponentially fast convergence of first ($k=1$) and second ($k=2$) moments, $\mu _{0,k}^{(n)} := n^{-1/6} \sum _{l=1}^n {{\hat{t}}}_l^k F_2'({{\hat{t}}}_l) \;\rightarrow \; \mu _{0,k}^{(\infty )} := \int _{-\infty }^\infty t^k F_2'(t)\,\textrm{d}t$

Remark 4.3

For the function $G_0(t) = F_2'(t)$, i.e., the density of the Tracy–Widom distribution, the asserted exponentially fast convergence (44) can be checked against numerical data which were obtained by applying the highly accurate numerical methods described in [9] to the representation (37b) of $F_2'(t)$; see Table 3.

We have thus derived from (41)—based on the assumptions of uniformity, exponential decay, and analytic continuation of the functions $G_j$ and their derivatives—the following expansion which adds three more terms to the expansion given in [4, Thm. 1.2]:

$$\begin{aligned} {\mathbb {E}}(L_n) = 2\sqrt{n} + \mu _0 n^{1/6} + \frac{1}{2} + \mu _1n^{-1/6} + \mu _2n^{-1/2} + O(n^{-5/6}); \end{aligned}$$

(45)

where we have, with the numerical value of $\mu _0$ taken from [9, Table 10], cf. also Table 3,

$$\begin{aligned} \mu _0&= \int _{-\infty }^\infty t \, F_2'(t)\,\textrm{d}t = -1.77108\,68074\,\cdots ,\\ \mu _1&= \int _{-\infty }^\infty t\Big (F_{2,1}'(t)+\frac{1}{24}F_2'''(t)\Big )\,\textrm{d}t = \int _{-\infty }^\infty t \, F_{2,1}'(t) \,\textrm{d}t,\\ \mu _2&= \int _{-\infty }^\infty t\Big (F_{2,2}'(t) + \frac{1}{24}F_{2,1}'''(t)+\frac{1}{1920}F_2^{(5)}(t)\Big )\,\textrm{d}t = \int _{-\infty }^\infty t \, F_{2,2}'(t) \,\textrm{d}t. \end{aligned}$$

We note that the higher derivatives do not contribute to the integral values, as can be shown using integration by parts and the assumed exponential decay to zero. Based on the polynomial approximations displayed in Figs. 4 and 6 we get the numerical estimates—comparing, in addition to $n=10^{10}$, with the analogous results for $n=10^9$ and $n=10^{11}$:

$$\begin{aligned} \mu _1 \approx \int _{-8}^{10} t {{\tilde{F}}}'_{2,1}(t)\,\textrm{d}t \approx 0.0659, \qquad \mu _2 \approx \int _{-7.5}^{9.5} t \tilde{F}'_{2,2}(t)\,\textrm{d}t \approx 0.25. \end{aligned}$$

(46)

Further evidence for the validity of the expansion (45) comes from looking at a least squares fit of the form³⁰ (where the upper bound $k=9$ has been chosen as to maximize the number of matching digits for the two data sets below)

$$\begin{aligned} {\mathbb {E}}(L_n) \approx 2\sqrt{n} + \frac{1}{2} + \sum _{k=0}^9 c_k n^{(1-2k)/6}, \end{aligned}$$

with ${\mathbb {E}}(L_n)$ obtained from the tabulated values of ${\mathbb {P}}(L_n=l)$ up to $n=1000$. If we do so in extended precision for two different data sets, first for n from 500 upwards and next for n from 600 upwards, we obtain as digits that are matching in both cases

$$\begin{aligned}{} & {} c_0 = -1.77108\,68074\cdots , \; c_1 = 0.06583\,238\cdots , \; c_2 = 0.26122\,27\cdots , \; \\{} & {} \quad c_3 = -0.11938\,4\cdots . \end{aligned}$$

Here the value of $c_0$ is in perfect agreement with the known value of $\mu _0$ and $c_1$, $c_2$ are consistent with the inaccuracies of the estimates in (46), cf. Remarks 4.1/4.2. Hence, the most accurate values that we can offer for $\mu _1$ and $\mu _2$ are those displayed in (12).

4.4 The Variance of $L_n$

By a shift and rescale, the variance of $L_n$ can we written as

$$\begin{aligned} \textrm{Var}(L_n)= & {} \sum _{l=1}^n l^2 \cdot {\mathbb {P}}(L_n=l) - {\mathbb {E}}(L_n)^2 = n^{1/3} \sum _{l=1}^n {{\hat{t}}}_l^2 \cdot {\mathbb {P}}(L_n=l) \\{} & {} - \left( {\mathbb {E}}(L_n)- 2\sqrt{n} -\frac{1}{2}\right) ^2. \end{aligned}$$

By inserting the expansions (43) and (45), and by arguing as in (44), we get the following expansion which adds two more terms to the leading order found in [4, Thm. 1.2]:

$$\begin{aligned} \textrm{Var}(L_n) = \nu _0 n^{1/3} + \nu _1 + \nu _2 n^{-1/3} + O(n^{-2/3}); \end{aligned}$$

(47)

where we have, with the numerical value of $\nu _0$ taken from [9, Table 10], cf. also Table 3,³¹

$$\begin{aligned} \nu _0= & {} \int _{-\infty }^\infty t^2 F_2'(t)\,\textrm{d}t - \mu _0^2 = 0.81319\,47928\,\cdots ,\\ \nu _1= & {} \int _{-\infty }^\infty t^2 F_{2,1}'(t)\,\textrm{d}t + \frac{1}{12} -2\mu _0\mu _1 ,\qquad \nu _2 = \int _{-\infty }^\infty t^2 F_{2,2}'(t)\,\textrm{d}t -\mu _1^2 - 2\mu _0\mu _2. \end{aligned}$$

By using the polynomial approximations ${{\tilde{F}}}_{2,1}$ and $\tilde{F}_{2,2}$ displayed in Figs. 4 and 6, as well as the numerical values of $\mu _0$, $\mu _1$, $\mu _2$ from (12), we get the estimates $\nu _1 \approx -1.2070$ and $\nu _2 \approx 0.57$.

Once again, further evidence and increased numerical accuracy comes from a least squares fit of the form³² (where the upper bound $k=8$ has been chosen as to maximize the number of matching digits for the two data sets below)

$$\begin{aligned} \textrm{Var}(L_n) \approx \sum _{k=0}^8 c_k n^{(1-k)/3}, \end{aligned}$$

with $\textrm{Var}(L_n)$ obtained from the tabulated values of ${\mathbb {P}}(L_n=l)$ up to $n=1000$. If we do so in extended precision for two different data sets, first for n from 500 upwards and then for n from 600 upwards, we obtain as digits that are matching in both cases

$$\begin{aligned}{} & {} c_0 = 0.81319\,47928\cdots ,\; c_1 = -1.20720\,507\cdots ,\\{} & {} c_2 = 0.56715\,6\cdots ,\; c_3 = 0.01669\cdots . \end{aligned}$$

Here the value of $c_0$ is in perfect agreement with the known value of $\nu _0$ and the values for $c_1$, $c_2$ are consistent with the inaccuracies of the estimates for $\nu _1$ and $\nu _2$ shown above, cf. Remarks 4.1/4.2. Hence, the most accurate values that we can offer for the coefficients $\nu _1$ and $\nu _2$ are those displayed in (13).

Remark 4.4

(added in proof) The conjectured functional form (11) of $F_{2,1}(t)$ gives numerical values of the coefficients $\mu _1$ and $\nu _1$ that agree with the values displayed in (12) and (13) to all digits shown.

A. Appendix: The Multidimensional Laplace Method

The classical one-dimensional method of Laplace can be generalized to provide the asymptotics, as $z\rightarrow \infty$, of multidimensional integrals of the form

$$\begin{aligned} \int _\Omega e^{-z S(x)} f(x)\,\textrm{d}x. \end{aligned}$$

Here, $\Omega \subset {\mathbb {R}}^n$ is a measurable set and $f, S : \Omega \rightarrow {\mathbb {R}}$ are subject to suitable assumptions. For instance, if we assume that f and S are sufficiently smooth and the phase function S(x) takes a unique minimum at an interior point $x_*$ of $\Omega$ with $\det S''(x_*)\ne 0$ then the standard result—going back to Hsu [34, Lemma 1]—states that for each $0<\delta \leqslant \frac{\pi }{2}$ as $z\rightarrow \infty$

$$\begin{aligned} \int _\Omega e^{-z S(x)} f(x)\,\textrm{d}x = \left( \frac{2\pi }{z}\right) ^{n/2} \frac{e^{-z S(x_*)}}{\sqrt{ \det S''(x_*) }} \left( f(x_*) + O\big (z^{-1}\big )\right) \quad (|\!\arg z| \leqslant \tfrac{1}{2}\pi -\delta ). \end{aligned}$$

(48)

Remark A.1

As is customary in asymptotic analysis in the complex plane, cf. [41, p. 7], we understand this asymptotic expansion (and similar expansions with o- or O-terms) to hold uniformly in $\{z\in {\mathbb {C}}: |z|\geqslant R_\delta , |\arg z| \leqslant \frac{\pi }{2}-\delta \}$ if $R_\delta >0$ has been chosen sufficiently large.

If $f(x_*)=0$, however, formula (48) fails to yield the precise leading order term of the expansion. On the other hand, it is known that there holds, for sufficiently smooth f and S, cf. [21, Eq. (1.27)], a general asymptotic expansion of the form

$$\begin{aligned} \int _\Omega e^{-z S(x)} f(x)\,\textrm{d}x \sim e^{-z S(x_*)}z^{-n/2}\sum _{k=0}^\infty c_k z^{-k} \qquad (z\rightarrow \infty ,\,|\!\arg z| \leqslant \tfrac{1}{2}\pi -\delta ). \end{aligned}$$

(49)

Still, it would be extremely awkward to determine the first nonzero coefficient from the standard proof³³ of (49), in which the $c_k$ depend on higher order derivatives of a nonlinear transform obtained from the Morse lemma, deforming S(x) to a quadratic form.

Building on a different technique introduced by Fulks and Sather [29], Kirwin [36] succeeded in establishing formulae for the higher order coefficients in terms of asymptotic expansions of f and S into homogeneous functions at $x_*$. However, these authors consider only the case of real $z\rightarrow \infty$, whereas we need uniformity, for arbitrary small $\delta >0$, as $z\rightarrow \infty$ with $|\arg z| \leqslant \frac{\pi }{2}-\delta$. Since the leading order term of their expansions suffices for the purposes of Sect. 2, we will give a much simplified version of their proof in this Appendix, explicitly tracing constants to establish the required uniformity.

1.1 Notation and Assumptions

By a simple transformation (see the proof of Theorem 2.2) we can restrict ourselves to

$$\begin{aligned} x_* = 0, \qquad S(x_*)=0. \end{aligned}$$

Writing $z = \sigma + i\tau$ we have

$$\begin{aligned} 0< \sigma \leqslant |z| \leqslant \sigma \csc \delta \qquad (|\!\arg z| \leqslant \tfrac{1}{2}\pi -\delta ). \end{aligned}$$

(50)

Let $\Omega \subset {\mathbb {R}}^n$ we a measurable set with 0 as an interior point, that is, there is $\epsilon _0 >0$ such that $B_\epsilon (0)\subset \Omega$ for all open balls of radius $0< \epsilon \leqslant \epsilon _0$ centered at 0. By denoting $|\cdot |$ the Euclidean norm on ${\mathbb {R}}^n$, we write ${\mathbb {R}}^n\setminus \{0\} \ni x = \rho \cdot \xi$ in spherical coordinates with $\rho =|x|$ and $\xi =x/\rho \in S^{n-1}$. We assume that $f, S: \Omega \rightarrow {\mathbb {R}}$ are measurable functions subject to the following conditions:

• (1) S(x) is positively bounded away from zero on $\Omega \setminus B_\epsilon (0)$ for each $\epsilon >0$;

• (2) there is a $\nu >0$ and a positive continuous function $S_0:S^{n-1}\rightarrow {\mathbb {R}}$ with

  $$\begin{aligned} S(x) = \rho ^\nu S_0(\xi ) + o(\rho ^\nu )\qquad (\rho \rightarrow 0), \end{aligned}$$

  uniformly in $\xi \in S^{n-1}$;

• (3) there is a $\lambda > 0$ and a bounded measurable function $f_0:S^{n-1}\rightarrow {\mathbb {R}}$ with

  $$\begin{aligned} f(x) = \rho ^{\lambda -n} f_0(\xi ) + o(\rho ^{\lambda -n})\qquad (\rho \rightarrow 0), \end{aligned}$$

  uniformly in $\xi \in S^{n-1}$;

• (4) integrability: there is $\sigma _0>0$ with $M = \int _\Omega e^{-\sigma _0 S(x)}|f(x)| \,dx <\infty$.

We extend $f_0$ and $S_0$ to all of ${\mathbb {R}}^n$ by homogeneity, that is, by

$$\begin{aligned} f_0(x) = \rho ^{\lambda -n} f_0(\xi ), \quad S_0(x) = \rho ^{\nu } S_0(\xi ) \qquad (x\ne 0) \end{aligned}$$

and, for definiteness, $f_0(0)=S_0(0)=0$. Finally, for purposes of reference we recall the following well known integral evaluation,

$$\begin{aligned} \int _0^\infty e^{-z \rho ^\nu } \rho ^{\lambda - 1}\,d\rho = \frac{\Gamma (\lambda /\nu )}{\nu } z^{-\lambda /\nu }\qquad (\sigma = \Re z > 0). \end{aligned}$$

(51)

The leading order result of [29, p. 186]³⁴ and [36, Thm. 1.1] is now as follows:

Theorem A.1

Under conditions (1)–(4) there holds for each $0< \delta \leqslant \frac{\pi }{2}$ as $z\rightarrow \infty$

$$\begin{aligned} \int _\Omega e^{-z S(x)} f(x)\,\textrm{d}x = z^{-\lambda /\nu } \int _{{\mathbb {R}}^n} e^{- S_0(x)} f_0(x)\,\textrm{d}x + o(z^{-\lambda /\nu })\qquad (|\!\arg z| \leqslant \tfrac{1}{2}\pi -\delta ). \end{aligned}$$

(52)

Denoting the surface measure on $S^{n-1}$ by $\omega$, the integral on the right evaluates to

$$\begin{aligned} \int _{{\mathbb {R}}^n} e^{- S_0(x)} f_0(x)\,\textrm{d}x = \frac{\Gamma (\lambda /\nu )}{\nu } \int _{S^{n-1}} \frac{f_0(\xi )}{S_0(\xi )^{\lambda /\nu }}\,\textrm{d}\omega (\xi ). \end{aligned}$$

We fix $\delta >0$ and break the proof of this theorem, based on the idea of “trading tails” (a notion popularized for Laplace’s method in [32, p. 466]), into some preparatory steps.

Lemma A.1

For $\sigma \geqslant \sigma _0$ and $0<\epsilon \leqslant \epsilon _0$, there is a constant $\alpha >0$ such that

$$\begin{aligned} \int _{\Omega \setminus B_\epsilon (0)} e^{-zS(x)} f(x)\,dx = O\left( e^{-(\sigma -\sigma _0)\alpha \epsilon ^\nu }\right) \end{aligned}$$

where the implied constant does not depend on z and $\epsilon$.

Proof

By Conditions (1) and (2) there is a constant $\alpha >0$ such that

$$\begin{aligned} S(x) \geqslant \alpha \epsilon ^\nu \qquad (x\in \Omega \setminus B_\epsilon (0)). \end{aligned}$$

Hence, Condition (4) yields straightforwardly

$$\begin{aligned} \left| \int _{\Omega \setminus B_\epsilon (0)} e^{-zS(x)} f(x)\,\textrm{d}x\right| \leqslant e^{-(\sigma - \sigma _0)\alpha \epsilon ^\nu } \int _\Omega e^{-\sigma _0 S(x)} |f(x)|\,\textrm{d}x = O\left( e^{-(\sigma -\sigma _0) \alpha \epsilon ^\nu }\right) , \end{aligned}$$

the implied constant being just M. $\square$

Lemma A.2

There is $\eta _\epsilon \rightarrow 0$ for $\epsilon \rightarrow 0$ such that for $\sigma \geqslant \sigma _0$ and $0<\epsilon \leqslant \epsilon _0$

$$\begin{aligned} \int _{B_\epsilon (0)} e^{-zS(x)} f(x)\,\textrm{d}x = \int _{B_\epsilon (0)} e^{-zS(x)} f_0(x)\,\textrm{d}x + \eta _\epsilon \cdot O\left( \sigma ^{-\lambda /\nu }\right) \end{aligned}$$

where the implied constant does not depend on z and $\epsilon$.

Proof

Conditions (1) and (2) give the existence of a constant $\gamma >0$ such that

$$\begin{aligned} S(x) \geqslant \gamma |x|^\nu \qquad (|x|\leqslant \epsilon _0) \end{aligned}$$

and Condition (3) gives $\eta _\epsilon \rightarrow 0$ for $\epsilon \rightarrow 0$ with

$$\begin{aligned} |f(x)-f_0(x)| \leqslant \eta _\epsilon |x|^{\lambda -n}\qquad (|x|\leqslant \epsilon ). \end{aligned}$$

Hence we get

$$\begin{aligned} \left| \int _{B_\epsilon (0)} e^{-zS(x)} (f(x)-f_0(x))\,\textrm{d}x \right|\leqslant & {} \eta _\epsilon \int _{{\mathbb {R}}^n} e^{-\sigma \gamma |x|^\nu } |x|^{\lambda -n}\,\textrm{d}x \\= & {} \eta _\epsilon \omega _{n-1} \int _0^\infty e^{-\sigma \gamma \rho ^\nu } \rho ^{\lambda -1}\,\textrm{d}\rho \\= & {} \eta _\epsilon \omega _{n-1} \nu ^{-1} \Gamma (\lambda /\nu ) (\sigma \gamma )^{-\lambda /\nu } \end{aligned}$$

where the integral over ${\mathbb {R}}^n$ was evaluated by spherical symmetry and the resulting gamma integral by Eq. (51); $\omega _{n-1}$ denotes the surface area of the sphere $S^{n-1}$. $\square$

Lemma A.3

There is $\eta '_\epsilon \rightarrow 0$ for $\epsilon \rightarrow 0$ such that, for $0<\epsilon \leqslant \epsilon _0$ with $\epsilon _0$ sufficiently small and $z\in {\mathbb {C}}$ with $|\!\arg z|\leqslant \frac{\pi }{2} - \delta$, $|z|\geqslant \sigma _0\cdot \csc \delta$,

$$\begin{aligned} \int _{B_\epsilon (0)} e^{-zS(x)} f_0(x)\,\textrm{d}x = \int _{B_\epsilon (0)} e^{-zS_0(x)} f_0(x)\,\textrm{d}x + \eta '_\epsilon \cdot O\left( z^{-\lambda /\nu }\right) \end{aligned}$$

where the implied constant does not depend on z and $\epsilon$.

Proof

Condition (2) gives $\eta ''_\epsilon \rightarrow 0$ for $\epsilon \rightarrow 0$ and $\alpha '>0$ with

$$\begin{aligned} |S(x)-S_0(x)| \leqslant \eta ''_\epsilon |x|^\nu \quad (|x|\leqslant \epsilon ) \quad \text {and}\quad S_0(x) \geqslant \alpha ' |x|^\nu \quad (x\in {\mathbb {R}}^n), \end{aligned}$$

(53)

whereas Condition (3) yields a constant $\gamma '>0$ such that

$$\begin{aligned} |f_0(x)| \leqslant \gamma ' |x|^{\lambda -n} \qquad (x\in {\mathbb {R}}^n). \end{aligned}$$

Because of $|e^w-1| \leqslant e^{|w|}-1$ for $w\in {\mathbb {C}}$ and by (50) we obtain for $x\in B_\epsilon (0)$

$$\begin{aligned} \left| e^{-zS(x)} -e^{-zS_0(x)}\right|\leqslant & {} e^{-\sigma S_0(x)} \left( e^{\eta ''_\epsilon |z| \cdot |x|^\nu } - 1 \right) \leqslant e^{-\sin (\delta ) \alpha ' |z| \cdot |x|^\nu } \left( e^{\eta ''_\epsilon |z| \cdot |x|^\nu } - 1 \right) \\\leqslant & {} e^{-\sin (\delta ) (\alpha '-\eta ''_\epsilon \csc \delta ) |z|\cdot |x|^\nu } - e^{-\sin (\delta ) \alpha ' |z| \cdot |x|^\nu } \end{aligned}$$

Thus, if $\epsilon _0$ is small enough to guarantee $\eta ''_\epsilon \csc \delta < \alpha '$, we obtain by the same calculations as previously in the proof of Lemma A.2

$$\begin{aligned}{} & {} \left| \int _{B_\epsilon (0)} \left( e^{-zS(x)}-e^{-zS_0(x)}\right) f_0(x)\,\textrm{d}x \right| \\{} & {} \quad \leqslant \gamma ' \omega _{n-1}\frac{\Gamma (\lambda /\nu )}{\nu } (\sin (\delta )|z|)^{-\lambda /\nu }\left( \frac{1}{(\alpha '-\eta ''_\epsilon \csc \delta )^{\lambda /\nu }}-\frac{1}{(\alpha ')^{\lambda /\nu }}\right) = \eta _\epsilon ' \cdot O(z^{-\lambda /\nu }) \end{aligned}$$

if we define the term given by the large bracket to be $\eta _\epsilon '$. $\square$

Proof of Theorem A.1

By splitting the integral as follows and by applying Lemmas A.1–A.3, we get for $0<\epsilon \leqslant \epsilon _0$ with $\epsilon _0$ sufficiently small and for $z\in {\mathbb {C}}$ with $|\!\arg z|\leqslant \frac{\pi }{2} - \delta$, $|z|\geqslant \sigma _0\cdot \csc \delta$

$$\begin{aligned} \int _\Omega e^{-zS(x)} f(x)\,\textrm{d}x= & {} \int _{B_\epsilon (0)} e^{-z S_0(x)} f_0(x)\,\textrm{d}x + \int _{B_\epsilon (0)} (e^{-z S(x)}-e^{-z S_0(x)}) f_0(x)\,\textrm{d}x\\{} & {} + \int _{B_\epsilon (0)} e^{-zS(x)} (f(x)-f_0(x))\,dx + \int _{\Omega \setminus B_\epsilon (0)} e^{-zS(x)} f(x)\,\textrm{d}x\\= & {} \int _{B_\epsilon (0)} e^{-zS_0(x)} f_0(x)\,\textrm{d}x + \eta _\epsilon ' \cdot O(z^{-\lambda /\nu }) \\{} & {} + \eta _\epsilon \cdot O(\sigma ^{-\lambda /\nu }) + O\left( e^{-(\sigma -\sigma _0)\alpha \epsilon ^\nu }\right) \end{aligned}$$

where the implied constants do not depend on z and $\epsilon$. By (53) and Lemma A.1, we get likewise

$$\begin{aligned} \int _{{\mathbb {R}}^n} e^{-z S_0(x)} f_0(x)\,\textrm{d}x = \int _{B_\epsilon (0)} e^{-zS_0(x)} f_0(x)\,\textrm{d}x + O\left( e^{-(\sigma -\sigma _0)\alpha '\epsilon ^\nu }\right) . \end{aligned}$$

Coupling $\epsilon \propto |z|^{-1/2\nu }$, we thus get some expression $\eta (|z|^{-1})\rightarrow 0$ for $z\rightarrow \infty$ with

$$\begin{aligned} \int _\Omega e^{-zS(x)} f(x)\,\textrm{d}x = \int _{{\mathbb {R}}^n} e^{-z S_0(x)} f_0(x)\,\textrm{d}x + \eta (|z|^{-1})\cdot O(z^{-\lambda /\nu }) \end{aligned}$$

where the implied constant does not depend on z. Finally, by transforming to spherical coordinates and using (51) for the inner integral once more, we calculate

$$\begin{aligned} \int _{{\mathbb {R}}^n} e^{-z S_0(x)} f_0(x)\,\textrm{d}x= & {} \int _{S^{n-1}}\int _0^\infty e^{- z S_0(\xi ) \rho ^\nu } f_0(\xi )\rho ^{\lambda -1}\,\textrm{d}\rho \,d\omega (\xi ) \\= & {} z^{-\lambda /\nu }\frac{\Gamma (\lambda /\nu )}{\nu } \int _{S^{n-1}} \frac{f_0(\xi )}{S_0(\xi )^{\lambda /\nu }}\,\textrm{d}\omega (\xi ) \end{aligned}$$

which finishes the proof. $\square$

1.2 Quadratic Leading Order Term in the Phase Function

It is straightforward from (52) to specialize Theorem A.1 to the case of a quadratic leading order term $S_0$ in the asymptotic expansion of the phase function S.

Corollary A.1

Under conditions (1)–(4) with a quadratic $S_0(x) = x^T H x/2$ defined by a symmetric positive definite matrix $H\in {\mathbb {R}}^{n\times n}$, there holds for each $0<\delta \leqslant \frac{\pi }{2}$ as $z\rightarrow \infty$

$$\begin{aligned} \int _\Omega e^{-zS(x)} f(x)\,\textrm{d}x = \frac{(2\pi )^{n/2} z^{-\lambda /2}}{\sqrt{\det H}}({\mathbb {E}}(f_0)+o(1)) \qquad (|\!\arg z| \leqslant \tfrac{1}{2}\pi -\delta ) \end{aligned}$$

where ${\mathbb {E}}$ denotes expectation with respect to the multivariate normal distribution with covariance matrix $H^{-1}$, namely

$$\begin{aligned} {\mathbb {E}}(f_0) := \frac{\sqrt{\det H}}{(2\pi )^{n/2}}\int _{{\mathbb {R}}^n} e^{-x^THx/2} f_0(x)\,\textrm{d}x. \end{aligned}$$

Remark A.2

Corollary A.1 is providing the precise leading order asymptotic of the integral only if the condition ${\mathbb {E}}(f_0) \ne 0$ is satisfied. In the case of a sufficiently smooth integrand the function $f_0$ is the first nonzero homogeneous polynomial appearing in the Taylor expansion of f at zero. For symmetry reasons, ${\mathbb {E}}(f_0) \ne 0$ implies that $\deg f_0$ must be even. Thus, if also S is sufficiently smooth, a comparison with (49) yields, if ${\mathbb {E}}(f_0)\ne 0$,

$$\begin{aligned} \int _\Omega e^{-zS(x)} f(x)\,\textrm{d}x = \frac{(2\pi )^{n/2} z^{-\frac{n+\deg f_0}{2}}}{\sqrt{\det H}}{\mathbb {E}}(f_0)\left( 1+O(z^{-1})\right) \qquad (|\!\arg z| \leqslant \tfrac{1}{2}\pi -\delta ). \end{aligned}$$

(54)

If $f(0)\ne 0$, this reproduces Hsu’s formula (48) since then $f_0\equiv f(0)$ and hence ${\mathbb {E}}(f_0)=f(0)$.