Computing marginal eigenvalue distributions
for the Gaussian and Laguerre orthogonal ensembles

Abstract

The Gaussian and Laguerre orthogonal ensembles are fundamental to random matrix theory, and the marginal eigenvalue distributions are basic observable quantities. Notwithstanding a long history, a formulation providing high precision numerical evaluations for $N$ large enough to probe asymptotic regimes, has not been provided. An exception is for the largest eigenvalue, where there is a formalism due to Chiani which uses a combination of the Pfaffian structure underlying the ensembles, and a recursive computation of the matrix elements. We augment this strategy by introducing a generating function for the conditioned gap probabilities. A finite Fourier series approach is then used to extract the sequence of marginal eigenvalue distributions as a linear combination of Pfaffians, with the latter then evaluated using an efficient numerical procedure available in the literature. Applications are given to illustrating various asymptotic formulas, local central limit theorems, and central limit theorems, as well as to probing finite size corrections. Further, our data indicates that the mean values of the marginal distributions interlace with the zeros of the Hermite polynomial (Gaussian ensemble) and a Laguerre polynomial (Laguerre ensemble).

⁽¹⁾ School of Mathematics and Statistics, The University of Melbourne, Victoria 3010, Australia.     Email: [email protected];   [email protected]

⁽²⁾ Department of Physics, Shiv Nadar Institution of Eminence, Gautam Buddha, Uttar Pradesh - 201314, India.    Email: [email protected]

1 Introduction

1.1 Background

Let $G_{n,N}$ denote an $n\times N$ rectangular standard real Gaussian matrix, meaning that each entry is chosen independently as a standard normal random variable. With $A_{N}=G_{N,N}$ and $B_{n,N}=G_{n,N}$ $(n\geq N)$ being particular square and rectangular real Gaussian matrices, introduce the real symmetric matrices

$$X_{N}={1\over 2}(A_{N}+A_{N}^{T}),\qquad Y_{n,N}=B_{n,N}^{T}B_{n,N}.$$ (1.1)

Random matrices with the construction of $X_{N}$ are said to form the Gaussian orthogonal ensemble (GOE) [28, §1.1]. Random matrices with the construction of $Y_{N}$ are examples of uncorrelated real Wishart matrices from multivariate statistics [28, §3.2], and through their eigenvalue PDF relate to the Laguerre orthogonal ensemble (LOE); see below.

Standard results in random matrix theory [28, Prop. 1.3.4 and 3.2.2] give that the eigenvalue probability density functions (PDFs) for the random matrices (1.1) are

$$\mathcal{P}^{\rm G}(x_{1},\dots,x_{N}) ={1\over Z_{N}^{\rm G}}\prod_{l=1}^{N}e^{-x_{l}^{2}/2}\prod_{1\leq j<k\leq N}|x_{k}-x_{j}|$$ (1.2)

$$\mathcal{P}^{\rm L}(x_{1},\dots,x_{N}) ={1\over Z_{N}^{\rm L}}\prod_{l=1}^{N}x_{l}^{a}e^{-x_{l}/2}\mathbbm{1}_{x_{l}>0}\prod_{1\leq j<k\leq N}|x_{k}-x_{j}|,\quad a=(n-N-1)/2,$$ (1.3)

with normalisations

$$Z_{N}^{\rm G}=\prod_{j=1}^{N}\Gamma(1+j/2),\qquad Z_{n,N}^{\rm L}=N!\pi^{-N/2}2^{Nn/2}\prod_{j=1}^{N}\Gamma(j/2)\Gamma((j+n-N)/2).$$ (1.4)

In the above, we have used the subscripts G (L) in relation to the GOE (LOE). One should mention at this stage that the terminology LOE applies more generally when the exponent $a=(n-N-1)/2$ in each factor $x_{l}^{a}$ of (1.3) is regarded as a continuous parameter $a>-1$. Random matrices that realise this more general eigenvalue PDF are known [23], but they cannot be constructed out of matrices with standard Gaussian entries.

Our interest in this paper is in the computation of the individual eigenvalue marginals corresponding to (1.2) and (1.3), and the underlying conditioned gap probabilities. The eigenvalue marginals are the PDFs $f_{N}^{(\cdot)}(k;s)$ for the $k$-th largest eigenvalue $(k=1,\dots,N)$. As is well known [28, §8.1], these PDFs can be expressed in terms of $\{E_{N}^{(\cdot)}(j;(s,\infty))\}_{j=0}^{k-1}$, where $E_{N}^{(\cdot)}(j;(s,\infty))$ denotes the conditioned gap probability that there are exactly $j$ eigenvalues in the interval $(s,\infty)$, and which themselves specify the distribution of the random variable $\sum_{l=1}^{N}\mathbbm{1}_{x_{l}>s}$ corresponding to the number of eigenvalues in the half interval $(s,\infty)$. Specifically, one has the recursive formula

$$f_{N}^{(\cdot)}(k;s)={d\over ds}E_{N}^{(\cdot)}(k-1;(s,\infty))+f_{N}^{(\cdot)}(k-1;s)$$ (1.5)

subject to the initial condition $f_{N}^{(\cdot)}(0;s)=0$, telling us that

$$f_{N}^{(\cdot)}(k;s)={d\over ds}\sum_{l=0}^{k-1}E_{N}^{(\cdot)}(l;(s,\infty)).$$ (1.6)

Regarding motivation, one notes that the individual eigenvalue marginals for the GOE appear explicitly in formulas relating to critical points of isotropic Gaussian random fields [2, 3]. Further, associated with individual eigenvalues, and for the conditioned gap probabilities, are various limit laws. While for the former these are best known at the spectrum edge [27, 57], there are also limit laws of individual eigenvalues in the bulk, and moving inwards from the edge with the matrix size [39, 48]. This is similarly true of the conditioned gap probabilities, or more specifically the random variable for the number of eigenvalues in an interval, in the setting that the number of eigenvalues in the interval tends to infinity [19, 54, 41].

Results in the literature have addressed the computation of (1.5) for the eigenvalue PDFs (1.2) and (1.3) in two distinct, fundamental ways. The first line of study aimed to give exact functional forms for small values of $N$. The tool for this was the discovery by Davis [20, 21] that associated with $\{f_{N}^{(\cdot)}(k;s)\}$ is an $(N+1)\times(N+1)$ first order matrix differential equation, which allows for a recursive computational scheme. In [21, Appendix B] there is a listing of exact functional forms of $\{f_{N}^{\rm L}(k;s)\}_{k=1}^{N}$ up to $N=5$. These ideas were developed in the Gaussian case (1.2) in the PhD thesis of Eckert [25] (this reference is available in electronic form on the internet), which contains the evaluations of $\{f_{N}^{\rm G}(k;s)\}_{k=1}^{N}$ up to $N=7$.

Another of the results of [25], was the identification of a transcendental basis for the PDF $f_{N}^{\rm G}(0;s)$ (a summary is given in [40]). To specify this basis, write $\phi(s)=e^{-s^{2}/2}$, $\Phi_{1}(s)=\int_{-\infty}^{s}\phi(y)\,dy$, $\Phi_{2}(s)=\int_{-\infty}^{s}(\phi(y))^{2}\,dy$, and define functions $\{\Omega_{j}(s)\}$

$$\Omega_{2k+1}(s)=\Phi_{1}(s)(\Phi_{2}(s))^{k}\>\>(k=0,1,\dots),\qquad\Omega_{2k}(s)=(\Phi_{2}(s))^{k}\>\>(k=0,1,\dots).$$ (1.7)

With this notation, it is established in [25, after the change of variables $s\mapsto\sqrt{2}s$] that

$$f_{N}^{\rm G}(1;s)=\sum_{j=0}^{N-1}\pi_{j}(s)(\phi(s))^{N-1-j}\Omega_{j}(s)$$ (1.8)

for some polynomials $\{\pi_{j}(s)\}$. Moreover, a family of linear operators $\{T_{p}\}_{p=1}^{N-1}$ was identified with an explicit action on any functional form with the structure of the right hand side of (1.8), and which allows for each member of $\{f_{N}^{\rm G}(n;s)\}_{n=2}^{N}$ to be deduced using knowledge of (1.8). For the Laguerre case (1.3) it was shown in [31] for $k=1$, and subsequently in [32] for general $k$ ($k=1,\dots,N$), that with $a:=(n-N-1)/2$ a non-negative integer, there are polynomials $q_{j}(s;k,N,a)$ of degree $j(a+N-j)$ such that

$$f_{N}^{\rm L}(k;s)=\sum_{j=n}^{N}q_{j}(s;k,N,a)e^{-js/2};$$ (1.9)

see too [45]. In addition, the references [31] and [32] implement in Mathematica notebooks the recursive method of Davis which provides for the exact evaluation of the polynomials. Regarding efficiency, in the case $k=1$ and $a=0$, $N=40$ the evaluation of $\{q_{j}(s;k,N,a)\}_{j=1}^{N}$ takes around 21 minutes using a 2017 model IMac machine. However the recurrence is such that to compute $\{f_{N}^{\rm L}(k;s)\}_{k=1}^{m}$ requires having first computed $\{f_{M}^{\rm L}(k;s)\}_{k=1}^{{\rm max}\{m,M\}}$ for $M=1,\dots,N-1$. While the algorithm is only cubic in $N$, with three main loops, there are overheads associated with the intermediate computer algebra specific operations. This limits the practical use of the code with respect to the size of $N$ and $M$; for example with $a=0$ and $N=M=11$ the run-time is around 93 minutes using the same machine.

A special functional form has also been identified in the case that $a+1/2$ is a non-negative integer [32]. This is most conveniently written in terms of the cumulative distribution $F_{N}^{\rm L}(1;x)=\int_{0}^{x}f_{N}^{\rm L}(1;s)\,ds$, and depends too on the parity of $N$. Setting $N$ even for convenience, it was exhibited in [32] for small values of $N$ and conjectured to hold in general that there are polynomials $p_{l,1}^{\rm e}(x),p_{l,2}^{\rm e}(x)$ such that

