Moments of random matrices and hypergeometric orthogonal polynomials
Abstract.
We establish a new connection between moments of random matrices and hypergeometric orthogonal polynomials. Specifically, we consider moments as a function of the complex variable , whose analytic structure we describe completely. We discover several remarkable features, including a reflection symmetry (or functional equation), zeros on a critical line in the complex plane, and orthogonality relations. An application of the theory resolves part of an integrality conjecture of Cunden et al. [F. D. Cunden, F. Mezzadri, N. J. Simm and P. Vivo, J. Math. Phys. 57 (2016)] on the timedelay matrix of chaotic cavities. In each of the classical ensembles of random matrix theory (Gaussian, Laguerre, Jacobi) we characterise the moments in terms of the Askey scheme of hypergeometric orthogonal polynomials. We also calculate the leading order asymptotics of the moments and discuss their symmetries and zeroes. We discuss aspects of these phenomena beyond the random matrix setting, including the Mellin transform of products and Wronskians of pairs of classical orthogonal polynomials. When the random matrix model has orthogonal or symplectic symmetry, we obtain a new duality formula relating their moments to hypergeometric orthogonal polynomials.
Contents
 1 Introduction
 2 Motivation and applications
 3 Notation and definitions
 4 Unitary ensembles
 5 Large asymptotics of the spectral zeta functions
 6 Beyond random matrices: Wronskians of orthogonal polynomials
 7 Higher order cumulants
 8 Orthogonal and Symplectic ensembles
 A Orthogonal and symplectic ensembles: duality
 B Mellin transform
 C Hypergeometric orthogonal polynomials
1. Introduction
In this paper we present a novel approach to the moments of the classical ensembles of random matrices. Much of random matrix theory is devoted to moments () of random matrices of finite or asymptotically large size . The Gaussian, Laguerre and Jacobi unitary ensembles have been extensively studied and virtually everything is known about the moments as functions of the matrix size . In particular, for the GUE, is a polynomial in . This fact is a consequence of Wick’s theorem, it is usually called ‘genus expansion’, and it is at the heart of several successful theories such as the topological recursion. For example, the th moment of GUE matrices of size is
In contrast to the wealth of results on moments as functions of the size , less attention has been devoted to them as functions of the order . One of the consequences is that some remarkable properties have been somehow missed. The theory described in this paper is intended to fill this gap. By looking at the moments as functions of , we gain access to additional structure. Several results contained in this paper are in fact facets of the same phenomenon, which appears to be a new observation: moments of classical matrix ensembles, if properly normalized, are hypergeometric orthogonal polynomials as functions of . For example, for a GUE random matrix of size
and this polynomial is actually a Meixner polynomial. In fact, the moments are essentially Meixner polynomials as functions of , too.
During our investigation it became natural to consider complex moments () or, equivalently, averages of spectral zeta functions of random matrices.
1.1. Spectral zeta functions of random matrices
There exist various generalizations of the Riemann function, associated with operator spectra and which are generically called spectral zeta functions.
Consider a compact operator on a separable Hilbert space. Then is a nonnegative operator, so that makes sense. The singular values of are defined as the (nonzero) eigenvalues of . If is selfadjoint with discrete spectrum , the singular values are . The Dirichlet series representation of the Riemann zeta function suggests to define the spectral zeta function of the operator as the maximal analytic continuation of the series
(this is also called Minakshisundaram–Pleijel [57] zeta function of .) In this sense, the Riemann is the spectral zeta of the integer spectrum . Several authors have posed the question of how the ‘spectral’ properties of Riemann’s zeta function carry over (or not) to various spectral zeta functions [68]; classical properties of the Riemann zeta function are:

Functional equation: the function satisfies ;

Meromorphic structure: is analytic in ;

Special values: trivial zeros for ;

