A Stirling-type formula for the distribution of the length of longest increasing subsequences, applied to finite size corrections to the random matrix limit

1. Introduction

As witnessed by a number of outstanding surveys and monographs (see, e.g., [MR1694204, MR3468920, MR3468738, MR2334203]

), a surprisingly rich topic in combinatorics and probability theory, deeply related to representation theory and to random matrix theory, is the study of the lengths

L_{n} (σ)

of longest increasing subsequences of permutations

σ

on the set

[n] = {1, 2, \dots, n}

and of the behavior of their distribution in the limit

n \to \infty

. Here,

L_{n} (σ)

is defined as the maximum of all

k

for which there are

1 ⩽ i_{1} < i_{2} < \dots < i_{k} ⩽ n

such that

σ_{i_{1}} < σ_{i_{2}} < \dots < σ_{i_{k}}

. Writing permutations in the form

σ = (σ_{1} σ_{2} \dots σ_{n})

we get, e.g.,

L_{9} (σ) = 5

for

σ = (412765893)

, where one of the longest increasing subsequences has been highlighted. Enumeration of the permutations with a given

L_{n}

can be encoded probabilistically: by equipping the symmetric group on

[n]

with the uniform distribution,

L_{n}

becomes a discrete random variable with cumulative probability distribution (CDF)

P (L_{n} ⩽ l)

and probability distribution (PDF)

P (L_{n} = l)

Constructive combinatorics.

Using the Robinson–Schensted correspondence [MR121305], one gets the distribution of $L_{n}$ in the following form, see, e.g., [MR843332, §§3.3–3.7]:

(1)

P (L_{n} ⩽ l) = \frac{1}{n!} \sum λ ⊢ n : l_{λ} ⩽ l d_{λ}^{2} .

Here $λ ⊢ n$ denotes an integer partition $λ_{1} ⩾ λ_{2} ⩾ \dots ⩾ λ_{l_{λ}} > 0$ of $n = \sum_{j = 1}^{l_{λ}} λ_{j}$ and $d_{λ}$ is the number of standard Young tableaux of shape $λ$ , as given by the hook length formula. By generating all partitions of order $n$ , in 1968 Baer and Brock [MR228216] computed tables of $P (L_{n} = l)$ up to $n = 36$ ; in 2000 Odlyzko and Rains [MR1771285] for $n = 15, 30, 60, 90, 120$ (the tables are online, see [OdlyzkoTable]), reporting a computing time for $n = 120$ of about 12 hours (here $p_{n}$ , the number of partitions, is of size $1.8 \times 10^{9}$ ). This quickly becomes infeasible,¹¹1See Remark 3.1 below for a method to compute the exact rational values of the distribution of $L_{n}$ based on random matrix theory, which has recently been used to tabulate $P (L_{n} = l)$ , $1 ⩽ l ⩽ n$ , up to $n = 1000$ . as $p_{n}$ is already as large as $2.3 \cdot 10^{14}$ for $n = 250$ . Another use of the combinatorial methods is the approximation of the distribution of $L_{n}$ by Monte Carlo simulations [MR228216, MR1771285]: one samples random permutations $σ$ and calculates $L_{n} (σ)$ by the Robinson–Schensted correspondence.²²2In Mathematica, a single trial is generated by the command

\tt Length@LongestOrderedSequence[PermutationList[RandomPermutation[n],% n]].

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

Analytic combinatorics and the random matrix limit.

For analytic methods of enumeration the starting points is a more or less explicit representation of a suitable generating function; here, the suitable one turns out to be the exponential generating function of the CDF $P (L_{n} ⩽ l)$ , when considered as a sequence of $n$ with the length $l$ fixed:

f_{l} (z) := \infty \sum n = 0 \frac{P (L_{n} ⩽ l)}{n!} z^{n} (z \in C, l \in N) .

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

(2)

f_{l} (z^{2}) = D_{l} (z), D_{l} (z) := l det j, k = 1 I_{j - k} (2 z),

in terms of a Toeplitz determinant of the modified Bessel functions $I_{m}$ , $m \in Z$ , which are entire functions of exponential type themselves. By relating, first, the Toeplitz determinant to the machinery of Riemann–Hilbert problems to study a double-scaling limit of the generating function and by using, next, a Tauberian theorem³³3The Tauberian part (“de-poissonization”) makes it hard to get more than the leading order of the asymptotics. to induce from that limit an asymptotics of the coefficients, Baik, Deift and Johansson [MR1682248] succeeded 1999 in establishing⁴⁴4In the first place, [MR1682248, Thm. 1.1] states the following limit to hold pointwise in $t \in R$ :

lim n \to \infty P (\frac{L_{n} - 2 \sqrt{n}}{n^{1 / 6}} ⩽ t) = F_{2} (t) .

However, since the limit distribution

F_{2}

is continuous, by a standard Tauberian follow-up [MR1652247, Lemma 2.1] of the Portmanteau theorem in probability theory, this convergence is known to hold, in fact, uniformly in

t

(3)

P (L_{n} ⩽ l) = F_{2} (\frac{l - 2 \sqrt{n}}{n^{1 / 6}}) + o (1) (n \to \infty),

uniformly in $l \in N$ . Here, $F_{2} (t)$

denotes the Tracy–Widom distribution for

β = 2