$$F_{N}^{\rm L}(1;x)=\sum_{l=1}^{N/2+1}\Big{(}e^{-(l-1)x}p_{l,1}^{\rm e}(x)+\sqrt{x}{\rm erf}(\sqrt{x/2})e^{-(l-1/2)x}p_{l,2}^{\rm e}(x)\Big{)}.$$ (1.10)

The second fundamental approach to the computation of $f_{N}^{(\cdot)}(k;s)$ available in the literature — albeit presented only in the case $k=1$ — is due to Chiani [18]. The starting point for this is the classical knowledge in random matrix theory due to Mehta [46], that computing averages with respect to (1.2) gives rise to a Pfaffian structure, and which moreover remains true in the Laguerre case (1.3). This is a consequence of the method of integration over alternate variables, which can be traced back to de Bruijn [13]. For the average corresponding to $E_{N}^{(\cdot)}(0;(s,\infty))$ in both the Gaussian and Laguerre cases it was shown in [18] that the matrix elements of the Pfaffian, which at first are given as double integrals, can be calculated recursively starting from a particular special function (error function in the Gaussian case, incomplete gamma function in the Laguerre case). This allows for an efficient computation of the matrix elements to high precision, with the later feature carrying over to the numerical evaluation of the Pfaffian, and leading via an application of (1.5) in the case $n=1$ to a high precision evaluation of $f_{N}^{(\cdot)}(1;s)$. The advantage of this approach relative to the one based on exact functional forms is that larger values of $N$ can be accessed. For example, in Chiani [18] a run-time of just 5 seconds is reported for the Laguerre case with $a=0$, $N=200$, albeit for a fixed value of $s$.

1.2 Specific aims and paper outline

As already stated, our primary concern in this paper is with the computation of the individual eigenvalue marginals for the joint PDFs (1.2) and (1.3), and the underlying conditioned gap probabilities. We will proceed using the (well known — recall the final paragraph above) Pfaffian structure, combined with the insight from [18] that the matrix elements can be computed by recurrence. The computer algebra system Mathematica allows for a high precision computation of these matrix elements. Essential too is use of the Mathematica code of Wimmer [59] for the efficient numerical evaluation of a Pfaffian.

The details of the recurrences for the Gaussian and Laguerre cases differ, necessitating that they be treated separately. Also, the feature of the GOE of being independent of a parameter beyond the matrix size gives rise to unique marginal probability density functions. For these, following on from work began in the thesis of Eckert [25], there is interest in providing high precision statistical data relating to the cumulants. Eckert provides such a tabulation up to and including $N=7$. Using our exact (internally stored within Mathematica) functional forms obtained from the Pfaffian formulation allows this tabulation to be extended up to and including $N=12$ (Appendix A). Unlike the LOE, the GOE is symmetrical about the origin, which distinguishes the conditioned probabilities $\{E_{N}^{(\rm G)}(k;(0,\infty))\}_{k=0}^{\infty}$, which determine the number of positive eigenvalues. References relating to the interest in this quantity are given in the first paragraph of §3.1. This provides motivation to investigate the large $N$ form of the variance of the number of positive eigenvalues, and also the accuracy of a known local central limit theorem. In the literature the case $k=0$ of these conditioned probabilities. In particular, there is a known large $N$ asymptotic formula which we are able to illustrate. As two extra applications of our numerics in the Gaussian case, we consider the rate of convergence to a known central limit law for the marginal probability density function near the centre of the spectrum, and also interlacing properties of the mean of the marginal distributions (for the latter see Appendix A). Analogous asymptotic questions are probed in the Laguerre case in Section 3.2, with an interlacing property of the means of the marginal densities considered in Appendix A.

2 Pfaffian formulation and exact evaluations

Fundamental to our working is the generating function

$$\Xi_{N}^{(\cdot)}((s,\infty);\zeta):=\sum_{k=0}^{N}\zeta^{k}E_{N}^{(\cdot)}(k;(s,\infty))=\bigg{\langle}\prod_{l=1}^{N}\Big{(}\mathbbm{1}_{x_{l}<s}+\zeta\mathbbm{1}_{x_{l}>s}\Big{)}\bigg{\rangle}^{(\cdot)},$$ (2.1)

where the average is with respect to one of the joint PDFs (1.2) or (1.3). A Pfaffian form of the average in (2.1) can be written down for any joint PDF structurally identical to (1.2) and (1.3), although the details depend on the parity of $N$ [28, Prop. 6.3.4 and Exercises 6.3 q.1]. To specify the Pfaffian, introduce

$$H_{0}^{\rm G}(j,k;s)={1\over\Gamma(j/2)\Gamma(k/2)}\int_{-\infty}^{s}dx\,x^{j-1}e^{-x^{2}/2}\int_{-\infty}^{s}dy\,y^{k-1}e^{-y^{2}/2}{\rm sgn}(y-x)$$

$$H_{1}^{\rm G}(j,k;s)={1\over\Gamma(j/2)\Gamma(k/2)}\int_{s}^{\infty}dx\,x^{j-1}e^{-x^{2}/2}\int_{-\infty}^{s}dy\,y^{k-1}e^{-y^{2}/2}-(j\leftrightarrow k)$$

$$H_{2}^{\rm G}(j,k;s)={1\over\Gamma(j/2)\Gamma(k/2)}\int_{s}^{\infty}dx\,x^{j-1}e^{-x^{2}/2}\int_{s}^{\infty}dy\,y^{k-1}e^{-y^{2}/2}{\rm sgn}(y-x)$$ (2.2)

and

$$H_{0}^{\rm L}(j,k;s)={1\over\Gamma(a+j)\Gamma(a+k)}\int_{0}^{s}dx\,x^{a+j-1}e^{-x}\int_{0}^{s}dy\,y^{a+k-1}e^{-y}{\rm sgn}(y-x)$$

$$H_{1}^{\rm L}(j,k;s)={1\over\Gamma(a+j)\Gamma(a+k)}\int_{s}^{\infty}dx\,x^{a+j-1}e^{-x}\int_{0}^{s}dy\,y^{a+k-1}e^{-y}-(j\leftrightarrow k)$$

$$H_{2}^{\rm L}(j,k;s)={1\over\Gamma(a+j)\Gamma(a+k)}\int_{s}^{\infty}dx\,x^{a+j-1}e^{-x}\int_{s}^{\infty}dy\,y^{a+k-1}e^{-y}{\rm sgn}(y-x).$$ (2.3)

Note that in our notation $H_{\mu}^{\rm L}(j,k;s)$ the dependence on the parameter $a$ has been suppressed. Further introduce the column vector $[\nu_{j}^{(\cdot)}(s)]_{j=1}^{N+1}$, where $\nu_{N+1}^{(\cdot)}(s)=0$ and

$$\nu_{j}^{\rm G}(s) ={1\over\Gamma(j/2)}\bigg{(}\int_{-\infty}^{s}dx\,x^{j-1}e^{-x^{2}/2}+\zeta\int_{s}^{\infty}dx\,x^{j-1}e^{-x^{2}/2}\bigg{)}$$ (2.4)

$$\nu_{j}^{\rm L}(s) ={1\over\Gamma(a+j)}\bigg{(}\int_{0}^{s}dx\,x^{a+j-1}e^{-x}+\zeta\int_{s}^{\infty}dx\,x^{a+j-1}e^{-x}\bigg{)},$$ (2.5)

for $j=1,\dots,N$, where here the dependence on $\zeta$ has been suppressed in the notation for both quantities, as well as the dependence on $a$ in $\nu_{j}^{\rm L}$.

Proposition 2.1.

For $N$ even

$$\Xi_{N}^{\rm G}((s,\infty);\zeta) =2^{N/2}{\rm Pf}\Big{[}H_{0}^{\rm G}(j,k;s)+\zeta H_{1}^{\rm G}(j,k;s)+\zeta^{2}H_{2}^{\rm G}(j,k;s)\Big{]}_{j,k=1}^{N},$$

$$\Xi_{N}^{\rm L}((s,\infty);\zeta) =\prod_{j=1}^{N}{\pi^{1/2}\Gamma(a+j)\over\Gamma(j/2)\Gamma(a+(j+1)/2)}$$

$$\qquad\times{\rm Pf}\Big{[}H_{0}^{\rm L}(j,k;s/2)+\zeta H_{1}^{\rm L}(j,k;s/2)+\zeta^{2}H_{2}^{\rm L}(j,k;s/2)\Big{]}_{j,k=1}^{N}.$$ (2.6)

In the case of $N$ odd, the size of the matrices that its Pfaffian is being computed in (2.1) must be increased from $N$ to $N+1$ by the bordering of the existing matrix by an additional column with entries $[\nu_{j}^{\rm G}(s)]_{j=1}^{N+1}$ and $[\nu_{j}^{\rm L}(s/2)]_{j=1}^{N+1}$ in the Gaussian and Laguerre cases respectively.

Remark 2.1.

For both (1.2) and (1.3) there are what may referred to as Fredholm Pfaffian formulas for $\Xi_{N}^{(\cdot)}((s,\infty);\zeta)$ [42, Eq. (180)], which relate to some specific $2\times 2$ anti-symmetric integral kernel; see also [7, §2.2.2 and §2.2.3]. In the circumstance that the integral kernel is an analytic function, such Fredholm expressions offer powerful computational properties [6, 7]. However, as noted in [7, §2.2.3], for the matrix integral kernel that results from the Pfaffian point processes with orthogonal symmetry such as (1.2) and (1.3) one of the matrix entries contains the non-analytic additive factor ${\rm sgn}(x-y)$, which nullifies its computational value.

