Asymptotic expansions relating to the distribution of the length of longest increasing subsequences

Abstract

We study the distribution of the length of longest increasing subsequences in random permutations of n integers as n grows large and establish an asymptotic expansion in powers of $n^{-1/3}$. Whilst the limit law was already shown by Baik, Deift and Johansson to be the GUE Tracy–Widom distribution F, we find explicit analytic expressions of the first few finite-size correction terms as linear combinations of higher order derivatives of F with rational polynomial coefficients. Our proof replaces Johansson’s de-Poissonization, which is based on monotonicity as a Tauberian condition, by analytic de-Poissonization of Jacquet and Szpankowski, which is based on growth conditions in the complex plane; it is subject to a tameness hypothesis concerning complex zeros of the analytically continued Poissonized length distribution. In a preparatory step an expansion of the hard-to-soft edge transition law of LUE is studied, which is lifted to an expansion of the Poissonized length distribution for large intensities. Finally, expansions of Stirling-type approximations and of the expected value and variance of the length distribution are given.

1. Introduction

The length $L_n(\sigma )$ of longest increasing subsequences ¹ of permutations $\sigma $ on $\{1,2,\ldots ,n\}$ becomes a discrete random variable when the permutations are drawn randomly with uniform distribution. This way, the problem of enumerating all permutations $\sigma $ that satisfy $L_n(\sigma )\leqslant l$ gets encoded in the discrete probability distribution ${\mathbb P}(L_n \leqslant l)$ . The present paper studies an asymptotic expansion of this distribution when n grows large. As there are relations to KPZ growth models (directly so for the PNG model with droplet initial condition, see [36, 68, 69] and [43, Chap. 10]), we expect our findings to have a bearing there, too.

1.1. Prior work

We start by recalling some fundamental results and notions. More details and references can be found in the outstanding surveys and monographs [3, 10, 71, 75].

1.1.1. Ulam’s problem

The study of the behavior as n grows large dates back to Ulam [80] in 1961, who mentioned that Monte-Carlo computations of E. Neighbor would indicate ${\mathbb E}(L_n)\approx 1.7\sqrt {n}$ . Ulam continued by stating, ‘Another question of interest would be to find the distribution of the length of the maximum monotone subsequence around this average’.

Refined numerical experiments by Baer and Brock [5] in 1968 suggested that ${\mathbb E}(L_n) \sim 2\sqrt {n}$ might be the precise leading order. In a 1970 lecture, Hammersley [50] presented a proof, based on subadditive ergodic theory, that the limit $c = \lim _{n\to \infty } {\mathbb E}(L_n)/\sqrt {n}$ exists. Finally, in 1977, Vershik and Kerov [82] as well as Logan and Shepp [59] succeeded in proving $c=2$ .

1.1.2. Poissonization

A major tool used by Hammersley was a random process that is, basically, equivalent to the following Poissonization of the random variable $L_n$ : by drawing from the different permutation groups independently and by taking $N_r \in \{0,1,2,\ldots \}$ to be a further independent random variable with a Poisson distribution of intensity $r>0$ , the combined random variable $L_{N_r}$ is distributed according to

$$\begin{align*}{\mathbb P}(L_{N_r} \leqslant l) = e^{-r} \sum_{n=0}^{\infty} {\mathbb P}(L_n\leqslant l)\frac{r^n}{n!} =: P(r;l). \end{align*}$$

The entire function ² $P(z;l)$ is the Poisson generating function of the sequence ${\mathbb P}(L_n \leqslant l)$ ( $n=0,1,\ldots $ ), and $f(z;l):=e^{z}P(z;l)$ is the corresponding exponential generating function. As it turns out, it is much simpler to analyze the Poissonized distribution of $L_{N_r}$ as the intensity r grows large than the original distribution of $L_n$ as n grows large.

There is, however, also a way back from $L_{N_r}$ to $L_n$ – namely, the expected value of the Poisson distribution being ${\mathbb E}(N_r) = r$ , combined with some level of concentration, suggests

$$\begin{align*}{\mathbb P}(L_n \leqslant l) \approx P(n;l) \end{align*}$$

when $n\to \infty $ while l is kept near the mode of the distribution. Being a Tauberian result, such a de-Poissonization is subject to additional conditions, which we will discuss in a moment.

Starting in the early 1990s, the Poisson generating function $P(z;l)$ (or the exponential one to the same end) has been represented in terms of one of the following interrelated forms:

• • a Toeplitz determinant in terms of modified Bessel functions [47],

• • Fredholm determinants of various (discrete) integral operators [8, 9, 21, 22, 56],

• • a unitary group integral [70].

A particular case of those representations plays a central role in our study – namely, ³

(1.1) $$ \begin{align} P(r;l) = E_2^{\text{hard}}(4r;l), \end{align} $$

where ⁴ $E_2^{\text {hard}}(s;\nu )$ denotes the probability that, in the hard-edge scaling limit, the scaled smallest eigenvalue of the Laguerre unitary ensemble (LUE) with real parameter $\nu> 0$ is bounded from below by $s\geqslant 0$ . This probability is known to be given in terms of a Fredholm determinant (see [42]):

(1.2) $$ \begin{align} E_2^{\text{hard}}(s;\nu) = \det(I - K_\nu^{\text{Bessel}})\big|_{L^2(0,s)}, \end{align} $$

where $K_\nu $ denotes the Bessel kernel in $x,y \geqslant 0$ (for the integral formula see [79, Eq. (2.2)]):

(1.3) $$ \begin{align} K_\nu^{\text{Bessel}}(x,y) := \frac{J_\nu(\sqrt{x}) \sqrt{y}J_\nu^{\prime}(\sqrt{y})-J_\nu(\sqrt{y}) \sqrt{x}J_\nu^{\prime}(\sqrt{x})}{2(x-y)} = \frac{1}{4}\int_0^1 J_\nu(\sqrt{\sigma x})J_\nu(\sqrt{\sigma y})\,d\sigma. \end{align} $$

Obviously, the singularities at the diagonal $x=y$ are removable.

The work of Tracy and Widom [79] establishes that the Fredholm determinant (1.2) can be expressed in terms of Painlevé III. Recently, based on Okamoto’s Hamiltonian $\sigma $ -PIII $'$ framework, Forrester and Mays [45] used that connection to compile a table of the exact rational values of ${\mathbb P}(L_n \leqslant l)$ for up to $n=700$ , ⁵ whereas in our work [19], based on an equivalent representation in terms of a Chazy I equation, we have compiled such a table ⁶ for up to $n=1000$ .

In their seminal 1999 work [7], by relating the representation of $P(z;l)$ in terms of the Toeplitz determinant to the machinery of Riemann–Hilbert problems and studying the underlying double-scaling limit by the Deift–Zhou method of steepest descent, Baik, Deift and Johansson answered Ulam’s question and proved that, for t being any fixed real number,

(1.4) $$ \begin{align} \lim_{r\to\infty}{\mathbb P}\left(\frac{L_{N_r}-2\sqrt{r}}{r^{1/6}}\leqslant t\right) = F(t), \end{align} $$

where F is the GUE Tracy–Widom distribution – that is, the distribution which expresses, among many other limit laws, 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. As for the Poissonized length distribution itself, the Tracy–Widom distribution can be represented in terms of a Fredholm determinant (see [42]) – namely,

(1.5) $$ \begin{align} F(t) = \det(I - K_0)|_{L^2(t,\infty)} \end{align} $$

where $K_0$ denotes the Airy kernel in $x,y\in {\mathbb R}$ (for the integral formula see [78, Eq. (4.5)]):

(1.6) $$ \begin{align} K_0(x,y) := \frac{\operatorname{Ai}(x)\operatorname{Ai}'(y)-\operatorname{Ai}'(x)\operatorname{Ai}(y)}{x-y} = \int_0^\infty \operatorname{Ai}(x+ \sigma) \operatorname{Ai}(y+\sigma)\,d \sigma. \end{align} $$

Obviously, also in this case, the singularities at the diagonal $x=y$ are removable. Since the limit distribution in (1.4) is continuous, by a standard Tauberian follow-up [81, Lemma 2.1] of the Portmanteau theorem, the limit law holds uniformly in t.

In 2003, Borodin and Forrester gave an alternative proof of (1.4) which is based on studying the hard-to-soft edge transition of LUE for $\nu \to \infty $ in form of the limit law [21, Thm. 1]

(1.7) $$ \begin{align} \lim_{\nu\to\infty} E_2^{\text{hard}}\left(\big(\nu-t (\nu/2)^{1/3}\big)^2;\nu\right) = F(t) \end{align} $$

(see also [43, §10.8.4]), which will be the starting point of our study. Still, there are other proofs of (1.4) based on representations in terms of Fredholm determinants of further (discrete) integral operators; for expositions and references, see the monographs [10, 71].

1.1.3. De-Poissonization

In the literature, the de-Poissonization of the limit law (1.4) has so far been based exclusively on variants of the following lemma (cf. [10, Cor. 2.5], originally stated as [55, Lemma 2.5]), which uses monotonicity as the underlying Tauberian condition.

Lemma 1.1 (Johansson’s de-Poissonization lemma [55]).

Suppose the sequence $a_n$ of probabilities $0\leqslant a_n \leqslant 1$ satisfies the monotonicity condition $a_{n+1} \leqslant a_n$ for all $n=0,1,2,\ldots $ and denote its Poisson generating function by

(1.8) $$ \begin{align} P(z) = e^{-z} \sum_{n=0}^\infty a_n \frac{z^n}{n!}. \end{align} $$

Then, for $s\geqslant 1$ and $n\geqslant 2\,$ , ⁷ $P\big (n + 2\sqrt {s n \log n} \big ) - n^{-s} \leqslant a_n \leqslant P \big (n - 2\sqrt {s n \log n} \big ) + n^{-s}$ .

After establishing the Tauberian condition of monotonicity and applying a variant of Lemma 1.1 to (1.4), Baik, Deift and Johansson [7, Thm. 1.1] got

(1.9) $$ \begin{align} \lim_{n\to\infty}{\mathbb P}\left(\frac{L_{n}-2\sqrt{n}}{n^{1/6}}\leqslant t\right) = F(t), \end{align} $$

which holds uniformly in t for the same reasons as given above. (The simple calculations based on Lemma 1.1 are given in [10, p. 239]; note that the uniformity of the limit law (1.4) is used there without explicitly saying so.) Adding tail estimates to the picture, those authors were also able to lift the limit law to the moments and got, expanding on Ulam’s problem, that the expected value satisfies [7, Thm. 1.2]

(1.10) $$ \begin{align} {\mathbb E}(L_n) = 2\sqrt{n} + M_1 n^{1/6} + o(n^{1/6}),\qquad M_1 = \int_{-\infty}^\infty t F'(t)\,dt \approx -1.7711. \end{align} $$

1.1.4. Expansions

To our knowledge, only for the Poissonized limit law (1.4) a finite-size correction term has been rigorously established prior to the present paper ⁸ – namely, as a by-product along the way of their study of the limiting distribution of maximal crossings and nestings of Poissonized random matchings, Baik and Jenkins [11, Thm. 1.3] obtained (using the machinery of Riemann–Hilbert problems and Painlevé representations of the Tracy–Widom distribution), as $r\to \infty $ with t being any fixed real number,

