Finite size corrections at the hard edge for the Laguerre \beta ensemble

# Finite size corrections at the hard edge for the Laguerre β ensemble

## Abstract.

A fundamental question in random matrix theory is to quantify the optimal rate of convergence to universal laws. We take up this problem for the Laguerre ensemble, characterised by the Dyson parameter , and the Laguerre weight , in the hard edge limit. The latter relates to the eigenvalues in the vicinity of the origin in the scaled variable . Previous work has established the corresponding functional form of various statistical quantities — for example the distribution of the smallest eigenvalue, provided that . We show, using the theory of multidimensional hypergeometric functions based on Jack polynomials, that with the modified hard edge scaling , the rate of convergence to the limiting distribution is , which is optimal. In the case , general the explicit functional form of the distribution of the smallest eigenvalue at this order can be computed, as it can for and general . An iterative scheme is presented to numerically approximate the functional form for general .

spacing=nonfrench

## 1. Introduction

Random matrix theory abounds in limit laws. As an example, consider an ensemble of complex Hermitian matrices, diagonal entries independently chosen from a distribution with mean and variance , and diagonal entries have mean zero and finite standard deviation. Assume too that all entries have finite -th moment for some fixed . Let denote the probability that the interval is free of eigenvalues. Then, independent of further details of the distributions, one has the limit theorem [34, 15, 1]

 limN→∞EN(0,(πa√2N,πb√2N))=det(IN−Ksine(0,b−a)), (1.1)

where is the integral operator on the interval with kernel

 K(x,y)=sinπ(x−y)π(x−y).

With such a limit law established, one can ask for the rate, as a function of , that the limiting distribution is approached. Such a question is core not only from a probabilistic viewpoint on random matrix theory, but also in the context of applications of random matrices, particularly to the Riemann zeros.

In this topic from number theory, the zeros of the Riemann zeta function in the upper half complex -plane are written as with in accordance with the Riemann hypothesis. For large a model put forward by Keating and Snaith [32] predicts that the zeros in intervals of size behave statistically like the eigenvalues of Haar distributed random matrices with . This is a refinement of the so-called Montgomery-Odlyzko law (see e.g. [31]) which asserts that for the statistical properties of the Riemann zeros, scaled to have mean spacing unity, coincide with the statistical properties of the bulk scaled eigenvalue of complex Hermitian matrices with . It immediately draws attention to finite corrections to limit laws such as (1.1). Indeed, using the model, the leading finite correction to various statistical quantities can be computed exactly and compared against Odlyzko’s high precision evaluation of large runs (over ) at very large heights () as announced in [35] and subsequently made available to interested researchers; see [3, 23, 8] for such studies.

The analysis of finite corrections to bulk scaled statistical quantities for Hermitian matrices, in contrast to matrices from , is complicated by these corrections not being translationally invariant. Thus for example one would expect that in (1.1) a dependence on both end points and at this order, rather than their difference as is the case of the limit law. This means that to obtain meaningful predictions the statistical quantity needs to be averaged over some mesoscopic interval. If one considers instead finite corrections at the spectrum edge, this complication is not present. Appreciation of this point, together with interest and earlier work within multivariate statistics [14, 30, 33], motivated our recent studies [26, 27] at the soft edge of various random matrix ensembles. One recalls that the defining feature of a soft edge of a random matrix spectrum is that at the spectrum edge the spectral density is non-zero on either side. An unexpected effect was observed: in all cases analysed a scaling and centring of eigenvalues could be found such that the leading correction to the statistical quantities being analysed occurred at order . One description of this is as a weak universality, since the rate of convergence to the universal laws is independent of the details of the random matrix model although the functional form of the leading correction term is ensemble dependent.