The computational scheme for the matrix elements depends on the details of the weight function and so differs in the Gaussian and Laguerre cases. Therefore each will be considered separately.

2.1 Gaussian case

Define

$$\Psi^{\rm G}(j;x):={1\over\Gamma(j/2)}\int_{-\infty}^{x}t^{j-1}e^{-t^{2}/2}\,dt,\quad I^{\rm G}(j,k;x):={1\over\Gamma(j/2)}\int_{-\infty}^{x}t^{j-1}e^{-t^{2}/2}\Psi^{\rm G}(k;t)\,dt$$

$$\tilde{\Psi}^{\rm G}(j;x):={1\over\Gamma(j/2)}\int_{x}^{\infty}t^{j-1}e^{-t^{2}/2}\,dt,\quad\tilde{I}^{\rm G}(j,k;x):={1\over\Gamma(j/2)}\int_{x}^{\infty}t^{j-1}e^{-t^{2}/2}\tilde{\Psi}^{\rm G}(k;t)\,dt.$$ (2.7)

Using this notation, the quantities in the Gaussian cases of (2) and (2.4) can be written

$$H_{0}^{\rm G}(j,k;s)=-I^{\rm G}(j,k;s)+I^{\rm G}(k,j;s),\quad H_{1}^{\rm G}(j,k;s)=-\tilde{\Psi}^{\rm G}(j;s){\Psi}^{\rm G}(k;s)+\tilde{\Psi}^{\rm G}(k;s){\Psi}^{\rm G}(j;s)$$

$$H_{2}^{\rm G}(j,k;s)=\tilde{I}^{\rm G}(j,k;s)-\tilde{I}^{\rm G}(k,j;s),\quad\nu_{j}^{\rm G}(s)=\Psi^{\rm G}(j;s)+\zeta\tilde{\Psi}^{\rm G}(j;s).$$ (2.8)

By writing $\tilde{\Psi}^{\rm G}$ in terms of ${\Psi}^{\rm G}$, the expression for $H_{1}^{\rm G}$ can be further simplified to read

$$H_{1}^{\rm G}(j,k;s)=-u_{j}{\Psi}^{\rm G}(k;s)+u_{k}{\Psi}^{\rm G}(j;s),$$ (2.9)

where $u_{l}=0$ for $l$ even and $u_{l}=2^{l/2}$ for $l$ odd. Following Chiani [18], using straightforward integration by parts, each of the functions in (2.1) can be determined by second order recurrences and appropriate initial conditions, thus allowing for the determination of the $H_{i}^{\rm G}$.

Proposition 2.2.

We have the recurrences

$$\Psi^{\rm G}(j;x)=-{x^{j-2}\over\Gamma(j/2)}e^{-x^{2}/2}+2\Psi^{\rm G}(j-2;x),\quad\tilde{\Psi}^{\rm G}(j;x)={x^{j-2}\over\Gamma(j/2)}e^{-x^{2}/2}+2\tilde{\Psi}^{\rm G}(j-2;x),$$

$$I^{\rm G}(j+2,k;x)=2I^{\rm G}(j,k;x)-{x^{j}e^{-x^{2}/2}\over\Gamma(j/2+1)}\Psi^{\rm G}(k;x)+{2^{-(j+k)/2}\Gamma((j+k)/2)\over\Gamma(j/2+1)\Gamma(k/2+1)}\Psi^{\rm G}(j+k;\sqrt{2}x),$$

$$\tilde{I}^{\rm G}(j+2,k;x)=2\tilde{I}^{\rm G}(j,k;x)+{x^{j}e^{-x^{2}/2}\over\Gamma(j/2+1)}\tilde{\Psi}^{\rm G}(k;x)-{2^{-(j+k)/2}\Gamma((j+k)/2)\over\Gamma(j/2+1)\Gamma(k/2+1)}\tilde{\Psi}^{\rm G}(j+k;\sqrt{2}x),$$ (2.10)

where in the first line $j\geq 2$, while in the second and third line $j\geq 0$. Associated initial conditions are

$$\Psi^{\rm G}(0;x)=0,\quad\Psi^{\rm G}(1;x)={1\over\sqrt{2}}\Big{(}1+{\rm erf}(x/\sqrt{2})\Big{)},$$

$$\tilde{\Psi}^{\rm G}(0;x)=0,\quad\tilde{\Psi}^{\rm G}(1;x)={1\over\sqrt{2}}\Big{(}1-{\rm erf}(x/\sqrt{2})\Big{)},$$

$$I^{\rm G}(0,k;x)=I^{\rm G}(j,0;x)=0,\quad\tilde{I}^{\rm G}(0,k;x)=\tilde{I}^{\rm G}(j,0;x)=0,\,j,k\geq 1,$$

$$I^{\rm G}(1,1;x)={1\over 2}\Big{(}\psi^{\rm G}(1;x)\Big{)}^{2},\quad\tilde{I}^{\rm G}(1,1;x)={1\over 2}\Big{(}\tilde{\psi}^{\rm G}(1;x)\Big{)}^{2}.$$ (2.11)

For a given (even) value of $N$, the recurrences in the first line of (2.2) are to be iterated for $j=2,3,\dots,2N$ using the initial conditions in the first line of (2.2). Thus we have available the values of

$$\{\Psi^{\rm G}(j;x),\tilde{\Psi}^{\rm G}(j;x)\}_{j=0}^{2N},$$ (2.12)

which are required to implement the recurrences for $I^{\rm G},\tilde{I}^{\rm G}$. In relation to the latter, as a start the initial conditions in the third line of (2.2) with $k=1$ and in the final line, are to be used in the recurrences of the second and third line of (2.2), allowing for the computation of

$$\{I^{\rm G}(j,1;x),\tilde{I}^{\rm G}(j,1;x)\}_{j=2}^{N}.$$ (2.13)

The quantities $I^{\rm G}$ and $\tilde{I}^{\rm G}$ also satisfy the symmetry relations

$$I^{\rm G}(j,k;x)=\psi^{\rm G}(j;x)\psi^{\rm G}(k;x)-I^{\rm G}(k,j;x),\quad\tilde{I}^{\rm G}(j,k;x)=\tilde{\psi}^{\rm G}(j;x)\tilde{\psi}^{\rm G}(k;x)-\tilde{I}^{\rm G}(k,j;x),$$ (2.14)

the derivation of which is an appropriate integration by parts (in the case of $I^{\rm G}(j,k;x)$ for example, the starting point is to note that $I^{\rm G}(j,k;x)=\int_{-\infty}^{x}\Big{(}{d\over dx}{\psi}^{\rm G}(j;x)\Big{)}{\psi}^{\rm G}(k;x)\,dx$). Setting $j=1$ and using knowledge of the values of (2.13) provides the values of

$$\{I^{\rm G}(1,k;x),\tilde{I}^{\rm G}(1,k;x)\}_{k=2}^{N}.$$ (2.15)

This, together with the $k$-dependent initial conditions in the third line of (2.2), provide initial conditions to use the recurrences of the second and third line of (2.2) to compute

$$\{I^{\rm G}(j,k;x),\tilde{I}^{\rm G}(j,k;x)\}_{j=2}^{N}$$ (2.16)

for each $k=2,\dots,N$. Combining then the evaluations in the final line of (2.2) with (2.15) and (2.16) we have the evaluation of the full $N\times N$ array

$$\{I^{\rm G}(j,k;x),\tilde{I}^{\rm G}(j,k;x)\}_{j,k=1}^{N}.$$ (2.17)

We see from (2.1) that knowledge of the values of (2.12) and (2.16) makes explicit all the matrix elements in the Pfaffian formula of Proposition 2.1 for $\Xi_{N}^{\rm G}((s,\infty);\zeta)$. It is well known that for $A$ a $2n\times 2n$ anti-symmetric matrix $[a_{j,k}]_{j,k=1}^{2n}$, one has for the corresponding Pfaffian the Laplace type expansion

$${\rm Pf}(A)=\sum_{r=1}^{2n-1}(-1)^{r+1}a_{r,2n}{\rm Pf}_{r,2n}(A),$$ (2.18)

where ${\rm Pf}_{r,2n}$ is the Pfaffian of the $(2n-2)\times(2n-2)$ anti-symmetric matrix obtained by deleting rows and columns $r,2n$; see e.g. [28, Exercises 6.1]. Iterative use of this allows for the $(2n-1)!!$ terms (each a degree $n$ monomial in the matrix elements $\{a_{j,k}\}$) in the fully expanded form of ${\rm Pf}(A)$ to be made explicit (a brief code in the computer algebra system Mathematica for this purpose is given in [55]). From (2.1), in relation to the first Pfaffian in (2.1) this is a polynomial in $\zeta$, the coefficients of which are the probabilities $\{E_{N}^{\rm G}(k;s)\}_{k=1}^{N}$. With knowledge of the latter, application of the formula (1.5) gives $\{f_{N}^{\rm G}(k;s)\}_{k=1}^{N}$. The symmetry of the GOE spectrum upon reflection about the origin implies

$$f_{N}^{\rm G}(k;s)=f_{N}^{\rm G}(N+1-k;-s),$$ (2.19)

and so it suffices to extract only $\{E_{N}^{\rm G}(k;s)\}_{k=1}^{\lceil N/2\rceil}$.

Code has been written in the computer algebra system Mathematica to carry out the above steps. Most time consuming is the need to use the command Simplify on the coefficients of $\zeta$ obtained by the expansion of the Pfaffian. Without applying this command, as $N$ increases it is not possible to make use of the output for purposes of characterising the statistical properties of the $f_{N}^{\rm G}(k;s)$, due to the subsequent failure of the command NIntegrate. Such overheads have limited our exact determination (albeit stored electronically) to $N$ no bigger than 12.