, that is, the probability that in the soft-edge scaling limit of the Gaussian unitary ensemble (GUE) the largest eigenvalue is bounded from above by

t

. This distribution can be evaluated numerically based on its representation either in terms of the Airy kernel determinant [MR1236195] or in terms of the Painlevé II transcendent [MR1257246]; see Remark 3.2 and [MR2895091] for details.

As impressive as the use of the limit (3) might look as a numerical approximation to the distribution of $L_{n}$ near its mode for larger $n$ , cf. Fig. 1, there are two notable deficiencies: first, since the error term in (3) is additive (i.e., w.r.t. absolute scale), the approximation is rather poor for $l ≪ 2 \sqrt{n}$ ; second, the convergence rate is rather slow, in fact conjectured to be just of the order $O (n^{- 1 / 3})$ , see [arxiv.2205.05257] and the discussion below. Both deficiencies are well illustrated in Table 1 for $n = 20$ and in Fig. 2 for $n = 1000$ .

$l$	$P (L_{20} ⩽ l)$	Stirling-type (5)	Monte-Carlo	random matrix (3)
1	$4.110 \cdot 10^{- 19}$	$4.119 \cdot 10^{- 19}$	$0.000$	$6.282 \cdot 10^{- 5}$
2	$2.698 \cdot 10^{- 9}$	$2.703 \cdot 10^{- 9}$	$0.000$	$1.422 \cdot 10^{- 3}$
3	$6.698 \cdot 10^{- 5}$	$6.710 \cdot 10^{- 5}$	$7.240 \cdot 10^{- 5}$	$1.485 \cdot 10^{- 2}$
4	$1.090 \cdot 10^{- 2}$	$1.092 \cdot 10^{- 2}$	$1.089 \cdot 10^{- 2}$	$8.014 \cdot 10^{- 2}$
5	$1.427 \cdot 10^{- 1}$	$1.429 \cdot 10^{- 1}$	$1.428 \cdot 10^{- 1}$	$2.503 \cdot 10^{- 1}$
6	$4.841 \cdot 10^{- 1}$	$4.846 \cdot 10^{- 1}$	$4.838 \cdot 10^{- 1}$	$5.079 \cdot 10^{- 1}$
7	$8.042 \cdot 10^{- 1}$	$8.064 \cdot 10^{- 1}$	$8.040 \cdot 10^{- 1}$	$7.513 \cdot 10^{- 1}$
8	$9.521 \cdot 10^{- 1}$	$9.581 \cdot 10^{- 1}$	$9.519 \cdot 10^{- 1}$	$9.041 \cdot 10^{- 1}$
9	$9.921 \cdot 10^{- 1}$	$9.996 \cdot 10^{- 1}$	$9.920 \cdot 10^{- 1}$	$9.716 \cdot 10^{- 1}$

Table 1. The values of

P (L_{n} ⩽ l)

and various of its approximations for

n = 20

. The Monte Carlo simulation was run with

T = 5 \cdot 10^{6}

samples. Clearly observable is a multiplicative (i.e., relative) error of already less than

1 %

for the Stirling-type formula (5) as well as an additive (i.e., absolute) error of order

4 \cdot 10^{- 4} \approx T^{- 1 / 2}

for Monte Carlo and an additive error of order

10^{- 1} \approx n^{- 1 / 3}

for the limit (3) from random matrix theory.

A Stirling-type formula

In this paper we suggest a different type of numerical approximation to the distribution of $L_{n}$ that enjoys the following advantages: (a) it has a small multiplicative (i.e., relative) error, (b) it has faster convergence rates, apparently even faster than (3) with its first finite size correction term added, and (c) it is much faster to compute than Monte Carlo simulations. In fact, the distribution of $L_{n}$ for $n = 10^{5}$ , as shown in Fig. 1, exhibits an estimated maximum additive error of less than $10^{- 5}$ and took just about five seconds to compute; whereas Forrester and Mays [arxiv.2205.05257] have recently reported a computing time of about 14 hours to generate $T = 5 \cdot 10^{6}$ Monte Carlo trials for this $n$ ; the error of such a simulation is expected to be of the order $1 / \sqrt{T} \approx 4 \cdot 10^{- 4}$ .

Specifically, we use Hayman’s generalization [Hayman56] of Stirling’s formula for $H$ -admissible functions; for expositions see [MR2483235, MR1373678, Wong89]. For simplicity, as is the case here for $f (z) = f_{l} (z)$ , assume that

f (z) = \infty \sum n = 0 a_{n} z^{n} (z \in C)

is an entire function with positive coefficients $a_{n}$ and consider the auxiliary functions

a (r) = r \frac{d}{d r} log f (r), b (r) = r \frac{d}{d r} a (r) (r > 0) .

If $f$ is $H$ -admissible, then for each $n \in N$ the equation $a (r_{n}) = n$ has a unique solution $r_{n} > 0$ such that $b (r_{n}) > 0$ and the following generalization of Stirling’s formula⁵⁵5A generalization of Stirling’s classical formula, indeed: for the $H$ -admissible function $f (z) = e^{z}$ , cf. Thm. 2.1 Criterion II.g, we have $a_{n} = 1 / n!$ , $b (r) = a (r) = r$ , $r_{n} = n$ and (4) specifies to