Complex zeros of and the Riemann hypothesis (RH): the complex Riemann zeros are in the critical strip and enjoy the reflection symmetries along the real axis and the line . It is conjectured (RH) that the nontrivial zeros all lie on the critical line .
Let be a random Hermitian matrix, and denote by its eigenvalues. Assume that, with full probability, is not in the spectrum of (this is certainly true for the classical ensembles of random matrices). We proceed to define the averaged spectral zeta function as the maximal analytic continuation of
Note that is not a random function. Much of the paper is devoted to pointing out the analytic structure of when comes from the Gaussian, Laguerre or Jacobi ensembles.
1.2. Timedelay matrix of chaotic cavities
Random matrix theory provides a mathematical framework to develop a statistical theory of quantum transport. This theory is believed to apply in particular to mesoscopic conductors confined in space, often referred as quantum dots, connected to the environment through ideal leads. For these systems, Brouwer, Frahm and Beenakker [16], showed that the proper delay times are distributed as the inverse of the eigenvalues of matrices in the Laguerre ensemble (the size being the number of scattering channels) with Dyson index labelling the classical symmetry classes, and parameter .
The moments of the proper delay times have been studied using both random matrix theory [23, 25, 26, 38, 51, 56, 64, 65] and semiclassical scattering orbit theory [62, 12, 45]. Of course, the subject of moments on rather general matrix ensembles has been extensively studied. There is, however, one important complication here: the moments of the timedelay matrix are singular spectral linear statistic on the Laguerre ensemble. (Invariant random matrix ensembles with singular potentials received considerable interest in recent years in mathematical physics, see, e.g. Refs. [3, 9, 13, 20, 54, 15].)
In [23, 25, 26], it was conjectured that the expansion of the cumulants of power traces for the timedelay matrix of quantum dots has positive integer coefficients. In this paper we prove that the conjecture is true for the first order cumulants, i.e. the moments, when (systems without broken time reversal symmetry).
1.3. Mellin transform of orthogonal polynomials
The averaged zeta function is related to the Mellin transform of the onepoint correlation function
where are the normalized Hermite, Laguerre or Jacobi wavefunctions with leading coefficient . The Mellin transform of a function is defined by the integral
when it exists. We set . Of course, . In all instances in this paper, the Mellin transforms have meromorphic extensions to all of with simple poles (see Appendix B.1).
Bump and Ng [18] and Bump, Choi, Kurlberg and Vaaler [19] made the remarkable observation that the Mellin transforms of Hermite and Laguerre functions form families of orthogonal polynomials (OP’s) and have zeros on the critical line . (Jacobi functions were not considered by them.) A few years later, Coffey [21] and Coffey and Lettington [22] pointed out that the polynomials described by Bump et al. were hypergeometric OP’s and investigated other families.
Indeed, we show that for the classical matrix ensembles, the Mellin transform of a Wronskian of two adjacent wavefunctions is a hypergeometric OP (up to a factor containing ratios of Gamma functions). We stress that the proof does not go along the lines of the method of Bump et al.. They started from the orthogonality of the classical wavefunctions which is preserved by the Mellin transform (a unitary operator in ). In our case, by explicit computations, we identify a discrete SturmLiouville (SL) problem satisfied by (as a function of ) and this turns out to be the same SL of the classical hypergeometric OP’s. We remark that the Mellin transforms of do have a probabilistic meaning (moments of random matrices). The Mellin transforms studied in [18, 19], while not unmotivated, do not have an obvious probabilistic interpretation.
Once the analytic structure of the Mellin transform of was established, it became natural for us to look for similar polynomial properties for Wronskians of nonadjacent wavefunctions , . Such Wronskians do not have a random matrix interpretation. Nevertheless, they have a certain interest in mathematical physics as they appear when applying DarbouxCrum [27] transformations on a Schrödinger operator to generate families of exceptional orthogonal polynomials [34, 44, 63, 36].
1.4. Orthogonal and symplectic ensembles
The theory developed for the classical ensembles of complex random matrices suggested to look for similar polynomial properties in the real and quaternionic cases (orthogonal and symplectic ensembles).
Now a fundamental insight came from recursion formulae satisfied by orthogonal / symplectic moments coupled with moments of the corresponding unitary ensembles (see [48] in the Gaussian case and [26] in the Laguerre ensemble). It turns out that for the classical ensembles of random matrices with orthogonal and symplectic symmetries, certain combinations of moments satisfy three term recursion formulae which, again, correspond to the SL equations defining families of hypergeometric OP’s. Therefore, this combination of moments plays the role of the single moments in the unitary case: they satisfy three term recursions, have hypergeometric OP factors, reflection symmetries, zeros on a vertical line, etc.
The (single) moments of the symplectic ensembles do have polynomial factors, but these do not belong to the Askey scheme. In the orthogonal cases, we use a novel duality formula (based on the results by Adler et al. [2]) to write the moments of real random matrices of odd dimension as quaternionic moments plus a remainder containing an orthogonal polynomial factor.
Coupling this result with a classical duality between orthogonal and symplectic ensembles, we discover a functional equation for moments of real random matrices.
1.5. Outline
The paper has the following structure:

In Section 2 the physics motivations and application to quantum transport in chaotic cavities are presented;

In Section 3 we set some notation and we recall the definition of the classical ensembles of random matrices and hypergeometric OP’s;

In Section 4 we present the main results along with their proofs for the Gaussian, Laguerre and Jacobi unitary ensembles;

In Section 5 we discuss the large asymptotics of the spectral zeta functions;

In Section 6 we discuss the relation of our findings with earlier works on the Mellin transform of classical orthogonal polynomials; then, we extend our results beyond random matrix theory by considering Mellin transforms of products and Wronskians of generic pairs of orthogonal polynomials;

In Section 7 we discuss the extension of duality formulae between moments of random matrices to higher order cumulants;

In Section 8 we consider the classical orthogonal and symplectic ensembles and, in particular, present a new duality formula relating their moments to hypergeometric orthogonal polynomials.
2. Motivation and applications
One of the original motivations for this work was to make some progress on an integrality conjecture for the expansion of the Laguerre ensemble put forward in [23, 25, 26]. The problem originated from a random matrix approach to quantum transport in chaotic cavities.
In this section, we will first present some of our findings on the Laguerre ensemble. Then we will briefly review the connection between the Laguerre ensemble and the timedelay matrix in chaotic cavities, and explain the applicability of our results to the integrality conjecture.
2.1. The Laguerre ensemble: reciprocity law and spectral zeta function
Let be a random matrix from the Laguerre Unitary Ensemble (LUE) with parameter . That is, is distributed according to
(2.1) 
on the space of nonnegative Hermitian matrices, where is Lebesgue measure on , the normalizing constant, and .
Consider for integer , the (inverse) moments
(2.2) 
where are the eigenvalues of . The above moments are finite if and only if [46]. We will prove the following remarkable property of these moments.
Proposition 2.1 (Reciprocity law for LUE).
(2.3) 
The above identity can be verified using the recurrence relation for moments proved by Haagerup and Thorbjørnsen [39] and extended to inverse moments in [26]. This gives an explicit formula for the inverse moments given known identities for positive moments. For instance,
It is natural to consider complex moments or, equivalently, the averaged LUE spectral zeta function defined as
(2.4) 
and by analytic continuation for other values of . We list below a few remarkable properties of the averaged LUE spectral zeta function.

Analytic structure: it turns out that can be analytically extended to the whole complex plane; In particular, is a polynomial of degree .

Special values: trivial zeros for ;