Denote by $\kappa_{p}(N,k)$ the $p$-th cumulant of $f_{N}^{\rm G}(k;s)$, where in particular $\kappa_{1}(N,k)=\mu$ is the mean and $\kappa_{2}(N,k)=\sigma^{2}$ is the variance. Starting from the third cumulant, define too their scaled form by introducing

$$\gamma_{p-2}=\kappa_{p}(N,k)/\sigma^{p};$$ (2.20)

here $\gamma_{1}$ is the skewness and $\gamma_{2}$ the kurtosis. In [25, Appendix 4], a table of $\{\mu,\sigma,\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4}\}$ accurate up to and including the 7th decimal place is presented for $N=2,\dots,7$ and $k=1,\dots,\lceil N/2\rceil$, although an adjustment is required due to a different scale. In Appendix A below, using the exact functional form obtained from our code, we extend this data by providing the same table, now for $N=8,\dots,12$.

2.2 Laguerre case

Relevant here are the quantities

$$\Psi^{\rm L}(\alpha;x):={1\over\Gamma(\alpha)}\int_{0}^{x}t^{\alpha-1}e^{-t}\,dt,\quad I^{\rm L}(\alpha,\beta;x):={1\over\Gamma(\alpha)}\int_{0}^{x}t^{\alpha-1}e^{-t}\Psi^{\rm L}(\beta;t)\,dt$$

$$\tilde{\Psi}^{\rm L}(\alpha;x):={1\over\Gamma(\alpha)}\int_{x}^{\infty}t^{\alpha-1}e^{-t}\,dt,\quad\tilde{I}^{\rm L}(\alpha,\beta;x):={1\over\Gamma(\alpha)}\int_{x}^{\infty}t^{\alpha-1}e^{-t}\tilde{\Psi}^{\rm L}(\beta;t)\,dt.$$ (2.21)

They allow (2) and (2.5) to be written

$$H_{0}^{\rm L}(j,k;s)=-I^{\rm L}(a+j,a+k;s)+I^{\rm L}(a+k,a+j;s),$$

$$H_{1}^{\rm L}(j,k;s)=-\tilde{\Psi}^{\rm L}(a+j;s){\Psi}^{\rm L}(a+k;s)+\tilde{\Psi}^{\rm L}(a+k;s){\Psi}^{\rm L}(a+j;s)$$

$$\hskip 49.79231pt=-{\Psi}^{\rm L}(a+k;s)+{\Psi}^{\rm L}(a+j;s)$$

$$H_{2}^{\rm L}(j,k;s)=\tilde{I}^{\rm L}(a+j,a+k;s)-\tilde{I}^{\rm L}(a+k,a+j;s),\quad\nu_{j}^{\rm L}(s)=\Psi^{\rm L}(a+j;s)+\zeta\tilde{\Psi}^{\rm L}(a+j;s).$$ (2.22)

Chiani’s [18] strategy of integration by parts used in relation to $\Psi^{\rm L}$ and $I^{\rm L}$ allows all the functions in (2.2) to be determined by first order recurrences and appropriate initial conditions.

Proposition 2.3.

We have the recurrences

$$\Psi^{\rm L}(\alpha;x)=-{x^{\alpha-1}\over\Gamma(\alpha)}e^{-x}+\Psi^{\rm L}(\alpha-1;x),\quad\tilde{\Psi}^{\rm L}(\alpha;x)={x^{\alpha-1}\over\Gamma(\alpha)}e^{-x}+\tilde{\Psi}^{\rm L}(\alpha-1;x),$$

$$I^{\rm L}(\alpha+1,\beta;x)=I^{\rm L}(\alpha,\beta;x)-{x^{\alpha}e^{-x}\over\Gamma(\alpha+1)}\Psi^{\rm L}(\beta;x)+{2^{-(\alpha+\beta)}\Gamma(\alpha+\beta)\over\Gamma(\alpha+1)\Gamma(\beta)}\Psi^{\rm L}(\alpha+\beta;{2}x),$$

$$\tilde{I}^{\rm L}(\alpha+1,\beta;x)=\tilde{I}^{\rm L}(\alpha,\beta;x)+{x^{\alpha}e^{-x}\over\Gamma(\alpha+1)}\tilde{\Psi}^{\rm L}(\beta;x)-{2^{-(\alpha+\beta)}\Gamma(\alpha+\beta)\over\Gamma(\alpha+1)\Gamma(\beta)}\tilde{\Psi}^{\rm L}(\alpha+\beta;2x).$$ (2.23)

Using the recurrences (2.3), together with the symmetry relations

$$I^{\rm L}(\alpha,\beta;x)=\psi^{\rm L}(\alpha;x)\psi^{\rm L}(\beta;x)-I^{\rm L}(\beta,\alpha;x),\quad\tilde{I}^{\rm L}(\alpha,k;x)=\tilde{\psi}^{\rm L}(\alpha;x)\tilde{\psi}^{\rm L}(\beta;x)-\tilde{I}^{\rm L}(\beta,\alpha;x)$$ (2.24)

(as with (2.14), these follow by an appropriate integration by parts), the formulas of (2.2) imply recurrences for $H_{0}^{\rm L}$ and $H_{2}^{\rm L}$, the first of which can be found in [18, Eq. (18)].

Corollary 2.1.

We have

$$H_{0}^{\rm L}(j,k+1;s)=H_{0}^{\rm L}(j,k;s)-{x^{a+k}\over\Gamma(a+k+1)}e^{-x}\Psi^{\rm L}(a+j;x)\\ +2^{-(2a+j+k-1)}{\Gamma(2a+j+k)\over\Gamma(a+k+1)\Gamma(a+j)}\Psi^{\rm L}(2a+j+k;2x)$$ (2.25)

and

$${H}_{2}^{\rm L}(j,k+1;s)={H}_{2}^{\rm L}(j,k;s)-{x^{a+k}\over\Gamma(a+k+1)}e^{-x}\tilde{\Psi}^{\rm L}(a+j;x)\\ +2^{-(2a+j+k-1)}{\Gamma(2a+j+k)\over\Gamma(a+k+1)\Gamma(a+j)}\tilde{\Psi}^{\rm L}(2a+j+k;2x).$$ (2.26)

To compute $\{H_{\mu}(j,k;s)\}_{j=k}^{N}$ ($\mu=0,1,2$) we see from (2.25), (2.26) and the formula for $H_{1}$ in (2.2) that it suffices to have knowledge of $\{\Psi^{\rm L}(a+j;x),\tilde{\Psi}^{\rm L}(a+j;x)\}_{j=1}^{N}$ and $\{\Psi^{\rm L}(2a+j;x),\tilde{\Psi}^{\rm L}(2a+j;x)\}_{j=1}^{2N-1}$ (note that the initial conditions for the recurrences (2.25) and (2.26) are $H_{0}^{\rm L}(j,j;s)=H_{2}^{\rm L}(j,j;s)=0$). If the aim is to obtain exact functional forms involving a minimal basis in the sense of (1.9) and (1.10), these functions are to be computed using the recurrences on the first line of Proposition 2.3, together with the initial conditions

$$\Psi^{\rm L}(0;x)=1,\quad\Psi^{\rm L}(1/2;x)={\rm erf}(\sqrt{x}),\quad\tilde{\Psi}^{\rm L}(0;x)=0,\quad\tilde{\Psi}^{\rm L}(1/2;x)=\Big{(}1-{\rm erf}(\sqrt{x})\Big{)}.$$ (2.27)

The latter are relevant to the setting of $\mathcal{P}^{\rm L}$ in (1.3) that the parameter $a$ is a non-negative integer, or a half integer no less than $-1/2$. But if the aim is numerical evaluation using Mathematica, the facts that

$$\Psi^{\rm L}(\alpha;x)=1-{\Gamma(\alpha;x)\over\Gamma(\alpha)},\quad\tilde{\Psi}^{\rm L}(\alpha;x)={\Gamma(\alpha;x)\over\Gamma(\alpha)},$$ (2.28)

where $\Gamma(\alpha;x):=\int_{x}^{\infty}t^{\alpha-1}e^{-t}\,dt$ is the incomplete gamma function, implies that there is no need for the use of a recurrence, as the special function $\Gamma(\alpha;x)$ is part of the package. Note too that $a$ in (1.3) can take continuous values from this viewpoint.

3 Pfaffian formulation and numerical evaluations

The use of the symbol $\zeta$ in (2.1) can be replaced by the use of complex roots of unity according to the Fourier sum formula