Random matrix ensembles can also exhibit a hard edge. Thus for example a positive definite matrix has its spectral density strictly zero for negative values, and so if the eigenvalue support borders the origin, this spectral boundary is a hard edge. Our interest in this paper is on the form of the optimal scaling of eigenvalues at the hard edge for various random matrix ensembles, chosen to make the leading correction term as small as possible. In fact there is some previous literature on this issue. Consider an complex Gaussian matrix , and form the product . Set . Let denote the probability that the domain contains no eigenvalues (the reason for the argument will become apparent later). Here is the so-called Dyson index, which for complex entries equals 2. Edelman, Guionnet and Péché [13] conjectured the large expansion

 EN,2(0;(0,s/4N);xae−x))=Ehard2(s;a)+a2NsddsEhard2(s;a)+O(1N2), (1.2)

where

 Ehard2(s;a):=limN→∞EN,2(0;(0,s/4N);xae−x)).

Proofs were subsequently provided by Perret and Schehr [36], and by Bornemann [7].

Furthermore Bornemann [7] observed that the hard edge scaling can be optimally tuned to

 s4N(1−a2N). (1.3)

A simple Taylor series expansion says that the term proportional to in (1.2) then cancels. Therefore the leading non-trivial large correction is .

One immediate question is to give the functional form of this correction. We give two such expressions in Section 2: one follows from an operator theoretic approach, and the other makes use of differential equations and Painlevé transcendents. In light of the weak universality observed at the soft edge, in the sense that optimal scaling gives the same order for the leading correction term in all cases where analysis has been possible (see [27] and references therein), it is natural to consider the large form of , and related statistical quantities, for other values of the Dyson index . Our working in Section 3, using the theory of multivariable hypergeometric functions based on Jack polynomials to analyse the hard edge asymptotics of for general and , shows that the extension of (1.3) to

 s4(N+a/β) (1.4)

gives an optimal convergence to the limiting distribution at rate . Hence, as already obtained at the soft edge, there is strong evidence for a weak universality of statistical quantities at the hard edge of random matrix ensembles: by tuning the scaling there is an optimal rate of convergence to the limiting distribution, which is the same for all ensembles considered. This is further confirmed in Section 5, where the methods of Section 3 are applied to analyse the hard edge asymptotics of the spectral density for even and general .

In the case , and for general , the results of Section 2 give the explicit functional form of the optimal correction term to the limiting hard edge scaled distribution of the smallest eigenvalue. A (scalar) hypergeometric representation gives the functional of the correction for general , in the case . This is given in Section 3.3. In Section 4, an iterative scheme is presented to numerically approximate the functional form of the correction for general .

## 2. Functional form of optimal correction term for complex Wishart matrices

### 2.1. Operator theoretic approach

The matrix product , where is an () matrix with complex standard Gaussian entries is referred to as a complex Wishart matrix. The eigenvalue PDF has the functional form (see e.g. [19, Ch. 3])

 1ZN,2N∏l=1w2(xl)∏1⩽j

with known as the Laguerre weight. The function if is true and otherwise, while denotes the normalisation. It is the exponent 2 on the product of differences which determines that the Dyson index here is .

For a general weight , a PDF of the form (2.1) has the special feature of its general -point correlation function being given by a determinant, fully determined by a single function of two variables . This kernel involves the orthogonal polynomials relating to the weight (see e.g. [19, Ch. 5]). In the present case

 ρk(x1,…,xk)=det[KN(xi,xj)]i,j=1,…,k (2.2)

where ,with denoting the Laguerre polynomials,

 KN(x,y)=(w2(x)w2(y))1/2N−1∑n=0L(a)n(x)L(a)n(y)hn,hn=n!Γ(a+n+1). (2.3)

Moreover the normalisation in (2.1) is given by . With use of the Christoffel-Darboux formula (see e.g. [19, Prop. 5.1.3]), and the recurrence , the correlation kernel (2.3) can be written

 KN(x,y)=N!e−(x+y)/2(xy)a/2Γ(a+N)L(a)N(x)L(a−1)N(y)−L(a−1)N(x)L(a)N(y)x−y. (2.4)