Complex zeros and Riemann hypothesis: as for the Riemann zeta function, the set of complex zeros is symmetric with respect to reflections along the real axis and the critical line . It is tempting to ask whether a RH holds true for the averaged LUE zeta function. Amusingly, the answer is ‘Yes’: the nontrivial zeros of all lie on the critical line .
These facts are an immediate consequences of the main results presented in Section 4. See Fig. 1 for a plot of the averaged LUE zeta function.
Remark.
For any fixed , the function is a finite sum of exponentials. Therefore, without taking the average, is a random analytic function in and never vanishes.
2.2. Application to quantum transport in chaotic cavities
In [23, 25], it was proposed that, for , the cumulants of the timedelay matrix of a ballistic chaotic cavity have a expansion with positive integer coefficients (similar to the genus expansion of Gaussian matrices). The precise conjectural statement is as follows.
Consider the measure (2.1) with , and the rescaled inverse power traces
It is known that the expectation of has a expansion
Conjecture ([26]).
.
The conjecture was supported by a systematic computation of certain generating functions, and it is in agreement with the diagrammatic expansions of scattering orbit theory. The integrality of the coefficients in the large expansion has been also conjectured in the real case (LOE) [26], and for higher order cumulants [23, 25].
The results reported in this paper resolve the conjecture in the complex case.
Theorem 2.2.
The above conjecture is true.
Proof.
To prove the Theorem we make advantage of the reciprocity law to use known results for positive moments of the Laguerre ensemble.
Let be in the LUE with parameter . For , from [40, Corollary 2.4] (see also [58, Exercise 12]) we read the formula
(2.7) 
where is the symmetric group, and for a permutation , denotes the number of cycles in . By we denote the cycle . If with , the above formula shows that is a polynomial in with positive coefficients (see Lemma 4.5 below).
Remark.
From (2.7) we see that, if , then ( integer) has a expansion with positive integer coefficients. The computation above shows that the integrality of the coefficients for the LUE negative moments also holds whenever . Therefore, is the only case when all moments () have integer coefficients in their large expansion.
3. Notation and definitions
3.1. Classical ensembles of random matrices
We will consider expectations with respect to the measures
(3.1) 
for finite and for any value of . The value of corresponds to ensembles of real symmetric (), complex hermitian () or quaternion selfdual matrices (). The function is the weight of the ensemble:
(3.2) 
for Gaussian, Laguerre, and Jacobi, respectively. is a normalization constant which may vary at each occurrence. For convenience we set in the Laguerre ensemble and in the Jacobi ensemble.
3.2. Hypergeometric orthogonal polynomials
We use the standard notation for hypergeometric functions
where . We need to introduce some families of hypergeometric OP’s. Recall that there are three types of hypergeometric OP’s:

Polynomials of the first type are solutions of usual SL problems for second order differential operators: Jacobi and its degenerations, Laguerre and Hermite . They have a hypergeometric representation, but they are perhaps better known by the Rodriguestype formulae
(3.5) (3.6) (3.7) Note that the Jacobi polynomials considered in this paper are orthogonal with respect to the measure on the unit interval ;

Polynomials of the second type are solutions of discrete SL problems (threeterms recurrence relations) with real coefficients: Racah , including its degenerations Hahn , dual Hahn , Meixner , etc. They can be represented as finite hypergeometric series
(3.8) (3.9) (3.10) (3.11) where . Note that some authors define the Meixner polynomials as ;