$$(N+1)E_{N}^{(\cdot)}(k;(s,\infty))=\sum_{l=0}^{N}e^{-2\pi ikl/(N+1)}\Xi_{N}^{(\cdot)}((s,\infty);e^{2\pi il/(N+1)})\\ =\begin{cases}(1+2{\rm Re}\,\sum_{l=1}^{N/2}e^{-2\pi ikl/(N+1)}\Xi_{N}^{(\cdot)}((s,\infty);e^{2\pi il/(N+1)}),&N\>{\rm even}\\ 1+(-1)^{k}\Xi_{N}^{(\cdot)}((s,\infty);-1)+2{\rm Re}\,\sum_{l=1}^{(N-1)/2}e^{-2\pi ikl/(N+1)}\Xi_{N}^{(\cdot)}((s,\infty);e^{2\pi il/(N+1)}),&N\>{\rm odd},\end{cases}$$ (3.1)

where to obtain the second line use has been made of the normalisation requirement that $\Xi_{N}^{(\cdot)}((s,\infty);1)=1$. The significance of this is that with $\zeta$ a specific (complex) value, for a given value of $s$, the matrix elements in the Pfaffian formulas for $\Xi_{N}^{(\cdot)}$ in Proposition 2.1 can all be calculated numerically rather than symbolically, thus allowing for a numerical determination of $E_{N}^{(\cdot)}(k;(s,\infty))$.

For $A$ a general $2n\times 2n$ anti-symmetric matrix, there is the standard result (see e.g. [28, Eq. (6.12)]) that

$$({\rm Pf}(A))^{2}=\det(A).$$ (3.2)

In the case that $A$ has numerical entries and one has prior knowledge that the Pfaffian is positive, by taking the square root this formula can be used to provide an efficient computation of ${\rm Pf}(A)$. However, when the Pfaffian is a general complex number as in (3.1) for $l\neq 0$, a decision has to be made about the correct branch of the square root. For the present setting, by noting from (2.1) that for large $s$ we have $\Xi_{N}^{(\cdot)}((s,\infty),\zeta)\sim E_{N}^{(\cdot)}(0;(s,\infty))\sim 1$ independent of $\zeta$, a possible way to proceed using (3.2) would be to start computing $\Xi_{N}^{(\cdot)}((s,\infty),\zeta)$ for a fixed complex $\zeta$ with $s$ large, for which its value is approximately 1, and then to reduce $s$ in small intervals. The branch of the square root required by (3.2) is to be chosen by the requirement of approximate (with respect to the small interval size) continuity in values. In Figure 1 this task is illustrated in a particular Gaussian case with $N=6,l=3$, by comparing the values in the complex plane of $\Xi_{N}^{\rm G}((s,\infty);e^{2\pi il/(N+1)})$ obtained by decreasing $s$ from $s=5$ to $s=-4$ with the values of $(\Xi_{N}^{\rm G}((s,\infty);e^{2\pi il/(N+1)}))^{2}$. Increasing the size of $N$, analogous plots are observed but with more and tighter rotations about the origin, making the task of choosing the correct branch of the square root for discrete increments in $s$ a delicate task.

Fortunately, since the 2012 work on Wimmer [59] and in particular the software provided as part of the corresponding arXiv posting, it is now possible to efficiently compute ${\rm Pf}(A)$ for $A$ a general $2n\times 2n$ anti-symmetric matrix with numerical entries directly, independent of (3.2) and the associated square root issue. This is based on unitary conjugations reducing $A$ to a skew symmetric tridiagonal form. For given $N,l,s$ it is the Mathematica code by Wimmer, applied to the Pfaffian formulas of Proposition 2.1, which we use to compute $\Xi_{N}^{(\cdot)}((s,\infty);e^{2\pi il/(N+1)})$ in (3.1).

There is an important point to make in relation to the loss of conditioning with respect to increasing values of $N$. This point is that the computation of determinants and Pfaffians is, in general, ill-conditioned with respect to truncations of the values of the entries [56]. Thus it is necessary to increase the number of digits in the floating point arithmetic with $N$, which is simple to do using Mathematica.

Figure 1: Plotted parametrically are the values in the complex plane of $\Xi_{N}^{\rm G}((s,\infty);e^{2\pi il/(N+1)})$ (left panel) and $(\Xi_{N}^{\rm G}((s,\infty);e^{2\pi il/(N+1)}))^{2}$ (right panel) for $N=6$ and $l=3$ The parameter is the variables $s$, which is taken in the range $s=5$ (at this value of $s$ the curve (graphically) is at the real value 1 in the complex plane) and decreased down to $s=-4$ (which corresponds to the other end of the curve).

3.1 Gaussian case

We will consider first the sequence of conditioned gap probabilities $\{E_{N}^{\rm G}(k;(s,\infty))\}_{k=0}^{N}\}$ determining the random variable ${\mathcal{N}}_{(s,\infty)}:=\sum_{l=1}^{N}\chi_{x_{l}\in(s,\infty)}$ for the number of eigenvalues in the semi-infinite interval $(s,\infty))$. The specific case $s=0$ corresponds to the number of positive eigenvalues, which has attracted particular interest for its relevance to stability questions in disordered systems [17], landscape based string theory [1], quantum cosmology [47], among other examples; see the introduction to [43] for more references.

On the basis of some approximate analysis, it was predicted in [17] that in a neighbourhood of the mean $\mu_{N}:=\langle{\mathcal{N}}_{(0,\infty)}\rangle=N/2$, and for large $N$, the random variable ${\mathcal{N}}_{(0,\infty)}$ satisfies the local central limit theorem

$${\rm Pr}\Big{(}({\mathcal{N}}_{(0,\infty)}-\mu_{N})=k\Big{)}\to{1\over\sqrt{2\pi\sigma_{N}^{2}}}e^{-k^{2}/(2\sigma_{N}^{2})},\quad\sigma_{N}^{2}={1\over\pi^{2}}\log N;$$ (3.3)

see also [43]. Note that for this to be a meaningful limit law one must have $k/\sigma_{N}$ to be of order unity as $N\to\infty$; also, for convenience, in (3.3) it is assumed that $N$ is even. A rigorous determination of the leading order variance $\sigma_{N}^{2}$ can be found in [52] (see the review [29, Remark 3.1.3] for related references). Our computational scheme applied to $\{E_{N}^{\rm G}(k;(0,\infty))\}_{k=0}^{N}\}$ can test the prediction (3.3) and moreover probe correction terms in a large $N$ expansion.

Before doing so, some remarks along these lines in relation to the random variable ${\mathcal{N}}_{(0,\infty)}$ for the Gaussian unitary ensemble (GUE) of complex Hermitian matrices (see e..g. [28, §1.3.1]) are in order. Up to the scaling $x_{l}\mapsto\sqrt{2}x_{l}$, the eigenvalue PDF for the GUE is given by (1.2) with the product therein now squared. Whereas the eigenvalue PDF for the GOE is an example of a Pfaffian point process, the eigenvalue PDF for the GUE is an example of a determinantal point process, which is a simpler structure with additional integrability properties. Leveraging the latter, it was proved in [60, Eq. (143)] (see also [43] for the first two leading orders)

$${\rm Var}\Big{(}{\mathcal{N}}_{(0,\infty)}^{\rm GUE}\Big{)}={\log 4N+\gamma+1\over 2\pi^{2}}-{\log 4N+\gamma\over 4\pi^{2}N}+{7\over 96\pi^{2}N^{2}}+{24\log 4N+24\gamma-41\over 192\pi^{2}N^{3}}\\ -{219\over 5120\pi^{2}N^{4}}+{\rm O}(N^{-5}\log N),$$ (3.4)

where $\gamma$ denotes Euler’s constant and for convenience $N$ is assumed even (the analogue for $N$ odd was also derived, as were the explicit form of the terms up to ${\rm O}(N^{-7}\log N)$). The determinantal structure implies that the generating function $\Xi_{N}^{\rm GUE}((0,\infty);\zeta)$, defined as in (2.1), has all its zeros on the negative real axis in the complex $\zeta$-plane. This, combined with knowledge of a central limit theorem for ${\mathcal{N}}^{\rm GUE}_{(0,\infty)}$ [19, 54], implies the validity of the local central limit theorem (3.3) in the GUE case, now with $\sigma_{N}^{2}={1\over 2\pi^{2}}\log N$ [33].

Returning now to the consideration of $\{E_{N}^{\rm G}(k;(0,\infty))\}_{k=0}^{N}\}$, we first list in Table 1 our computed values of

$${\rm Var}\Big{(}{\mathcal{N}}_{(0,\infty)}^{\rm GOE}\Big{)}:=\sum_{k=0}^{N}(k-N/2)^{2}E_{N}^{\rm G}(k;(0,\infty))$$ (3.5)

for $N$ starting at 10 and finishing at 100, with increments of 10.

$N$ ${\rm Var}({\mathcal{N}}_{(0,\infty)}^{\rm GOE})$ 10 0.4735719613 20 0.5464313214 30 0.5883873972 40 0.6179704941 50 0.6408395677

$N$ ${\rm Var}({\mathcal{N}}_{(0,\infty)}^{\rm GOE})$ 60 0.6594853769 70 0.6752272387 80 0.6888489801 90 0.7008544880 100 0.7115869419

Table 1: 10 decimal place accurate variances (3.5)

Guided by (3.4) and knowledge of the leading term as given by $\sigma_{N}^{2}$ in (3.3), we make the ansatz

$$\pi^{2}{\rm Var}({\mathcal{N}}_{(0,\infty)}^{\rm GOE})-\log N={c_{1}}+{c_{2}\log N\over N}+{c_{3}\over N}+\cdots$$ (3.6)

to fit to the data of Table 1. We do this by choosing various combinations of 3 rows of the tables, which results in the values $c_{1}\approx 2.4229$, $c_{2}\approx 0.006$, $c_{3}\approx-0.52$. The small value of $c_{2}$ relative to $c_{3}$ puts in doubt the correctness of (3.6), which may then be distinct to (3.4) in the second leading term. On the other hand, the value of $c_{1}$ is available in the literature as $3\log 2+\gamma+1-{\pi^{2}\over 8}=2.42295\dots$ [53, Eq. (40) with $a=0$, $\beta=1$], thus providing evidence for the accuracy of our numerical values.

Next we make some remarks relating to correction terms to the local central limit law (3.3). As already remarked, the appropriate scaling variable is $k/\sigma_{N}$. But with $k$ taking on integer values only, and $\sigma_{N}={1\over\pi}(\log N)^{1/2}$, numerical tabulation for moderate values of $N$ (say between 50 and 100) are aways probing only the tails of the limiting distribution. Even a value of $N=1000$ (which is out of reach for the practical implementation of our numerical methods due to the ill-conditioning), the value of $\sigma_{N}$ is $0.91...$, which implies varying $k$ by 1 is more than one unit of standard deviation. The expansion for the variance (3.5) suggests the scaled large $N$ expansion