\frac{1}{n!} = \frac{e^{n}}{n^{n} \sqrt{2 π n}} (1 + o (1)) (n \to \infty) .

As in Table 1, the error is already below

1 %

for

n

as small as

n = 20

. holds true:

(4)

a_{n} = \frac{f (r_{n})}{r_{n}^{n} \sqrt{2 π b (r_{n})}} (1 + o (1)) (n \to \infty) .

We observe that the error is multiplicative here. In Thm. 2.2 we will prove, using some theory of entire functions, the $H$ -admissibility of the generating functions $f_{l}$ . Hence the Stirling-type formula (4) applies without further ado to their coefficients $P (L_{n} ⩽ l) / n!$ . Since the error is multiplicative, nothing changes if we multiply the approximation by $n!$ and we get

(5)

P (L_{n} ⩽ l) = \frac{n! \cdot f_{l} (r_{l, n})}{r_{l, n}^{n} \sqrt{2 π b_{l} (r_{l, n})}} (1 + o (1)) (n \to \infty),

where we have labeled all quantities when applied to $f = f_{l}$ by an additional index $l$ . An approximation to $P (L_{n} = l)$ is then obtained simply by taking differences. The power of these approximations, if used as a numerical tool even for $n$ as small as $n = 20$ , is illustrated in Table 1 and Figs. 1–3, as well as in Table 2 below.

Figure 2. Discrete distribution of $L_{n}$ for $n = 1000$ near its mode vs. the random matrix limit given by the leading order terms in (8) and (33a) (solid red line). Left: CDF $P (L_{n} ⩽ l)$ ; right: PDF $P (L_{n} = l)$ . The exact values of the distribution of $L_{n}$ and their approximation by the Stirling-type formula (5) differ just by additive errors of the order $10^{- 4}$ , see Fig. 3, which is well below plotting accuracy.

Numerical evaluation of the generating function

For the Stirling-type formula (5) to be easily accessible in practice, we require an expression for $f_{l} (r)$ that can be numerically evaluated, for $r > 0$ , in a stable, accurate, and efficient fashion. Since the direct evaluation of the Toeplitz determinant (2) is numerically highly unstable, and has a rather unfavorable complexity of $O (l^{3})$ for larger $l$ , we look for alternative representations. One option—used in [MR2754188] to numerically extract $P (L_{n} ⩽ l)$ from $f_{l} (z)$ by Cauchy integrals over circles in the complex plane that are centered at the origin with the same radius $r_{l, n}$ as in (5)—is the machinery, cf., e.g., [MR1935759], to transform Toeplitz determinants into Fredholm determinants which are then amenable for the numerical method developed in [MR2600548]. However, since we need the values of the generating function $f_{l} (r)$ for real $r > 0$ only, there is a much more efficient option, which comes from yet another connection to random matrix theory.

To establish this connection we first note that an exponentially generating function $f (r)$ of a sequence of probability distributions has a probabilistic meaning if multiplied by $e^{- r}$ : a process called poissonization. Namely, if the draws from the permutation groups are independent for all orders and if we take $N \in N_{0} := {0, 1, 2, 3, \dots}$

to be a further independent random variable with Poisson distribution of intensity

r ⩾ 0

, we see that

P (L_{N} ⩽ l) = \infty \sum n = 0 \frac{e^{- r} r^{n}}{n!} P (L_{n} ⩽ t) = e^{- r} f_{l} (r)

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

(6)

E_{2}^{(hard)} (0; [0, 4 r], l) = e^{- r} f_{l} (r) .

Here, $E_{2}^{(hard)} (0; [0, t], l)$ denotes the probability that, in the hard-edge scaling limit of the Laguerre unitary ensemble (LUE) with parameter $l$ , the smallest eigenvalue is bounded from below by $t ⩾ 0$ . Now, the point here is that this distribution can be evaluated numerically, stable and accurate with a complexity that is largely independent of $l$ , based on two alternative representations: either in terms of the Bessel kernel determinant [MR1236195] or in terms of the Jimbo–Miwa–Okamoto $σ$ -form of the Painlevé III ( $σ$ -PIII) transcendent [MR1266485]; see [MR2895091] for details. We will show in Sect. 3 that the auxiliary functions $a_{l} (r)$ and $b_{l} (r)$ fit into both frameworks, too.

Figure 3. Maximum absolute (i.e., additive) errors of various approximations to (left panel) the CDF $P (L_{n} ⩽ l)$ and (right panel) the PDF $P (L_{n} = l)$ in double logarithmic scaling of the axes, based on tabulated exact values up to $n = 1000$ , cf. Remark 3.1; solid lines are fits of the form $c n^{- α}$ for $n ⩾ 100$ . Red $+$ : error of the leading order terms (random matrix limit) in (8), (33a), with $α = \frac{1}{3}$ (CDF) and $α = \frac{1}{3} + \frac{1}{6} = \frac{1}{2}$ (PDF). Green $\circ$ : error of the expansions (8), (33a) truncated after the first finite size correction term, with $α = \frac{2}{3}$ (CDF) and $α = \frac{2}{3} + \frac{1}{6} = \frac{5}{6}$ (PDF), where $F_{2, 1}$ has been approximated as in Fig. 4. Blue $∙$ : error of the Stirling-type formula (5), with $α = \frac{3}{4}$ (CDF) and $α = \frac{3}{4} + \frac{1}{6} = \frac{11}{12}$ (PDF).