Polynomials of the third type are solutions of discrete SL problems with complex coefficients: Wilson including its degenerations, continuous dual Hahn , continuous Hahn , MeixnerPollaczek , etc. They have the following hypergeometric representations:
(3.12) (3.13) (3.14) (3.15)
At the top of the hierarchy of hypergeometric orthogonal polynomials are the Wilson and Racah polynomials. In this paper, the parameter ranges are such that it is most natural to consider the polynomials which appear as Wilson polynomials and their degenerations, specifically continuous dual Hahn and MeixnerPollaczek polynomials. The reader can find the common notation and the most important properties of hypergeometric OP’s in [43, Section 9].
4. Unitary ensembles
It is known that the th moments of the classical unitary invariant ensembles of random matrices of dimension are polynomials in (or in after rescaling). Here we show that the (completed) moments can also be seen as polynomials in the parameter . These polynomials are hypergeometric orthogonal polynomials OP’s belonging to the Askey scheme [7]. The polynomial property suggests to consider complex moments or, equivalently, averages of spectral zeta functions. For the three classical ensembles, define
and
when the expectations exist (, , and for , , and , respectively) and by analytic continuation otherwise (see Appendix B.1).
We can now state the first result.
Theorem 4.1.
For all , is a hypergeometric orthogonal polynomial:
where . In particular, satisfies the functional equation in the LUE and JUE cases, and for the GUE. Moreover, all its zeros lie on the critical line .
The weights of the OP’s in Theorem 4.1 are
For the GUE and LUE, the following orthogonality conditions hold
(4.1) 
where is an explicit constant depending on the ensemble (see Appendix C).
Remark.
The orthogonality in the JUE is slightly different to the GUE and LUE cases. First of all, the fourth parameter of the Wilson polynomial is negative. An orthogonality relation in this case is far from obvious and was established by Neretin [61], see also Appendix C. This fourth parameter also depends on , therefore each belongs to a distinct family of orthogonal polynomials obtained by fixing the fourth parameter. As before, this orthogonality implies that the zeros lie on the line .
4.1. Gaussian Unitary Ensemble
The GUE is a classical orthogonal polynomial ensemble. In particular, the correlation functions can be compactly and conveniently written in terms of Hermite polynomials. It turns out that the complex moments are essentially a MeixnerPollaczek polynomial. The moments for GOE and GSE can be expressed using known formulae relating the onepoint correlation functions of the three Gaussian ensembles.
Let be a GUE random matrix of dimension . Define for all for which the expectation exists. It is known [41] that, for ,
(4.2) 
where (this is equal to for integer). For each , the moment is a polynomial in with positive integer coefficients:
This fact is the wellknown genus expansion for Gaussian complex matrices.
As observed in [69, Theorem 8], (4.2) can be written in terms of a (terminating) hypergeometric series:
(4.3) 
From this hypergeometric representation, we see that the moment , if properly normalised, is a MeixnerPollaczek polynomial in of degree .
Theorem 4.2.
If we write , then for
(4.4) 
In particular, can be extended to an analytic function in (a polynomial), invariant up to a change of sign under reflection , with complex zeros on the vertical line .
Proof.
Consider the polynomials
From the definition of MeixnerPollaczek polynomials (3.15) of and the hypergeometric representation (4.3) we see that when is a nonnegative integer. In order to prove the general complex case, we use a procedure of analytic continuation from integer points to a complex domain via Carlson’s theorem [6, Theorem 2.8.1]. A standard calculation in random matrix theory shows that, for ,
where denotes the Hermite polynomial (3.5) of degree . This shows that is analytic in the halfplane . If we write for some constants , then
We use now the elementary inequality for . Setting (with ) we have then
Therefore as with . In conclusion, and the polynomial coincide on nonnegative integers and their difference is for any . By Carlson’s theorem the two functions coincide in the whole domain .
The polynomials enjoy the symmetries
Therefore thus explaining the reflection symmetry . Recall that the MeixnerPollaczek polynomial form an orthogonal family with respect to a positive weight on the real line. Therefore, the zeros of are purely imaginary; they occur in conjugate pairs, with zero included if is odd. This proves that all the zeros of lie on the line . ∎
Remark.
The polynomials satisfy the difference equation
and the threeterm recurrence
Since , these are in fact equivalent.
The MeixnerPollaczek polynomials can be thought as continuous version of the Meixer polynomials [8]:
In fact, the normalized moment (4.3) is a Meixner polynomial (see (3.11)) in of degree or, by symmetry, a Meixner polynomial in of degree . The alternative form of Theorem 4.2 using Meixner polynomials is the following.
Theorem
(4.7) 
∎
The first polynomials are