A corollary of the determinantal structure (2.2) is that the gap probability (the extra argument on top of the arguments used in (1.2) can be thought of as a thinning parameter: each eigenvalue is deleted independently with probability and therefore ; see e.g. the recent works [20, 23, 8, 10, 2, 9] as well as the pioneering paper [4]) has the operator theoretic form

 EN,2(0;J;xae−x;ξ))=det(I−ξKN,J), (2.5)

where is the identity operator and is the integral operator with the kernel (2.4) supported on . As a consequence, the large expansion of (2.5) for  — referred to as a hard edge scaling — can be deduced from knowledge of the corresponding large asymptotics of (2.4).

###### Lemma 1.

For general fixed, the large expansion of the scaled Laguerre polynomials is given by

 (2.6)

where the are Bessel functions of the first kind. This holds uniformly in in a compact set on the positive half line.

###### Proof.

Begin with the integral representation for the Laguerre polynomials,

 2aL(a)N(x/4N)=∫Cdz2πie−xz/8N(z+2)N+azN+1

where is a positively oriented Hankel loop contour that encloses the origin but does not contain . Change variables and use the elementary expansion

 (1+yN)N=ey(1−y22N+3y4+8y324N2+⋯).

Then (2.6) follows by recognising

 Ja(x)=∫Cdz2πiex2(z−1/z)z−a−1.

###### Proposition 2.

For large , and valid uniformly for and contained in a compact set on the positive half line,

 14NKN(x4N,y4N)=K(a)∞(x,y)+1NL(a)1(x,y)+1N2L(a)2(x,y)+O(1N3), (2.7)

where

 K(a)∞(x,y) =√yJa(√x)J′a(√y)−√xJ′a(√x)Ja(√y)2(x−y), (2.8) L(a)1(x,y) =a8Ja(√x)Ja(√y), (2.9)

and

 L(a)2(x,y)=−1192[(a2+x+y)Ja(√x)Ja(√y)+(xy)1/2J′a(√x)J′a(√y)−(3a2−2)(√yJa(√x)J′a(√y)+√xJ′a(√x)Ja(√y))]. (2.10)

In the above .

###### Proof.

The expansion (2.6) can readily be substituted into (2.4). With use of Stirling’s formula for the prefactor terms, and after simplification using the Bessel recurrences

 2αuJα(u)=Jα−1(u)+Jα+1(u),2J′α(u)=Jα−1(u)−Jα+1(u),

(2.7) follows. ∎

It can be observed that (2.9) is related to (2.8) by the simple derivative operation given by

 L(a)1(x,y)=a2(x∂∂x+y∂∂y+1)K(a)∞(x,y). (2.11)

The quantity is well known as the limiting hard edge kernel [17]. The task now is to appropriately centre the hard edge scaling as in (1.3) so that the term in (2.7) proportional to cancels and explicitly express the leading correction term now proportional to . It should be noted that since the scaling factor is multiplicative, the term (2.10) is dependent on the scaling variable up to . For this reason and for the ease of calculations, we work with the rescaled hard edge variable

 x4N+2a=x4N(1−a2N+a24N2−⋯) (2.12)

which agrees with (1.3) in the first two terms for large . Making the replacement , then performing a Taylor series expansion of (2.7) making use of the relation (2.11) eliminates the term in (2.7).

###### Proposition 3.

For large ,

 14N+2aKN(x4N+2a,y4N+2a)=K(a)∞(x,y)+1N2^L(a)2(x,y)+O(1N3), (2.13)

where

 ^L(a)2(x,y)=−1192((a2+x+y)Ja(√x)Ja(√y)+(xy)1/2J′a(√x)J′a(√y)+2√yJa(√x)J′a(√y)+2√xJ′a(√x)Ja(√y)). (2.14)

Immediate from (2.13) is that the spectral density (or the one-point correlation function) at the hard edge has the large expansion

 14N+2aρN(x4N+2a)=ρ∞,0(x)+1N2^ρ∞,2(x)+O(1N3), (2.15)