$${\rm Pr}\Big{(}({\mathcal{N}}_{(0,\infty)}^{\rm GOE}-\mu_{N})=k\Big{)}\sim{1\over\sqrt{2\pi\sigma_{N}^{2}}}e^{-k^{2}/(2\sigma_{N}^{2})}+{1\over\sigma_{N}^{2}}P_{1}(k/\sigma_{N})+\cdots$$ (3.7)

for some functional form $P_{1}$. A question in keeping with recent literature [38, 10, 15] relates to (3.7) identifying the optimal rate of convergence, meaning can $\sigma_{N}^{2}$ herein be replace by a function of $N$ with leading large $N$ form $\sigma_{N}^{2}$ and with the effect of eliminating the correction term in (3.7)? This will happen when $P_{1}$ is related to the leading term by a derivative operation [50, 37].

As far as illustrating (3.3) goes (where we replace $\sigma_{N}^{2}$ by ${\rm Var}({\mathcal{N}}_{(0,\infty)}^{\rm GOE})$ as is usual in stating a local central limit theorem), in Table 2, where we compare $p_{N}^{\rm exact}(k):={\rm Pr}({\mathcal{N}}_{(0,\infty)}^{\rm GOE}-\mu_{N}=k)$ against

$$p_{N}^{\rm approx}(k):={1\over\sqrt{2\pi{\rm Var}({\mathcal{N}}_{(0,\infty)}^{\rm GOE})}}e^{-k^{2}/(2{\rm Var}({\mathcal{N}}_{(0,\infty)}^{\rm GOE}))}.$$ (3.8)

We remark that a theoretical justification of such a relation can, following [33], be undertaken by knowledge of the location of the zeros of the generating function (2.1); see the numerically evidenced discussion in [32, §3.5].

$k$ $p_{70}^{\rm exact}(k)$ $p_{70}^{\rm approx}(k)$ $\delta_{70}(k)$ 0 0.4838115 0.4854953 $-$0.0016838 1 0.2325255 0.2315227 0.0010028 2 0.0250092 0.0251084 $-$0.0000992 3 0.0005570 0.0006192 $-$0.0000622

$k$ $p_{100}^{\rm exact}(k)$ $p_{100}^{\rm approx}(k)$ $\delta_{100}(k)$ 0 0.4714611 0.4729291 $-$0.0014681 1 0.2350801 0.2342270 0.0008531 2 0.0284044 0.0284550 $-$0.0000506 3 0.0007803 0.0008479 $-$0.0000676

Table 2: Here $\delta_{N}(k):=p_{N}^{\rm exact}(k)-p_{N}^{\rm approx}(k)$

In applications to the stability questions cited in the beginning paragraph of this section, there is specific interest in the large $N$ form of ${\rm Pr}({\mathcal{N}}_{(0,\infty)}=0)$, which is in the large deviation regime with respect to (3.3). For the GUE, it has been proved [12, 22] that

$$\log{\rm Pr}\Big{(}{\mathcal{N}}_{(0,\infty)}^{\rm GUE}=0\Big{)}=c_{1}N^{2}+c_{2}\log N+c_{3}+\cdots$$ (3.9)

with $c_{1}=-{1\over 2}\log 3$, $c_{2}=-{1\over 12}$, $c_{3}={1\over 8}\log 3-{1\over 6}\log 2+\zeta^{\prime}(-1)$. In fact the results of [12] apply to the GOE case well, telling us that

$$\log{\rm Pr}\Big{(}{\mathcal{N}}_{(0,\infty)}^{\rm GOE}=0\Big{)}=\tilde{c}_{1}N^{2}+\tilde{c}_{2}N+\tilde{c}_{3}\log N+\tilde{c}_{4}+\cdots$$ (3.10)

with $\tilde{c}_{1}=-{1\over 4}\log 3$, $\tilde{c}_{2}=-{1\over 2}\log(1+{2\over\sqrt{3}})$, $\tilde{c}_{3}=-{1\over 24}$, $\tilde{c}_{4}=-{1\over 12}\log 2-{1\over 16}\log 3+{1\over 4}\log(1+{2\over\sqrt{3}})+{1\over 2}\zeta^{\prime}(1)$. Our high precision evaluation of $E_{N}^{\rm G}(0;(0,\infty))$ (which by definition is equal to ${\rm Pr}({\mathcal{N}}_{(0,\infty)}^{\rm GOE}=0)$) for various $N$ is used in Table 3 to illustrate the accuracy of the asymptotic expansion (3.10).

$N$	$\log E_{N}^{\rm G}(0;(0,\infty))$	$\log\tilde{E}_{N}^{\rm G}(0;(0,\infty))$	$\delta_{N}$
10	$-$31.4183282	$-$31.4167301	$-$0.0015980
20	$-$117.6805735	$-$117.6797917	$-$0.0007817
30	$-$258.8619980	$-$258.8614809	$-$0.0005170

Table 3: Here $\log\tilde{E}_{N}^{\rm G}(0;(0,\infty))$ is defined as the asymptotic expansion on the right hand side of (3.10), and $\delta_{N}:=\log{E}_{N}^{\rm G}(0;(0,\infty))-\log\tilde{E}_{N}^{\rm G}(0;(0,\infty))$. Notice that $\delta_{N}$ is decreasing at the rate $1/N$ as $N$ increases.

We turn our attention now to the marginal eigenvalue PDFs in the bulk, choosing for simplicity $f_{N}^{\rm G}((N+1)/2;s)$ with $N$ odd, this (by symmetry) being an even function of $s$. We use our computational scheme to compute the sum in (1.5) for discrete values of $s$, small with respect to the scale of $\sigma_{N}$. An interpolating function connecting these values is then formed within the Mathematica software, which allows for the derivative operation required in (1.5) to be carried out. As mentioned in the Introduction, there are limit laws available for the individual eigenvalues, moving inwards from the edge with the matrix size. Such limit laws were first derived in the GUE case by Gustavsson [39], with the main tool being the determinantal structure and its knowledge of a central limit theorem for the counting statistic ${\mathcal{N}}_{(s,\infty)}$ from [54]. Subsequently, it was realised by O’Rourke [48] that the use of the inter-relation even(GOE${}_{N}\cup{\rm GOE}_{N+1}$) = GUE_N from [36] (here the operation GOE${}_{N}\cup{\rm GOE}_{N+1}$ denotes the random superposition of the spectrum of a GOE ensemble with $N$ eigenvalues, and a GOE ensemble of $N+1$ eigenvalues, while the operation even$(\cdot)$ denotes observing only the even labelled eigenvalues when reading from either edge with the eigenvalues labelled successively) allows the results of [39] to be extended to the GOE case. Specifically, for a GOE eigenvalue $x_{k}$ near or at the centre of the spectrum (for us $x_{(N-1)/2}$ with $N$ odd), we have from [48, Remark 8] that

$$\lim_{N\to\infty}\Big{(}{2N\over\log N}\Big{)}^{1/2}x_{k}\mathop{=}\limits^{\rm d}{\tt N}[0,1],$$ (3.11)

where ${\tt N}[0,1]$ denotes the standard normal distribution. (An alternative derivation of this result, extended to the Gaussian $\beta$ ensemble, has recently been given in [26].) One remarks that in this formula the factor of $\sqrt{N}$ can be interpreted as resulting from the value of the density of the GOE at the origin, which to leading order is a scaled semi-circle; see [28, Eq. (1.52)].

In keeping with (3.11), introduce the scaled eigenvalues

$$X_{k}=\Big{(}{2N\over\log N}\Big{)}^{1/2}x_{k}.$$ (3.12)

Denote by $f_{N}^{\rm G}(k;X)$ the PDF of $X_{k}$, as is consistent with the notation used in (1.6). Our interest is in using our ability to compute this PDF, with $k=(N-1)/2$ for a sequence of $N$ values, to probe the leading rate of convergence to the limit law (3.11). Thus we seek the function of $N$, $\alpha(N)$ say, with $\alpha(N)\to 0$ as $N\to\infty$, and the function of $X$, $h(X)$ say, such that for large $N$ one has the asymptotic expansion

$$\Big{(}{2N\over\log N}\Big{)}^{1/2}f_{N}^{\rm G}((N-1)/2;X)\sim{1\over\sqrt{2\pi}}e^{-X^{2}/2}+\alpha(N)h(X)+\cdots,$$ (3.13)

where terms not written decay at a rate fast than $\alpha(N)$. A numerical approach can access the difference

$$\Big{(}{2N\over\log N}\Big{)}^{1/2}f_{N}^{\rm G}((N-1)/2;X)-{1\over\sqrt{2\pi}}e^{-X^{2}/2},$$ (3.14)

which is only an approximation to $\alpha(N)h(X)$ as it contains all the higher order terms in the expansion (3.13) as well. Nonetheless, one is lead to the prediction that $\alpha(N)={1\over\log N}$, and to (an approximation of) the graphical form of $h(X)$ after computing (3.14) for just the two values of $N$, $N=21$ and $N=41$. Moreover, inspection of the graphical form shows that it closely resembles (but is not equal to)

$$c\Big{(}1-X{d\over dX}\Big{)}e^{-X^{2}/2}$$ (3.15)

for a certain constant $c$; see Figure 2.

Suppose now that instead of the scaling (3.12) one was to introduce $\tilde{X}_{k}:=\Big{(}{\log N\over 2N}\Big{)}^{1/2}(1+\alpha(N))x_{k}$. Taylor expanding with respect to the small variable $\alpha(N)$ to reclaim the variable (3.12) shows

$$\Big{(}{2N\over\log N}\Big{)}^{1/2}(1+\tilde{c}\alpha(N))f_{N}^{\rm G}((N-1)/2;\tilde{X})\\ \sim{1\over\sqrt{2\pi}}e^{-X^{2}/2}+\alpha(N)\bigg{(}\tilde{c}\Big{(}1-X{d\over dX}\Big{)}e^{-X^{2}/2}+h(X)\bigg{)}+\cdots,$$ (3.16)