(1.11a) $$ \begin{align} {\mathbb P}\left(\frac{L_{N_r} - 2\sqrt{r}}{r^{1/6}}\leqslant t\right) = F\big(t^{(r)}\big) - \frac{1}{10}\Big(F"(t) + \frac{t^2}{6} F'(t)\Big)r^{-1/3} + O(r^{-1/2}), \end{align} $$

where (with $\lfloor \cdot \rfloor $ denoting the Gauss bracket)

(1.11b) $$ \begin{align} t^{(r)} := \frac{\lfloor 2\sqrt{r} + t r^{1/6}\rfloor - 2\sqrt{r}}{r^{1/6}}. \end{align} $$

However, even if there is enough uniformity in this result and the option to Taylor expand the Poisson generating function $P(r)$ at n with a uniform bound while l is kept near the mode of the distributions (see Section 5.1 for details on this option), the sandwiching in Johansson’s de-Poissonization Lemma 1.1 does not allow us to obtain a result better than (cf. [11, §9])

(1.12). $$ \begin{align} {\mathbb P}\left(\frac{L_n - 2\sqrt{n}}{n^{1/6}}\leqslant t\right) = F\big(t^{(n)}\big) + O\big(n^{-1/6} \sqrt{\log n}\,\big). \end{align} $$

In their recent study of finite-size effects, Forrester and Mays [45, Prop. 1.1] gave a different proof of (1.11) based on the Bessel kernel determinant (1.2). Moreover, suggested by exact data for $n=700$ and a Monte-Carlo simulation for $n=20\,000$ , they were led to conjecture [45, Conj. 4.2]

(1.13) $$ \begin{align} {\mathbb P}\left(\frac{L_n - 2\sqrt{n}}{n^{1/6}}\leqslant t\right) = F\big(t^{(n)}\big) + F^D_1(t)n^{-1/3} + \cdots \end{align} $$

with the approximate graphical form of $F^D_1(t)$ displayed in [45, Fig. 7].

The presence of the Gauss bracket in $t^{(r)}$ and $t^{(n)}$ , while keeping t at other places of the expansions (1.11) and (1.13), causes undesirable effects in the error terms (see Remark 4.2 below for a detailed discussion). Therefore, in our work [19] on a Stirling-type formula approximating the distribution ${\mathbb P}(L_n \leqslant l)$ , we suggested to use the integer l in the continuous expansion terms instead of introducing the continuous variable t into the discrete distribution in the first place, with the latter variant turning the discrete distribution into a piecewise constant function of t. By introducing the scaling

$$\begin{align*}t_\nu(r) := \frac{\nu-2\sqrt{r}}{r^{1/6}} \qquad (r>0), \end{align*}$$

we were led (based on numerical experiments using the Stirling-type approximation for n getting as large as $10^{10}$ ) to conjecture the expansion

$$\begin{align*}{\mathbb P}(L_n \leqslant l) = F(t) + F_{1}^D(t)\, n^{-1/3} + F_{2}^D(t)\, n^{-2/3} + O(n^{-1})\,\Big|_{t=t_l(n)}, \end{align*}$$

displaying the graphical form of the functions $F^D_1(t)$ , $F^D_2(t)$ in the left panels of [19, Figs. 4/6]. Moreover, as a note added in proof (see [19, Eq. (11)]), we announced that inserting the Baik–Jenkins expansion (1.11) into the Stirling-type formula and using its (numerically observed) apparent order $O(n^{-2/3})$ of approximation would yield the functional form of $F_1^D$ to be

(1.14) $$ \begin{align} F_1^D(t) = -\frac{1}{10}\left(6 F"(t) + \frac{t^2}{6} F'(t) \right). \end{align} $$

The quest for a proof, and for a similar expression for $F_2^D$ , motivated our present work.

1.2. The new findings of the paper

In the analysis of algorithms in theoretical computer science, or the enumeration of combinatorial structures to the same end, the original enumeration problem is often represented in form of recurrences or functional/differential equations. For instance, this situation arises in a large class of algorithms involving a splitting process, trees or hashing. Embedding such processes into a Poisson process ⁹ often leads to more tractable equations, so that sharp tools for a subsequent de-Poissonization were developed in the 1990s; for references and details, see [53, 54, 77] and Appendix A.1. In particular, if the Poisson generating function $P(z)$ of a sequence of real $a_n>0$ , as defined in (1.8), is an entire function, an application of the saddle point method to the Cauchy integral

$$\begin{align*}a_n = \frac{n!}{2\pi i} \oint P(z) e^z \frac{dz}{z^{n+1}} \end{align*}$$

yields, under suitable growth conditions on $P(z)$ as $z\to \infty $ in the complex plane, the Jasz ¹⁰ expansion

$$\begin{align*}a_n \sim P(n) + \sum_{j=2}^\infty b_j(n) P^{(j)}(n), \end{align*}$$

where the polynomial coefficients $b_j(n)$ are the diagonal Poisson–Charlier polynomials (that is, with intensity $r=n$ ). In Appendix A.1, we give a heuristic derivation of that expansion and recall, in the detailed form of Theorem A.2, a specific analytic de-Poissonization result from the comprehensive memoir [53] of Jacquet and Szpankowski – a result which applies to a family of Poisson generating functions at once, providing uniform error bounds.

Now, the difficult part of applying Theorem A.2 is checking the Tauberian growth conditions in the complex plane, which are required to hold uniformly for the family of Poisson generating functions (recall that, in the case of the longest increasing subsequence problem, $P(z)=P(z;l)$ depends on the integer l near the mode of the length distribution). After observing a striking similarity of those growth conditions with the notion of H-admissibility for the corresponding exponential generating function (as introduced by Hayman in his memoir [51] on the generalization of Stirling’s formula), a closer look at the proof of Hayman’s [51, Thm. XI] revealed the following result (see Theorem A.5 for a precise quantitative statement):

The family of all entire functions of genus zero which have, for some $\epsilon>0,$ no zeros in the sector $|\!\arg z| \leqslant \pi /2+ \epsilon $ satisfies a universal bound that implies Tauberian growth conditions suitable for analytic de-Poissonization.

However, in our work [19, Thm. 2.2] on Stirling-type formulae for the problem of longest increasing subsequences, when proving the H-admissibility of the exponential generating functions $f(z;l)=e^z P(z;l)$ (for each l), we had established, based on the representation [70] of $P(z;l)$ as a group integral, the following:

For any integer $l\geqslant 0$ and any $\delta>0$ , the exponential generating function $f(z;l)$ is an entire function of genus zero having at most finitely many zeros in the sector $|\!\arg z| \leqslant \pi - \delta $ , none of them being real.

Under the reasonable assumption (supported by numerical experiments) that those finitely many complex zeros do not come too close to the real axis and do not grow too fast as $n\to \infty $ while l stays near the mode of the length distribution, the uniformity of the Tauberian growth conditions can be preserved (see Corollary A.7 for the technical details). We call this assumption the tameness hypothesis ¹¹ concerning the zeros of the family of $P(z;l)$ .

Subject to the tameness hypothesis, the main result of the present paper, Theorem 5.1, gives the asymptotic expansion

$$\begin{align*}{\mathbb P}(L_n \leqslant l) = F(t) + \sum_{j=1}^m F_{j}^D(t)\, n^{-j/3} + O\big(n^{-(m+1)/3}\big)\,\bigg|_{t=t_l(n)}, \end{align*}$$

which is uniformly valid when $n,l \to \infty $ while $t_l(n)$ stays bounded ¹² and the $F_j^D$ are certain smooth functions.

Finally, now without any detour via the Stirling-type formula, Theorem 5.1 confirms that the expansion term $F_1^D$ is given by (1.14), indeed, and yields the striking formula

$$\begin{align*}F_2^D(t) = \Big(-\frac{139}{350} + \frac{2t^3}{1575}\Big) F'(t) + \Big(-\frac{43t}{350} + \frac{t^4}{7200}\Big) F"(t) + \frac{t^2}{100} F"'(t) + \frac{9}{50} F^{(4)}(t); \end{align*}$$

see (3.4) for a display of a similarly structured expression for $F_3^D(t)$ .

Put to the extreme, with the help of a CAS such as Mathematica, the methods of the present paper can be used to calculate the concrete functional form of the expansion terms for up to $m=10$ and larger. ¹³ In all cases inspected, we observe that the expansion terms take the form of a linear combination of higher order derivatives of the limit law (that is, the Tracy–Widom distribution F) with certain rational polynomials as coefficients; we conjecture that this is generally true for the problem at hand.

1.3. Generalization: involutions, orthogonal and symplectic ensembles

In our subsequent work [18], we present a similar structure for the expansion terms relating to longest monotone subsequences in (fixed-point free) random involutions, F then being one of the Tracy–Widom distributions for $\beta =1$ or $\beta =4$ . The limit laws were first obtained by Baik and Rains [12], using the machinery of Riemann–Hilbert problems, and later reclaimed by Borodin and Forrester [21] through establishing hard-to-soft edge transition laws for LOE and LSE similar to (1.7). In [18], we derive the asymptotic expansions by using determinantal formulae [29, 37] of the hard and soft edge limits for $\beta =1,4$ while taking advantage of their algebraic interrelations with the $\beta =2$ case studied in the present paper (basically, the expansions of the hard-to-soft edge limit laws for $\beta =1,4$ turn out to correspond to a certain factorization of the $\beta =2$ case). This is then followed by applying analytic de-Poissonization, once again subject to a tameness hypothesis.

1.4. Organization of the paper

The paper splits into two parts: a first one, where all results and proofs are unconditional, addressing asymptotic expansions of Fredholm determinants, of the hard-to-soft edge transition law and the Poissonized length distribution; and a second one, where we restrict ourselves to assuming the tameness hypothesis when addressing analytic de-Poissonization and its various consequences.

Part I: Unconditional results. In Section 2, we start with a careful discussion of expansions of perturbed Airy kernel determinants. We stress the importance of such kernel expansions to be differentiable (i.e., one can differentiate into the error term) to easily lift the error bounds to trace norms. The subtle, but fundamental difficulty of such a lift seems to have been missed, more often than not, in the existing literature on convergence rates and expansions of limit laws in random matrix theory (notable exceptions are, for example, [33, 57, 58]).

In the rather lengthy Section 3, we study the asymptotic expansion of the Borodin–Forrester hard-to-soft edge transition law (1.7). It is based on a uniform version of Olver’s asymptotic expansion of Bessel functions of large order in the transition region, which we discuss in Appendix A.3. In Section 3, we lay the foundational work for the concrete functional form of all subsequent finite-size correction terms. We reduce the complexity of computing these terms by using a coordinate transform on the level of kernels to simplify the kernel expansion – a coordinate transform which gets subsequently reversed on the level of the distributions.¹⁴ As yet another application of that technique, we simplify in Section 3.4 the finite-size correction terms of Choup [24] to the soft-edge limit law of GUE and LUE.

In Section 4, we expand the Poissonized length distribution, thereby generalizing the result (1.12) of Baik and Jenkins. In Section 4, we also discuss the potentially detrimental effect of using Gauss brackets, as in (1.12), alongside with the continuous variable t in the expansion terms.

Part II: Results based on the tameness hypothesis. In Section 5, we state and prove the main result of the paper: the expansion of the Baik–Deift–Johansson limit law (1.9) of the length distribution. Here we use the Jasz expansion of analytic de-Poissonization (as detailed in Appendix A.1). The universal bounds for entire functions of genus zero, which are used to prove the Tauberian growth conditions in the complex plane, and their relation to the theory of H-admissibility are prepared for in Appendix A.2. Additionally, in Section 5, we discuss the modifications that apply to the discrete density ${\mathbb P}(L_n=l)$ (that is, to the PDF of the length distribution).

In Section 6, we apply our findings to the asymptotic expansion of the Stirling-type formula which we introduced in [19] as an accurate tool for the numerical approximation of the length distribution. Subject to the tameness hypothesis, we prove the observation [19, Eq. (8b)] about a leading $O(n^{-2/3})$ error of that formula.

Finally, in Section 7, we study the asymptotic expansion of the expected value and variance. Based on a reasonable hypothesis about some uniformity in the tail bounds, we add several more concrete expansion terms to the Baik–Deift–Johansson solution (1.10) of Ulam’s problem.

Part I: Unconditional results

2. Expansions of perturbed Airy kernel determinants

In Section 3, we will get, with m being some non-negative integer and $t_0$ some real number, kernel expansions of the form

(2.1) $$ \begin{align} K_{(h)}(x,y) = K_0(x,y) + \sum_{j=1}^m h^j K_j(x,y) + h^{m+1} O\big(e^{-(x+y)}\big), \end{align} $$

which are

• • uniformly valid for $t_0 \leqslant x,y < c h^{-1}$ as $h\to 0^+$ , where $c>0$ is some constant;

• • repeatedly differentiable w.r.t. x, y as uniform expansions under the same conditions.

Here, $K_{(h)}$ is a family of smooth kernels, $K_0$ denotes the Airy kernel (1.6) and the $K_j(x,y)$ are finite sums of rank one kernels $u(x)v(y)$ , with factors $u(\xi )$ , $v(\xi )$ of the functional form

(2.2) $$ \begin{align} p(\xi) \operatorname{Ai}(\xi)\quad \text{or}\quad p(\xi) \operatorname{Ai}'(\xi), \end{align} $$

where $p(\xi )$ is a polynomial in $\xi $ . Since the existing literature tends to neglect the issue of estimating trace norms in terms of kernel bounds, this section aims at establishing a relatively easy framework for lifting such an expansion to one of the Fredholm determinant.

2.1. Bounds on the kernels and induced trace norms

Bounds on the kernels $K_j$ , and on the trace norms of the induced integral operators, can be deduced from the estimates, ¹⁵ with p being any polynomial,

$$\begin{align*}|p(\xi)| \cdot \max\big(|\operatorname{Ai}(\xi)|,|\operatorname{Ai}'(\xi)|\big) \leqslant a_p e^{-\xi}\quad (\xi\in {\mathbb R}), \end{align*}$$

where the constant $a_p$ does only depend on p (note that we can take $a_p=1$ when $p(\xi )\equiv 1$ ). This way, we get from (1.6)

$$\begin{align*}|K_0(x,y)| \leqslant \int_0^\infty |\operatorname{Ai}(x+\sigma) \operatorname{Ai}(y+\sigma)|\,d\sigma \leqslant e^{-(x+y)} \int_0^\infty e^{-2\sigma}\,d\sigma = \tfrac{1}{2}e^{-(x+y)} \quad (x,y \in {\mathbb R}) \end{align*}$$

and, for $1\leqslant j\leqslant m$ , constants $c_j$ such that

$$\begin{align*}|K_j(x,y)| \leqslant c_j e^{-(x+y)} \quad (x,y \in {\mathbb R}). \end{align*}$$

For a given continuous kernel $K(x,y)$ , we denote the induced integral operator on $L^2(t,c h^{-1})$ by $\bar {\mathbf K}$ and the one on $L^2(t,\infty )$ , if defined, by ${\mathbf K}$ (suppressing the dependence on t in both cases). The spaces of trace class and Hilbert–Schmidt operators acting on $L^2(t,s)$ are written as ${\mathcal J}^p(t,s)$ with $p=1$ and $p=2$ . By using the orthogonal projection of $L^2(t,\infty )$ onto the subspace $L^2(t,ch^{-1})$ , we see that

(2.3) $$ \begin{align} \|\bar {\mathbf K}\|_{{\mathcal J}^1(t,ch^{-1})} \leqslant \|{\mathbf K}\|_{{\mathcal J}^1(t,\infty)}. \end{align} $$

The Airy operator ${\mathbf K}_0$ (being by (1.6) the square of the Hilbert–Schmidt operator ${\mathbf A}_t$ with kernel $\operatorname {Ai}(x+y-t)$ ) and the expansion operators ${\mathbf K}_j$ ( $j\geqslant 1$ ) (being finite rank operators) are trace class on the space $L^2(t,\infty )$ . Their trace norms are bounded by

$$ \begin{align*} \|{\mathbf K}_0\|_{{\mathcal J}^1(t,\infty)} &\leqslant \|{\mathbf A}_t\|_{{\mathcal J}^2(t,\infty)}^2 = \int_t^\infty \int_t^\infty |\operatorname{Ai}(x+y-t)|^2 \,dx\,dy \\ &= \int_0^\infty \int_0^\infty |\operatorname{Ai}(x+y+t)|^2 \,dx\,dy \leqslant e^{-2t} \int_0^\infty\int_0^\infty e^{-2(x+y)}\,dx\,dy = \tfrac{1}{4} e^{-2t} \quad (t\in{\mathbb R}), \end{align*} $$

and, likewise with some constants $c_j^*$ , by

$$\begin{align*}\|{\mathbf K}_j\|_{{\mathcal J}^1(t,\infty)} \leqslant c_j^* e^{-2t}\qquad (1\leqslant j\leqslant m, \; t\in {\mathbb R}). \end{align*}$$

Here, we have used

$$\begin{align*}\|u \otimes v\|_{{\mathcal J}^1(t,\infty)} \leqslant \|u\|_{L^2(t,\infty)}\|v\|_{L^2(t,\infty)} \leqslant \frac{a_p a_q}{2} e^{-2t} \quad (t\in {\mathbb R}) \end{align*}$$

for factors $u(\xi )$ , $v(\xi )$ of the form (2.2) with some polynomials p and q.

However, there is in general no direct relation between kernel bounds and bounds of the trace norm of induced integral operators (see [74, p. 25]). Writing the kernel of the error term in (2.1), and its bound, in the form

(2.4) $$ \begin{align} h^{m+1} R_{m+1,h}(x,y) = h^{m+1} O\big(e^{-(x+y)}\big) \end{align} $$

does therefore not offer, as it stands, any direct method of lifting the bound to trace norm. ¹⁶

By taking, however, the explicitly assumed differentiability of the kernel expansion into account, there is a constant $c_\sharp $ such that

$$\begin{align*}S_{m+1,h}(x,y):=\partial_y R_{m+1,h}(x,y),\quad |S_{m+1,h}(x,y)| \leqslant c_\sharp e^{-(x+y)} \end{align*}$$

holds true for all $t_0 \leqslant x,y < c h^{-1}$ and $0<h\leqslant h_0$ ( $h_0$ chosen sufficiently small). If we denote an indicator function by $\chi $ and choose some $c_* < c$ close to c, integration gives

$$\begin{align*}R_{m+1,h}(x,y) = R_{m+1}(x,c_*h^{-1}) -\int_t^{c_* h^{-1}} S_{m+1,h}(x,\sigma)e^{\sigma/2} \cdot e^{-\sigma/2} \chi_{[t,\sigma]}(y)\,d\sigma \end{align*}$$

if $t_0 \leqslant x,y \leqslant c_* h^{-1}$ . This shows that the induced integral operator on $L^2(t,c_* h^{-1})$ , briefly written accordingly as

$$\begin{align*}\bar{\mathbf R}_{m+1,h} = \bar{\mathbf R}_{m+1,h}^{\flat} + \bar{\mathbf R}_{m+1,h}^{\sharp}, \end{align*}$$

is trace class since it is the sum of a rank-one operator and a product of two Hilbert–Schmidt operators. In fact, for $t_0 \leqslant t \leqslant c_* h^{-1}$ , the trace norms of those terms are bounded by (denoting the implied constant in (2.4) by $c_\flat $ )

$$ \begin{align*} \| \bar{\mathbf R}_{m+1,h}^{\flat}\|_{{\mathcal J}^1(t,c_* h^{-1})}^2 &\leqslant \|1\|_{L^2(t,c_*h^{-1})}^2 \cdot \|R_{m+1,h}(\,\cdot\,,c_*h^{-1})\|_{L^2(t,c_*h^{-1})}^2 \\ &\leqslant (ch^{-1} - t_0) \cdot c_{\flat}^2 e^{-2c_*h^{-1}} \int_t^\infty e^{-2x}\,dx = h^{-1}e^{-2c_*h^{-1}} O(e^{-2t}) \end{align*} $$

and

$$ \begin{align*} \| \bar{\mathbf R}_{m+1,h}^{\sharp}\|_{{\mathcal J}^1(t,c_* h^{-1})}^2 &\leqslant \left(\int_t^{c_*h^{-1}}\int_t^{c_*h^{-1}} S_{m+1,h}(x,y)^2e^{y}\,dx\,dy\right)\cdot \left(\int_t^{c_*h^{-1}}\int_t^x e^{-x}\,dx\,dy\right)\\ &\leqslant c_\sharp^2 \left(\int_t^\infty\int_t^\infty e^{-2x-y}\,dx\,dy\right)\cdot \left(\int_t^\infty\int_t^x e^{-x}\,dx\,dy\right) = O(e^{-4t}). \end{align*} $$

Here, the implied constants are independent of the particular choice of $c_*$ . By noting that the orthogonal projection of $L^2(t,ch^{-1})$ onto $L^2(t,c_*h^{-1})$ converges, as $c_* \to c$ , to the identity in the strong operator topology, we obtain by a continuity theorem ¹⁷ of Grümm [49, Thm. 1]

$$ \begin{align*} \| \bar{\mathbf R}_{m+1,h}\|_{{\mathcal J}^1(t,c h^{-1})} & = \lim_{c_*\to c} \| \bar{\mathbf R}_{m+1,h}\|_{{\mathcal J}^1(t,c_* h^{-1})} \\ &\leq \limsup_{c_*\to c}\| \bar{\mathbf R}_{m+1,h}^{\flat}\|_{{\mathcal J}^1(t,c_* h^{-1})} + \limsup_{c_*\to c}\| \bar{\mathbf R}_{m+1,h}^{\sharp}\|_{{\mathcal J}^1(t,c_* h^{-1})}\\ &= h^{-1/2} e^{-c h^{-1}} O(e^{-t}) + O(e^{-2t}). \end{align*} $$

We have thus lifted the kernel expansion (2.1) to an operator expansion in ${\mathcal J}^1(t,ch^{-1})$ – namely,

(2.5a) $$ \begin{align} \bar{\mathbf K}_{(h)} = \bar{\mathbf K}_0 + \sum_{j=1}^m h^j \bar{\mathbf K}_j + h^{m+1} \bar{\mathbf R}_{m+1,h}, \end{align} $$

with the bounds (recall (2.3) and observe that we can absorb $h^{m+1/2} e^{-c h^{-1}}$ into $O(e^{-c h^{-1}})$ )

(2.5b) $$ \begin{align} \|{\mathbf K}_j \|_{{\mathcal J}^1(t,\infty)} &= O(e^{-2t}), \end{align} $$

(2.5c) $$ \begin{align} h^{m+1}\|\bar{\mathbf R}_{m+1,h}\|_{{\mathcal J}^1(t,ch^{-1})} &= h^{m+1} O(e^{-2t}) + e ^{-c h^{-1}} O(e^{-t}), \end{align} $$

uniformly valid for $t_0 \leqslant t < c h^{-1}$ as $h\to 0^+$ .

2.2. Fredholm determinants

Given a continuous kernel $K(x,y)$ on the bounded rectangle $t_0 \leqslant x,y \leqslant t_1$ , the Fredholm determinant (cf. [4, §3.4])

(2.6) $$ \begin{align} \det(I - K)|_{L^2(t,s)} := \sum_{m=0}^\infty \frac{(-1)^m}{m!} \int_t^s \cdots\int_t^s \det_{j,k=1}^m K(x_j,x_k) \,dx_1\cdots\, dx_m \end{align} $$

is well-defined for $t_0\leqslant t < s \leqslant t_1$ . If, as is the case for the kernels $K_j$ (see Section 2.1), the kernel is also continuous for all $x,y\geqslant t_0$ and there is a weighted uniform bound of the form

$$\begin{align*}\sup_{x,y \geqslant t_0} e^{x+y} \,|K(x,y)| \leqslant M < \infty, \end{align*}$$

the Fredholm determinant (2.6) is well-defined even if we choose $s=\infty $ (this can be seen by writing the integrals in terms of the weighted measure $e^{-x}\,dx$ ; cf. [4, §3.4]).

Now, if the induced integral operator ${\mathbf K}|_{L^2(t,s)}$ on $L^2(t,s)$ is trace class, the Fredholm determinant can be expressed in terms of the operator determinant (cf. [74, Chap. 3]):

(2.7) $$ \begin{align} \det(I - K)|_{L^2(t,s)} = \det\big({\mathbf I}- {\mathbf K}|_{L^2(t,s)}\big). \end{align} $$

Using the orthogonal decomposition

(2.8) $$ \begin{align} L^2(t,\infty) = L^2(t,ch^{-1}) \oplus L^2(ch^{-1},\infty) \end{align} $$

and writing the operator expansion (2.5) in the form

$$\begin{align*}\bar{\mathbf K}_{(h)} = \bar{\mathbf K}_{m,h} + h^{m+1} \bar{\mathbf R}_{m+1,h}, \qquad \bar{\mathbf K}_{m,h}:=\bar{\mathbf K}_0 + \sum_{j=1}^m h^j \bar{\mathbf K}_j, \end{align*}$$

we get from the error estimate in (2.5c), by the local Lipschitz continuity of operator determinants w.r.t. trace norm (cf. [74, Thm. 3.4]), that

$$ \begin{align*} \det\big(I - K_{(h)}\big)|_{L^2(t,ch^{-1})} &= \det({\mathbf I}- \bar{\mathbf K}_{(h)}) = \det\big({\mathbf I}- \bar{\mathbf K}_{m,h}\big) + h^{m+1} O(e^{-2t}) + e ^{-c h^{-1}} O(e^{-t}) \\ &= \det\big({\mathbf I}- {\mathbf K}_{m,h}\big) + h^{m+1} O(e^{-2t}) + e ^{-c h^{-1}} O(e^{-t}), \end{align*} $$

uniformly for $t_0 \leqslant t < ch^{-1}$ as $h\to 0^+$ . Here, the last estimate comes from a block decomposition of $\bar {\mathbf K}_{(h)}$ according to (2.8) and applying the trace norm bounds of ${\mathbf K}_j$ in Section 2.1 to the boundary case $t=ch^{-1}$ , which yields

$$\begin{align*}\det\big({\mathbf I}- \bar{\mathbf K}_{m,h}\big) = \det\big({\mathbf I}- {\mathbf K}_{m,h}\big) + O(e^{-2ch^{-1}}); \end{align*}$$

note that the error term of order $O(e^{-2ch^{-1}})$ can be absorbed into $e ^{-c h^{-1}} O(e^{-t})$ if $t < c h^{-1}$ .

2.3. Expansions of operator determinants

Plemelj’s formula gives, for trace class perturbations $\mathbf E$ of the identity on $L^2(t,\infty )$ bounded by $\|{\mathbf E}\|_{{\mathcal J}^1(t,\infty )} < 1$ , the convergent series expansion (cf. [74, Eq. (5.12)])

(2.9) $$ \begin{align} \det({\mathbf I} - {\mathbf E}) &= \exp\left(-\sum_{n=1}^\infty n^{-1} \operatorname{tr}({\mathbf E}^n)\right) \\ = 1- \operatorname{tr} {\mathbf E} + \frac{1}{2}\left((\operatorname{tr} {\mathbf E})^2-\operatorname{tr}({\mathbf E}^2)\right) + \frac{1}{6}&\left(-(\operatorname{tr} {\mathbf E})^3 + 3 \operatorname{tr}{\mathbf E}\operatorname{tr} {\mathbf E}^2 - 2 \operatorname{tr} {\mathbf E}^3 \right)+ O\Big(\|{\mathbf E}\|_{{\mathcal J}^1(t,\infty)}^4\Big).\notag \end{align} $$

Thus, since ${\mathbf I}-{\mathbf K_0}$ is invertible with a uniformly bounded inverse as $t\geqslant t_0$ , ¹⁸ we have

$$\begin{align*}\det\big({\mathbf I}-{\mathbf K}_{m,h}\big) = \det({\mathbf I}-{\mathbf K_0}) \det\big({\mathbf I}-{\mathbf E}_{m,h} \big), \qquad {\mathbf E}_{m,h} := \sum_{j=1}^m h^j ({\mathbf I}-{\mathbf K_0})^{-1}{\mathbf K_j}, \end{align*}$$

with the trace norm of ${\mathbf E}_{m,h} $ being bounded as follows (using the results of Section 2.1 and observing that the trace class forms an ideal within the algebra of bounded operators):

$$\begin{align*}\|{\mathbf E}_{m,h} \|_{{\mathcal J}^1(t,\infty)} \leqslant \|({\mathbf I}-{\mathbf K_0})^{-1}\|\, \frac{ h}{1-h}\cdot O\big(e^{-2t}\big) = h\cdot O(e^{-2t}), \end{align*}$$

uniformly for $t\geqslant t_0$ as $h\to 0^+$ . By (2.9), this implies, uniformly under the same conditions,

$$\begin{align*}\det\big({\mathbf I}-{\mathbf E}_{m,h} \big) = 1 + \sum_{j=1}^m{d_j(t) h^j} + h^{m+1}\cdot O(e^{-2t}). \end{align*}$$

Here, $d_j(t)$ depends smoothly on t and satisfies the (right) tail bound $d_j(t) = O(e^{-2t})$ . If we write briefly

$$\begin{align*}{\mathbf E}_j = ({\mathbf I}-{\mathbf K_0})^{-1}{\mathbf K_j}, \end{align*}$$

the first few cases are given by the expressions (because the traces are taken for trace class operators acting on $L^2(t,\infty )$ they depend on t)

(2.10a) $$ \begin{align} d_1(t) &= - \operatorname{tr} {\mathbf E}_1, \end{align} $$

(2.10b) $$ \begin{align} d_2(t) &= \frac{1}{2}(\operatorname{tr}{\mathbf E_1})^2 - \frac{1}{2}\operatorname{tr}{\mathbf E}_1^2 - \operatorname{tr}{\mathbf E_2}, \end{align} $$

(2.10c) $$ \begin{align} d_3(t) &= -\frac{1}{6}(\operatorname{tr}{\mathbf E_1})^3 +\frac{1}{2}\operatorname{tr} {\mathbf E}_1\operatorname{tr} {\mathbf E}_1^2 -\frac{1}{2}\operatorname{tr}({\mathbf E}_1{\mathbf E}_2 + {\mathbf E}_2 {\mathbf E}_1) \end{align} $$

$$\begin{align*} &\qquad - \frac{1}{3}\operatorname{tr} {\mathbf E}_1^3 + \operatorname{tr}{\mathbf E}_1 \operatorname{tr}{\mathbf E}_2 - \operatorname{tr}{\mathbf E}_3.\end{align*}$$

Since the ${\mathbf K}_{j}$ are trace class, those trace expressions can be recast in terms of resolvent kernels and integral traces (cf. [74, Thm. (3.9)]).

Taking the bound $0\leqslant F(t)\leqslant 1$ of the Tracy–Widom distribution (being a probability distribution) into account, the results of Section 2 can be summarized in form of the following:

Theorem 2.1. Let $K_{(h)}(x,y)$ be a continuous kernel, $K_0$ the Airy kernel (1.6), and let the $K_j(x,y)$ be finite sums of rank one kernels with factors of the form (2.2). If, for some fixed non-negative integer m and some real number $t_0$ , there is a kernel expansion of the form

$$\begin{align*}K_{(h)}(x,y) = K_0(x,y) + \sum_{j=1}^m h^j K_j(x,y) + h^{m+1} \cdot O\big(e^{-(x+y)}\big), \end{align*}$$

which, for some constant $c>0$ , holds uniformly in $t_0\leqslant x,y < ch^{-1}$ as $h\to 0^+$ and which can be repeatedly differentiated w.r.t. x and y as uniform expansions, then the Fredholm determinant of $K_{(h)}$ on $(t,ch^{-1})$ satisfies

(2.11) $$ \begin{align} \det\big(I-K_{(h)}\big)|_{L^2(t,ch^{-1})} = F(t) + \sum_{j=1}^m G_j(t) h^j + h^{m+1} O(e^{-2t}) + e ^{-c h^{-1}} O(e^{-t}), \end{align} $$

uniformly for $t_0\leqslant t < ch^{-1}$ as $h\to 0^+$ . Here, F denotes the Tracy–Widom distribution (1.5) and the $G_j(t)$ are smooth functions depending on the kernels $K_1,\ldots ,K_j$ , satisfying the (right) tail bounds $G_j(t) = O(e^{-2t})$ . The first two are

$$ \begin{align*} G_1(t) &= - F(t) \cdot \left.\operatorname{tr}\big((I-K_0)^{-1}K_1\big)\right|{}_{L^2(t,\infty)}, \\ G_2(t) &= F(t)\cdot \bigg(\frac{1}{2} \Big(\!\left.\operatorname{tr}\big((I-K_0)^{-1}K_1\big)\right|{}_{L^2(t,\infty)}\Big)^2\\ &\qquad- \frac{1}{2}\left.\operatorname{tr}\big(((I-{K_0})^{-1}{K_1})^2\big)\right|{}_{L^2(t,\infty)} - \left.\operatorname{tr}\big(({I}-{K_0})^{-1}{K_2}\big)\right|{}_{L^2(t,\infty)}\bigg), \end{align*} $$

where $(I-K_0)^{-1}$ is understood as a resolvent kernel and the traces as integrals. The determinantal expansion (2.11) can repeatedly be differentiated w.r.t. t, preserving uniformity.

3. Expansion of the hard-to-soft edge transition

In this section, we prove an expansion for the hard-to-soft edge transition limit (1.7). To avoid notational clutter, we use the quantity

(3.1) $$ \begin{align} h_\nu := 2^{-1/3} \nu^{-2/3} \end{align} $$

and study expansions in powers of $h_\nu $ as $h_\nu \to 0^+$ . The transform $s=\phi _\nu (t)$ used in the transition limit can briefly be written as

(3.2) $$ \begin{align} \phi_\nu(t) = \nu^2(1-h_\nu t)^2. \end{align} $$

Theorem 3.1. There holds the hard-to-soft edge transition expansion

(3.3) $$ \begin{align} E_2^{\text{hard}}(\phi_\nu(t);\nu) = F(t) + \sum_{j=1}^m F_{j}(t) h_\nu^j + h_\nu^{m+1}\cdot O(e^{-3t/2}), \end{align} $$

which is uniformly valid when $t_0\leqslant t < h_\nu ^{-1}$ as $h_\nu \to 0^+$ , m being any fixed non-negative integer and $t_0$ any fixed real number. Preserving uniformity, the expansion can be repeatedly differentiated w.r.t. the variable t. Here, the $F_{j}$ are certain smooth functions starting with

(3.4a) $$ \begin{align} F_{1}(t) &= \frac{3t^2}{10} F'(t) - \frac{1}{5} F"(t), \end{align} $$

(3.4b) $$ \begin{align} F_{2}(t) &=\Big(\frac{2}{175} + \frac{32t^3}{175}\Big) F'(t) + \Big(-\frac{16t}{175} + \frac{9t^4}{200}\Big) F"(t) - \frac{3t^2}{50} F"'(t) + \frac{1}{50} F^{(4)}(t), \end{align} $$

(3.4c) $$ \begin{align} F_3(t) &= \Big(\frac{64 t}{7875} + \frac{1037 t^4}{7875}\Big) F'(t) +\Big( -\frac{9t^2}{175} + \frac{48t^5}{875} \Big) F"(t) \end{align} $$

$$\begin{align*} &\qquad +\Big(-\frac{122}{7875} -\frac{8 t^3}{125} + \frac{9 t^6}{2000}\Big) F"'(t) +\Big(\frac{16t}{875} -\frac{9 t^4}{1000} \Big)F^{(4)}(t)\notag\\ &\qquad +\frac{3 t^2 }{500}F^{(5)}(t) -\frac{1}{750} F^{(6)}(t).\end{align*}$$

It is rewarding to validate intriguing formulae such as (3.4a–c) by numerical methods: Figure 1 plots the functions $F_{1}(t)$ , $F_{2}(t)$ , $F_3(t)$ next to the approximation

(3.5) $$ \begin{align} F_{3}(t) \approx h_\nu^{-3} \cdot \big( E_2^{\text{hard}}(\phi_\nu(t);\nu) - F(t) - F_{1}(t)h_\nu - F_{2}(t) h_\nu^2\big) \end{align} $$

for $\nu =100$ and $\nu =800$ : the close matching with $F_3(t)$ as displayed by the latter is a very strong testament of the correctness of (3.4a–c) (in fact, some slips in preliminary calculations have been caught looking at plots which exhibited mismatches).

Figure 1 Plots of $F_{1}(t)$ (left panel) and $F_{2}(t)$ (middle panel) as in (3.4a/b). The right panel shows $F_{3}(t)$ as in (3.4c) (black solid line) with the approximations (3.5) for $\nu =100$ (red dotted line) and $\nu =800$ (green dashed line): the close agreement validates the functional forms displayed in (3.4). Details about the numerical method can be found in [14, 15, 19, 20].

The proof of Theorem 3.1 is split into several steps and will be concluded in Sections 3.2 and 3.3.

3.1. Kernel expansions

We start with an auxiliary result.

Lemma 3.2. Define for $h>0$ and $x,y<h^{-1}$ the function

(3.6) $$ \begin{align} \Phi(x,y;h):=\frac{\sqrt{(1- hx)(1-hy)}} {1-h(x+y)/2}. \end{align} $$

This function $\Phi $ satisfies the bound

(3.7) $$ \begin{align} 0 < \Phi(x,y;h) \leqslant 1 \end{align} $$

and has the convergent power series expansion

(3.8a) $$ \begin{align} \Phi(x,y;h) = 1 - (x-y)^2 \sum_{k=2}^\infty r_k(x,y) h^k, \end{align} $$

where the $r_k(x,y)$ are certain homogeneous symmetric rational ¹⁹ polynomials of degree $k-2$ , the first few of them being

(3.8b) $$ \begin{align} r_2(x,y) = \frac{1}{8},\quad r_3(x,y) = \frac{1}{8}(x+y), \quad r_4(x,y) = \frac{1}{128}\big(13(x^2+y^2)+22xy\big). \end{align} $$

The series converges uniformly for $x,y<(1-\delta )h^{-1}$ , $\delta $ being any fixed real positive number.

Proof. The bound (3.7) is the inequality of arithmetic and geometric means for the two positive real quantities $1-hx$ and $1-hy$ . By using

$$\begin{align*}\lim_{y\to x} \frac{1}{(x-y)^2}\Big(\Phi(x,y;h) -1\Big) = -\frac{1}{8}\left(\frac{h}{1-hx}\right)^2, \end{align*}$$

the analyticity of $\Phi (x,y;h)$ w.r.t. h, and the scaling law

$$\begin{align*}\Phi(\lambda^{-1}x,\lambda^{-1}y;\lambda h) = \Phi(x,y;h) \quad (\lambda>0), \end{align*}$$

we deduce the claims about the form and uniformity of the power series expansion (3.8).

Because of the representation (1.2) of $E_2^{\text {hard}}(s;\nu )$ in terms of a Fredholm determinant of the Bessel kernel (1.3), we have to expand the induced transformation of that kernel.

Lemma 3.3. The change of variables $s=\phi _\nu (t)$ , mapping $t< h_\nu ^{-1}$ monotonically decreasing to $s> 0$ , induces the symmetrically transformed Bessel kernel

(3.9a) $$ \begin{align} \hat K_\nu^{\mathrm{Bessel}}(x,y) &:=\sqrt{\phi_\nu^{\prime}(x)\phi_\nu^{\prime}(y)} \, K_\nu^{\mathrm{Bessel}}(\phi_\nu(x),\phi_\nu(y)). \end{align} $$

There holds the kernel expansion

(3.9b) $$ \begin{align} \hat K_\nu^{\mathrm{Bessel}}(x,y) &= K_0(x,y) + \sum_{j=1}^m K_j(x,y) h_\nu^j + h_\nu^{m+1}\cdot O\big(e^{-(x+y)}\big), \end{align} $$

which is uniformly valid when $t_0 \leqslant x,y < h_\nu ^{-1}$ as $h_\nu \to 0^+$ , m being any fixed non-negative integer and $t_0$ any fixed real number. Here, the $K_j$ are certain finite rank kernels of the form

$$\begin{align*}K_j(x,y) = \sum_{\kappa,\lambda\in\{0,1\}} p_{j,\kappa\lambda}(x,y) \operatorname{Ai}^{(\kappa)}(x)\operatorname{Ai}^{(\lambda)}(y), \end{align*}$$

where $p_{j,\kappa \lambda }(x,y)$ are rational polynomials; the first two kernels are

(3.10a) $$ \begin{align} K_1(x,y) &= \frac{1}{10} \Big( -3\big(x^2+xy+y^2\big)\operatorname{Ai}(x)\operatorname{Ai}(y)\\&\qquad + 2 \big(\!\operatorname{Ai}(x)\operatorname{Ai}'(y) + \operatorname{Ai}'(x)\operatorname{Ai}(y)\big)+ 3(x+y)\operatorname{Ai}'(x)\operatorname{Ai}'(y) \Big),\notag \end{align} $$

(3.10b) $$ \begin{align}K_2(x,y) &= \frac{1}{1400}\Big(\big(-235 \left(x^3+y^3\right)-319 x y (x+y)+56\big) \operatorname{Ai}(x)\operatorname{Ai}(y)\\&\qquad +\big(63(x^4+x^3 y-x^2 y^2-x y^3-y^4)-55 x+239 y\big) \operatorname{Ai}(x)\operatorname{Ai}'(y) \notag\\ &\qquad\quad +\big(63(-x^4-x^3 y-x^2 y^2+x y^3+y^4)+239 x-55 y\big) \operatorname{Ai}'(x)\operatorname{Ai}(y) \notag\\ &\qquad\qquad +\big(340(x^2+y^2)+256 x y\big) \operatorname{Ai}'(x)\operatorname{Ai}'(y)\Big). \notag \end{align} $$

Preserving uniformity, the kernel expansion (3.9) can repeatedly be differentiated w.r.t. x, y.

Proof. We have to prove, analytically, the claim about the domain of uniformity of the error of the expansion and, algebraically, the finite-rank structure of the expansion kernels $K_j$ . In our original proof, ²⁰ presented below, we use Olver’s expansion of Bessel functions of large order in the transition region (see Appendix A.3), and the finite-rank structure is obtained by explicitly inspecting (with Mathematica) a certain algebraic condition (see (3.14)) for the first instances $j=1,\ldots ,m_*$ – we choose to stop at $m_*=100$ . However, using the machinery of Riemann–Hilbert problem, this restriction was recently removed by Yao and Zhang [85] (their proof extends over 18 pages); we will comment on their work at the end.

The original proof, requiring $m\leqslant m_*$ . By using $\Phi (x,y;h)$ as defined in (3.6) and writing

$$\begin{align*}\begin{aligned} \phi_\nu(t) = \omega_\nu(t)^2,\quad \omega_\nu(t)=\nu(1-h_\nu t),\\ a_\nu(x,y)= (1-h_\nu y) \cdot \frac{1}{\sqrt{2h_\nu}}J_\nu\big(\omega_\nu(x)\big)\cdot \frac{d}{dy} \frac{1}{\sqrt{2h_\nu}}J_\nu\big(\omega_\nu(y)\big), \end{aligned} \end{align*}$$

we can factor the transformed Bessel kernel in the simple form

$$\begin{align*}\hat K_\nu^{\text{Bessel}}(x,y) = \Phi(x,y;h_\nu) \cdot \frac{a_\nu(x,y)-a_\nu(y,x)}{x-y}, \end{align*}$$

noting, by symmetry, the removability of the singularities at $x=y$ of the second factor.

First, if x or y is between $\tfrac 34 \cdot h_\nu ^{-1}$ and $h_\nu ^{-1}$ , using the bound $0<\Phi \leqslant 1$ (see (3.7)), one can argue as in the proof of Lemma A.9: since at least one of the Bessel factors is of the form $J_\nu ^{(\kappa )}(z)$ with $0\leqslant z \leqslant 1/4$ , which plainly falls into the superexponentially decaying region as $\nu \to \infty $ , and since at least one of the Airy factors of each term of the expansion is also superexponentially decaying as $\nu \to \infty $ , the transformed Bessel kernel and the expansion terms in (3.9) get completely absorbed into the error term (bounding the other factors as in Section 2.1)

$$\begin{align*}h_\nu^{m+1}\cdot O\big(e^{-(x+y)}\big). \end{align*}$$

Here, the removable singularities at $x=y$ are dealt with by using the differentiability of the corresponding bounds (or by extending to the complex domain and using Cauchy’s integral formula as in the proof of [21, Prop. 8]).

Therefore, we may suppose from now on that $t_0 \leqslant x,y \leqslant \tfrac 34 \cdot h_\nu ^{-1}$ . By Lemma 3.2, in this range of x and y, the power series expansion

(3.11) $$ \begin{align} \Phi(x,y;h_\nu) = 1 - (x-y)^2 \sum_{k=2}^\infty r_k(x,y) h_\nu^k \end{align} $$

converges uniformly. Here, the $r_k(x,y)$ are certain homogeneous symmetric rational polynomials of degree $k-2$ , the first of them being $r_2(x,y)=1/8$ .

Next, we rewrite the uniform version of the large order expansion of Bessel functions in the transition region, as given in Lemma A.9, in the form

(3.12) $$ \begin{align} \frac{1}{\sqrt{2h_\nu}} J_\nu\big(\omega_\nu(t)\big) =\big(1 + h_\nu p_m(t;h_\nu)\big) \operatorname{Ai}(t) + h_\nu q_m(t;h_\nu) \operatorname{Ai}'(t) + h_\nu^{m+1} O(e^{-t}), \end{align} $$

where the estimate of the remainder is uniform for $t_0 \leqslant t < h_\nu ^{-1}$ as $h_\nu \to 0^+$ and

$$\begin{align*}p_m(t;h) = 2^{1/3}\sum_{k=0}^{m-1} A_{k+1}(-t/2^{1/3}) \,(2^{1/3}h)^{k}, \quad q_m(t;h) = 2^{2/3}\sum_{k=0}^{m-1} B_{k+1}(-t/2^{1/3}) \,(2^{1/3}h)^{k}, \end{align*}$$

with the polynomials $A_k(\tau )$ and $B_k(\tau )$ from (A.16). It follows from Remark A.8 that $p_m(t;h)$ and $q_m(t;h)$ are rational polynomials in t and h, starting with

$$\begin{align*}p_2(t;h) = \frac{2}{10} t + \frac{h}{1400}(63t^5+120t^2),\quad q_2(t;h) = \frac{3}{10} t^2 + \frac{h}{1400}(340t^3+40). \end{align*}$$

Also given in Lemma A.9, under the same conditions, the expansion (3.12) can be repeatedly differentiated w.r.t t while preserving uniformity. From this, we obtain, using the Airy differential equation $\operatorname {Ai}"(\xi ) = \xi \operatorname {Ai}(\xi )$ , that uniformly (given the range x and y) ²¹

$$\begin{align*}a_\nu(x,y) = \sum_{\kappa,\lambda\in\{0,1\}} p_{\kappa\lambda}^m(x,y;h_\nu) \operatorname{Ai}^{(\kappa)}(x)\operatorname{Ai}^{(\lambda)}(y) + h_\nu^{m+1} O(e^{-(x+y)}), \end{align*}$$

where

$$ \begin{align*} p^m_{00}(x,y;h) &= h(1- hy)(1 + hp_m(x;h))\big(yq_m(y;h)+\partial_y p_m(y;h)\big), \\ p^m_{01}(x,y;h) &= (1- hy)(1 + hp_m(x;h))\big(1+h(p_m(y;h)+\partial_y q_m(y;h))\big),\\ p^m_{10}(x,y;h) &= h^2(1- hy)q_m(x;h)\big(y q_m(y;h)+\partial_y p_m(y;h)\big),\\ p^m_{11}(x,y;h) &= h(1- hy)q_m(x;h)\big(1+h(p_m(y;h)+\partial_y q_m(y;h))\big) \end{align*} $$

are rational polynomials in x, y and h. In particular, those factorizations show

$$\begin{align*}\begin{aligned} p^m_{00}(x,y;h) = O(h), \qquad p^m_{11}(x,y;h) = O(h),\\ p^m_{01}(x,y;h) = 1 + O(h), \qquad p^m_{10}(x,y;h) = O(h^2). \end{aligned} \end{align*}$$

If we denote by $\hat p^m_{\kappa \lambda }(x,y;h)$ the polynomials obtained from $p^m_{\kappa \lambda }(x,y;h)$ after dropping all powers of h that have an exponent larger than m (thus contributing terms to the expansion that get absorbed in the error term), we obtain

$$\begin{align*}a_\nu(x,y) = \operatorname{Ai}(x)\operatorname{Ai}'(y) + \sum_{\kappa,\lambda\in\{0,1\}} \hat p^m_{\kappa\lambda}(x,y;h_\nu) \operatorname{Ai}^{(\kappa)}(x)\operatorname{Ai}^{(\lambda)}(y) + h_\nu^{m+1} O(e^{-(x+y)}), \end{align*}$$

with a polynomial expansion

$$\begin{align*}\hat p^m_{\kappa\lambda}(x,y;h) = \sum_{j=1}^m \hat p_{j,\kappa\lambda}(x,y) h^j, \end{align*}$$

whose coefficient polynomials $\hat p_{j,\kappa \lambda }(x,y)$ , being the unique expansion coefficients as $h_\nu \to 0^+$ , are now independent of m. Hence, the anti-symmetrization of $a_\nu (x,y)$ satisfies the uniform expansion (given the range of x and y)

(3.13) $$ \begin{align} a_\nu(x,y)-a_\nu(y,x) = \operatorname{Ai}(x)\operatorname{Ai}'&(y)- \operatorname{Ai}'(x)\operatorname{Ai}(y) \notag \\ &+ \sum_{\kappa,\lambda\in\{0,1\}} \hat q^m_{\kappa\lambda}(x,y;h_\nu) \operatorname{Ai}^{(\kappa)}(x)\operatorname{Ai}^{(\lambda)}(y) + h_\nu^{m+1} O(e^{-(x+y)}) \end{align} $$

with the polynomial expansion

$$\begin{align*}\begin{aligned} \hat q^m_{\kappa\lambda}(x,y;h) &= \sum_{j=1}^m \hat q_{j,\kappa\lambda}(x,y) h^j,\\ \hat q_{j,00}(x,y) = \hat p_{j,00}(x,y) - \hat p_{j,00}(y,x),&\qquad \hat q_{j,11}(x,y) = \hat p_{j,11}(x,y) - \hat p_{j,11}(y,x)\\ \hat q_{j,01}(x,y) = \hat p_{j,01}(x,y) - \hat p_{j,10}(y,x),&\qquad \hat q_{j,10}(x,y) = - \hat q_{j,01}(y,x). \end{aligned} \end{align*}$$

Since the rational polynomials $\hat q_{j,00}(x,y)$ and $\hat q_{j,11}(x,y)$ are anti-symmetric in x, y, they factor in the form

(3.14) $$ \begin{align} (x-y) \times (\text{symmetric rational polynomial in } x \text{ and } y); \end{align} $$

the first few cases are

$$ \begin{align*} \hat q_{1,00}(x,y) &= -\frac{3}{10}(x-y)(x^2+xy+y^2),\\ \hat q_{1,11}(x,y) &= \frac{3}{10}(x-y)(x+y),\\ \hat q_{2,00}(x,y) &= \frac{1}{1400} (x-y)\Big(-235\big(x^3+y^3\big)-319xy(x+y)+56\Big),\\ \hat q_{2,11}(x,y) &= \frac{1}{1400} (x-y)\Big(340\big(x^2+y^2\big)+256xy\Big). \end{align*} $$

Even though there is no straightforward structural reason for the rational polynomials $\hat q_{j,01}$ (and thus $\hat q_{j,10}$ ) to be divisible by $x-y$ as well, an inspection ²² of the first cases reveals this to be true for at least $j=1,\ldots ,m_*$ , the first two of them being

$$ \begin{align*} \hat q_{1,01}(x,y) &= \frac{2}{10}(x-y),\\ \hat q_{2,01}(x,y) &= \frac{1}{1400}(x-y) \Big(63\big(x^4+x^3 y-x^2 y^2-x y^3-y^4\big) +120 x + 64y\Big). \end{align*} $$

Now, by restricting ourselves to the explicitly checked cases $m\leqslant m_*$ , we denote by $q^m_{\kappa \lambda }(x,y;h)$ the polynomials obtained from $\hat q^m_{\kappa \lambda }(x,y;h)$ after division by the factor $x-y$ . Since (3.13) is an expansion of an anti-symmetric function with anti-symmetric remainder which can repeatedly be differentiated, division by $x-y$ yields removable singularities at $x=y$ and does not change the character of the expansion (see also the argument given in the proof of [21, Prop. 8]):

$$\begin{align*}\frac{a_\nu(x,y)-a_\nu(y,x)}{x-y} = K_0(x,y)+ \sum_{\kappa,\lambda\in\{0,1\}} q^m_{\kappa\lambda}(x,y;h_\nu) \operatorname{Ai}^{(\kappa)}(x)\operatorname{Ai}^{(\lambda)}(y) + h_\nu^{m+1} O\big(e^{-(x+y)}\big). \end{align*}$$

The lemma now follows by multiplying this expansion with (3.11), noting that the terms

$$\begin{align*}-r_k(x,y) (x-y)^2 K_0(x,y) = -r_k(x,y) (x-y) \big(\!\operatorname{Ai}(x)\operatorname{Ai}'(y) - \operatorname{Ai}'(x)\operatorname{Ai}(y) \big) \end{align*}$$

also take the form asserted for the terms in the kernels $K_j$ ( $j\geqslant 1$ ).

Finally, since all the expansions can repeatedly be differentiated under the same conditions while preserving their uniformity, the same holds for the resulting expansion of the kernel.

Comments on the unconditional proof of Yao and Zhang [85]. Instead of using expansions of the Bessel functions of large order in the transition region, Yao and Zhang address expanding the integrable kernel

$$\begin{align*}\frac{a_\nu(x,y)-a_\nu(y,x)}{x-y} \end{align*}$$

directly by the machinery of Riemann–Hilbert problems (see [27] for a general discussion of kernels of that form and their induced integral operators). The advantage of such a direct approach is a better understanding of the polynomial coefficients in (3.13), and the divisibility by $x-y$ follows from an interesting algebraic structure of the Airy functions; namely, by the Airy differential equation, there are rational polynomials $\alpha _k, \beta _k \in {\mathbb Q}[\xi ]$ such that

$$\begin{align*}\operatorname{Ai}^{(k)}(\xi) = \alpha_k(\xi) \operatorname{Ai}(\xi) + \beta_k(\xi) \operatorname{Ai}'(\xi), \end{align*}$$

and the divisibility (3.14) turns out to be equivalent [85, pp. 18–20] to the relation [85, Lemma 4.1]

$$\begin{align*}\sum_{j=0}^m \binom{m}{j} \begin{vmatrix} \alpha_j & \alpha_{m+1-j} \\ \beta_j & \beta_{m+1-j} \end{vmatrix} = 0 \qquad (m\geqslant 1), \end{align*}$$

which can be proved by induction.

Remark 3.4. The case $m=0$ of Lemma 3.3,

$$\begin{align*}\hat K_\nu^{\text{Bessel}}(x,y) = K_0(x,y) + h_\nu\cdot O\big(e^{-(x+y)}\big), \end{align*}$$

is the $\beta =2$ case of [21, Eq. (4.8)] in the work of Borodin and Forrester. There, in [21, Prop. 8], it is stated that this expansion would be uniformly valid for $x,y \geqslant t_0$ . However, stated in such a generality, it is not correct (see Footnote A.10) and, in fact, similar to our proof given above, their proof is restricted to the range $t_0\leqslant x,y < h_\nu ^{-1}$ , which completely suffices to address the hard-to-soft edge transition. (See Remark A.10 for yet another issue with [21, Prop. 8].)

To reduce the complexity of calculating the functional form of the first three finite-size correction terms in the hard-to-soft edge transition (3.3), we consider a second kernel transform.

Lemma 3.5. For $h>0$ , the Airy kernel $K_0$ and the first expansion kernels $K_1$ , $K_2$ , $K_3$ from Lemma 3.3 we consider

$$\begin{align*}K_{(h)}(x,y) := K_0(x,y) + K_1(x,y) h + K_2(x,y)h^2 + K_3(x,y)h^3 \end{align*}$$

and the transformation, ²³ where $\zeta (z)$ is defined as in Section A.3,

(3.15) $$ \begin{align} s = \psi_h^{-1}(t):=2^{-1/3} h^{-1} \zeta(1-h t). \end{align} $$

Then $t=\psi _h(s)$ maps $s \in {\mathbb R}$ monotonically increasing to $-\infty < t < h^{-1}$ , with $t\leqslant \mu h^{-1}$ , $\mu = 0.94884\cdots $ , when $s\leqslant 2h^{-1}$ , and induces the symmetrically transformed kernel

(3.16a) $$ \begin{align} \tilde K_{(h)}(x,y) &:= \sqrt{\psi_h^{\prime}(x)\psi_h^{\prime}(y)} \, K_{(h)}(\psi_h(x),\psi_h(y)) \end{align} $$

which expands as

(3.16b) $$ \begin{align} \tilde K_{(h)}(x,y) &= K_0(x,y) + \tilde K_1(x,y)h + \tilde K_2(x,y)h^2 + \tilde K_3(x,y)h^3 + h^4 \cdot O\big(e^{-(x+y)}\big), \end{align} $$

uniformly in $s_0 \leqslant x,y \leqslant 2h^{-1}$ as $h\to 0^+$ , $s_0$ being a fixed real number. The three kernels are

(3.16c) $$ \begin{align} \tilde K_1 &= \frac{1}{5}(\operatorname{Ai}\otimes \operatorname{Ai}' + \operatorname{Ai}'\otimes \operatorname{Ai}), \end{align} $$

(3.16d) $$ \begin{align} \tilde K_2 &= -\frac{48}{175}\operatorname{Ai}\otimes\operatorname{Ai} + \frac{11}{70}(\operatorname{Ai}\otimes\operatorname{Ai}"'+\operatorname{Ai}"'\otimes\operatorname{Ai}) - \frac{51}{350}(\operatorname{Ai}'\otimes\operatorname{Ai}"+\operatorname{Ai}"\otimes\operatorname{Ai}'), \end{align} $$

(3.16e) $$ \begin{align} \tilde K_3 &= -\frac{176}{1125} (\operatorname{Ai}\otimes \operatorname{Ai}" + \operatorname{Ai}"\otimes \operatorname{Ai}) +\frac{13}{450} (\operatorname{Ai}\otimes \operatorname{Ai}^{\mathrm V} + \operatorname{Ai}^{\mathrm V}\otimes \operatorname{Ai}) + \frac{3728}{7875}\operatorname{Ai}'\otimes \operatorname{Ai}' \end{align} $$

$$\begin{align*} & \qquad -\frac{583}{5250} (\operatorname{Ai}'\otimes \operatorname{Ai}^{\mathrm IV} + \operatorname{Ai}^{\mathrm IV}\otimes \operatorname{Ai}') + \frac{13}{225} (\operatorname{Ai}"\otimes \operatorname{Ai}"' + \operatorname{Ai}"'\otimes \operatorname{Ai}").\notag \end{align*}$$

Preserving uniformity, the kernel expansion can repeatedly be differentiated w.r.t. x, y.

Proof. Reversing the power series (A.19) gives

$$\begin{align*}t = \psi_h(s) = s - \frac{3s^2}{10} h - \frac{s^3}{350} h^2 +\frac{479 s^4}{63000}h^3 + \cdots, \end{align*}$$

which is uniformly convergent for $s_0 \leqslant s \leqslant 2h^{-1}$ since $|ht| \leqslant \mu < 1$ . Taking the expressions for $K_1$ , $K_2$ given in (3.10), and for $K_3$ from the supplementary Mathematica notebook referred to in Footnote 13, a routine calculation with truncated power series gives formula (3.16c) for $\tilde K_1$ and

$$\begin{align*}\begin{aligned} \tilde K_2(x,y) = \frac{1}{350} \Big(14\operatorname{Ai}(x)\operatorname{Ai}(y) + (-51 x + 55 y)\operatorname{Ai}(x)\operatorname{Ai}'(y) + (55x -51 y) \operatorname{Ai}'(x)\operatorname{Ai}(y) \Big),\\ \tilde K_3(x,y) = \frac{1}{15750}\Big(266(x+y)\operatorname{Ai}(x)\operatorname{Ai}(y) + (-1749 x^2+910 x y+455 y^2) \operatorname{Ai}(x)\operatorname{Ai}'(y)\\ \qquad + (455 x^2+910 x y-1749 y^2) \operatorname{Ai}'(x)\operatorname{Ai}(y) + 460 \operatorname{Ai}'(x)\operatorname{Ai}'(y)\Big). \end{aligned} \end{align*}$$

The Airy differential equation $\xi \operatorname {Ai}(\xi ) = \operatorname {Ai}"(\xi )$ implies the replacement rule

(3.17) $$ \begin{align} \xi^j \operatorname{Ai}^{(k)}(\xi) = \xi^{j-1} \operatorname{Ai}^{(k+2)}(\xi) - k \xi^{j-1} \operatorname{Ai}^{(k-1)}(\xi) \qquad (j\geqslant 1, \; k\geqslant 0) \end{align} $$

which, if repeatedly applied to a kernel of the given structure, allows us to absorb any powers of x and y into higher order derivatives of $\operatorname {Ai}$ . This process yields the asserted form of $\tilde K_2$ and $\tilde K_3$ , which will be the preferred form in course of the calculations in Section 3.3.

Since we stay within the range of uniformity of the power series expansions and calculations with truncated powers series are amenable to repeated differentiation, the result now follows from the bounds given in Section 2.1.

3.2. Proof of the general form of the expansion

By Lemma 3.3 and Theorem 2.1, we get (the Fredholm determinants are seen to be equal by transforming the integrals)

$$ \begin{align*} E_2^{\text{hard}}(\phi_\nu(t);\nu) = \det(I - K_\nu^{\text{Bessel}})|_{L^2(0,\phi_\nu(t))} &= \det(I - \hat K_\nu^{\text{Bessel}})|_{L^2(t,h_\nu^{-1})} \\ &= F(t) + \sum_{j=1}^m F_{j}(t) h_\nu^j + h_\nu^{m+1} O(e^{-2t}) + e^{-h_\nu^{-1}} O(e^{-t}), \end{align*} $$

uniformly for $t_0\leqslant t < h_\nu ^{-1}$ as $h_\nu \to 0^+$ ; preserving uniformity, this expansion can be repeatedly differentiated w.r.t. the variable t. By Theorem 2.1, the $F_{j}(t)$ are certain smooth functions that can be expressed in terms of traces of integral operators of the form given in Theorem 2.1. Observing

$$\begin{align*}e^{-h_\nu^{-1}} < e^{-h_\nu^{-1}/2} e^{-t/2} = h_\nu^{m+1} O(e^{-t/2})\qquad (h_\nu\to0^+), \end{align*}$$

we can combine the two error terms as $h_\nu ^{m+1} O(e^{-3t/2})$ . This finishes the proof of (3.3).

3.3. Functional form of $F_{1}(t)$ , $F_{2}(t)$ and $F_3(t)$

Instead of calculating $F_{1}$ , $F_{2}$ , $F_3$ directly from the formulae in Theorem 2.1 applied to the kernels $K_1$ , $K_2$ in (3.10) (and to the unwieldy expression for $K_3$ obtained in the supplementary Mathematica notebook referred to in Footnote 13), we will reduce them to the corresponding functions $\tilde F_1, \tilde F_2, \tilde F_3$ induced by the much simpler kernels $\tilde K_1$ , $\tilde K_2$ , $\tilde K_3$ in (3.16).

3.3.1. Functional form of $\tilde F_{1}(t)$ and $\tilde F_{2}(t)$

Upon writing

$$\begin{align*}u_{jk}(t) = \operatorname{tr}\big((I-K_0)^{-1} \operatorname{Ai}^{(j)} \otimes \operatorname{Ai}^{(k)} \big)\big|_{L^2(t,\infty)} \end{align*}$$

and observing (the symmetry of the resolvent kernel implies the symmetry $u_{jk}(t)= u_{kj}(t)$ )

$$ \begin{align*} \operatorname{tr}\big((I-K_0)^{-1} \tilde K_1\big)\big|_{L^2(t,\infty)} &= \frac{2}{5} u_{10}(t),\\ \operatorname{tr}\big(((I-K_0)^{-1} \tilde K_1)^2\big)\big|_{L^2(t,\infty)} &= \frac{2}{25} \big(u_{00}(t)u_{11}(t) + u_{10}(t)^2\big),\\ \operatorname{tr}\big((I-K_0)^{-1} \tilde K_2\big)\big|_{L^2(t,\infty)} &= \frac{1}{175}\big(-48u_{00}(t)-51u_{21}(t) +55u_{30}(t)\big), \end{align*} $$

the formulae of Theorem 2.1 applied to $\tilde K_1$ and $\tilde K_2$ give

(3.18a) $$ \begin{align} \tilde F_{1}(t) &= -\frac{2}{5} F(t) u_{10}(t), \end{align} $$

(3.18b) $$ \begin{align} \tilde F_{2}(t) &= F(t)\left(\frac{48}{175} u_{00}(t) - \frac{1}{25} \begin{vmatrix} u_{00}(t) & u_{01}(t) \\ u_{10}(t) & u_{11}(t) \end{vmatrix} +\frac{51}{175} u_{21}(t) - \frac{11}{35} u_{30}(t) \right). \end{align} $$

We recall from [19, Remark 3.1] that the simple recursion

$$\begin{align*}u_{jk}'(t) = u_{j+1,k}(t) + u_{j,k+1}(t) - u_{j0}(t) u_{k0}(t) \end{align*}$$

yields similar formulae for the first few derivatives of the distribution $F(t)$ :

(3.19) $$ \begin{align} \begin{aligned} F'(t) = F(t)\cdot &u_{00}(t),\quad F"(t) = 2F(t)\cdot u_{10}(t),\quad F"'(t) = 2F(t)\cdot \big(u_{11}(t)+u_{20}(t)\big),\\ &F^{(4)}(t) = 2F(t)\cdot\left( \begin{vmatrix} u_{00}(t) & u_{01}(t) \\ u_{10}(t) & u_{11}(t) \end{vmatrix} + 3u_{21}(t) + u_{30}(t) \right). \end{aligned} \end{align} $$

By a linear elimination of the terms

$$\begin{align*}u_{00}(t), \quad u_{10}(t),\quad \begin{vmatrix} u_{00}(t) & u_{01}(t) \\ u_{10}(t) & u_{11}(t) \end{vmatrix}, \end{align*}$$

we obtain, as an intermediate step,

(3.20a) $$ \begin{align} \tilde F_{1}(t) = -\frac{1}{5} F"(t),\quad \tilde F_{2}(t) = \frac{48}{175} F'(t) - \frac{1}{50} F^{(4)}(t) + \frac{24}{175}F(t)\big(3u_{21}(t) - 2 u_{30}(t)\big). \end{align} $$

To simplify even further, we have to refer to the full power of the general Tracy–Widom theory (i.e., representing F in terms of Painlevé II): by advancing its set of formulae, Shinault and Tracy [73, p. 68] showed, through an explicit inspection of each single case, that the functions $F(t)\cdot u_{jk}(t)$ in the range $0 \leqslant j+k\leqslant 8$ are linear combinations of the form

(3.20b) $$ \begin{align} p_1(t) F'(t) + p_2(t) F"(t) + \cdots + p_{j+k+1}(t) F^{(j+k+1)}(t) \end{align} $$

with rational polynomials $p_\kappa (t)$ (depending, of course, on j, k). They conjectured this structure to be true for all j, k. In particular, their table [73, p. 68] has the entries

$$\begin{align*}F(t) \cdot u_{21}(t) = -\frac{1}{4}F'(t) +\frac{1}{8}F^{(4)}(t), \quad F(t)\cdot u_{30}(t) = \frac{7}{12} F'(t) + \frac{t}{3} F"(t) + \frac{1}{24} F^{(4)}(t). \end{align*}$$

This way, we get, rather unexpectedly, the simple and short form ²⁴

(3.20c) $$ \begin{align} \tilde F_2(t) = \frac{2}{175} F'(t) - \frac{16t}{175}F"(t) + \frac{1}{50} F^{(4)}(t). \end{align} $$

3.3.2. Functional form of $\tilde F_{3}(t)$

For the $\tilde F_j(t)$ with $j\geqslant 3$ , calculating such functional forms requires a more systematic, algorithmic approach. By ‘reverse engineering’ the remarks of Shinault and Tracy about validating their table [73, p. 68], we have actually found an algorithm to compile such a table; see Appendix B. This algorithm can also be applied to nonlinear rational polynomials of the terms $u_{jk}$ , resulting either in an expression of the desired form (3.20b) (with $j+k+1$ replaced by some n), or a message that such a form does not exist (for the given n).

Now, if we evaluate the expansion function $G_3(t)=F(t)d_3(t)$ of Theorem 2.1 by using (2.10c) and rewrite the traces in terms of the $u_{jk}$ , we obtain

(3.21a) $$ \begin{align} \tilde F_3(t) &= F(t) \left(-\frac{3728 }{7875}u_{11}(t) +\frac{352 }{1125}u_{20}(t) -\frac{51}{875} \begin{vmatrix} u_{10}(t) & u_{11}(t) \\ u_{20}(t) & u_{21}(t) \end{vmatrix} \right. \end{align} $$

$$\begin{align*} &\hspace{50pt} \left. -\frac{11}{175} \begin{vmatrix} u_{00}(t) & u_{01}(t) \\ u_{30}(t) & u_{31}(t) \end{vmatrix}-\frac{26}{225} u_{32}(t)+\frac{583 }{2625}u_{41}(t)-\frac{13}{225} u_{50}(t)\right)\notag, \end{align*}$$

which evaluates by the algorithm of Appendix B further to

(3.21b) $$ \begin{align} &= \frac{64 t }{7875}F'(t) -\frac{24 t^2 }{875}F"(t) -\frac{122}{7875}F"'(t) +\frac{16 t}{875} F^{(4)}(t) -\frac{1}{750} F^{(6)}(t). \end{align} $$

Alternatively, we arrive here by combining the tabulated expressions for $u_{11}$ , $u_{20}$ , $u_{32}$ , $u_{41}$ , $u_{50}$ as displayed in [73, p. 68] with those for both of the $2\times 2$ determinants in (B.8).

3.3.3. Lifting to the functional form of $F_{1}(t)$ , $F_{2}(t)$ , $F_{3}(t)$

The relation between $F_{1}(t)$ , $F_{2}(t)$ , $F_3(t)$ and their counterparts with a tilde is established by Lemma 3.5. By using the notation introduced there, with t being any fixed real number, the expansion parameter h sufficiently small and $s=\psi _{h}^{-1}(t)$ , Theorem 2.1 yields (the Fredholm determinants are seen to be equal by transforming the integrals)

(3.22) $$ \begin{align} \begin{aligned} & \det\big(I-K_{(h)}\big)|_{L^2(t, \,\mu h^{-1})} = F(t) + F_{1}(t) h + F_{2}(t)h^2 + F_{3}(t)h^3 + O(h^4) \\ = &\det\big(I-\tilde K_{(h)}\big)|_{L^2(s, \, 2 h^{-1})} = F(s) + \tilde F_{1}(s) h + \tilde F_{2}(s)h^2 + \tilde F_{3}(s)h^3 + O(h^4), \end{aligned} \end{align} $$

where we have absorbed the exponentially small contributions $e^{-\mu h^{-1}} O(e^{-t})$ and $e^{-2h^{-1}} O(e^{-s})$ into the $O(h^4)$ error term. Using the power series (A.19),

$$\begin{align*}s = 2^{-1/3} h^{-1} \zeta(1-h t) = t+ \frac{3t^2}{10}h + \frac{32t^3}{175} h^2 + \frac{1037t^4}{7875} h^3 + \cdots, \end{align*}$$

we get by Taylor expansion, for any smooth function $G(s)$ ,

$$ \begin{align*} G(s) &= G(t) + \frac{3t^2}{10} G'(t) h + \left(\frac{32t^3}{175} G'(t) + \frac{9t^4}{200} G"(t) \right) h^2 \\ &\qquad +\left(\frac{1037 t^4}{7875} G'(t)+\frac{48t^5}{875} G"(t)+ \frac{9 t^6 }{2000}G"'(t)\right)h^3 + O(h^4). \end{align*} $$

By plugging this into (3.22) and comparing coefficients, we obtain

$$ \begin{align*} F_{1}(t) &= \tilde F_{1}(t) + \frac{3t^2}{10} F'(t),\\ F_{2}(t) &=\tilde F_{2}(t) + \frac{3t^2}{10} \tilde F_{1}^{\prime}(t) + \frac{32t^3}{175} F'(t) + \frac{9t^4}{200} F"(t),\\ F_{3}(t) &= \tilde F_3(t) + \frac{32t^3}{175} \tilde F_1^{\prime}(t)+\frac{9t^4}{200} \tilde F_1^{\prime\prime}(t)+\frac{3t^2}{10} \tilde F_2^{\prime}(t) +\frac{1037 t^4 }{7875}F'(t)+\frac{48t^5 }{875} F"(t)+ \frac{9 t^6 F"'(t)}{2000}. \end{align*} $$

Combined with (3.20), this finishes the proof of (3.4).

3.4. Simplifying the form of Choup’s Edgeworth expansions

When, instead of the detour via $\tilde K_1$ , Theorem 2.1 is directly applied to the kernel $K_1$ in (3.10a), we get

$$\begin{align*}F_{1}(t) = - \frac{1}{5}F"(t) + \frac{3}{10} F(t) \operatorname{tr}\big((I-K_0)^{-1} L\big)\big|_{L^2(t,\infty)}, \end{align*}$$

where

(3.23a) $$ \begin{align} L(x,y) = (x^2+xy+y^2)\operatorname{Ai}(x)\operatorname{Ai}(y) - (x+y)\operatorname{Ai}'(x)\operatorname{Ai}'(y). \end{align} $$

Now, a comparison with (3.4a) proves the useful formula ²⁵

(3.23b) $$ \begin{align} F(t) \operatorname{tr}\big((I-K_0)^{-1} L\big)\big|_{L^2(t,\infty)} = t^2 F'(t). \end{align} $$

As an application to the existing literature, this formula helps us to simplify the results obtained by Choup for the soft-edge limit expansions of GUE and LUE – that is, when studying the distribution of the largest eigenvalue distribution function in $\text {GUE}_n$ and $\text {LUE}_{n,\nu }$ (dimension n, parameter $\nu $ ) as $n\to \infty $ . In fact, since the kernel L appears in the first finite-size correction term of a corresponding kernel expansion [24, Thm. 1.2/1.3], lifting that expansion to the Fredholm determinant by Theorem 2.1 allows us to recast [24, Thm. 1.4] in a simplified form; namely, denoting the maximum eigenvalues by $\lambda _{n}^{\text {G}}$ and $\lambda _{n,\nu }^{\text {L}}$ , we obtain, locally uniform in t as $n\to \infty $ ,

(3.24) $$ \begin{align} {\mathbb P}\Big( \lambda_{n}^{\text{G}} \leqslant \sqrt{2n} + t \cdot 2^{-1/2} n^{-1/6} \Big) &= F(t) + \frac{n^{-2/3}}{40}\big(2t^2 F'(t) - 3F"(t)\big) + O(n^{-1}),\end{align} $$

a result, which answers a question suggested by Baik and Jenkins [11, p. 4367].

4. Expansion of the Poissonized length distribution

The Poissonization of the length distribution requires the hard-to-soft edge transition of Theorem 3.1 to be applied to the probability distribution $E_2^{\text {hard}}(4r;\nu )$ (for integer $\nu =l$ , but we consider the case of general $\nu>0$ first). For large intensities r, the mode of this distribution is located in the range of those parameters $\nu $ for which the scaled variable

(4.1a) $$ \begin{align} t_\nu(r) := \frac{\nu-2\sqrt{r}}{r^{1/6}} \qquad (r>0) \end{align} $$

stays bounded. It is convenient to note that $t_\nu (r)$ satisfies the differential equation

(4.1b) $$ \begin{align} t_\nu^{\prime}(r) = -r^{-2/3} - \frac{r^{-1}}{6} t_\nu(r)\qquad (r>0). \end{align} $$

In these terms, we get the following theorem.

Theorem 4.1. There holds the expansion

(4.2) $$ \begin{align} E_2^{\text{hard}}(4r;\nu) = F(t) + \sum_{j=1}^m F_{j}^P(t)\, r^{-j/3} + r^{-(m+1)/3} \cdot O\big(e^{-t}\big)\bigg|_{t=t_\nu(r)}, \end{align} $$

which is uniformly valid when $r,\nu \to \infty $ subject to

$$\begin{align*}t_0 \leqslant t_\nu(r) \leqslant r^{1/3}, \end{align*}$$

with m being any fixed non-negative integer and $t_0$ any fixed real number. Preserving uniformity, the expansion can be repeatedly differentiated w.r.t. the variable r. Here, the $F_{j}^P$ are certain smooth functions; the first three are ²⁶

(4.4a) $$ \begin{align} F_{1}^P(t) &= -\frac{t^2}{60} F'(t) - \frac{1}{10} F"(t),\end{align} $$

(4.4b) $$ \begin{align} F_{2}^P(t) &= \Big(\frac{1}{350} + \frac{2t^3}{1575}\Big) F'(t) + \Big(\frac{11t}{1050} + \frac{t^4}{7200}\Big) F"(t) + \frac{t^2}{600} F"'(t) + \frac{1}{200} F^{(4)}(t), \end{align} $$

(4.4c) $$ \begin{align} F_3^P(t) &= -\Big(\frac{t}{1125}+\frac{41t^4}{283500}\Big) F'(t) -\Big(\frac{11 t^2}{6300}+ \frac{t^5}{47250}\Big) F"(t) \end{align} $$

$$\begin{align*} &\qquad-\Big(\frac{61}{31500}+\frac{19 t^3}{63000}+\frac{t^6}{1296000} \Big) F"'(t) -\Big(\frac{11 t}{10500} + \frac{t^4}{72000}\Big) F^{(4)}(t) \notag \\ &\qquad -\frac{t^2}{12000}F^{(5)}(t) -\frac{1}{6000}F^{(6)}(t). \notag \end{align*}$$

Figure 2 Plots of $F_{1}^P(t)$ (left panel) and $F_{2}^P(t)$ (middle panel) as in (4.4a/b). The right panel shows $F_3^P(t)$ as in (4.4c) (black solid line) with the approximations (4.3) for $r=250$ (red dotted line) and $r=2000$ (green dashed line); the parameter $\nu $ has been varied such that $t_\nu (r)$ covers the range of t on display. Note that the functions $F_{j}^P(t)$ ( $j=1,2,3$ ) are about two orders of magnitude smaller in scale than their counterparts in Figure 1.

Proof. For $r, \nu> 0$ (i.e., equivalently, $t>-2r^{1/3}$ and $s < h_\nu ^{-1}$ ), the transformations

$$\begin{align*}4r = \phi_\nu(s), \quad t = t_\nu(r) \end{align*}$$

are inverted by the expressions

(4.5) $$ \begin{align} s = \frac{t}{\big(1+\tfrac{t}2 r^{-1/3}\big)^{1/3}},\quad h_\nu = \frac{r^{-1/3}}{2\big(1+\tfrac{t}2 r^{-1/3}\big)^{2/3}}. \end{align} $$

For $t_0\leqslant t \leqslant r^{1/3}$ , we get

$$\begin{align*}s_0:= (\tfrac{2}{3})^{1/3} t_0 \leqslant (\tfrac{2}{3})^{1/3} t \leqslant s < h_\nu^{-1} \end{align*}$$

and observe that in this range of t the expressions in (4.5) expand as uniformly convergent power series in powers of $r^{-1/3}$ , starting with

$$\begin{align*}s = t - \frac{t^2}{6}r^{-1/3} + \frac{t^3}{18}r^{-2/3} -\frac{7 t^4}{324}r^{-1} + \cdots, \quad h_\nu = \frac{1}{2}r^{-1/3} - \frac{t}{6}r^{-2/3} +\frac{5 t^2}{72}r^{-1} + \cdots. \end{align*}$$

If we plug these uniformly convergent power series into the uniform expansion of Theorem 3.1,

$$\begin{align*}E_2^{\text{hard}}(4r;\nu) = E_2^{\text{hard}}(\phi_\nu(s);\nu) = F(s) + \sum_{j=1}^m F_{1}(s) h_\nu^j + h_\nu^{m+1} O(e^{-3s/2}), \end{align*}$$

we obtain the asserted form of the expansion (4.2) (as well as the claim about the repeated differentiability), simplifying the exponential error term by observing that $(3/2)^{2/3}> 1$ . In particular, the first three correction terms in (4.2) are thus

$$ \begin{align*} F_{1}^P(t) &= \frac{1}{2} F_{1}(t) - \frac{t^2}{6} F'(t),\\ F_{2}^P(t) &= \frac{1}{4} F_{2}(t) - \frac{t}{6} F_{1}(t) - \frac{t^2}{12} F_{1}^{\prime}(t) + \frac{t^3}{18} F'(t) + \frac{t^4}{72}F"(t),\\ F_{3}^P(t) &= \frac{1}{8}F_3(t) -\frac{t}{6} F_2(t)-\frac{t^2}{24} F_2^{\prime}(t)+\frac{5t^2}{72} F_1(t)+\frac{ t^3}{18} F_1^{\prime}(t)+\frac{t^4}{144} F_1^{\prime\prime}(t)\\ &\qquad -\frac{7t^4}{324} F'(t) -\frac{t^5}{108} F"(t)-\frac{t^6 }{1296}F"'(t). \end{align*} $$

Together with the expressions given in (3.4), this yields the functional form asserted in (4.4).

By (1.1), specializing Theorem 4.1 to the case of integer parameter $\nu = l$ yields the expansion

(4.6) $$ \begin{align} {\mathbb P}\big(L_{N_r} \leqslant l\big) = F(t) + \sum_{j=1}^m F_{j}^P(t)\, r^{-j/3} + r^{-(m+1)/3} \cdot O\big(e^{-t}\big) \bigg|_{t=t_l(r)}, \end{align} $$

which is uniformly valid under the conditions stated there.

Remark 4.2. In the literature, scalings are often applied to the probability distribution rather than to the expansion terms. Since $L_{N_r}$ is a an integer-valued random variable, one has to exercise some care with the scaled distribution function being piecewise constant. Namely, for $t\in {\mathbb R}$ being any fixed number, one has

$$\begin{align*}{\mathbb P}\left(\frac{L_{N_r} - 2\sqrt{r}}{r^{1/6}}\leqslant t\right) = {\mathbb P}\big(L_{N_r} \leqslant l\big),\qquad l = \big\lfloor2\sqrt{r} + tr^{1/6}\big\rfloor, \end{align*}$$

where $\lfloor \cdot \rfloor $ denotes the Gauss bracket. Thus, by defining

$$\begin{align*}t^{(r)} = \frac{\big\lfloor2\sqrt{r} + tr^{1/6}\big\rfloor - 2\sqrt{r}}{r^{1/6}} \end{align*}$$

and noting that $t^{(r)}$ stays bounded when $r\to \infty $ while t is fixed, (4.6) takes the form

(4.7) $$ \begin{align} {\mathbb P}\left(\frac{L_{N_r} - 2\sqrt{r}}{r^{1/6}}\leqslant t\right) = F\big(t^{(r)}\big) + \sum_{j=1}^m F_{j}^P\big(t^{(r)}\big)\, r^{-j/3} + O\big(r^{-(m+1)/3}\big)\qquad (r\to\infty). \end{align} $$

If one chooses to re-introduce the continuous variable t in (parts of) the expansion terms, one has to take into account that

(4.8) $$ \begin{align} t^{(r)} = t + O(r^{-1/6})\qquad (r\to\infty), \end{align} $$

where the exponent $-1/6$ in the error term is sharp. For example, this gives (as previously obtained by Baik and Jenkins [11, Thm. 1.3] using the technology of Riemann–Hilbert problems to prove the expansion and Painlevé representations to put $F_{1}^P$ into the simple functional form (4.4a))

(4.9) $$ \begin{align} {\mathbb P}\left(\frac{L_{N_r} - 2\sqrt{r}}{r^{1/6}}\leqslant t\right) = F\big(t^{(r)}\big) + F_{1}^P(t)\, r^{-1/3} + O(r^{-1/2})\qquad (r\to\infty), \end{align} $$

where the $O(r^{-1/2})$ error term is governed by the Gauss bracket in (4.8) and cannot be improved upon – completely dominating the order $O(r^{-2/3})$ correction term in (4.7). Therefore, claiming an $O(r^{-2/3})$ error term to hold in (4.9) as stated in [45, Prop. 1.1] neglects the effect of the Gauss bracket. ²⁷

Part II: Results based on the tameness hypothesis

5. De-Poissonization and the expansion of the length distribution

5.1. Expansion of the CDF

In this section, we prove (subject to a tameness hypothesis on the zeros of the generating functions in a sector of the complex plane) an expansion of the CDF ${\mathbb P}(L_n \leqslant l)$ of the length distribution near its mode. The general form of such an expansion was conjectured in the recent papers [19, 45] where approximations of the graphical form of the first few terms were provided (see [19, Figs. 4/6] and [45, Fig. 7]). Here, for the first time, we give the functional form of these terms. The underlying tool is analytic de-Poissonization, a technique that was developed in the 1990s in theoretical computer science and analytic combinatorics.

To prepare for the application of analytic de-Poissonization in the form of the Jacquet–Szpankowski Theorem A.2, we consider any fixed compact interval $[t_0,t_1]$ and a sequence of integers $l_n \to \infty $ such that

(5.1) $$ \begin{align} t_0 \leqslant t_n^*:= t_{l_n}(n)\leqslant t_1 \qquad (n=1,2,3,\ldots). \end{align} $$

When $n-n^{3/5} \leqslant r \leqslant n+n^{3/5}$ and $n\geqslant n_0$ with $n_0$ large enough (depending only on $t_0, t_1$ ), we thus get the uniform bounds ²⁸

$$\begin{align*}2\sqrt{r} + (t_0-1) r^{1/6} \leqslant l_n \leqslant 2\sqrt{r} + (t_1+1) r^{1/6}. \end{align*}$$

We write the induced Poisson generating function, and exponential generating function, of the length distribution as

(5.2) $$ \begin{align} P_k(z) := P(z; l_k) = e^{-z} \sum_{n=0}^\infty {\mathbb P}(L_n \leqslant l_k) \frac{z^n}{n!},\qquad f_k(z):= e^z P_k(z). \end{align} $$

By (1.1), we have $P_n(r) = E_2^{\text {hard}}(4r;l_n)$ for real $r>0$ , so that Theorem 4.1 (see also (4.6)) gives the expansion

(5.3) $$ \begin{align} P_n(r) = F(t) + \sum_{j=1}^m F_{j}^P(t) \, r^{-j/3} + O(r^{-(m+1)/3})\,\bigg|_{t=t_{l_n}(r)}, \end{align} $$

uniformly valid when $n - n^{3/5} \leqslant r \leqslant n + n^{3/5}$ as $n\to \infty $ , m being any fixed non-negative integer. Here, the implied constant in the error term depends only on $t_0, t_1$ , but not on the specific sequence $l_n$ . Preserving uniformity, the expansion can be repeatedly differentiated w.r.t. the variable r. In particular, using the differential equation (4.1b), we get that $P_n^{(j)}(n)$ expands in powers of $n^{-1/3}$ , starting with a leading order term of the form

(5.4a) $$ \begin{align} P_n^{(j)}(n)= (-1)^j F^{(j)} (t^*_n) n^{-2j/3} + O(n^{-(2j+1)/3})\qquad (n\to\infty); \end{align} $$

the specific cases to be used below are (see (6.7) for $P_n^{\prime }(n)$ )

(5.4b) $$ \begin{align} P_n(n) &= F(t) + F_{1}^P(t) n^{-1/3} + F_{2}^P(t) n^{-2/3} + F_{3}^P(t) n^{-1} + O(n^{-4/3})\,\bigg|_{t=t_n^*},\end{align} $$

(5.4c) $$ \begin{align} P_n^{\prime\prime}(n) &= F"(t) n^{-4/3} + \left(F_{1}^{P}\!\!\phantom{|}"(t) + \frac{5}{6} F'(t) + \frac{t}{3} F"(t) \right) n^{-5/3} \end{align} $$

(5.4d) $$ \begin{align} &+ \left(\frac{3}{2}F_{1}^{P}\!\!\phantom{|}'(t) + \frac{t}{3} F_{1}^{P}\!\!\phantom{|}"(t) + F_{2}^{P}\!\!\phantom{|}"(t)+ \frac{7 t}{36} F'(t)+\frac{t^2}{36} F"(t) \right) n^{-2}+ O(n^{-7/3})\,\bigg|_{t=t_n^*},\notag\\ P_n^{\prime\prime\prime}(n) &= -F"'(t) n^{-2} + O(n^{-7/3}) \,\bigg|_{t=t_n^*}, \end{align} $$

(5.4e) $$ \begin{align} P_n^{(4)}(n) &= F^{(4)}(t) n^{-8/3} + \left(F_{1}^{P}\!\!\phantom{|}^{(4)}(t)+5 F"'(t)+ \frac{2 t}{3} F^{(4)}(t)\right) n^{-3} + O(n^{-10/3})\,\bigg|_{t=t_n^*}, \end{align} $$

(5.4f) $$ \begin{align} P_n^{(6)}(n) &= F^{(6)}(t) n^{-4} + O(n^{-13/3})\,\bigg|_{t=t_n^*}, \end{align} $$

where the implied constants in the error terms depend only on $t_0, t_1$ .

We recall from the results of [19, Sect. 2] (note the slight differences in notation), and the proofs given there, that the exponential generating functions $f_n(z)$ are entire functions of genus zero having, for each $0<\epsilon <\pi /2$ , only finitely many zeros ²⁹ in the sector $|\!\arg z| \leqslant \pi /2 + \epsilon $ . If we denote the real auxiliary functions (cf. Definition A.1) of $f_n(r) = e^r P_n(r)$ by $a_n(r)$ and $b_n(r)$ , the expansion (5.3), and its derivatives based on (4.1b), give (cf. also (6.8))

(5.5) $$ \begin{align} a_n(r) = r + O(r^{1/3}),\qquad b_n(r) = r + O(r^{2/3}), \end{align} $$

uniformly valid when $n - n^{3/5} \leqslant r \leqslant n + n^{3/5}$ as $n\to \infty $ ; the implied constants in the error terms depend only on $t_0, t_1$ .

Analytically, we lack the tools to study the asymptotic distribution of the finitely many zeros of $f_n(z)$ in the sector $|\!\arg z| \leqslant \pi /2 + \epsilon $ as $n\to \infty $ . Numerically, we proceed as follows. The meromorphic logarithmic derivative of $f_n(z)$ takes the form [19, §3.1]

$$\begin{align*}\frac{f_n^{\prime}(z)}{f_n(z)} = 1 - \frac{v_{l_n}(z)}{z}, \end{align*}$$

where $v_l$ satisfies a Jimbo–Miwa–Okamoto $\sigma $ -form of the Painlevé III equation [19, Eq. (31)], or alternatively, a certain Chazy I equation [19, Eq. (34)]. Because of $f_n(0)=1$ , the zeros of $f_n$ are in a one-to-one correspondence to the pole field of the meromorphic function $v_{l_n}$ . Fornberg and Weideman [39] developed a numerical method, the pole field solver, specifically for the task of numerically studying the pole fields of equations of the Painlevé class. They documented results for Painlevé I [39], Painlevé II [40], its imaginary variant [41] and, together with Fasondini, for (multivalued) variants of the Painlevé III, V and VI equations [34, 35].

Now, extensive numerical experiments with the pole field solver applied to $v_{l_n}$ (which will be documented in a separate publication) strongly hint at the property that the zeros of the exponential generating functions $f_n(z)$ in the sectors $|\!\arg z| \leqslant \pi /2 + \epsilon $ satisfy a uniform tameness condition as in Definition A.2 (see also Remark A.6): the zeros are neither coming too close to the positive real axis nor are they getting too large. Given this state of affairs, the results on the expansions of the length distribution will be subject to the following:

Figure 3 Plots of $F_{1}^D(t)$ (left panel), $F_{2}^D(t)$ (middle panel) as in (5.8); both agree with the numerical prediction of their graphical form given in the left panels of [19, Figs. 4/6]. The right panel shows $F_3^D(t)$ as in (5.8) (black solid line) with the approximations (5.7) for $n=250$ (red $+$ ), $n=500$ (green $\circ $ ) and $n=1000$ (blue $\bullet $ ); the integer l has been varied such that $t_l(n)$ spreads over the range of t displayed here. Evaluation of (5.7) uses the table of exact values of ${\mathbb P}(L_n\leqslant l)$ up to $n=1000$ that was compiled in [19].

Tameness hypothesis. For any real $t_0<t_1$ and sequence of integers $l_n\to \infty $ satisfying (5.1), the zeros of the induced family $f_n(z)$ of exponential generating functions (5.2) are uniformly tame (see Definition A.2), with parameters and implied constants only depending on $t_0$ and $t_1$ .

Theorem 5.1. Let $t_0<t_1$ be any real numbers and assume the tameness hypothesis. Then there holds the expansion

(5.6) $$ \begin{align} {\mathbb P}(L_n \leqslant l) = F(t) + \sum_{j=1}^m F_{j}^D(t)\, n^{-j/3} + O(n^{-(m+1)/3})\,\bigg|_{t=t_l(n)}, \end{align} $$

which is uniformly valid when $n,l \to \infty $ subject to $t_0 \leqslant t_l(n) \leqslant t_1$ with m being any fixed non-negative integer. Here, the $F_{j}^D$ are certain smooth functions; the first three are ³⁰

(5.8a) $$ \begin{align} F_{1}^D(t) &= -\frac{t^2}{60} F'(t) - \frac{3}{5} F"(t), \end{align} $$

(5.8b) $$ \begin{align} F_{2}^D(t) &= \Big(-\frac{139}{350} + \frac{2t^3}{1575}\Big) F'(t) + \Big(-\frac{43t}{350} + \frac{t^4}{7200}\Big) F"(t) + \frac{t^2}{100} F"'(t) + \frac{9}{50} F^{(4)}(t), \end{align} $$

(5.8c) $$ \begin{align}F_3^D(t) &= -\Big(\frac{562 t }{7875} +\frac{41t^4}{283500}\Big) F'(t) +\Big(\frac{t^2}{300}-\frac{t^5}{47250}\Big) F"(t)\\&\qquad+\Big(\frac{5137}{15750}+\frac{9 t^3}{7000}-\frac{t^6}{1296000}\Big) F"'(t) +\Big(\frac{129 t}{1750}-\frac{t^4}{12000}\Big) F^{(4)}(t)\notag\\ &\qquad-\frac{3 t^2}{1000}F^{(5)}(t) -\frac{9}{250}F^{(6)}(t).\notag \end{align} $$

Proof. Following up the preparations preceding the formulation of the theorem, the tameness hypothesis allows us to apply Corollary A.7, bounding $f_n(z)=e^z P_n(z)$ by

(5.9) $$ \begin{align} \big| f_n(re^{i\theta}) \big| \leqslant \begin{cases} 2f_n(r) e^{-\frac{1}{2}\theta^2 r}, &\qquad 0\leqslant |\theta| \leqslant r^{-2/5},\\ 2f_n(r) e^{-\frac{1}{2}r^{1/5}}, &\qquad r^{-2/5}\leqslant |\theta| \leqslant \pi, \end{cases} \end{align} $$

for $n - n^{3/5} \leqslant r \leqslant n + n^{3/5}$ and $n\geqslant n_0$ when $n_0$ is sufficiently large (depending on $t_0$ , $t_1$ ). Using the trivial bounds (for $r>0$ and $|\theta |\leqslant \pi $ )

$$\begin{align*}0\leqslant f_n(r)\leqslant e^r, \qquad 0\leqslant P_n(r) \leqslant 1, \qquad 1 - \tfrac{1}{2}\theta^2 \leqslant \cos\theta, \end{align*}$$

the first case in (5.9) can be recast in form of the bound

$$\begin{align*}\big|P_n(re^{i\theta})\big| \leqslant 2 P_n(r) e^{r \left(1-\cos\theta - \tfrac12 \theta^2\right)} \leqslant 2, \end{align*}$$

which proves condition (I) of Theorem A.2 with $B=2$ , $D=1$ , $\beta =0$ and $\delta =2/5$ , whereas the second case implies

$$\begin{align*}|f_n(n e^{i\theta})| \leqslant 2 f_n(n) e^{-\tfrac12 n^{1/5}} \leqslant 2 \exp\big(n - \tfrac12 n^{1/5}\big), \end{align*}$$

which proves condition (O) of Theorem A.2 with $A=1/2$ , $C=0$ , $\alpha =1/5$ and $\gamma =0$ . Hence, there holds the Jasz expansion (A.7) – namely,

$$\begin{align*}{\mathbb P}(L_n \leqslant l_n) = P_n(n) + \sum_{j=2}^M b_j(n) P_n^{(j)}(n) + O(n^{-(M+1)/5}) \end{align*}$$

for any $M=0,1,2,\ldots $ as $n\geqslant n_1$ ; here, $n_1$ and the implied constant depend on $t_0$ , $t_1$ . By noting that the diagonal Poisson–Charlier polynomials $b_j$ have degree $\leqslant \lfloor j/2\rfloor $ and by choosing M large enough, the expansions (5.4) of $P_n^{(j)}(n)$ in terms of powers of $n^{-1/3}$ yield that there are smooth functions $F_{j}^D$ such that

$$\begin{align*}{\mathbb P}(L_n \leqslant l_n) = F(t) + \sum_{j=1}^m F_{j}^D(t)\, n^{-j/3} + O(n^{-(m+1)/3})\,\bigg|_{t=t_n^*} \end{align*}$$

as $n\to \infty $ , with m being any fixed non-negative integer. Given the uniformity of the bound for fixed $t_0$ and $t_1$ , we can replace $l_n$ by l and $t_n^*$ by $t_l(n)$ as long as we respect $t_0 \leqslant t_l(n) \leqslant t_1$ . This finishes the proof of (5.6).

The first three functions $F_{1}^D$ , $F_{2}^D(t)$ , $F_{3}^D(t)$ can be determined using the particular case (A.8) of the Jasz expansion from Example A.3 (which applies here because of (5.4a)) – namely,

$$\begin{align*}{\mathbb P}(L_n \leqslant l_n) =P_n(n) - \frac{n}{2} P_n^{\prime\prime}(n) + \frac{n}{3} P_n^{\prime\prime\prime}(n) + \frac{n^2}{8} P_n^{(4)}(n) - \frac{n^3}{48}P_n^{(6)}(n) + O(n^{-4/3}). \end{align*}$$

Inserting the formulae displayed in (5.4), we thus obtain

(5.10a) $$ \begin{align} F_{1}^D(t) &= F_{1}^P(t) - \frac{1}{2} F"(t), \end{align} $$

(5.10b) $$ \begin{align} F_{2}^D(t) &= F_{2}^P(t) - \frac{1}{2} F_{1}^{P}\!\!\phantom{|}"(t) - \frac{5}{12}F'(t) - \frac{t}{6}F"(t) + \frac{1}{8} F^{(4)}(t), \end{align} $$

(5.10c) $$ \begin{align} F_3^D(t) &= F_3^P(t)-\frac{1}{2} F_{2}^{P}\!\!\phantom{|}"(t)-\frac{3 }{4} F_{1}^{P}\!\!\phantom{|}'(t)-\frac{t}{6} F_{1}^{P}\!\!\phantom{|}"(t)+\frac{1}{8} F_{1}^{P}\!\!\phantom{|}^{(4)}(t) \end{align} $$

$$\begin{align*} &\qquad -\frac{7t}{72} F'(t)-\frac{t^2}{72} F"(t)+\frac{7}{24} F^{(3)}(t)+\frac{t}{12} F^{(4)}(t)-\frac{1}{48} F^{(6)}(t).\notag \end{align*}$$

Together with the expressions given in (4.4), this yields the functional form asserted in (5.8).

5.2. Expansion of the PDF

Subject to its assumptions, Theorem 5.1 implies for the PDF of the length distribution that

$$ \begin{align*} {\mathbb P}(L_n=l) &= {\mathbb P}(L_n \leqslant l) - {\mathbb P}(L_n \leqslant l-1) \\ &\qquad = \big(F(t_l(n)) - F(t_{l-1}(n))\big) + \sum_{j=1}^m \big(F_{j}^D(t_l(n))-F_{j}^D(t_{l-1}(n))\big)\, n^{-j/3} + O(n^{-(m+1)/3}). \end{align*} $$

Applying the central differencing formula (which is, basically, just a Taylor expansion for smooth G centered at the midpoint)

$$\begin{align*}G(t+h) - G(t) = h G'\big(t+\tfrac{h}{2} \big) + \frac{h^3}{24}G"'\big(t+\tfrac{h}{2} \big) + \frac{h^5}{1920}G^{(5)} \big(t+\tfrac{h}{2} \big) +\frac{h^7 }{322560} G^{(7)} \big(t+\tfrac{h}{2} \big)+ \cdots , \end{align*}$$

with increment $h=n^{-1/6}$ , we immediately get the following corollary of Theorem 5.1.

Corollary 5.2. Let $t_0<t_1$ be any real numbers, and assume the tameness hypothesis. Then there holds the expansion

(5.11) $$ \begin{align} n^{1/6}\, {\mathbb P}(L_n=l) = F' (t) + \sum_{j=1}^m F_{j}^*(t) n^{-j/3} + O(n^{-(m+1)/3})\,\bigg|_{t=t_{l-1/2}(n)}, \end{align} $$

which is uniformly valid when $n,l \to \infty $ subject to the constraint $t_0 \leqslant t_{l-1/2}(n) \leqslant t_1$ with m being any fixed non-negative integer. Here, the $F_{j}^*$ are certain smooth functions; in particular, ³¹

(5.13a) $$ \begin{align} F_{1}^*(t) &= -\frac{t}{30} F'(t) - \frac{t^2}{60}F"(t) - \frac{67}{120} F"'(t),\end{align} $$

(5.13b) $$ \begin{align} F_{2}^*(t) &= \frac{2t^2}{525}F'(t) + \Big(-\frac{629}{1200}+\frac{23t^3}{12600}\Big)F"(t) + \Big(-\frac{899t}{8400}+\frac{t^4}{7200}\Big)F"'(t) \end{align} $$

$$\begin{align*} &\qquad\quad + \frac{67t^2}{7200}F^{(4)}(t) + \frac{1493}{9600}F^{(5)}(t), \notag \end{align*}$$

(5.13c) $$ \begin{align}F_3^*(t) &= -\Big(\frac{373}{5250} + \frac{41 t^3}{70875} \Big) F'(t) -\Big(\frac{1781 t}{28000} + \frac{71 t^4 }{283500}\Big) F"(t)\\&\qquad+\Big(\frac{63 t^2}{8000}-\frac{13 t^5}{504000}\Big) F"'(t) +\Big(\frac{41473 }{112000}+\frac{13 t^3}{12096} -\frac{t^6}{1296000} \Big) F^{(4)}(t)\notag\\ &\qquad+\Big(\frac{131057 t}{2016000} - \frac{67 t^4}{864000}\Big) F^{(5)}(t) -\frac{1493 t^2}{576000} F^{(6)}(t) -\frac{232319}{8064000}F^{(7)}(t).\notag \end{align} $$

Remark 5.3. The case $m=0$ of Corollary 5.2 gives

$$\begin{align*}n^{1/6} \, {\mathbb P}(L_n=l) = F'(t_{l-1/2}(n)) + O(n^{-1/3}), \end{align*}$$

where the exponent in the error term cannot be improved. By noting

$$\begin{align*}t_{l-1/2}(n) = t_l(n) - \tfrac{1}{2}n^{-1/6}, \end{align*}$$

we understand that, for fixed large n, visualizing the discrete length distribution near its mode by plotting the points

$$\begin{align*}\big(t_l(n),n^{1/6}\,{\mathbb P}(L_n=l)\big) \end{align*}$$

next to the graph $(t,F'(t))$ introduces a perceivable bias; namely, all points are shifted by an amount of $n^{-1/6}/2$ to the right of the graph. Exactly such a bias can be observed in the first ever published plot of the PDF vs. the density of the Tracy–Widom distribution by Odlyzko and Rains in [62, Fig. 1]: the Monte-Carlo data for $n=10^6$ display a consistent shift by $0.05$ .

A bias free plot is shown in Figure 5, which, in addition, displays an improvement of the error of the limit law by a factor of $O(n^{-2/3})$ that is obtained by adding the first two finite-size correction terms.

Figure 4 Plots of $F_{1}^*(t)$ (left panel) and $F_{2}^*(t)$ (middle panel) as in (5.13a/b); both agree with the numerical prediction of their graphical form given in the right panels of [19, Figs. 4/6]. The right panel shows $F_3^*$ as in (5.13c) (black solid line) with the approximations (5.12) for $n=250$ (red $+$ ), $n=500$ (green $\circ $ ) and $n=1000$ (blue $\bullet $ ); the integer l has been varied such that $t_{l-1/2}(n)$ spreads over the range of t displayed here. Evaluation of (5.12) uses the table of exact values of ${\mathbb P}(L_n= l)$ up to $n=1000$ that was compiled in [19].

Figure 5 The exact discrete length distribution ${\mathbb P}(L_{n}=l)$ (blue bars centered at the integers l) vs. the asymptotic expansion (5.11) for $m=0$ (the Baik–Deift–Johansson limit, dotted line) and for $m=2$ (the limit with the first two finite-size correction terms added, solid line). Left: $n=100$ ; right: $n=1000$ . The expansions are displayed as functions of the continuous variable $\nu $ , evaluating the right-hand side of (5.11) in $t=t_{\nu -1/2}(n)$ . The exact values are from the table compiled in [19]. Note that a graphically accurate continuous approximation of the discrete distribution must intersect the bars right in the middle of their top sides: this is, indeed, the case for $m=2$ (except at the left tail for $n=100$ ). In contrast, the uncorrected limit law ( $m=0$ ) is noticeable inaccurate for this range of n.

6. Expansions of Stirling-type formulae

In our work [19], we advocated the use of a Stirling-type formula to approximate the length distribution for larger n (because of being much more efficient and accurate than Monte-Carlo simulations). To recall some of our findings there, let us denote the exponential generating function and its Poisson counterpart simply by

$$\begin{align*}f(z) = \sum_{n=0}^\infty {\mathbb P}(L_n\leqslant l) \frac{z^n}{n!},\qquad P(z) = e^{-z} f(z), \end{align*}$$

suppressing the dependence on the integer parameter l from the notation for the sake of brevity. It was shown in [19, Thm. 2.2] that the entire function f is H-admissible so that there is the normal approximation (see Definition A.1 and Theorem A.4)

(6.1) $$ \begin{align} {\mathbb P}(L_n\leqslant l) = \frac{n! f(r)}{r^n \sqrt{2\pi b(r)}} \left(\exp\left(-\frac{(n-a(r))^2}{2b(r)}\right)+ o(1)\right) \qquad (r\to \infty) \end{align} $$

uniformly in $n=0,1,2,\ldots $ while l is any fixed integer. Here, $a(r)$ and $b(r)$ are the real auxiliary functions

$$\begin{align*}a(r) = r \frac{f'(r)}{f(r)},\qquad b(r) = r a'(r). \end{align*}$$

We consider the two cases $r=r_n$ , $a(r_n)=n$ and $r=n$ . After dividing (6.1) by the classical Stirling factor (which does not change anything of substance; see Remark 6.2)

(6.2) $$ \begin{align} \tau_n := \frac{n!}{\sqrt{2\pi n}} \left(\frac{e}{n}\right)^n \sim 1+\frac{n^{-1}}{12} + \frac{n^{-2}}{288} + \cdots, \end{align} $$

we get after some rearranging of terms the Stirling-type formula $S_{n,l}$ ( $r=r_n$ ) and the simplified Stirling-type formula $\tilde S_{n,l}$ ( $r=n$ ):

(6.3a) $$ \begin{align} S_{n,l} &:= \frac{P(r_n)}{\sqrt{b(r_n)/n}} \exp\left(n\,\Lambda\left(\frac{r_n-n}n\right)\right),\quad \Lambda(h) = h - \log(1+h), \end{align} $$

(6.3b) $$ \begin{align} \tilde S_{n,l} &:= \frac{P(n)}{\sqrt{b(n)/n}} \exp\left(-\frac{(n-a(n))^2}{2b(n)}\right). \end{align} $$

As shown in [19], both approximations are amenable for a straightforward numerical evaluation using the tools developed in [14, 15, 20]. For fixed l, the normal approximation (6.1) implies

$$\begin{align*}{\mathbb P}(L_n\leqslant l) = S_{n,l} \cdot (1 + o(1))\qquad (n\to\infty); \end{align*}$$

but numerical experiments reported in [19, Fig. 3, Eq. (8b)] suggest that there holds

$$\begin{align*}{\mathbb P}(L_n\leqslant l) = S_{n,l} + O(n^{-2/3}) \end{align*}$$

uniformly when $n,l\to \infty $ while $t_l(n)$ stays bounded. Subject to the tameness hypothesis of Theorem 5.1, we prove this observation as well as its counterpart for the simplified Stirling-type formula, thereby unveiling the functional form of the error term $O(n^{-2/3})$ : ³²

Theorem 6.1. Let $t_0<t_1$ be any real numbers, and assume the tameness hypothesis. Then, for the Stirling-type formula $S_{n,l}$ and its simplification $\tilde S_{n,l}$ , there hold the expansions (note that both are starting at $j=2$ )

(6.4a) $$ \begin{align} {\mathbb P}(L_n \leqslant l) &= S_{n,l} + \sum_{j=2}^m F_{j}^S\big(t_l(n)\big)\, n^{-j/3} + O(n^{-(m+1)/3}), \end{align} $$

(6.4b) $$ \begin{align} {\mathbb P}(L_n \leqslant l) &= \tilde S_{n,l} + \sum_{j=2}^m \tilde F_{j}^S\big(t_l(n)\big)\, n^{-j/3} + O(n^{-(m+1)/3}), \end{align} $$

which are uniformly valid when $n,l \to \infty $ subject to $t_0 \leqslant t_l(n) \leqslant t_1$ with m being any fixed non-negative integer. Here, the $F_{j}^S$ and $\tilde F_{j}^S$ are certain smooth functions – the first being ³³

(6.5a) $$ \begin{align} F_{2}^S(t) &= -\frac{3}{4}\frac{F'(t)^4}{F(t)^3} + \frac{3}{2}\frac{F'(t)^2F"(t)}{F(t)^2}- \frac{3}{8}\frac{F"(t)^2}{F(t)} - \frac{1}{2}\frac{F'(t) F"'(t)}{F(t)} + \frac{1}{8} F^{(4)}_2(t) , \end{align} $$

(6.5b) $$ \begin{align} \tilde F_{2}^S(t) &= -\frac{1}{2}F'(t) + \frac{1}{4}\frac{F'(t)^4}{F(t)^3}- \frac{3}{8}\frac{F"(t)^2}{F(t)} + \frac{1}{8} F^{(4)}_2(t). \end{align} $$

The solution $r_n$ of the equation $a(r_n)=n$ , required to evaluate $S_{n,l}$ , satisfies the expansion ³⁴

$$\begin{align*}r_n = n + \frac{F'\big(t_l(n) \big)}{F \big(t_l(n) \big)} n^{1/3} + O(1), \end{align*}$$

which is uniformly valid under the same conditions.

Proof. We restrict ourselves to the case $m=2$ , focusing on the concrete functional form of the expansion terms; nevertheless, the general form of the expansions (6.4) should become clear along the way.

Preparatory steps. Because of $P(r) = E_2^{\text {hard}}(4r;l)$ (using the notation preceding Theorem 6.1), Theorem 4.1 gives that

(6.6) $$ \begin{align} P(r) = {\left.F (t) + F_{1}^P (t) \cdot r^{-1/3} + F_{2}^P (t) \cdot r^{-2/3} + O(r^{-1})\,\right|}_{t=t_l(r)}, \end{align} $$

which is uniformly valid when $r,l \to \infty $ subject to the constraint $t_0 \leqslant t_l(r) \leqslant t_1$ (the same constraint applies to the expansions to follow). Preserving uniformity, the expansion can be repeatedly differentiated w.r.t. the variable r, which yields by using the differential equation (4.1b) satisfied by $t_l(r)$ (cf. also (5.4a))

(6.7) $$ \begin{align} P'(r) = {\left. -F'(t) r^{-2/3} - \Big(F_{1}^P\!\!\phantom{|}'(t) + \frac{t}{6} F'(t)\Big)r^{-1} + O(r^{-4/3})\, \right|}_{t=t_l(r)}. \end{align} $$

Recalling $f(r) = e^r P(r)$ , we thus get

(6.8a) $$ \begin{align} a(r) = r + r \frac{P'(r)}{P(r)} = {\left. r + a_1(t) r^{1/3} + a_2(t) + O(r^{-1/3})\, \right|}_{t=t_l(r)} \end{align} $$

with the coefficient functions

(6.8b) $$ \begin{align} a_1(t) = -\frac{F'(t)}{F(t)},\qquad a_2(t) = -\frac{1}{F(t)}\Big(F_{1}^P\!\!\phantom{|}'(t) + \frac{t}{6} F'(t) \Big) + \frac{F_{1}^P(t)F'(t)}{F(t)^2}; \end{align} $$

a further differentiation yields

(6.8c) $$ \begin{align} b(r) = r a'(r) = {\left. r- a_1'(t) r^{2/3} + \Big(\frac{1}{3} a_1(t) - \frac{t}{6}a_1^{\prime}(t) - a_2^{\prime}(t) \Big) r^{1/3} + O(1)\, \right|}_{t=t_l(r)}. \end{align} $$

The simplified Stirling-type formula. Here, we have $r=n$ , and we write $t_* := t_l(n)$ to be brief. By inserting the expansions (6.6) and (6.8) into the expression (6.3b), we obtain after a routine calculation with truncated power series and collecting terms as in (5.10) that

(6.9) $$ \begin{align} \tilde S_{n,l} = F(t_*) + F_{1}^D(t_*) n^{-1/3} + \Big(F_{2}^D(t_*) - \tilde F_{2}^S(t_*) \Big) n^{-2/3} + O(n^{-1}), \end{align} $$

where the remaining $\tilde F_{2}^S(t)$ is given by (6.5b); a subtraction from (5.6) yields (6.4b).

The Stirling-type formula. Here, we have $r=r_n$ , and we have to distinguish between $t_*$ and

$$\begin{align*}t_l(r) = t_* \cdot (n/r)^{1/6} + 2 \frac{\sqrt{n}-\sqrt{r}}{r^{1/6}}. \end{align*}$$

By inserting the expansion

(6.10a) $$ \begin{align} r_n = n + r_1(t_*) n^{1/3} + r_2(t_*) + O(n^{-1/3}) \end{align} $$

into $t_n(r_n)$ and $a(r_n)$ , we obtain

(6.10b) $$ \begin{align} t_l(r_n) &= t_* - r_1(t_*) n^{-1/3} - \Big(r_2(t_*) + \frac{t_*}{6} r_1(t_*)\Big) n^{-2/3} + O(n^{-1}) \end{align} $$

(6.10c) $$ \begin{align} a(r_n) &= n +(a_1(t_*) + r_1(t_*)) n^{1/3} + \Big(a_2(t_*) + r_2(t_*) -r_1(t_*) a_1^{\prime}(t_*) \Big) + O(n^{-1/3}). \end{align} $$

Thus, the solution of $a(r_n) = n$ , which by Theorem A.4 is unique, leads to the relations

(6.10d) $$ \begin{align} r_1(t) = -a_1(t),\qquad r_2(t) = - a_2(t) -a_1(t)a_1^{\prime}(t). \end{align} $$

By inserting, first, the expansions (6.10) into the expansions (6.6) and (6.8) for the particular choice $r=r_n$ and, next, the thus obtained results into the expression (6.3a), we obtain after a routine calculation with truncated power series and collecting terms as in (5.10) that

(6.11) $$ \begin{align} S_{n,l} = F(t_*) + F_{1}^D(t_*) n^{-1/3} + \Big(F_{2}^D(t_*) - F_{2}^S(t_*) \Big) n^{-2/3} + O(n^{-1}), \end{align} $$

where the remaining $F_{2}^S(t)$ is given by (6.5a); a subtraction from (5.6) yields (6.4a).

To validate the expansions (6.4) and the formulae (6.5a/b), Figure 6 plots the approximations

Figure 6 Left panel: plots of $\tilde F_{2}^S(t)$ (solid line) and $\tilde F_{2}^S(t)$ (dash-dotted line) as in (6.5). The middle and right panel show the approximations of $F_{3}^S(t)$ and $\tilde F_{3}^S(t)$ in (6.12) for $n=250$ (red $+$ ), $n=500$ (green $\circ $ ) and $n=1000$ (blue $\bullet $ ); the integer l has been varied such that $t_l(n)$ spreads over the range of the variable t displayed here; the dotted line displays a polynomial fit to the data points of degree $30$ to help visualizing their joint graphical form. Evaluation of (6.12) uses the table of exact values of ${\mathbb P}(L_n\leqslant l)$ up to $n=1000$ that was compiled in [19].

(6.12a) $$ \begin{align} F_{3}^S\big(t_l(n)\big) &\approx n^{-1} \Big( {\mathbb P}(L_n \leqslant l) - S_{n,l} - F_{2}^S\big(t_l(n)\big) n^{-2/3}\Big), \end{align} $$

(6.12b) $$ \begin{align} \tilde F_{3}^S\big(t_l(n)\big) &\approx n^{-1} \Big( {\mathbb P}(L_n \leqslant l) - \tilde S_{n,l} - \tilde F_{2}^S\big(t_l(n)\big) n^{-2/3}\Big) \end{align} $$

for $n=250$ , $n=500$ and $n=1000$ , varying the integer l in such a way that $t=t_l(n)$ spreads over $[-6,3]$ . The plot suggests the following observations:

• • Apparently there holds $\tilde F_{2}^S(t) < F_{2}^S(t) < 0$ for $t\in [-6,3]$ , which if generally true would imply

  $$\begin{align*}{\mathbb P}(L_n \leqslant l) < S_{n,l} < \tilde S_{n,l} \end{align*}$$
  for n being sufficiently large and l near the mode of the distribution. This one-sided approximation of the length distribution by the Stirling formula $S_{n,l}$ from above is also clearly visible in [19, Tables 1/2].

• • Comparing $F_{2}^S(t)$ in Figure 6 to $F_{2}^D(t)$ in Figure 3 shows that the maximum error

  $$\begin{align*}\qquad\quad \max_{l=1,\ldots,n} \left| {\mathbb P}(L_n \leqslant l) - \big(F(t) + F_{1}^D(t) n^{-1/3}\big)\big|_{t=t_l(n)}\right| \approx n^{-2/3} \|F_{2}^D\|_\infty \approx 0.25n^{-2/3} \end{align*}$$
  of approximating the length distribution by the first finite-size correction in Theorem 5.1 is about an order of magnitude larger than the maximum error of the Stirling-type formula,
  $$\begin{align*}\max_{l=1,\ldots,n} \left| {\mathbb P}(L_n \leqslant l) - S_{n,l}\right| \approx n^{-2/3} \|F_{2}^S\|_\infty \approx 0.031n^{-2/3}. \end{align*}$$
  This property of the Stirling-type formula was already observed in [19, Fig. 3] and was used there to approximate the graphical form of $F_{2}^D(t)$ (see [19, Fig. 6]).

Remark 6.2. If one includes the classical Stirling factor (6.2) into the Stirling-type formula by replacing (6.3a) with the unmodified normal approximation (6.1), that is, with

$$\begin{align*}S_{n,l}^* := \tau_n\frac{P(r_n)}{\sqrt{b(r_n)/n}} \exp\left(n\,\Lambda\left(\frac{r_n-n}n\right)\right) = \tau_n S_{n,l}, \end{align*}$$

Theorem 6.1 would remain valid: in fact, multiplication of (6.4a) by the expansion (6.2) of $\tau _n$ in powers of $n^{-1}$ gives, by taking (5.6) into account,

$$\begin{align*}{\mathbb P}(L_n \leqslant l) = S^*_{n,l} + \sum_{j=2}^m F_{j}^{S^*}\big(t_l(n)\big)\, n^{-j/3} + O(n^{-(m+1)/3}), \end{align*}$$

where the first two coefficient functions are

$$\begin{align*}F_{2}^{S^*}(t)= F_{2}^{S}(t), \quad F_{3}^{S^*}(t)= F_{3}^{S}(t) - \frac{1}{12} F(t), \quad\ldots \;. \end{align*}$$

Because of $\lim _{t\to +\infty }F(t) = 1$ , we would loose the decay of $F_{3}^S(t)$ for large t, leaving us with a nonzero residual value coming from the classical Stirling factor. For this reason, we recommend dropping the factor $\tau _n$ , thereby resolving an ambiguity expressed in [19, Fn. 28].

7. Expansions of expected value and variance

Lifting the expansion (5.11) of the PDF of the length distribution to one of the expected value and variance requires a control of the tails (of the distribution itself and of the expansion terms) which, at least right now, we can only conjecture to hold true.

To get to a reasonable conjecture, we recall the tail estimates for the discrete distribution (see [7, Eqn. (9.6/9.12)]),

$$ \begin{align*} {\mathbb P}(L_n= l) &\leqslant C e^{- c |t_l(n)|^3} & (t_l(n) \leqslant t_0 < 0),\\ {\mathbb P}(L_n= l) &\leqslant C e^{-c |t_l(n)|^{3/5}} & ( 0< t_1 \leqslant t_l(n)), \end{align*} $$

when n is large enough with $c>0$ being some absolute constant and $C>0$ a constant that depends on $t_0, t_1$ .

However, from Theorem 5.1 and its proof we see that the $F_{j}^*(t)$ take the form

(7.1) $$ \begin{align} F(t) \cdot \big(\text{rational polynomial in terms of the form } \operatorname{tr}((I-K_0)^{-1}K)|_{L^2(t,\infty)}\big), \end{align} $$

where the kernels K are finite sums of rank one kernels with factors of the form (2.2). The results of Section 2.3 thus show that the $F_{j}^*(t)$ are exponentially decaying when $t\to \infty $ . Now, looking at the left tail, the (heuristic) estimate of the largest eigenvalue of the Airy operator ${\mathbf K_0}$ on $L^2(t,\infty )$ as given in the work of Tracy and Widom [78, Eq. (1.23)] shows a superexponential growth bound of the operator norm

$$\begin{align*}\|({\mathbf I}-{\mathbf K_0})^{-1}\| \leqslant C |t|^{-3/4} e^{c |t|^{3/2}} \qquad (t\leqslant t_0). \end{align*}$$

This, together with the superexponential decay (see [6, Cor. 1.3] and [28, Thm. 1] for the specific constants)

$$\begin{align*}0\leqslant F(t) \leqslant C |t|^{-1/8} e^{-c |t|^3} \qquad (t\leqslant t_0) \end{align*}$$

of the Tracy–Widom distribution itself, and with an at most polynomial growth of the trace norms $\|{\mathbf K}\|_{{\mathcal J}^1(t,\infty )}$ as $t\to -\infty $ , shows that the bounds for the discrete distribution find a counterpart for the expansion terms $F_{j}^*(t)$ :

(7.2) $$ \begin{align} \begin{aligned} |F_{j}^*(t)| &\leqslant C e^{- c |t|^3} &\quad (t\leqslant t_0 < 0),\\ |F_{j}^*(t)| &\leqslant C e^{-c |t|} &\quad ( 0< t_1 \leqslant t). \end{aligned} \end{align} $$

Thus, assuming an additional amount of uniformity that would allow us to absorb the exponentially small tails in the error term of Corollary 5.2, we conjecture the following:

Uniform Tails Hypothesis. The expansion (5.11) can be sharpened to include the tails in the form

(7.3) $$ \begin{align} n^{1/6}\,{\mathbb P}(L_n=l) = F' (t) + \sum_{j=1}^m F_{j}^*(t) n^{-j/3} + n^{-(m+1)/3} \cdot O\Big( e^{-c |t|^{3/5}}\Big)\,\bigg|_{t=t_{l-1/2}(n)}, \end{align} $$

uniformly valid in $l=1,\ldots ,n$ as $n\to \infty $ .

We now follow the ideas sketched in our work [19, §4.3] on the Stirling-type formula. By shift and rescale, the expected value of $L_n$ can be written in the form

$$\begin{align*}{\mathbb E}(L_n) = \sum_{l=1}^n l\cdot {\mathbb P}(L_n=l) = 2\sqrt{n} +\frac{1}{2} + \sum_{l=1}^n t_{l-1/2}(n)\cdot n^{1/6}\, {\mathbb P}(L_n=l). \end{align*}$$

Inserting the expansion (7.3) of the uniform tail hypothesis gives, since its error term is uniformly summable,

(7.4) $$ \begin{align} {\mathbb E}(L_n) = 2\sqrt{n} +\frac{1}{2} + \sum_{j=0}^m \mu_j^{(n)} n^{1/6 - j/3} + O(n^{-1/6-m/3}), \end{align} $$

with coefficients (still depending on n, though), writing $F_0^* := F'$ ,

$$\begin{align*}\mu^{(n)}_j := n^{-1/6} \sum_{l=1}^n t_{l-1/2}(n) F_j^*\big(t_{l-1/2}(n)\big). \end{align*}$$

By the tail estimates (7.2), we have, writing $a :=- n^{-1/6}(2\sqrt {n} + \tfrac 12)$ and $h:=n^{-1/6}$ ,

$$\begin{align*}\mu^{(n)}_j = h \sum_{l=-\infty}^\infty (a+lh) F_j^*(a+l h) + O(e^{-c n^{1/3}}). \end{align*}$$

Now, based on a precise description of its pole field in [52], it is known that the Hastings–McLeod solution of Painlevé II and a fortiori, by the Tracy–Widom theory [78], also F and its derivatives can be continued analytically to the strip $|\Im z| < 2.9$ . Therefore, we assume the following:

Uniform strip hypothesis. The $F_j^*$ ( $j=1,\ldots ,m$ ) extend analytically to a strip $|\Im z| \leqslant s$ of the complex z-plane, uniformly converging to $0$ as $z\to \infty $ in that strip.

Under that hypothesis, a classical result about the rectangular rule in quadrature theory (see, for example, [26, Eq. (3.4.14)]) gives

$$\begin{align*}h \sum_{l=-\infty}^\infty (a+lh) F_j^*(a+l h) = \int_{-\infty}^\infty F_j^*(t)\,dt + O(e^{-\pi s /h}). \end{align*}$$

Thus, the $\mu _j^{(n)}$ and their limit quantities

$$\begin{align*}\mu_j = \int_{-\infty}^\infty t F_j^*(t)\,dt \end{align*}$$

differ only by an exponential small error of at most $O(e^{-c n^{1/6}})$ , which can be absorbed in the error term of (7.4); an illustration of such a rapid convergence is given in [19, Table 3] for the case $j=0$ .

The functional form of $F_0^* = F'$ , $F_1^*$ , $F_2^*$ and $F_3^*$ – namely, being a linear combination of higher order derivatives of F with polynomial coefficients (see (5.13)) – allows us to express $\mu _0,\mu _1,\mu _2, \mu _3$ in terms of the moments

$$\begin{align*}M_j := \int_{-\infty}^\infty t^j F'(t)\,dt \end{align*}$$

of the Tracy–Widom distribution F. In fact, repeated integration by parts yields the simplifying rule (where $k\geqslant 1$ )

$$\begin{align*}\int_{-\infty}^\infty t^j\, F^{(k)}(t) \,dt = \begin{cases} \displaystyle\frac{(-1)^{k-1} j!}{(j-k+1)!} M_{j-k+1} &\quad k\leqslant j+1,\\ 0 &\quad\text{otherwise}. \end{cases} \end{align*}$$

Repeated application of that rule proves, in summary, the following contribution to Ulam’s problem about the expected value when n grows large.

Theorem 7.1. Let m be any fixed non-negative integer. Then, subject to the tameness, the uniform tails and the uniform strip hypotheses there holds, as $n\to \infty $ ,

(7.5) $$ \begin{align} {\mathbb E}(L_n) = 2\sqrt{n} +\frac{1}{2} + \sum_{j=0}^m \mu_j n^{1/6 - j/3} + O(n^{-1/6-m/3}), \end{align} $$

where the constants $\mu _j$ are given by

$$\begin{align*}\mu_0 = \int_{-\infty}^\infty t F'(t)\,dt, \qquad \mu_j = \int_{-\infty}^\infty t F_j^*(t)\,dt \quad (j=1,2,\ldots). \end{align*}$$

The first few cases can be expressed in terms of the moments of the Tracy–Widom distribution:

(7.6) $$ \begin{align} \begin{aligned} \mu_0 = M_1, \quad \mu_1 = \frac{1}{60}M_2,\quad \mu_2 = \frac{89}{350} - \frac{1}{1400}M_3,\quad \mu_3 = \frac{538}{7875} M_1+\frac{281}{4536000} M_4; \end{aligned} \end{align} $$

highly accurate numerical values are listed in Table 1.

Table 1 Highly accurate values of $\mu _0,\ldots ,\mu _3$ and $\nu _0,\ldots ,\nu _3$ as computed from (7.6), (7.8) based on values for $M_j$ obtained as in [19, Table 3] (cf. Prähofer’s values for $M_1,\ldots ,M_4$ , published in [73, p. 70]). For the values of $\mu _4,\mu _5$ and $\nu _4,\nu _5$ , see the supplementary material mentioned in Footnote 13

Likewise, by shift and rescale, the variance of $L_n$ can be written in the form

$$\begin{align*}\begin{aligned} {\mathrm 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 t_{l-1/2}(n)^2 \cdot {\mathbb P}(L_n=l) - \Big( {\mathbb E}(L_n) - 2\sqrt{n} - \frac{1}{2}\Big)^2. \end{aligned} \end{align*}$$

By inserting the expansions (7.3), (7.5) and arguing as for Theorem 7.1 we get the following:

Corollary 7.2. Let m be any fixed non-negative integer. Then, subject to the tameness, the uniform tails and the uniform strip hypotheses there holds, as $n\to \infty $ ,

(7.7) $$ \begin{align} {\mathrm Var}(L_n) = \sum_{j=0}^m \nu_j n^{1/3 - j/3} + O(n^{-m/3}), \end{align} $$

with certain constants $\nu _j$ . The first few cases can be expressed in terms of the moments of the Tracy–Widom distribution

(7.8) $$ \begin{align} \begin{aligned} \nu_0 &= -M_1^2 + M_2, \quad \nu_1 = -\frac{67}{60} + \frac{1}{30} \big(-M_1 M_2 + M_3\big),\\ \nu_2 &= -\frac{57}{175} M_1 + \frac{1}{700} M_1 M_3 - \frac{1}{3600}M_2^2 - \frac{29}{25200}M_4,\\ \nu_3 &= -\frac{1076}{7875}M_1^2 -\frac{281}{2268000}M_1M_4 +\frac{893}{7875}M_2 +\frac{1}{42000}M_2M_3 + \frac{227}{2268000} M_5; \end{aligned} \end{align} $$

highly accurate numerical values are listed in Table 1.

The expansions of expected value and variance can be cross-validated by looking at the numerical values for the coefficients $\mu _1,\mu _2,\mu _3$ and $\nu _1,\nu _2,\nu _3$ that we predicted in [19, §4.3]: those values were computed by fitting, in high precision arithmetic, expansions (back then only conjectured) of the form (7.5) with $m=9$ and (7.7) with $m=8$ to the exact tabulated data for $n=500,\ldots ,1000$ . A decision about which digits were to be considered correct was made by comparing the result against a similar computation for $n=600,\ldots ,1000$ . As it turns out, the predictions of [19, §4.3] agree to all the decimal places shown there (that is, to $7$ , $7$ , $6$ and $9$ , $6$ , $4$ places) with the theory-based, highly accurate values given in Table 1.