with

 ρ∞,0(x) =14(Ja(√x)2−Ja+1(√x)Ja−1(√x)) (2.16) ^ρ∞,2(x) =−1192((2x+a2)Ja(√x)2+4√xJa(√x)J′a(√x)+xJ′a(√x)2). (2.17)
###### Remark 4.

The functional form (2.17) appears as a sum of products of a special function (i.e. the Bessel function) and its derivative. This is analogous to the known density in the soft edge for both Gaussian and Laguerre cases where the special function of interest is the Airy function (see e.g. [19, Ch. 7]).

As shown in [8], the integral operator formula (1.2) for a bounded interval can be expanded in large up to its first correction. Let denote the integral operator on with kernel (2.8) and denote the integral operator on with kernel (2.14). In keeping with the notation of [8], define

 Ω(K∞,(0,s)):^L2,(0,s)=−det(I−K∞,(0,s))Tr((I−K∞,(0,s))−1^L2,(0,s)). (2.18)

Then we have

 det(I−ξKN,(0,s/(4N+2a)))=det(I−ξK∞,(0,s))+1N2Ω(ξK∞,(0,s)):ξ^L2,(0,s)+O(1N3). (2.19)

This framework has the numerical advantage, in the methods demonstrated by Bornemann [5, 6], in that the evaluation of converges exponentially fast if the kernel is analytic in a neighbourhood of . Explicit methodology for the numerical computation of has been presented in [8] and later used in [26].

### 2.2. Painlevé theory approach

The operator theoretic formula (2.5) can be reformulated in terms of the -function expression [38, 28]

 det(I−ξKN,J)=exp(∫s0U(a)N(x;ξ)dxx), (2.20)

where is a solution to the -form Painlevé V differential equation

 (xσ′′)2−(σ−xσ′+2(σ′)2+(a+2N)σ′)2+4(σ′)2(σ′+N)(σ′+a+N)=0 (2.21)

subject to the boundary condition

 U(a)N(x;ξ)∼x→0+−ξxKN(x,x). (2.22)

Under the hard edge scaled variable, the expansion

 U(a)N(x4N;ξ)=σ0(x;ξ)+1Nσ1(x;ξ)+1N2σ2(x;ξ)+O(1N3) (2.23)

is required for consistency with the results outlined in the operator theoretic approach.

###### Proposition 5.

Let , be specified by (2.8), (2.9). Then satisfies the -form Painlevé III equation [37, 28]

 (xσ′′)2+σ′(1+4σ′)(xσ′−σ)=(aσ′)2 (2.24)

subject to the boundary condition

 σ0(x;ξ)∼x→0+−ξxK(a)∞(x,x) (2.25)

and satisfies the second order linear differential equation

 A(x)σ′′(x)+B(x)σ′(x)+C(x)σ(x)=D(x) (2.26)

where

 A(x) =4x2σ′′0(x) B(x) =24xσ′0(x)2−16σ0(x)σ′0(x)+4(x−a2)σ′0(x)−2σ0(x) C(x) =−2σ′0(x)(4σ′0(x)+1) D(x) =−aσ′0(x)(xσ′0(x)−σ0(x)) (2.27)

subject to the boundary condition

 σ1(x;ξ)∼x→0+−ξxL(a)1(x,x). (2.28)
###### Proof.

Substitute the proposed expansion given by (2.23) and replace with the usual hard edge scaled variable into the -form Painlevé equation (2.21). Comparing the coefficients proportional to and gives the equations (2.24) and (2.26) respectively. The boundary condition follows directly from (2.22). ∎

It can be verified with use of (2.24) that

 σ1(x;ξ)=a2xσ′0(x;ξ) (2.29)