Thus, if it were to be that (3.15) was proportional to $h(X)$ it would be possible to choose $\tilde{c}$ to improve the rate of convergence by eliminating the term proportional to $\alpha(N)$; see [50, 38] for some examples. However, in the present situation, while the functional forms are very similar, they are not exactly the same, so such an improvement is not possible. On the other hand, from the viewpoint of numerical values rather than the rate of convergence, it follows that the PDF for the scaled eigenvalue $X_{(N-1)/2}/{\rm Var}(X_{(N-1)/2})$, rather than for (3.12), more accurately follows ${\tt N}[0,1]$ in distribution for finite $N$ (cf. (3.8)).

Remark 3.1.

A question of much interest (see e.g. [11] for motivation), but not accessible via our present results, is the finite size corrections for the scaled spacing distribution of two GOE eigenvalues near the centre of the spectrum (say $x_{N/2}$ and $x_{(N/2)+1}$ for $N$ even). In the case of the circular version of the GOE — the circular orthogonal ensemble (COE) (see [28, §2.2.2]) — the leading correction term to the large $N$ limit law is proportional to $1/N^{2}$ [34].

Figure 2: [color online] Plot of the difference (3.14) for $N=41$ (blue solid curve), superimposed with the approximation (3.15) (red dashed curve) for $c$ chosen for a matching at the origin. The graphs are shown for $X\geq 0$ only, as both are symmetrical about the origin.

3.2 Laguerre case

Relevant to the consideration of large $N$ limit laws in the Laguerre case is the limiting eigenvalue density (the limiting eigenvalue density is also relevant in the Gaussian case, but didn’t appear in our discussion since our consideration of the bulk was restricted to the neighbourhood of the origin for simplicity of presentation — in the Laguerre case such a simplification is not possible as the spectrum has no such point of symmetry). First, scale the LOE eigenvalues by writing $\Lambda_{j}=\lambda_{j}/N$. Set ${n\over N}={1\over c}$, where $n\geq N$ is defined as in (1.1). The density of $\{\Lambda_{j}\}$, normalised to integrate to unity, is given by the Marchenko-Pastur law [49]

$$\rho^{\rm MP}(x)={1\over 2\pi}{\sqrt{(c_{+}-x)(x-c_{-})}\over cx}\mathbbm{1}_{x\in(c_{-},c_{+})},$$ (3.17)

with $c_{\pm}=(1\pm\sqrt{c})^{2}$. By definition of the density, for large $N$ the expected number of the scaled eigenvalues in the interval $(\tilde{s},c_{+})$ is $N\mu_{\tilde{s}}$, $\mu_{\tilde{s}}=\int_{\tilde{s}}^{c_{+}}\rho^{\rm MP}(x)\,dx$. With $\mathcal{N}_{(\tilde{s},c_{+})}^{\rm L}(\{\Lambda_{j}\})$ the random variable for the number of scaled eigenvalues in $(\tilde{s},c_{+})$, considerations from the topic of principal component analysis motivated a study in [44] of the fluctuation ${\rm Var}\,\mathcal{N}_{(\tilde{s},c_{+})}^{\rm L}(\{\Lambda_{j}\})$. This quantity was shown to exhibit the asymptotic form ${1\over\pi^{2}}\log N$ (which is the same as $\sigma_{N}^{2}$ in (3.3)) independent of the value of $\tilde{s}$ ($c_{-}<\tilde{s}<c_{+}$), and moreover it was argued that a local central limit theorem quantitatively the same as (3.3) holds true.

Our ability to compute $\{E_{N}^{\rm L}(k;(s,\infty))\}_{k=0}^{N}$ from §2.2 allows for the above predictions to be illustrated and checked. First we consider