Finite size corrections to the random matrix limit

A double logarithmic plot of the additive errors (taking the maximum w.r.t. $l \in {1, 2, \dots, n}$ ) in approximating the distribution $P (L_{n} ⩽ l)$ by either the random matrix limit (3) or by the Stirling-type formula (5) exhibits nearly straight lines, see Fig. 3 for $n$ chosen between $10$ and $1000$ . Fitting the slopes to simple rationals strongly suggests that, as $n \to \infty$ , uniformly in $l \in {1, 2, \dots, n}$ ,


(7a)	$P (L_{n} ⩽ l)$	$= F_{2} (\frac{l - 2 \sqrt{n}}{n^{1 / 6}}) + O (n^{- 1 / 3}),$
(7b)	$P (L_{n} ⩽ l)$	$= \frac{n! \cdot f_{l} (r_{l, n})}{r_{l, n}^{n} \sqrt{2 π b_{l} (r_{l, n})}} + O (n^{- 3 / 4}) .$

The approximation order (and the size of the implied constant) in (7b) is much better than the one in (7a) so that the Stirling-type formula can be used to reveal the structure of the $O (n^{- 1 / 3})$ term in the random matrix limit. In fact, as the error plot in Fig. 3 suggests and we will more carefully argue in Sect. 4.1, this can even be iterated yet another step and we are led to the specific conjecture⁶⁶6Note that $l \in N$ is always a discrete variable in this paper, so there is no need for taking integer parts here.

(8)

P (L_{n} ⩽ l) = F_{2} (t_{l}) + n^{- 1 / 3} F_{2, 1} (t_{l}) + n^{- 2 / 3} F_{2, 2} (t_{l}) + O (n^{- 1}), t_{l} := \frac{l - 2 \sqrt{n}}{n^{1 / 6}},

as $n \to \infty$ , uniformly in $l \in N$ . Compelling evidence for the existence of the functions $F_{2, 1}$ and $F_{2, 2}$ is given in the left panels of Figs. 4 and 6.⁷⁷7Based on Monte-Carlo simulations, Forrester and Mays [arxiv.2205.05257, Eq. (1.10) and Fig. 7] were recently also led to conjecture an expansion of the form

P (L_{n} ⩽ l) = F_{2} (t_{l}) + n^{- 1 / 3} F_{2, 1} (t_{l}) + \dots .

However, the substantially larger errors of Monte Carlo simulations as compared to the Stirling-type formula (5) would inhibit them from getting, in reasonable time, much evidence about the next finite size correction term.