is a solution to (2.26). This is in keeping with the feature that the rescaled hard edge variable (2.12) gives a correction term of . Therefore consider the new expansion

 U(a)N(x4N+2a;ξ)=σ0(x;ξ)+1N2^σ2(x;ξ)+O(1N3) (2.30)

and repeat the exercise outlined below Proposition 5 by substituting the above equation into (2.21) and using the variable (2.12). This shows that too satisfies a second order linear differential equation.

###### Proposition 6.

The correction satisfies the second order linear differential equation (2.26) but with replaced with

 ^D2(x)=18(xσ′0−σ0)2, (2.31)

subject to the boundary condition

 ^σ2(x)∼x→0+−ξx^L(a)2(x,x) (2.32)

with specified by (2.14). Furthermore,

 EN,2(0;(0,s4N+2a);xae−x;ξ)=exp(∫s0σ0(x;ξ)dxx)+1N2(∫s0^σ2(x;ξ)dxx)exp(∫s0σ0(x;ξ)dxx)+O(1N3). (2.33)
###### Remark 7.

Results analogous to the above are known for the soft edge scaling of the Laguerre and Gaussian unitary ensembles [26], although as emphasised in the Introduction, upon optimal centring and scaling, the leading correction term is of order . An equation with a structure analogous to (2.26) has also appeared in a recent study of the two-point diagonal spin-spin correlation function of the two-dimensional Ising model [24].

Let denote the probability density function for the smallest eigenvalue in the ensemble (2.1) with . As noted below (2.4) the parameter is a thinning parameter: each eigenvalue is deleted independently with probability . This quantity relates to the gap probability according to

 pN,β(s;xae−x;ξ)=−ddsEN,β(0;(0,s);xae−x;ξ). (2.34)

According to the above working, this probability density for has a well defined hard edge scaling limit, obtained by replacing by as specified by (2.12), and the use of this variable gives an optimal rate of convergence at . Moreover the above working gives formulas implying the functional form of the correction term at this order.