$$\Big{\langle}\mathcal{N}_{(\tilde{s},\infty)}(\{\Lambda_{j}\}\Big{\rangle}=\sum_{k=1}^{\infty}kE_{N}^{\rm L}(k;(s,\infty))$$ (3.18)

in the specific case that the Laguerre parameter $a$ in (1.3) is equal to $1$ and $\tilde{s}=1$ (which implies $\tilde{\mu}=0.391\dots$). in Table 4 we compare our high precision numerical evaluations of (3.18) against the leading order theoretical value $N\tilde{\mu}$, for a sequence of $N$ values, tabulating too their difference. The latter shows evidence of converging to a constant value at a rate proportional to $1/N$.

$N$	$\Big{\langle}\mathcal{N}_{(\tilde{s},\infty)}(\{\Lambda_{j}\}\Big{\rangle}$	$N\tilde{\mu}$	$\delta_{N}$
20	8.711490	7.8200	0.8914
40	16.543998	15.6401	0.9039
60	24.368259	23.4601	0.9081
80	32.190428	31.2802	0.9102

Table 4: 8 digit accurate means (3.18), their leading order theoretical value $N\tilde{\mu}$, and the difference $\delta_{N}$.

For definiteness we take the Laguerre parameter $a$ in (1.3) equal to $1$, which implies $c=1,c_{-}=0$ and $c_{+}=4$, and we take $s=N$ (i.e. $\tilde{s}=1$ which implies $\mu_{\tilde{s}}={2\over\pi}\int^{1}_{1/4}x^{-1/2}(1-x)^{1/2}\,dx$. The variance can be computed from the formula

$${\rm Var}\Big{(}{\mathcal{N}}_{(\tilde{s},c_{+})}^{\rm L}(\{\Lambda_{j}\}\Big{)}:=\sum_{k=0}^{N}(k-M_{\tilde{s}})^{2}E_{N}^{\rm L}(k;(N\tilde{s},\infty)),\quad M_{\tilde{s}}:=\Big{\langle}\mathcal{N}_{(\tilde{s},\infty)}(\{\Lambda_{j}\}\Big{\rangle}$$ (3.19)

(cf. (3.5)), which we do in Table 5 for the same parameters as in Table 4. We compare these values against the leading order theoretical prediction ${1\over\pi^{2}}\log N$, and display too the difference, which shows evidence of tending to a constant.

$N$	${\rm Var}({\mathcal{N}}_{(\tilde{s},c_{+})}^{\rm L})$	${1\over\pi^{2}}\log N$	$\delta_{N}$
20	0.5096148748	0.3035	0.2060
40	0.5779495211	0.3737	0.2041
60	0.6182432846	0.4148	0.2033
80	0.6469688861	0.4439	0.2029

Table 5: 10 decimal place accurate variances (3.5) against their leading order theoretical value, and the difference $\delta_{N}$

It is easy to check from the definitions that ${\rm Var}\Big{(}{\mathcal{N}}_{(\tilde{s},c_{+})}^{\rm L}(\{\Lambda_{j}\}\Big{)}={\rm Var}\Big{(}{\mathcal{N}}_{(c_{-},\tilde{s})}^{\rm L}(\{\Lambda_{j}\}\Big{)}$. The significance of this is the prediction that in the case $c_{-}=0$ [53, Eq. (48) with $\beta=1$, $\tilde{a}=\sqrt{\tilde{s}/4}$, $\tilde{\gamma}=0$],

$$\lim_{N\to\infty}\Big{(}{\rm Var}\Big{(}{\mathcal{N}}_{(c_{-},\tilde{s})}^{\rm L}(\{\Lambda_{j}\}\Big{)}-{1\over\pi^{2}}\log N\Big{)}={1\over\pi^{2}}\Big{(}3\log 2+\gamma-{\pi^{2}\over 8}+1+\log(\tilde{s}^{1/2}(1-\tilde{s}/4)^{3/2})\Big{)}.$$ (3.20)

Substituting $\tilde{s}=1$ gives the numerical value $0.20177\dots$ which is consistent with the values of $\delta_{N}$ in Table 5.

In addition we have carried out numerical computations comparing

$$p_{N}^{\rm exact,L}(k):={\rm Pr}\Big{(}({\mathcal{N}}_{(\tilde{s},c_{+})}^{\rm L}(\{\Lambda_{j}\}-\lfloor\langle{\mathcal{N}}_{(\tilde{s},c_{+})}^{\rm L}(\{\Lambda_{j}\}\rangle\rfloor)=k\Big{)}$$ (3.21)

against the prediction from the validity of a local central limit theorem

$$p_{N}^{\rm approx,L}(k):={1\over\sqrt{2\pi{\rm Var}({\mathcal{N}}_{(\tilde{s},c_{+})}^{\rm L}(\{\Lambda_{j}\})}}e^{-k^{2}/(2{\rm Var}({\mathcal{N}}_{(\tilde{s},c_{+})}^{\rm L}(\{\Lambda_{j}\})}.$$ (3.22)

Note that both these quantities are functions of $\tilde{s}$, which has been suppressed in our notation. Some results of our computations, which were carried out with $\tilde{s}=1$, are presented in Table 6. The accuracy of (3.22) is evident.

$k$ $p_{60}^{\rm exact,L}(k)$ $p_{60}^{\rm approx,L}(k)$ $\delta_{60}^{L}(k)$ 0 0.4536364 0.4546717 $-$0.0010352 1 0.3676810 0.3674131 0.00026781 2 0.0591954 0.0589041 $-$0.00029126 3 0.0017303 0.0018735 $-$0.00014319

$k$ $p_{90}^{\rm exact,L}(k)$ $p_{90}^{\rm approx,L}(k)$ $\delta_{90}^{L}(k)$ 0 0.4860010 0.4877209 $-$0.0017198 1 0.2671577 0.2662165 0.0009412 2 0.0318105 0.0318451 $-$0.0000346 3 0.0007564 0.0008348 $-$0.0000784

Table 6: Here $\delta_{N}^{L}(k):=p_{N}^{\rm exact,L}(k)-p_{N}^{\rm approx,L}(k)$

As for the GOE, and after fixing a value of the Laguerre parameter $a$ as well as $N$, we have access to a tabulation of the marginal eigenvalue PDFs in the bulk. On the theory side, using a very different set of ideas than those used to prove (3.11), now based on the tridiagonal model of [23] and martingale arguments, it is established in [51] that

$$\lim_{N,l\to\infty}{\pi\rho^{\rm MP}(\gamma_{l})\over\log N}\Big{(}x_{l}-N\gamma_{l}\Big{)}\mathop{=}\limits^{\rm d}{\tt N}[0,1],$$ (3.23)

where $\gamma_{l}$ is such that $\int_{\gamma_{l}}^{c_{+}}\rho^{\rm MP}(x)\,dx=l/N$. We recognise $N\gamma_{l}$ as the leading order of the mean (3.18), and $\pi^{2}/\log N$ as being the reciprocal of the leading order of the variance (3.19). We can therefore substitute the mean and variance into (3.23). With this done, and $l={2\over 5}N$ (then $\gamma_{l}=0.9677\dots$) we have computed the PDF for the scaled random variable, $s$ say, on the LHS of (3.23) for various values of $N$. Accurate agreement with the PDF of the standard normal is found in all cases. For example, with $N=60$, the absolute difference is no greater that $3\times 10^{-3}$ in the range of $s$ between $-4$ and 4.

Ancillary Mathematica files

The various tables and graphs in the text were produced by implementing the specified computational schemes as Mathematica notebooks. The symbolic computations of the marginal distributions for the GOE reported in Appendix A use GOEsymbolic.nb. Numerical calculations for the probabilities $\{E_{N}^{\rm G}(k;(s,\infty)\}_{k=0}^{N}$ with $s$ fixed as required to produce Tables 1–3 use GOEnumeric1.nb, while the numerical computation of the marginal distributions in the Gaussian case, which were required to produce Figure 2 use GOEnumeric2.nb. The Tables 4–6 in relation to the Laguerre case use LOEnumeric.nb. Table 7 in Appendix A relates to the LOE, and uses LOEsymbolic.nb.

Acknowledgements

PJF and BJS are supported by the Australian Research Council Discovery Project grant DP210102887, and a 2024 University of Melbourne Science Faculty small grant. SK acknowledges the support provided by SERB, DST, Government of India, via Grant No. CRG/2022/001751.

Appendix A

In Section 2.1 a formalism to compute the exact functional forms of the independent marginal individual eigenvalue PDFs $\{f_{N}^{\rm G}(k;s)\}_{k=1}^{\lfloor N/2\rfloor}$ has been presented. Here, extending results in [25, Appendix 4], the results of using these exact functional forms to compute the corresponding statistical quantities $\{\mu,\sigma,\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4}\}$ (recall (2.20)) is presented in tabular form, where in keeping with [25, Appendix 4] the decimals are truncated after the 7th place.

$N$	$k$	$\mu$	$\sigma$	$\gamma_{1}$	$\gamma_{2}$	$\gamma_{3}$	$\gamma_{4}$
8	1	3.2451029	0.6431706	0.2175978	0.0866120	0.0161454	$-$0.0388095
8	2	2.1372504	0.5423684	0.0972884	0.0126373	$-$0.0059766	$-$0.00572405
8	3	1.2367953	0.5064137	0.0468383	$-$0.0040771	$-$0.0046966	0.0017507
8	4	0.4060264	0.4926050	0.0142245	$-$0.0089303	$-$0.0015978	0.0039260
9	1	3.5094346	0.6296768	0.2232643	0.0913421	0.0201166	$-$0.0370827
9	2	2.4308448	0.5285622	0.1033451	0.0148436	$-$0.0054812	$-$0.0062513
9	3	1.5608761	0.4908939	0.0542856	$-$0.0028530	$-$0.0051481	0.0011737
9	4	0.7660494	0.4741134	0.0241204	$-$0.0085843	$-$0.0026359	0.0036123
9	5	0	0.4691792	0	$-$0.0100241	0	0.0042217
10	1	3.7575287	0.6179386	0.2279701	0.0953997	0.0236601	$-$0.0352918
10	2	2.7037563	0.5168106	0.1082262	0.0167597	$-$0.0049572	$-$0.0066215
10	3	1.8588219	0.4780208	0.0600734	$-$0.0017051	$-$0.0054038	0.0006890
10	4	1.0923667	0.4592937	0.2109507	0.0314529	$-$0.0080085	$-$0.0033338
10	5	0.3607214	0.4513082	0.0098626	$-$0.0102013	$-$0.0011011	0.0041399
11	1	3.9920188	0.6075750	0.2319524	0.0989270	0.0268408	$-$0.0335013
11	2	2.9597112	0.5066156	0.11225752	0.0184367	$-$0.0044339	$-$0.0068800
11	3	2.1358866	0.4670738	0.0647201	$-$0.0006566	$-$0.0055397	0.0002862
11	4	1.3926567	0.4470071	0.0371329	$-$0.0073705	$-$0.0038169	0.0029094
11	5	0.6881186	0.4369907	0.0171500	$-$0.0100138	$-$0.0018800	0.0039388
11	6	0	0.4339162	0	$-$0.0107352	0	0.0042171
12	1	4.2148992	0.5983148	0.2353747	0.1020278	0.0297125	$-$0.0317470
12	2	3.2014467	0.4976369	0.1156525	0.0199163	$-$0.0039259	$-$0.0070582
12	3	2.3957884	0.4575857	0.0685453	0.0002926	$-$0.0056002	$-$0.0000490
12	4	1.6720878	0.4365640	0.0416807	$-$0.0067366	$-$0.0041593	0.0026001
12	5	0.9897011	0.4251242	0.0227820	$-$0.0096782	$-$0.0024465	0.0037082
12	6	0.3278102	0.4199860	0.0072728	$-$0.0108093	$-$0.0008020	0.0041287

We take this opportunity to draw attention to some inequalities associated with the means $\mu=\mu_{N,k}$. In relation to this, we require the fact that upon appropriate Householder similarity transformations, GOE matrices can be demonstrated to be similar to certain symmetric random tridiagonal matrices $T_{N}$ [58]. The $T_{N}$ have independent standard Gaussian entries on the diagonal, and leading off diagonal entries given by the random variables $\{\tilde{\chi}_{N-j}\}_{j=1,\dots,N-j}$, where $\tilde{\chi}_{k}$ is the square root of the usual $\chi_{k}^{2}$ random variable, scaled by $1/\sqrt{2}$. In particular $\mathbb{E}\,T_{N}$ is the symmetric tridiagonal matrix with zero on its diagonal entries and leading off diagonal entries equal to $\{{1\over\sqrt{2}}\sqrt{N-j}\}_{j=1,\dots,N-1}$. From the implied three term recurrence one can check that $\det(x\mathbb{I}_{N}-\mathbb{E}\,T_{N})\propto H_{N}(x)$, where $H_{N}(x)$ denotes the Hermite polynomial of degree $N$, telling us that eigenvalues of $\mathbb{E}\,T_{N}$ are equal to the zeros of this Hermite polynomial.

On the other hand, for general random real symmetric matrices $T_{N}$ a theorem of Cacoullos and Olkin [16, Corollary 2.8] gives that

$$\mathbb{E}\,\lambda_{1}(T_{N})+\cdots+\mathbb{E}\,\lambda_{k}(T_{N})\geq\lambda_{1}(\mathbb{E}\,\mathbb{T}_{N})+\cdots+\lambda_{k}(\mathbb{E}\,T_{N}).$$ (A.1)

By computing the zeros of the Hermite polynomials, and comparing with the $\mu_{N,k}$ in the table, we can readily verify (A.1). In fact in doing so, one observes the stronger interlacing inequality

$$\mathbb{E}\,\lambda_{k}(T_{N})>\lambda_{k}(\mathbb{E}\,\mathbb{T}_{N})>\mathbb{E}\,\lambda_{k-1}(T_{N}),\qquad(k=1,\dots,\lfloor N/2\rfloor)$$ (A.2)

exhibited by the listed $\mu_{N,k}$.

An analogous property is observed in the Laguerre case. The tridiagonal matrix $T_{N}$ similar to a LOE matrix has the property that $\det(x\mathbb{I}_{N}-\mathbb{E}\,\mathbb{T}_{N})\propto L_{N}^{2a-1}(x)$, where $L_{N}^{2a-1}(x)$ denotes the Laguerre polynomial of degree $N$ with parameter $2a-1$ [24, Table 2(a) with $\gamma\mapsto 2a$]. The following table presents the means of the marginal eigenvalue PDFs $\{f_{N}^{\rm L}(k;s)\}_{k=1}^{N}$ obtained by numerical integration of its exact functional form, which is computed according to the formalism in Section 2.2. By comparing the means with the zeros of $L_{N}^{2a-1}(x)$, one can verify the validity of the interlacing (A.2) for each $k=1,\dots,N$.

	$N$	$2$	$3$	$4$	$5$	$6$
$a$	$k$
4	1	15.0634920	19.4986342	23.7003816	27.7874859	31.8100369
4	2	6.9365079	10.9999999	14.7379489	18.3726752	21.9688258
4	3		5.5013657	8.9489690	12.2197905	15.4601777
4	4			4.6127004	7.6273247	10.5502764
4	5				3.9927234	6.6808343
4	6					3.5298487
7/2	1	13.8656315	18.1848768	22.3053570	26.3308207	30.3044360
7/2	2	6.1343684	10.0000000	13.6136374	17.1577682	20.6833885
7/2	3		4.8151231	8.0729350	11.2189206	14.3657513
7/2	4			4.0080704	6.8422317	9.6424231
7/2	5				3.4502585	5.9670831
7/2	6					3.0369177

Table 7: Means of the LOE marginal density $f_{N}^{L}(k;s)$