Figure 4. Rescaled differences between the distributions of $L_{n}$ and their expansions truncated after the leading order term (i.e., the random matrix limit), see (8) for the CDF resp. (33a) for the PDF; data points (to avoid clutter just every $5^{th}$ is displayed) have been calculated using the Stirling-type formula (5) for $n = 10^{6}$ (red $+$ ), $n = 10^{8}$ (green $\circ$ ), $n = 10^{10}$ (blue $∙$ ). Left: CDF errors rescaled by $n^{1 / 3}$ , horizontal axis is $t = (l - 2 \sqrt{n}) / n^{1 / 6}$ , cf. [arxiv.2205.05257, Fig. 7] for a similar figure with data from Monte Carlo simulations for $n = 2 \cdot 10^{4}$ and $n = 10^{5}$ . The solid line is a polynomial ${~ F}_{2, 1} (t)$ of degree $64$ fitted to the $836$ data points for $n = 10^{10}$ with $- 8 ⩽ t ⩽ 10$ ; it approximates $F_{2, 1} (t)$ in that interval. Right: PDF errors rescaled by $n^{1 / 2}$ , horizontal axis is $t = (l - \frac{1}{2} - 2 \sqrt{n}) / n^{1 / 6}$ . The solid line displays the function ${~ F}_{2, 1}^{'} (t) + F_{2}^{'''} (t) / 24$ as an approximation of $F_{2, 1}^{'} (t) + F_{2}^{'''} (t) / 24$ , with the polynomial ${~ F}_{2, 1} (t)$ taken from the left panel. The dotted line shows the term ${~ F}_{2, 1}^{'} (t)$ only.

Additionally, we will argue in Sect. 4.3 that (8) allows us to derive an expansion of the expected value of $L_{n}$ , specifically

(9)		$E (L_{n}) = 2 \sqrt{n} + μ_{0} n^{1 / 6} + \frac{1}{2} + μ_{1} n^{- 1 / 6} + μ_{2} n^{- 1 / 2} + O (n^{- 5 / 6}),$
	$μ_{0} = \int_{- \infty}^{\infty} t F_{2}^{'} (t) d t = - 1.77108 68074 \dots,$
	$μ_{1} = \int_{- \infty}^{\infty} t F_{2, 1}^{'} (t) d t = 0.06583 238 \dots, μ_{2} = \int_{- \infty}^{\infty} t F_{2, 2}^{'} (t) d t = 0.2 612227 \dots .$

Similarly, we will derive in Sect. 4.4 an expansion of the variance of $L_{n}$ of the form

(10)		$V a r (L_{n})$	$= ν_{0} n^{1 / 3} + ν_{1} + ν_{2} n^{- 1 / 3} + O (n^{- 2 / 3}),$
	$ν_{0}$	$= \int_{- \infty}^{\infty} t^{2} F_{2}^{'} (t) d t - μ_{0}^{2} = 0. 8131947928 \dots,$
	$ν_{1}$	$= \int_{- \infty}^{\infty} t^{2} F_{2, 1}^{'} (t) d t - 2 μ_{0} μ_{1} = - 1.20720 507 \dots,$
	$ν_{2}$	$= \int_{- \infty}^{\infty} t^{2} F_{2, 2}^{'} (t) d t - μ_{1}^{2} - 2 μ_{0} μ_{2} = 0.56715 6 \dots .$

The values of $μ_{0}$ and $ν_{0}$ are the known values of mean and variance of the Tracy–Widom distribution $F_{2}$ , cf. [MR2895091, Table 10]. The leading parts of (9) up to $μ_{0} n^{1 / 6}$ and of (10) up to $ν_{0} n^{1 / 3}$ had been established previously by Baik, Deift and Johansson [MR1682248, Thm. 1.2].

2. $H$ -Admissibility of the Generating Function and its Implications

2.1. $H$ -admissible functions

For simplicity we restrict ourselves to entire functions. We refrain from displaying the rather lengthy technical definition of $H$ -admissibility,⁸⁸8Since we consider entire functions only, $H$ -admissibility is here understood to hold in all of $C$ . which is difficult to be verified in practice and therefore seldomly directly used. Instead, we start by collecting some usefuls facts and criteria from Hayman’s original paper [Hayman56]:⁹⁹9Interestingly, the powerful criterion in part III (which is [Hayman56, Thm. XI]) is missing from the otherwise excellent expositions [MR2483235, MR1373678, Wong89] of $H$ -admissibility.

Theorem 2.1 (Hayman 1956).

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

I. If $f$ is $H$ -admissible, then:

$f (r) > 0$ for large $r > 0$ , so that in particular the auxiliary functions¹⁰¹⁰10In terms of differential operators we have $r \frac{d}{d r} = \frac{d}{d log r}$ .

(11) $a (r) = r \frac{d}{d r} log f (r), b (r) = r \frac{d}{d r} a (r)$

are well defined there;
for large $r > 0$ there is $log f (r)$ strictly convex in $log r$ , $a (r)$ strictly monotonically increasing, and $b (r) > 0$ such that $a (r), b (r) \to \infty$ as $r \to \infty$ ; in particular, for large integers $n$ there is a unique $r_{n} > 0$ that solves $a (r_{n}) = n$ , it is $r_{n} \to \infty$ as $n \to \infty$ ;
if the coefficients $a_{n}$ of $f$ are all positive, then I.b holds for all $r > 0$ ;
as $r \to \infty$ , uniformly in $n \in N_{0}$ ,

(12) $\frac{a_{n} r^{n}}{f (r)} = \frac{1}{\sqrt{2 π b (r)}} (exp (- \frac{(n - a (r))^{2}}{2 b (r)}) + o (1)) .$

II. If $f$ and $g$ are $H$ -admissible, then:

$f (z) g (z)$ , $e^{f (z)}$ and $f (z) + p (z)$ are $H$ -admissible;
if the leading coefficient of $p$ is positive, $f (z) p (z)$ and $p (f (z))$ are $H$ -admissible;
if the Taylor coefficients of $e^{p (z)}$ are eventually positive, $e^{p (z)}$ is $H$ -admissible.

III. If $f$ has genus zero¹¹¹¹11By definition, an entire function $f$ has genus zero if it is a polynomial or if it can be represented as a convergent infinite product of the form

f (z) = c z^{m} \infty \prod n = 1 (1 - \frac{z}{z_{n}}) (z \in C),

where

c \in C

is a constant,

m \in N_{0}

is the order of the zero at

z = 0

, and

z_{1}, z_{2}, \dots

is the sequence of the non-zero zeros, where each one is listed as often as multiplicity requires. with, for some

δ > 0

, at most finitely many zeros in the sector

| arg z | ⩽ \frac{π}{2} + δ

and satisfies I.a such that

b (r) \to \infty

r \to \infty

, then

f

H

-admissible.

Obviously, the Stirling-type formula (4) is obtained from the approximation result (12) by just inserting the particular choice $r = r_{n}$ .

Remark 2.1.

We observe that (12) has an interesting probabilistic content:¹²¹²12Note that the particular case $f (z) = e^{z}$ (which is $H$ -admissible by Thm. 2.1

.II.g) specifies to the well-known normal approximation of the Poisson distribution for large intensities—which is a simple consequence of the central limit theorem if we observe that the sum of

k

independent Poisson random variates of intensity

ρ

is one of intensity

r = ρ k

. as a distribution in the discrete variable

n \in N_{0}

, the Boltzmann¹³¹³13We follow the terminology in the theory of Boltzmann samplers [MR2095975], a framework for the random generation of combinatorial structures. probabilities

a_{n} r^{n} / f (r)

associated with an

H

-admissible entire function

f

are, for large intensities

r > 0

, approximately normal with mean

a (r)

and variance

b (r)

; see the right panel of Fig. 5 for an illustrative example using the generating function

f_{5} (z)

. The additional freedom that is provided in the normal approximation (12) by the uniformity w.r.t.

n

will be put to good use in Sect. 2.3.

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

The classification of entire functions (by quantities such as genus, order, type, etc.) and their distribution of zeros is deeply related to the analysis of their essential singularity at $z = \infty$ . For the purposes of this paper, the following simple criterion is actually all we need. The proof uses some theory of entire functions, which can be found, e.g., in [MR589888].

Lemma 2.1.

Let $f (z)$ be an entire function of exponential type with positive Maclaurin coefficients. If there are constants $c, τ, ν > 0$ such that there holds, for the principal branch of the power function and for each $0 < δ ⩽ \frac{π}{2}$ , the asymptotic expansion¹⁴¹⁴14See Remark A.1 for the uniformity implied with the notation $O (z^{- 1})$ as $z \to \infty$ while $| arg z | ⩽ \frac{π}{2} - δ$ .

(13)

D (z) := f (z^{2}) = c z^{- ν} e^{τ z} (1 + O (z^{- 1})) (z \to \infty, | arg z | ⩽ \frac{π}{2} - δ),

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

(14)

a (r) = \frac{τ}{2} r^{1 / 2} - \frac{ν}{2} + O (r^{- 1 / 2}), b (r) = \frac{τ}{4} r^{1 / 2} + O (r^{- 1 / 2}) .

Proof.

The expansion (13) is equivalent to

(15)

f (z) = c z^{- ν / 2} e^{τ z^{1 / 2}} (1 + O (z^{- 1 / 2})) (z \to \infty, | arg z | ⩽ π - 2 δ),

which readily implies:

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

$\operator@fontlimsupr→∞log|f(reiθ|r1/2=τcos(θ/2)(−π<θ<π),$

so that $f$ has order $\frac{1}{2}$ and type $τ$ , hence genus zero;
for sufficiently large $R = R_{δ} > 0$ , there are no zeros $z$ of $f$ with

$| z | ⩾ R, | arg z | ⩽ π - 2 δ .$

Since the Maclaurin coefficients of $f$ are positive we have $f (r) > 0$ for $r > 0$ and the auxiliary functions $a (r)$ , $b (r)$ in (11) are well-defined for $r > 0$ . In fact, both functions can be analytically continued into the domain of uniformity of the expansion (15) and by differentiating this expansion (which is, because of analyticity, legitimate by a theorem of Ritt, cf. [MR0435697]) we obtain (14); this implies, in particular, $b (r) \to \infty$ as $r \to \infty$ . Thus, all the assumptions of Thm. 2.1.III are satisfied and $f$ is shown to be $H$ -admissible. ∎

2.2. Singularity analysis of the generating function at $z = \infty$

Establishing an expansion of the form (13) for $D_{l} (z) = f_{l} (z^{2})$ as given by (2), that is to say, for the Toeplitz determinant

(16)

D_{l} (z) = l det j, k = 1 I_{j - k} (2 z),

suggests to start with the expansions (valid for all $0 < δ ⩽ \frac{π}{2}$ , see [MR0435697, p. 251])


(17a)	$I_{m} (z)$	$\sim \frac{e^{z}}{(2 π z)^{1 / 2}} \infty \sum n = 0 (- 1)^{n} \frac{A_{n} (m)}{z^{n}} (z \to \infty, \| arg z \| ⩽ \frac{π}{2} - δ),$
(17b)	$A_{n} (m)$	$= \frac{(4 m^{2} - 1^{2}) (4 m^{2} - 3^{2}) \dots (4 m^{2} - (2 n - 1)^{2})}{n! 8^{n}} .$

This does not yield (13) at once, as there could be, however unlikely it would be, eventually a catastrophic cancellation of all of the expansion terms when being inserted into the determinant expression defining $D_{l} (z)$ . For the specific cases $l = 1, 2, \dots, 8$ a computer algebra system shows that exactly the first $l - 1$ terms of the expansion (17) mutually cancel each other in forming the determinant, and we get by this approach¹⁵¹⁵15Odlyzko [MR1373678, Ex. 10.9] reports that he and Wilf had used this approach, before 1995, for small $l$ in the framework of the method of “subtraction of singularities” in asymptotic enumeration. He states the expansions for $D_{4} (z)$ and $D_{5} (z)$ , cf. [MR1373678, Eqs. (10.30)/(10.39)], with a misrepresented constant factor in $D_{5} (z)$ , though. No attempt, however, was made back then to guess the general form. the expansions

$D_{1} (z)$	$= \frac{e^{2 z}}{2 π^{1 / 2} z^{1 / 2}} (1 + O (z^{- 1})),$	$D_{2} (z)$	$= \frac{e^{4 z}}{8 π z^{2}} (1 + O (z^{- 1})),$
$D_{3} (z)$	$= \frac{e^{6 z}}{32 π^{3 / 2} z^{9 / 2}} (1 + O (z^{- 1})),$	$D_{4} (z)$	$= \frac{3 e^{8 z}}{256 π^{2} z^{8}} (1 + O (z^{- 1})),$
$D_{5} (z)$	$= \frac{9 e^{10 z}}{1024 π^{5 / 2} z^{25 / 2}} (1 + O (z^{- 1})),$	$D_{6} (z)$	$= \frac{135 e^{12 z}}{8192 π^{3} z^{18}} (1 + O (z^{- 1})),$
$D_{7} (z)$	$= \frac{6075 e^{14 z}}{65536 π^{7 / 2} z^{49 / 2}} (1 + O (z^{- 1})),$	$D_{8} (z)$	$= \frac{1913625 e^{16 z}}{1048576 π^{4} z^{32}} (1 + O (z^{- 1})) .$

All of them, inherited from (17), are valid as $z \to \infty$ while $| arg z | ⩽ \frac{π}{2} - δ$ with the uniformity content implied by the symbol $O (z^{- 1})$ . From these instances, in view of (17) and the multilinearity of the determinant, we guess that

D_{l} (z) = c_{l} \frac{e^{2 l z}}{(4 π z)^{l / 2} (2 z)^{l (l - 1) / 2}} (1 + O (z^{- 1})) (z \to \infty, | arg z | ⩽ \frac{π}{2} - δ)

and observe

c_{1}, c_{2}, c_{3}, c_{4}, c_{5}, c_{6}, c_{7}, c_{8}, \dots = 1, 1, 2, 12, 288, 34560, 24883200, 125411328000, \dots .

Consulting the OEIS¹⁶¹⁶16https://oeis.org/A000178 (On-Line Encyclopedia of Integer Sequences) suggests the coefficients to be generally of the form

c_{l} = 0! \cdot 1! \cdot 2! \dots (l - 1)! .

Though this is very likely to hold for all $l \in N$ —a fact that would at once yield the $H$ -admissiblity of all the generating function $f_{l}$ by Lemma 2.1—a proof seems to be elusive along these lines, but see Remark 2.3 for a remedy.

Inspired by the fact that the one-dimensional Laplace’s method easily gives the leading order term in (17) when applied to the Fourier representation

(18)

I_{m} (2 z) = \frac{1}{2 π} \int_{- π}^{π} e^{2 z cos θ} e^{- i m θ} d θ (z \in C, m \in Z),

we represent the Toeplitz determinant $D_{l} (z)$ in terms of a multidimensional integral and study the limit $z \to \infty$ by the multidimensional Laplace method discussed in the Appendix. In fact, (18) shows that the symbol of the Toeplitz determinant $D_{l} (z)$ is $exp (2 z cos θ)$ and a classical formula of Szegő’s [Sz1915, p. 493] from 1915, thus gives, without further calculation, the integral representation¹⁷¹⁷17This induces, see (6), an integral representation of the distribution $E_{2}^{(hard)} (0; [0, s], l)$ which has been derived in 1994 by Forrester [MR1271945] using generalized hypergeometric functions defined in terms of Jack polynomials.

(19)

D_{l} (z) = \frac{1}{(2 π)^{l} l!} \int_{- π}^{π} \dots \int_{- π}^{π} e^{2 z \sum_{j = 1}^{l} cos θ_{j}} \cdot ∣ ∣ Δ (e^{i θ_{1}}, \dots, e^{i θ_{l}}) {∣ ∣}^{2} d θ_{1} \dots d θ_{l},

where

Δ (w_{1}, \dots, w_{l}) := \prod j > k (w_{j} - w_{k})

denotes the Vandermonde determinant of the complex numbers $w_{1}, \dots, w_{l}$ .

Remark 2.2.

By Weyl’s integration formula on the unitary group $U (l)$ , cf. [MR2105088, Eq. 1.5.89], the integral (19) can be written as

Dn(z)=EU∈U(l)ez\operator@fonttr(U+U∗),

where the expectation $E$ is taken with respect to the Haar measure. Without any reference to (2), this form was derived in 1998 by Rains [Rains98, Cor. 4.1] directly from the identity

P(Ln⩽l)=EU∈U(l)(|\operator@fonttrU|2n),

which he had obtained most elegantly from the representation theory of the symmetric group.

We are now able to prove our main theorem.

Theorem 2.2.

For each $0 < δ ⩽ π / 2$ and $l \in N$ there holds the asymptotic expansion

D_{l} (z) = f_{l} (z^{2}) = \frac{0! \cdot 1! \cdot 2! \dots (l - 1)! \cdot e^{2 l z}}{(2 π)^{l / 2} (2 z)^{l^{2} / 2}} (1 + O (z^{- 1})) (z \to \infty, | arg z | ⩽ \frac{π}{2} - δ) .

Thus, by Lemma 2.1, the generating functions $f_{l} (z)$ are $H$ -admissible and their auxiliary functions satisfy, as $r \to \infty$ ,

(20)

a_{l} (r) = l r^{1 / 2} - \frac{1}{4} l^{2} + O (r^{- 1 / 2}), b_{l} (r) = \frac{1}{2} l r^{1 / 2} + O (r^{- 1 / 2}) .

Proof.

We write (19) in the form

e^{- l z} D_{l} (z / 2) = \frac{1}{(2 π)^{l} l!} \int_{- π}^{π} \dots \int_{- π}^{π} e^{- z \sum_{j = 1}^{l} (1 - cos θ_{j})} \cdot ∣ ∣ Δ (e^{i θ_{1}}, \dots, e^{i θ_{l}}) {∣ ∣}^{2} d θ_{1} \dots d θ_{l} .

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

S (θ_{1}, \dots, θ_{l}) := l \sum j = 1 (1 - cos θ_{j}),

takes it minimum $S (θ_{*}) = 0$ at $θ_{*} = 0$ with the expansion $S (θ_{1}, \dots, θ_{l}) = \frac{1}{2} θ^{T} θ + O (| θ |^{4})$ as $θ \to 0$ , where $| \cdot |$ denotes Euclidian length. Likewise we get for the non-exponential factor¹⁸¹⁸18This factor is zero at $θ_{*} = 0$ , so that Hsu’s variant (39) of Laplace’s method, which is the one predominantly found in the literature, does not yield the leading term of $D_{l} (z)$ . That we have to expand the non-exponential factor up to degree $l (l - 1)$ for the first non-zero contribution to show up, corresponds to the mutual cancellation of the leading terms of (17) when being inserted into the Toeplitz determinant that defines $D_{l} (z)$ .

	$∣ ∣ Δ (e^{i θ_{1}}, \dots, e^{i θ_{l}}) {∣ ∣}^{2}$	$= \prod j > k \| e^{i θ_{j}} - e^{i θ_{j}} \|^{2} = \prod j > k {∣ ∣ i θ_{j} - i θ_{j} + O (\| θ \|^{2}) ∣ ∣}_{j}^{2}$
		$= Δ (θ_{1}, \dots, θ_{l})^{2} + O (\| θ \|^{l (l - 1) + 1}) (θ \to 0),$

where the degree of the homogeneous polynomial $Δ (θ_{1}, \dots, θ_{l})^{2}$ is $l (l - 1)$ .

Therefore, by the multidimensional Laplace method as given in Corollary A.1, see also formula (45) following it, we obtain immediately

e^{- l z} D_{l} (z / 2) = c_{l} \frac{(2 π)^{l / 2} z^{- \frac{l + l (l - 1)}{2}}}{(2 π)^{l}} (1 + O (z^{- 1})) (z \to \infty, | arg z | ⩽ \frac{π}{2} - δ)

with

(21)

c_{l} := \frac{1}{l!} \frac{1}{(2 π)^{l / 2}} \int_{R^{l}} e^{- θ^{T} θ / 2} \cdot ∣ ∣ Δ (θ_{1}, \dots, θ_{l}) {∣ ∣}^{2} d θ = 0! \cdot 1! \cdot 2! \dots (l - 1)!,

where the evaluation of this multiple integral is well-known in random matrix theory, e.g., as a consequence of Selberg’s integral formula, cf. [MR2760897, Eq. (2.5.11)]. ∎

Remark 2.3.

By (6), Thm. 2.2 implies

E_{2}^{(hard)} (0; [0, s], l) = \frac{0! \cdot 1! \cdot 2! \dots (l - 1)!}{(2 π)^{l / 2}} \cdot s^{- l^{2} / 4} e^{- s / 4 + l s^{1 / 2}} (

A Stirling-type formula for the distribution of the length of longest increasing subsequences, applied to finite size corrections to the random matrix limit

On the convergence rate of some nonlocal energies

Convergence Rate of Empirical Spectral Distribution of Random Matrices from Linear Codes

Convergence Analysis of the Randomized Newton Method with Determinantal Sampling

The volume-of-tube method for Gaussian random fields with inhomogeneous variance

Computation of the Expected Euler Characteristic for the Largest Eigenvalue of a Real Non-central Wishart Matrix

Economic Design of Memory-Type Control Charts: The Fallacy of the Formula Proposed by Lorenzen and Vance (1986)

Criteria for the numerical constant recognition

1. Introduction

Constructive combinatorics.

Analytic combinatorics and the random matrix limit.

A Stirling-type formula

Numerical evaluation of the generating function

Finite size corrections to the random matrix limit

2. $H$ -Admissibility of the Generating Function and its Implications

2.1. $H$ -admissible functions

Theorem 2.1 (Hayman 1956).

Remark 2.1.

Lemma 2.1.

Proof.

2.2. Singularity analysis of the generating function at $z = \infty$

Remark 2.2.

Theorem 2.2.

Proof.

Remark 2.3.

A Stirling-type formula for the distribution of the length of longest increasing subsequences, applied to finite size corrections to the random matrix limit

Related Research

On the convergence rate of some nonlocal energies

Convergence Rate of Empirical Spectral Distribution of Random Matrices from Linear Codes

Convergence Analysis of the Randomized Newton Method with Determinantal Sampling

The volume-of-tube method for Gaussian random fields with inhomogeneous variance

Computation of the Expected Euler Characteristic for the Largest Eigenvalue of a Real Non-central Wishart Matrix

Economic Design of Memory-Type Control Charts: The Fallacy of the Formula Proposed by Lorenzen and Vance (1986)

Criteria for the numerical constant recognition

1. Introduction

Constructive combinatorics.

Analytic combinatorics and the random matrix limit.

A Stirling-type formula

Numerical evaluation of the generating function

Finite size corrections to the random matrix limit

2. H-Admissibility of the Generating Function and its Implications

2.1. H-admissible functions

Theorem 2.1 (Hayman 1956).

Remark 2.1.

Lemma 2.1.

Proof.

2.2. Singularity analysis of the generating function at z=∞

Remark 2.2.

Theorem 2.2.

Proof.

Remark 2.3.

2. $H$ -Admissibility of the Generating Function and its Implications

2.1. $H$ -admissible functions

2.2. Singularity analysis of the generating function at $z = \infty$