Figure 2.1 is a numerical demonstration of the theoretical correction after rescaling the hard edge variables. Let as abbreviate the notation in the case , and as , and write . The figure given by the blue steps is the difference

 N2(p#N(sN,2)−p∞,0(s)). (2.35)

Here, in keeping with (2.34) and the leading term in (2.19),

 p∞,0(s)=−ddsdet(I−K∞,(0,s))∣∣a=1

is the limiting PDF of the smallest hard edge scaled eigenvalue and is the normalised histogram of the smallest eigenvalue from numerically generated complex Wishart matrices scaled by the factor with , . The red solid line is the predicted leading correction term of implied by either (2.19) or (2.33) substituted in (2.34).

## 3. General β case

### 3.1. Background theory

It was noted above that the matrix product , where is an , , matrix with complex standard Gaussian entries has the eigenvalue PDF (2.1) (set ). In the same setting except that the entries are now real standard normals, (2.1) requires modification to now read (see e.g. [19, Ch. 3])

 1ZN,1N∏l=1w1(xl)∏1⩽j

with . Notice that here the Dyson index is . A natural interpolation between (2.1) and (3.1) is the PDF

 1ZN,βN∏l=1wβ(xl)∏1⩽j

with the weight , referred to as the Laguerre ensemble. For general it can in fact be realised as specifying the PDF of the squared singular values of certain bi-diagonal matrices with independent entries [12].

The exact evaluation of for has been studied using methods of Selberg integral theory, and associated special functions (see [19, Ch. 13]). (Note that in distinction to the results of the previous section, it is not possible to include a general thinning parameter  — the results below are restricted to the case , which corresponds to no thinning.) To summarise the findings, introduce the generalised hypergeometric function

 pF(α)q(a1,…,ap;b1,…,bq;x1,…,xm)=∑k=(k1,…,km)1|k|![a1](α)k⋯[ap](α)k[b1](α)k⋯[bq](α)kC(α)k(x1,…,xm) (3.3)

where denotes a partition with

 |k|=m∑j=1kj,[a](α)k=m∏j=1Γ(a−(j−1)/α+kj)Γ(a−(j−1)/α)

and denotes the renormalised symmetric Jack polynomials (see [19, Eq. 13.1]). In particular when , (3.3) recovers the usual definition of the classical hypergeometric function.

It was shown [18] (see also [19, Proposition 13.2.6]) that given ,

 EN,β(0,(0,s);xae−βx/2)=e−βNs/21F(β/2)1(−N;2a/β;(−s)a), (3.4)

where has the multi-dimensional integral representation

 1F(β/2)1(−N;2a/β;(−s)a)=1Ma(2/β−1,N,2/β)∫1/2−1/2dx1⋯∫1/2−1/2dxa×a∏l=1e2πixl(2/β−1)(1+e−2πixl)N+(2/β−1)ese2πixl∏1⩽j

with

 Mn(A,B,C)=n∏j=1Γ(1+A+B−C+jC)Γ(1+jC)Γ(1+A−C+jC)Γ(1+B−C+jC)Γ(1+C). (3.6)

Here the notation appearing in (3.4) denotes the -tuple .

### 3.2. Hard edge scaling

We seek the large asymptotic expansion of (3.5), after replacing with the hard edge scaling , and under the assumption that is a positive integer. Regarding the leading term, it is immediate from the definition (3.3), and the fact that is a homogeneous function of its arguments, that (see [19, Prop. 13.2.7])

 Ehardβ(s;a):=limN→∞EN,β(0,(0,s/4N);xae−βx/2)=e−βs/80F(β/2)1(2a/β;(s/4)a). (3.7)

Moreover, beginning with (3.5) it is possible to express this function as an -dimensional integral,

 0F(β/2)1(2a/β;(s/4)a)=Γ(2/β)a(2π)aΓ(1+a)∫π−πdθ1⋯∫π−πdθaa∏l=1eiθl(2/β−1)eseiθl/4+e−iθl×∏1⩽j

This integral formula is a variant of the one given in earlier studies [18], [19]. In the latter the integrand involves the variable instead of as above. The form given in (3.8) has the advantage of being better suited to the analysis of subleading terms in the large expansion.

To see how to obtain (3.8) from (3.5), and to furthermore extend the expansion beyond leading order, begin by changing variables and then scaling . Following this prescription, and using Cauchy’s theorem to argue that the contour return the unit circle (see [19, Exercise 13.1 q4(ii)]) gives

 1F1(−N;2a/β;(−s/4N)a)=Na(2/β−1)Ma(2/β−1,N,2/β)1(2πi)a∫γdz1z1⋯∫γdzaza×a∏l=1z2/β−1le(N+2/β−1)log(1+1/(Nzl))eszl/4∏1⩽j

Expanding the logarithm then gives

 1F1(−N;2a/β;(−s/4N)a)=Na(2/β−1)Ma(2/β−1,N,2/β)1(2π)a∫π−πdθ1⋯∫π−πdθa×a∏l=1eiθl(2/β−1)eseiθl/4+e−iθl(1+(2β−1)e−iθlN−12e−2iθlN+O(1N2))∏1⩽j

Additionally the normalisation appearing in (3.5) can be checked via Stirling’s formula to have the large form

 1Ma(2/β−1,N,2/β)=a∏j=1Γ(2j/β)Γ(1+N+2(j−1)/β)Γ(1+2/β)Γ(N+2j/β)Γ(1+2j/β)=Na−2a/βΓ(2/β)aΓ(1+a)(1+a2βN(1−2β)+O(1N2)). (3.11)

Substituting shows (3.5) exhibits the large form

 1F1(−N;2a/β;(−s/4N)a)=Γ(2/β)a(2π)aΓ(1+a)∫π−πdθ1⋯∫π−πdθa×(1+a2βN(1−2β)−a∑l=1(1−2β)e−iθlN−12a∑l=1e−2iθlN+O(1N2))×