Moments of Random Matrices and Hypergeometric Orthogonal Polynomials

We establish a new connection between moments of n×n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${n \times n}$$\end{document} random matrices Xn and hypergeometric orthogonal polynomials. Specifically, we consider moments ETrXn-s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb{E}{\rm Tr} X_n^{-s}}$$\end{document} as a function of the complex variable s∈C\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${s \in \mathbb{C}}$$\end{document} , 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. (J Math Phys 57:111901, 2016) on the time-delay 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 n→∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${n \rightarrow \infty}$$\end{document} 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.


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 E Tr X k n (k ∈ N) of random matrices of finite or asymptotically large size n. 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 n. In particular, for the GUE, E Tr X k n is a polynomial in n. 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 [5,33]. For example, the 4-th moment of GUE matrices of size n is 1 n E Tr X 8 n = 14n 4 + 70n 2 + 21.
In contrast to the wealth of results on moments as functions of the size n, less attention has been devoted to them as functions of the order k. 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 k, 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 E Tr X k n of classical matrix ensembles, if properly normalized, are hypergeometric orthogonal polynomials as functions of k. For example, for a GUE random matrix of size n = 4 1 (2k − 1)!! E Tr X 2k 4 = 4 3 k 3 + 4k 2 + 20 3 k + 4, and this polynomial is actually a Meixner polynomial. In fact, the moments are essentially Meixner polynomials as functions of (n − 1), too. During our investigation it became natural to consider complex moments E Tr X k n (k ∈ C) or, equivalently, averages of spectral zeta functions of random matrices.

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 A on a separable Hilbert space. Then A A * is a nonnegative operator, so that |A| = A A * makes sense. The singular values of A are defined as the (nonzero) eigenvalues of |A|. If A is self-adjoint with discrete spectrum λ 1 , λ 2 , . . . , the singular values are |λ 1 |, |λ 2 |, . . . . The Dirichlet series representation of the Riemann zeta function ζ(s) suggests to define the spectral zeta function ζ A (s) of the operator A as the maximal analytic continuation of the series j≥1 |λ j | −s (this is also called Minakshisundaram-Pleijel [67] zeta function of A). In this sense, the Riemann ζ(s) is the spectral zeta of the integer spectrum λ j = j. 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 [79]; classical properties of the Riemann zeta function are: Let X n be a n ×n random Hermitian matrix, and denote by λ 1 , . . . , λ n its eigenvalues. Assume that, with full probability, 0 is not in the spectrum of X n (this is certainly true for the classical ensembles of random matrices). We proceed to define the averaged spectral zeta function ζ X n (s) as the maximal analytic continuation of E ζ X n (s) = E Tr |X n | −s = E n j=1 |λ j | −s .
Note that E ζ X n (s) is not a random function. Much of the paper is devoted to pointing out the analytic structure of E ζ X n (s) when X n comes from the Gaussian, Laguerre or Jacobi ensembles.

Time-delay 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 [17], showed that the proper delay times are distributed as the inverse of the eigenvalues of matrices X n in the Laguerre ensemble (the size n being the number of scattering channels) with Dyson index β ∈ {1, 2, 4} labelling the classical symmetry classes, and parameter α = n.
The moments of the proper delay times have been studied using both random matrix theory [25,27,28,43,44,58,[64][65][66]75,76] and semiclassical scattering orbit theory [13,52,72,73]. 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 time-delay 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,14,16,21,62].) In [25,27,28], it was conjectured that the 1/n-expansion of the cumulants of power traces for the time-delay 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 β = 2 (systems without broken time reversal symmetry).

Mellin transform of orthogonal polynomials.
The averaged zeta function is related to the Mellin transform of the one-point correlation function 1 . In the classical unitary invariant ensembles, by using the well-known determinantal formulae and Christoffel-Darboux formula, the one-point correlation function can be written as a Wronskian ρ (2) n (x) = where ψ j (x) are the normalized Hermite, Laguerre or Jacobi wavefunctions with leading coefficient k j . (The superscript stands for β = 2.) The Mellin transform of a function f (x) is defined by the integral In all instances in this paper, the Mellin transforms have meromorphic extensions to all of C with simple poles (see Appendix B.1).
Bump and Ng [19] and Bump, Choi, Kurlberg, and Vaaler [20] made the remarkable observation that the Mellin transforms of Hermite and Laguerre functions ψ j (x) form families of orthogonal polynomials (OP's) and have zeros on the critical line Re(s) = 1/2. (Jacobi functions were not considered by them.) A few years later, Coffey [22,23] and Coffey and Lettington [24] 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 ρ (2), * n (s) 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 L 2 ). In our case, by explicit computations, we identify a discrete Sturm-Liouville (S-L) problem satisfied by ρ (2), * n (s) (as a function of s) and this turns out to be the same S-L of the classical hypergeometric OP's. We remark that the Mellin transforms ρ (2), * n (s) 1 See Section 3 for the definition of the one-point function and the important identity (3.4). of ρ (2) n (x) do have a probabilistic meaning (moments of random matrices). The Mellin transforms ψ * j (s) studied in [19,20], while not unmotivated, do not have an obvious probabilistic interpretation.
Once the analytic structure of the Mellin transform of Wr(ψ n (x), ψ n+1 (x)) was established, it became natural for us to look for similar polynomial properties for Wronskians of nonadjacent wavefunctions Wr(ψ n (x), ψ n+k (x)), k > 1. Such Wronskians do not have a random matrix interpretation. Nevertheless, they have a certain interest in mathematical physics as they appear when applying Darboux and Crum [29] transformations on a Schrödinger operator to generate families of exceptional orthogonal polynomials [39,41,51,74].

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 [55] in the Gaussian case and [28] 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 S-L 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-n 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.

Motivation and Applications
One of the original motivations for this work was to make some progress on an integrality conjecture for the 1/n-expansion of the Laguerre ensemble put forward in [25,27,28]. 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 time-delay 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 X n be a random matrix from the Laguerre Unitary Ensemble (LUE) with parameter m ≥ n. That is, X n is distributed according to on the space P n of nonnegative Hermitian matrices, where d X is Lebesgue measure on P n R n 2 , Z the normalizing constant, and α = m − n.
Consider for integer k ∈ N, the (inverse) moments where λ 1 , . . . , λ n are the eigenvalues of X n . The above moments are finite if and only if k ≤ α [53]. We will prove the following remarkable property of these moments.

Proposition 2.1 (Reciprocity law for LUE).
E Tr X −(k+1) The above identity can be verified using the recurrence relation for moments proved by Haagerup and Thorbjørnsen [45] and extended to inverse moments in [28]. 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 so that (2.3) becomes the functional equation (2) Analytic structure: it turns out that E ζ X n (s) can be analytically extended to the whole complex plane; In particular, ξ n (s) is a polynomial of degree 2(n − 1). (3) Special values: trivial zeros E ζ X n (1 + α + j) = 0 for j = 1, 2, . . . ; (4) 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 Re(s) = 1/2. 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 E ζ X n (s) all lie on the critical line Re(s) = 1/2.
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 and Fig. 2 for an illustration of the zeros.
Remark. For any fixed n, the function ζ X n (s) is a finite sum of exponentials. Therefore, without taking the average, ζ X n (s) is a random analytic function in C and never vanishes.

Application to quantum transport in chaotic cavities.
In [25,27], it was proposed that, for β ∈ {1, 2}, the cumulants of the time-delay matrix of a ballistic chaotic cavity have a 1/n-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 α = n, and the rescaled inverse power traces It is known that the expectation of τ k (n) has a 1/n-expansion 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-n expansion has been also conjectured in the real case (LOE) [28], and for higher order cumulants [25,27].
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 take advantage of the reciprocity law to use known results for positive moments of the Laguerre ensemble. Let X n be in the LUE with parameter α = m − n. For k ≥ 0, from [46, Corollary 2.4] (see also [68,Exercise 12]) we read the formula where S k is the symmetric group, and for a permutation σ ∈ S k , #(σ ) denotes the number of cycles in σ . By γ k we denote the k-cycle (1 2 3 . . . k). If m = cn with c > 0, the above formula shows that 1 n k+1 E Tr X k n is a polynomial in n −2 with positive coefficients (see Lemma 4.5 below).
By the reciprocity law (2.3) with α = n (or, equivalently, c = 2) The factor k j=1 1 − j 2 n 2 −1 in (2.8) is a product of geometric series. Therefore, we have (2.9) and this readily prove that Eτ k has an expansion in n −2 with positive integer coefficients.
Remark. From (2.7) we see that, if c ∈ N, then 1 n k+1 E Tr X k n (k ≥ 0 integer) has a 1/n-expansion with positive integer coefficients. The computation above shows that the integrality of the coefficients for the LUE negative moments also holds whenever c/(c − 1) ∈ N. Therefore, c = 2 is the only case when all moments 1 n k+1 E Tr X k n (k ∈ Z) have integer coefficients in their large-n expansion.

Notation and Definitions
3.1. Classical ensembles of random matrices. We will consider expectations with respect to the measures for finite n and for any value of β ∈ {1, 2, 4}. The value of β corresponds to ensembles of real symmetric (β = 1), complex hermitian (β = 2) or quaternion self-dual matrices (β = 4). The function w β (x) is the weight of the ensemble: for Gaussian, Laguerre, and Jacobi, respectively. C n,β is a normalization constant which depends on the ensemble and is known explicitly [36]. For convenience we set α = m − n > 0 in the Laguerre ensemble and α 1 = m 1 − n > 0, α 2 = m 2 − n > 0 in the Jacobi ensemble. We define the one-point eigenvalue density ρ (β) n (x) corresponding to (3.1) by We will call (3.3) the eigenvalue density corresponding to the weight w β (x) defining the expectation over (3.1). The following identity easily follows from the definitions (3.3) and (3.1):

3.2.
Hypergeometric orthogonal polynomials. We use the standard notation for hypergeometric functions where (q) n = (q + n)/ (q). We need to introduce some families of hypergeometric OP's [50]. Recall that there are three types of hypergeometric OP's: (1) Polynomials of the first type are solutions of usual S-L problems for second order differential operators: Jacobi P (α 1 ,α 2 ) n (x) and its degenerations, Laguerre L (α) n (x) and Hermite H n (x). They have a hypergeometric representation, but they are perhaps better known by the Rodrigues-type formulae Note that the Jacobi polynomials considered in this paper are orthogonal with respect to the measure (1 − x) α 1 x α 2 dx on the unit interval [0, 1]; (2) Polynomials of the second type are solutions of discrete S-L problems (three-terms recurrence relations) with real coefficients: Racah R n (λ(x); α, β, γ , δ), including its degenerations Hahn Q n (x; α, β, N ), dual Hahn R n (λ(x); γ, δ, N ), Meixner M n (x; β, c), etc. They can be represented as finite hypergeometric series 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 Meixner-Pollaczek polynomials. The reader can find the common notation and the most important properties of hypergeometric OP's in [50,Section 9].

Unitary Ensembles
It is known that the k-th moments of the classical unitary invariant ensembles of random matrices X n of dimension n are polynomials in n (or in 1/n after rescaling). Here we show that the (completed) moments can also be seen as polynomials in the parameter k. 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 when the expectations exist (s < 1/4, s < α + 1, and s < α 2 + 1 for GUE, LUE, and JUE, respectively) and by analytic continuation otherwise (see Appendix B.1).
We can now state the first result.
For the GUE and LUE, the following orthogonality conditions hold 1 2πi where h m is an explicit constant depending on the ensemble (see Appendix C).
For an illustration of the zeros on the critical line, see Figs. 2 and 3. For the reader's convenience, the relation between moments of the unitary ensembles and hypergeometric OP's is summarised in Table 1.
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 [71], see also Appendix C. This fourth parameter also depends on n, therefore each ξ n (s) 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 Re(s) = 1/2.

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 Meixner-Pollaczek polynomial. The moments for GOE and GSE can be expressed using known formulae relating the one-point correlation functions of the three Gaussian ensembles.
Let X n be a GUE random matrix of dimension n. Define Q C k (n) = E Tr X 2k n for all k ∈ C for which the expectation exists. It is known [47] that, for k ∈ N, is a polynomial in n with positive integer coefficients: This fact is the well-known genus expansion for Gaussian complex matrices. As observed in [81,Theorem 8], (4.2) can be written in terms of a (terminating) hypergeometric series: From this hypergeometric representation, we see that the moment Q C k (n), if properly normalised, is a Meixner-Pollaczek polynomial in i(k + 1) of degree n − 1.
In particular, i n−1 Q C k (n)/(2k − 1)!! can be extended to an analytic function in C (a polynomial), invariant up to a change of sign under reflection k → −2 − k, with complex zeros on the vertical line Re(k) = −1.

Proof. Consider the polynomials
From the definition of Meixner-Pollaczek polynomials (3.15) of P (1) n−1 (x; π/2) and the hypergeometric representation (4.3) we see that Q C k (n) = n (2k − 1)!! q k (n) when k 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 Re(k) > −1/2, where H j denotes the Hermite polynomial (3.5) of degree j. This shows that We use now the elementary inequality a +by +cy 2 ≤ (a +b)+(b+c)y 2 for a, b, c, y ≥ 0.
; π/2) coincide on nonnegative integers and their difference is O(e c|k| ) for any c > 0. By Carlson's theorem the two functions coincide in the whole domain Re(k) > −1/2.
The polynomials q r (s) enjoy the symmetries Recall that the Meixner-Pollaczek polynomial form an orthogonal family with respect to a positive weight on the real line. Therefore, the zeros of q r (s) are purely imaginary; they occur in conjugate pairs, with zero included if r is odd. This proves that all the zeros of i n−1 Q C k (n)/(2k − 1)!! lie on the line Re(k) = −1.
Remark. The polynomials q r (s) satisfy the difference equation and the three-term recurrence Since q r −1 (s) = q s−1 (r ), these are in fact equivalent.
Recalling (4.5) and the recursion in n The Meixner-Pollaczek 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 n − 1 of degree k or, by symmetry, a Meixner polynomial in k of degree n − 1. The alternative form of Theorem 4.2 using Meixner polynomials is the following.
The first polynomials are Remark. The normalised GUE moments can be written as products of moments (2k−1)!! of a standard Gaussian, times a Meixner polynomial It is natural to ask whether the Meixner polynomials form moment sequences of some random variables, so that one can 'decompose' the GUE one-point function as multiplicative convolution of a standard Gaussian and another probability distribution (product of two independent random variables). In fact, Ismail and Stanton [48,49] considered the problem of orthogonal polynomials as moments. It turns out that the Meixner polynomials are moments of translated Beta random variables for Re(x) < 0 and Re(x + β) > 0. Note however that, in our setting, the Meixner polynomials have nonnegative argument x = n − 1, so that this representation of the one-point function as a 'convolution' is purely formal.
As a corollary of Theorem 4.2 we have the following two identities.
Proof. The reflection formula follows from the hypergeometric representation in (4.3).
To prove (4.9) we start from the convolution property of the Meixner-Pollaczek poly- It follows that To complete the proof we use the Forward Shift Operator for the Meixner-Pollaczek polynomials [50, Sec. 9, Eq. (9.7.6)], the reflection formula (4.8), and the recursion (4.6).
Remark. The reflection formula . Let X n be a LUE random matrix with parameter m. Denote α = m − n ≥ 0. Set Q C k (m, n) = E Tr X k n for all k ∈ C for which the expectation exists. Then Q C 0 (m, n) = n and, for k ∈ N, it is known that [46] is a symmetric polynomial in m, n of degree k + 1, with positive integer coefficients: In fact, for each positive integer n, Q C k (m, n) is a polynomial in k of degree 2(n−1). After some manipulations, the moments (4.10) can be expressed in terms of a hypergeometric function: This formula can be extended for k ∈ C and satisfies Q C 0 (m, n) = n.
where S n−1 denotes the continuous dual Hahn polynomial of degree n − 1. In particular this shows that Q C k (m, n)/(k + α)! can be extended to a polynomial invariant under the reflection k → −1 − k (reciprocity law) and, moreover, its complex zeros lie on the critical line Re(k) = −1/2.
Proof. Comparing (4.11) with the hypergeometric representation of continuous dual Hahn polynomials (3.13) we get the result for k integer. The extension to complex k is again an application of Carlson's theorem.
An alternative formulation in terms of Hahn and dual Hahn polynomials (3.9)-(3.10) is as follows. (4.14) Remark. The difference equations / three-term recurrence relations for these polynomials (see [50,Sections 9.5 and 9.6]) yield the Haagerup-Thorbjørnsen recursion [28,45] ) If k is a positive integer, and we treat α as a parameter, then Q C k (m, n) is a polynomial in n of degree k + 1. Moreover, we can write where p k−1 denotes the continuous Hahn polynomial (3.14) of degree k − 1.
Note the alternative formula: If α is fixed and we write Q C k (m, n) = a n R n , where a n = n(n + α), then (k − 1)(k + 2)R n = a n+1 R n+1 − (a n+1 + a n−1 )R n + a n−1 R n−1 .
It is very natural to consider m dependent on n. An interesting situation (from the point of view of the large-n limit) is the case m = cn with c > 0 and fixed. The next result shows that the moments, as functions of n, are polynomials with all zeros on a vertical line in the complex plane.
For any permutation σ ∈ S k , one has Hence and so Theorem 4.6. Fix c > 0. The zeros of the polynomials Q C k (cn, n) as a function of n are purely imaginary and satisfy the interlacing property.
Proof. Let q k (n) = Q C k (cn, n)/n. Then for each k, q k (n) is a polynomial of degree k, with positive coefficients, and only powers n k , n k−2 , . . . (see Lemma 4.5). It follows that if we define p k (x) = i k q k (−i x), then p k (x) is a polynomial of degree k, with alternating signs, and satisfies the (Haagerup-Thorbjørnsen) recursion It now follows from [57, Corollary 2.4] that { p k (x)} is a 'Sturm sequence' of polynomials. Hence the p k 's have only real zeros, and they satisfy the interlacing property.

Jacobi unitary ensemble.
Let X n be a JUE random matrix of size n with parameters α 1 = m 1 − n, α 2 = m 2 − n. It turns out that the suitable statistics in this ensemble are differences of consecutive moments Q C k (α 1 , α 2 , n), defined as for all k ∈ C for which the expectations exist.
In this case, our strategy is to look for a three-term recursion for Q C k (α 1 , α 2 , n) when k is an integer. In fact, adapting a method due to Ledoux [54, Eq. (30)-(31)], we find the following recurrence relation for the JUE which is the analogue of the Harer-Zagier and Haagerup-Thorbjørnsen recursions.
The coefficient R k , S k , and T k are given by Proof of Theorem 4.7. The proof is immediate when k is a nonnegative integer by observing that (4.20) is the discrete S-L problem for Wilson polynomials, and by checking the initial conditions. For k complex we can use the same method of Theorem 4.2.
Remark. Using the hypergeometric representation of Wilson polynomials (3.12) we have the explicit formula To our knowledge, this hypergeometric representation is new.
Remark. Ledoux [54] obtained a fourth order recursion for moment differences of the Jacobi ensemble, but the ensemble he considers is shifted compared to ours. In our notation, the moments L(k) considered in [54] can be written as for which it is shown that L(k) − L(k + 2) satisfies a fourth order recursion. Using It follows that Ledoux's moment differences can be expressed as a linear combination of hypergeometric functions.
The difference of moments Q C k (α 1 , α 2 , n) can alternatively be written in terms of Racah polynomials.
Using the recurrence relation (4.20), one can verify the following identity between positive and negative moments (this is the analogue of (2.3)).

Generating functions.
It is sometimes convenient to define the moments of random matrices in terms of their generating function. The first example of such a generating function was constructed by Harer and Zagier for the GUE of fixed size n (see Eq. (4.25) below). This convergent series is a rational function. As emphasised by Morozov and Shakirov [69], from the point of view of random matrices and enumeration problems, this is a highly non-trivial result: a generating function for moments at all genera appears to be rational. The generating function of covariances of the GUE computed in [69] turns out to be again an elementary function. The generating function of higher order cumulants of the GUE have been studied recently by Dubrovin and Yang [31] who expressed them in terms of traces of 2 × 2 matrix-valued series.
One of the advantages of the representation of the moments in terms of hypergeometric OP's discussed in the present work, is that we can write explicit formulae for the generating functions of the moments of GUE and LUE for fixed n and/or k. Remarkably, these closed expressions are elementary functions. .
The formula for generating function (4.25) follows from the identity The generating function (4.26) for fixed k, is a direct consequence of the representation of the moments in terms of Meixner-Pollaczek polynomials (4.4) and their generating function [50,Eq. (9.7.11)].
Finally, the joint generating function (4.27) is the resummation in n of (4.25) (or the resummation in k of (4.26)).
For the LUE we use the generating series of continuous (dual) Hahn polynomials. From the representation of the LUE moments as continuous Hahn (4.16), using [50, Eq. (9.5.11)] we get Note that the hypergeometric functions on the right-hand side are terminating series. In fact, they are Laguerre polynomials [50, Eq. (9.12.1)], thus proving (4.28). For (4.29) we use the representation in terms of continuous dual Hahn polynomials (4.12) and the formula of the generating series [50, Eq. (9.3.11)]. We have which is equal to (4.29) by Euler's tranformation. Note that the hypergeometric series is terminating. In fact, one could also write it in terms of Jacobi polynomials. The joint generating series in (4.30) is a resummation of (4.28) over n and m using the known formula for the generating function of Laguerre polynomials [50, Eq. (9.12.10)].
Remark. The series (4.25) was computed by Harer and Zagier [47] using different methods. It is surprising that, although the three-term recurrence in k and the generating function were known, nobody recognized the moments of the GUE as Meixner polynomials. The generating function of the GUE for fixed k (Eq. (4.26)) does not seem to appear in the previous literature. The joint series (4.27) appears in the work of Morozov and Shakirov [69] who stressed the nontrivial fact that it is a rational function in both variables.
The generating functions (4.28)-(4.29)-(4.30) for the LUE seem to be new. It is remarkable, again, that these series sum to elementary functions.

Large-n Asymptotics of the Spectral Zeta Functions
It is a classical result that, after rescaling, the one-point function ρ (β) n (x) of the random matrix ensembles considered in this paper weakly converges to a compactly supported probability measure, as n goes to infinity. The limit ρ ∞ (x) is known as equilibrium measure of the ensemble and does not depend on the Dyson index β. In formulae, for all k ∈ N, This suggests to define the limit zeta function ζ ∞ (s) as the analytic continuation of |x| −s ρ ∞ (x)dx. The limit zeta functions for the classical ensembles turn out to be meromorphic functions. For the LUE and JUE, ζ ∞ (s) has infinitely many nontrivial zeros, and they all lie on a critical line.
We discuss the three classical ensembles separately. For notational convenience we consider the GUE, LUE and JUE, although the results hold true for any β-ensemble.

Gaussian ensemble.
The equilibrium measure is given by the semicircular law After a suitable rescaling, in the large-n limit, the integer moments converge to the Catalan numbers: This formula can be analytically continued and suggests to define the limit GUE zeta function as the analytic continuation of |x| −s ρ ∞ (x)dx: This function has alternating simple poles and zeros on the positive integers, with no other zeros in the rest of the complex plane. The large-n limit is more interesting for matrices in Laguerre and Jacobi ensembles.

Laguerre ensemble.
Let X n be in the Laguerre unitary ensemble. Set α = m − n = (c − 1)n, with c ≥ 1. Define the equilibrium measure Define the limit LUE zeta function as . From the finite-n functional equation ξ n (s) = ξ n (1 − s), we see that a good definition for the limit n → ∞ is Then apply Euler's transformation formula to show the functional equation.
To show that the zeros of ξ ∞ (s) are on the critical line we use an argument based on Sturm-Liouville theory (we borrowed this argument from a similar problem in a paper by Biane [15]).

Jacobi ensemble.
Let X n be in the JUE. Set α 1,2 = (c 1,2 − 1)n, with c 1,2 ≥ 1. Then, the equilibrium measure is Define the limit JUE zeta function as Proof. Using Euler's integral formula, we have The proof of Proposition 5.1 is easily adapted. [19] and Bump, Choi, Kurlberg and Vaaler [20] made the remarkable discovery that the Mellin transforms of Hermite and Laguerre functions have zeros on the critical line Re(s) = 1/2. Their proof is based on the observation that the Mellin transform preserves orthogonality. Hence, Mellin transforms of orthogonal polynomials are themselves orthogonal with respect to some inner product. Later [22][23][24] it was noticed that those orthogonal functions are hypergeometric OP's (multiplied by some nonnegative integrable weight). Consider, for concreteness, the Hermite polynomials H n (x) and the normalised Hermite wavefunctions φ n (x) = (2 n n! √ π) −1/2 e −x 2 /2 H n (x). The following proposition follows from the result of Bump et al. [19,20].

Proposition 6.1 (Mellin transform of Hermite functions). For all integers n
The averaged spectral zeta function of unitary invariant ensembles of random matrices can be interpreted as Mellin transform of Wronskians of adjacent Hermite, Laguerre or Jacobi wavefunctions. Given our results, it is natural to ask whether more general Wronskians have the property that their Mellin transforms can be written in terms of hypergeometric OP's.

Mellin transforms of products and Wronskians of classical orthogonal polynomials. In this section we will use repeatedly the properties (B.2)-(B.3)-(B.4) of the Mellin transform.
The Wronskian of smooth functions f 1 (x) . . . , f m (x) is defined as .
The fact that M xω n, (x); s = −sω * n, (s) follows from integration by parts (or property (B.4) of the Mellin transform, with m = 1). We have proved (6.4). To complete the proof of (6.2) we compute the initial conditions where we have used the identity (6.6). Taking the Mellin transform of both sides, substituting the identity (6.2), and using the Forward Shift Operator for the Meixner-Pollaczek polynomials [50, Sec. 9, Eq. (9.7.6)] we complete the proof of (6.3).
Using the very same method one can show that the Mellin transform of products and Wronskians of two Jacobi wavefunctions is essentially a Wilson polynomial. The proof follows the same lines as the Hermite and Laguerre cases and is omitted. (α 1 + α 2 + 2n + 1) (α 1 + α 2 + n + 1) n! (α 1 + n + 1) (α 2 + n + 1) × (α 1 + α 2 + 2(n + ) + 1) (α 1 + α 2 + n + + 1) (n + )! (α 1 + n + + 1) (α 2 + n + + 1) × (−1) n (α 1 + n + + 1) (s) (α 2 + s) (s − ) (α 1 + α 2 + 2n + + 1) We can also calculate the Mellin transform of a single Jacobi wavefunction in terms of a continuous Hahn polynomial. The interesting case turns out to be for weights with slightly shifted parameters. Theorem 6.6. Consider the functions (6.14) Then the Mellin transform can be written in terms of continuous Hahn polynomials: Proof. To identify the continuous Hahn polynomial, one employs the standard expansion of Jacobi polynomials and calculates the Mellin transform integrating term by term. The result is a hypergeometric sum which can be identified with definition 3.13, leading to (6.15). The difficulty in establishing the conclusion is that these polynomials are only known to be orthogonal when the parameters have positive real part. We proceed by expressing the right-hand side of (6.15) in terms of Wilson polynomials and apply a result of Neretin [71]. We claim that for generic parameters a, b ∈ C, we have up to a constant independent of x. If a and b have positive real part, this follows by writing the orthogonality condition for the Wilson polynomials with the parameters given in (6.16) or (6.17). Use of the duplication formula for the Gamma function shows that the weight reduces to | (a + i x)| 2 | (b + i x)| 2 which is the weight function for continuous Hahn polynomials. If a or b have negative real part the identity follows by analytic continuation. Now inserting (6.16) and (6.17) into (6.15) and applying the orthogonality (C.2) completes the proof.
Going back to random matrices, we obtain exact ordinary differential equations for the one-point functions of the classical ensembles. Denote by ρ (2) n (x) the one-point function of the ensembles GUE, LUE, or JUE. As already discussed in the Introducion, ρ (2) n (x) can be represented as the Wronskian of two adjacent wavefunctions.
The following proposition is a corollary of the previous theorems on Mellin transforms of Wronskians. For GUE and LUE we recover a result obtained and used earlier by Götze (2) n (x) satisfies the differential equation Dρ (2) Proof. We present the proof for the GUE. The one-point function ρ (2) n (x) of the GUE is proportional to the Wronskian of consecutive Hermite functions W n−1, ). To prove the Theorem for the GUE it is therefore sufficient to show that the function W n,1 (x) satisfies the third-order differential equation Remark. It is natural to ask whether Wronskians of nonadjacent wavefunctions satisfy similar differential equations. The Mellin transforms of those Wronskians are essentially hypergeometric OP's (see Theorems 6.2, 6.4, and 6.5). Hence, they satisfy a discrete Sturm-Liouville problem. Upon inversion of the Mellin transform this discrete problem correspond to a differential equation.
We discuss, for concreteness the case of Wronskian of Hermite wavefunctions W n, (x) = Wr(φ n (x), φ n+ (x)). Its Mellin transform is given in (6.2)-(6.3). From the difference equation of Meixner-Pollaczek polynomials the Mellin transform of the Wronskian satisfies the difference equation This implies that, for generic n and , W n, (x) satisfies a fifth-order differential equation. When when = 1, formula (6.19) simplifies as (6.18), which corresponds to the thirdorder equation of Proposition 6.7.  (6.20) with c in the fundamental strip of convergence of the Mellin transform. Note that φ n (x) is in L 2 (R + ) so that the fundamental strip always contains the line 1 2 + iR. Given that φ * n , φ * n+ and ω * n, have hypergeometric OP's factors, the above formula is a 'convolution formula' for hypergeoemetric OP's. Note that this is different from the usual (discrete) convolutions formulas of orthogonal polynomials.

Convolution of hypergeometric
When φ n (x) is a Hermite wavefunction, φ * n (s) has a Meixner-Pollaczeck polynomial factor whose parameter depends on the parity of n. Using the explicit expressions in Proposition 6.1 and Theorem 6.2 we can write the special cases of (6.20): with r = min(m, n) and = |m − n|.
In a similar way, from Proposition 6.3 and Theorem 6.4, the convolution property (6.20) gives the identity 1 2π with r = min(m, n) and = |m − n|.

Higher Order Cumulants
It is tempting to look for reciprocity formulae for cumulants of higher order (covariances, etc.). Write the moments as Q C k (m, n) = E Tr X k n , second order moments as Q C k,l (m, n) = E Tr X k n Tr X l n , and covariances as C C k,l = Q C k,l (m, n) − Q C k (m, n)Q C l (m, n). Positive and negative moments of LUE matrices satisfy recursion relations known as loop equations. The following lemma can be proved using standard methods in random matrix theory (see, e.g. [33] for similar loop equations for positive moments of the GUE).
Proof Note that Q C 0 (m, n) = n, Q C 1 (m, n) = mn and Q C −1 (m, n) = n/α. The proof of (7.3) uses the loop equations (7.1)-(7.2) together with the reciprocity law (2.3) for moments, as follows. By the loop equations, and The first gives Using Q C −1 (m, n) = n/α, the second gives as required.
There exists also a precise reflection symmetry for the covariances of one-cut βensembles at leading order in n. Suppose that the eigenvalues x 1 , . . . , x n of X n have a joint probability density proportional to 1≤ j<k≤n on the real line. Special cases of one-cut β-ensembles are the Laguerre and Jacobi ensembles defined by the weights (3.2) with α = (c − 1)n and α 1,2 = (c 1,2 − 1)n, with c, c 1 , c 2 > 0 (see discussion in Section 5). Denote the covariances by C k,l = E Tr X k n Tr X l n − E Tr X k n E Tr X l n . (7.5) The two-point connected correlator is the generating function of covariances of positive and negative moments For one-cut β-ensembles the large n limit of G 2 (z, w) exists and depends only on the edges of the cut [5,10,12,26,27], see Eq. (7.13) below. (On the other hand, for multicut ensembles the asymptotics of G 2 (z, w) is more delicate due to the presence of oscillating terms [1,11]). In the Laguerre ensemble set α = m − n = (c − 1)n, with c > 1. The edges x ± of the cut are strictly positive, see (5.4). In the Jacobi ensemble set α 1,2 = (c 1,2 − 1), with c 2 > 1, so that the cut [x − , x + ] is contained in the interval (0, 1], see (5.8).

Theorem 7.3 (Covariances of Laguerre and Jacobi ensembles at leading order in n)
. Let X n be in the Laguerre (resp. Jacobi) ensemble with c > 1 (resp. c 2 > 1). Denote the edges of the cut by 0 < x − < x + . Then, for all β > 0, More explicitly: Proof By the one-cut property, the limit is given by the explicit formula Moreover, since the cut does not contain zero, the negative covariances lim n→∞ C −k,−l exist. From (7.13) it is easy to verify the following functional equation 0 (z, w). (7.14) Using (7.7)-(7.8), the claim follows.

Orthogonal and Symplectic Ensembles
We will now discuss analogous results for the orthogonal and symplectic ensembles of random matrices, corresponding to averages over the density (3.1) with β = 1 or β = 4 respectively. Our aim will be to isolate the polynomial factors of the moments (Mellin transforms), again for the Gaussian, Laguerre and Jacobi ensembles. These ensembles are characterized by their joint eigenvalue distribution as in (3.1) with corresponding weight functions (3.2). We will consider various expectation values of power of traces with respect to (3.1) with β = 1 and 4. We will use the shorthand GSE and GOE to mean Gaussian Symplectic Ensemble, GSE (β = 4), Gaussian Orthogonal Ensemble, GOE (β = 1) and similarly for the Laguerre (LSE / LOE) and Jacobi (JSE / JOE) cases. As in the complex case, we will denote by α = m − n in the Laguerre case, and α 1 = m 1 − n and α 2 = m 2 − n in the Jacobi case, which we treat as fixed n-independent parameters. Moments of real and quaternionic Gaussian ensembles have already received some attention in the literature, however much less is known compared with the complex case β = 2. One of the first explicit formulas was derived by Goulden and Jackson [40] who were motivated by the fact that, for β = 1, moments of Gaussian matrices describe the genus expansion of non-orientable surfaces. An important development was achieved in the work of Ledoux [55] who discovered recursion relations for moments of the GOE and GSE (which can be viewed as real and quaternionic analogues of the Harer-Zagier recurrence relations). This was extended to the Laguerre ensemble in [27]. Results holding for complex moments were obtained in [64][65][66].

Recurrence relations and hypergeometric representations. We define
for all k ∈ C for which the expectations exist. Moments of the classical orthogonal ensembles satisfy recursions similar to those of the unitary ensembles. To our knowledge this was first noticed by Ledoux for the GOE [55] and extended to the LOE in [28]. The first question is whether moments of orthogonal / symplectic ensembles enjoy reflection symmetries and have orthogonal polynomial factors as in the unitary case. This is not the case as can be ascertained from the following observation. The Harer-Zagier recursion for moments Q C k (n) of the GUE is a three terms recursion in k which can be interpreted as the discrete S-L problem of some families of hypergeometric (Meixner / Meixner-Pollaczek) polynomials. Moments of the classical orthogonal ensembles satisfy recursion formulae too. For the GOE, Ledoux [55] discovered that Q R k (n) satisfy a five term recurrence relation which cannot be interpreted as a S-L equation (a second order difference equation). In fact, Ledoux also found an alternative inhomogeneous recursion formula for Q R k (n) coupled with the moments of the GUE that is somewhat more convenient for some application. An analogue of this coupled recursion was obtained later for the LOE [28]. These results suggest that suitable combinations of moments, rather than the moments themselves, have nice hypergeometric polynomial factors similar to the unitary cases.
Moments of the symplectic ensembles can be analysed similarly given the duality relations between moments of GSE, LSE and JSE of size n and the (formal) moments of the GOE, LOE and JSE of size −2n For the GOE/GSE the duality was put forward by Mulase and Waldron in terms diagrammatic expansion of Gaussian integrals [70]. See also [18,55]. This duality between orthogonal and symplectic Laguerre ensembles appeared in the paper of Hanlon, Stanley and Stembridge [46,Corollary 4.2]. The duality in the Jacobi ensembles has been observed by Forrester, Rahman and Witte [37,Eq. (4.15)]. See also [32] and [38,Appendix B].

Theorem 8.1 The combinations of GOE and GSE moments
have Meixner polynomial factors: In particular, for any integer n, S R k (n)/(2k − 1)!! and S H k (n)/(2k − 1)!! are Meixner-Pollaczek polynomials in x = −i(k + 3/2) invariant up to a change of sign under the reflection k → −3 − k, with complex zeros on the vertical line Re(k) = −3/2.
The S-L problem satisfied by the Meixner polynomials is a three term recursion formula for S R k (n) and S H k (n), These recursions, which are very similar to the Harer-Zagier formula, become five term recurrences for the moments Q R k (n) (this is Ledoux recursion [55, Theorem 2]) and Q H k (n). Corollary 8.2 Then, S R k (m, n and S H k (m, n) can be written in terms of dual Hahn polynomials (as functions of k) or Hahn polynomials (as functions of n): that can be expressed in terms of continuous dual Hahn polynomials . This proves (8.21). We use the identity (3 − 2α) k−2 Q k−2 (−2n − 2; 2, 2, −3 + 2α) = (−1) k (3+2α) k−2 Q k−2 (2n − 1; 2, 2, −3 − 2α) and find (8.21). Now, this can be written as a polynomial in k which can be cast as (8.20).
Remark As in the Gaussian case, there is a reflection formula under the transformation 2m + 1 → −2m and 2n + 1 → −2n S R (2m + 1, 2n + 1) = (−1) k S R (−2m, −2n). (8.22) The S-L problem of continuous dual Hahn polynomials corresponds to the three term recursions (in k) which are similar to the Haagerup-Thorbjørnsen formula for the moments of the LUE. Writing S R k (m, n) and S H k (m, n) in terms of the moments we obtain the following five term recursions. To our knowledge these recursion formulae are new.
Using methods similar to those in [28] it is possible to write a recursion formula for the moments of the JOE. Denote Proposition 8.5 (Recursion for moments of the JOE) Set α 1 = m 1 − n and α 2 = m 2 − n. Then the differences of adjacent moments Q R k (α 1 , α 2 , n) of the JOE satisfy the following inhomogeneous three term recursion As in the Gaussian and Laguerre cases, the above recursion formula suggests that S R k (α 1 , α 2 , n), and not the moments themselves, have a nice polynomial property. This is the content of the next theorem whose proof goes along the same lines as the Gaussian and Laguerre cases. (8.28) In particular, S R k (α 1 , α 2 , n)((α 1 + α 2 + k + 2n − 1)!/(α 2 + k)!) is a polynomial of degree 2(n − 2) in k, invariant under the reflection k → −1 − k, with zeros on the vertical line Re(k) = −1/2.

Symplectic ensembles.
The goal of this section is to establish the following polynomial property for the moments of the symplectic ensembles.
Proof We will discuss the Gaussian case in detail, as the Laguerre and Jacobi cases follow a similar pattern. By [64][65][66]Eq. (33)], we have the explicit formula where Q C k (2n) denotes the moments of the GUE (see Section 4.1) and a n = (n + 1) (n) √ π (2n)4 1−n .
Although formula (8.30) was only stated in [64][65][66] for k ∈ N, it naturally defines a meromorphic continuation to k ∈ C, as follows. As a function of k ∈ N, we have that Q C 2k (2n)/( (k + 1/2)) is a polynomial of degree 2n − 1 in k (see equation (4.4)) and hence is defined for any k ∈ C. It remains to study the second term in (8.30). Note that is a polynomial of degree n − i − j, while k i and k i+ j are polynomials of degree i and i + j respectively. Hence Q H k (n)/ (k + 1/2) is a finite sum of polynomials in k and is therefore a polynomial. To compute degrees, notice that the highest degree term in the summand of (8.30) occurs when 2i + j + n − i − j = n + i is maximal, namely when i = n −1 implying a degree of 2n −1. That the degree of the combined polynomials (i.e. unitary plus symplectic contribution) is really 2n − 2 is a consequence of the following cancellation. Setting j = 1 and i = n − 1 in the summand of (8.30) and dividing by (k + 1/2) gives the polynomial a n k n − 1 k n .
Then Stirling's formula gives the estimate a n k n − 1 Similarly, consider the complex moments which has the same leading coefficient (setting j = 2n − 1) as Hence, the terms of order k 2n−1 in (8.30) cancel, yielding a polynomial of degree 2n −2.
To compute the normalizing factor c n requires studying terms of order k 2n−2 . This is a straightforward but tedious task and we omit the details. The only contributions to the monomial k 2n−2 come from (8.31) when j = 2n − 1 and j = 2n − 2, and from the double sum in (8.30) with indices (i, j) = (n − 1, 1) and (i, j) = (n − 2, 1), (n − 2, 2). Then studying the asymptotics of these five terms as k → ∞ with Stirling's formula gives the result. For the Laguerre and Jacobi cases, this computation can be repeated with the formulae [64][65][66]Eq. (89) and Eq. (98)] which have an identical structure to (8.30) and is therefore omitted.
Below are the first few polynomials p H n (k) for the GSE, whose zeros appear to settle onto an explicit contour in the complex plane as n becomes larger (see Fig. 4). n (x) with β = 1. One can expect this case to be more complicated in general, since now (3.1) contains a non-analytic term (the absolute value of the Vandermonde determinant, which happened to be a polynomial in the cases β = 2 and β = 4). In the case of n odd we are saved by a remarkable duality principle for the Mellin transform, relating the orthogonal and symplectic ensembles. This duality involves a simple correction term which is a single hypergeometric OP.
The case of n even has a different analytic structure, evident already at n = 2. Indeed, it was known since the beginnings of random matrix theory that the parity n plays an important role for ensembles with orthogonal symmetry (see [60,Chapter 6] for example or more recently [35]), with most authors assuming n to be even for simplicity. Here it is the converse, we describe the analytic structure for n odd and give an explicit analytic continuation. First we need a proposition relating the orthogonal and symplectic ensembles. Proposition 8.8 Given the notation of Section 3, let p n (x) denote the degree n monic polynomial orthogonal with respect to the weight w 2 (x) on the interval I . Then the one-point eigenvalue density (3.3) satisfies the following duality ρ (1) 2n+1 (x) = 2ρ (4) n (x) + whereρ (4) (x) is the β = 4 eigenvalue density with respect to the modified weights Remark We give a complete proof of Propostition 8.8 in Appendix A, which is based on the skew-orthogonal polynomial formalism developed in [2]. In the specific case of the GOE, formula (8.31) was mentioned in [55]. Actually, the statement of Proposition 8.8 is implicit in Forrester's book [36, (6.120)-(6.122)]. It is worth emphasizing that this duality goes beyond the one-point function and can be formulated as a duality between the correlation kernels of n-odd orthogonal and symplectic ensembles. This suggests a possibly simpler route to studying correlation functions of n-odd orthogonal ensembles, but this lies beyond the scope of the current investigation.
We now study the consequences of the duality (8.31) for the Mellin transforms of the orthogonal ensembles. Theorem 8.9 (Duality in the n-odd orthogonal ensemble) In the three orthogonal ensembles the following identity holds for all k ∈ C and n ∈ N:
The correction ψ n (k) is a weighted Mellin transform of the corresponding orthogonal polynomial. In the JOE and LOE this takes the form while for the GOE x k is replaced with |x| 2k . The integral (8.40) can be computed explicitly by expanding p 2n (x) as a sum and integrating term by term. This expansion turns out to be a terminating hypergeometric series which can be identified as one of the hypergeometric polynomials appearing in the claimed result. In fact, for the Gaussian and Laguerre cases, precisely this calculation is carried out in a different context in [22,23], so let us just explain the Jacobi case. Then the monic polynomials p 2n (x) are proportional to the usual Jacobi polynomials which can be written down explicitly (see e.g. [78,Eq. 4.32]). Integrating term by term in (8.40) gives The sum can be matched with a hypergeometric function and we obtain ψ n (k) = η n,α 1 ,α 2 α 2 +1 2 + k 1 + α 1 +α 2 2 + k + 2n 3 F 2 k + (α 2 + 1)/2, −2n, −α 1 − 2n Finally, comparing (8.41) with the definition of the continuous Hahn polynomial in (3.14) gives the result.

Remark
In the Gaussian and Laguerre cases, the evaluation of the integral (8.40) already appeared in the literature on special functions, see the work of Bump et al. [20], Coffey et al. [22][23][24], though no connection to random matrix theory is made. These works show that the quantity (8.40) satisfies a functional equation and Riemann hypothesis with critical line Re(k) = −1/2 (Hermite polynomials) and Re(k) = 0 (Laguerre polynomials) in our notation. We believe it is new that precisely these Mellin transforms should appear in the context of random matrices. The last and most complicated case of Jacobi appears to be absent from the literature. This turns out to be a continuous Hahn polynomial p 2n (−ik; a, b, c, d) with a = c > 0 and b = d < 0. The analogous properties in this case are most easily proved by noticing that the continuous Hahn polynomial can be represented in terms of the Wilson polynomial with a negative fourth parameter. Explicitly, one has This identity demonstrates that the Mellin transforms satisfy a symmetry on the line Re(k) = 0 (the polynomials are invariant under k → −k). Furthermore, by the orthogonality property (C.2), we can deduce that the zeros all lie on the imaginary axis (this does not seem to be obvious from the Hahn polynomial representation).
We now study the orthogonal ensemble with n even. In this case the analytic structure of the Mellin transform seems to be more complicated and remains somewhat mysterious to us. For this reason we restrict ourselves to the Gaussian case, though analogous results for Laguerre and Jacobi could easily be derived. We are able to prove an analytic continuation of Q R k (2n) to an entire function of k as in the previous sections, but with a more complicated structure. We first consider the simplest case n = 2 where this structure already appears. Directly integrating |x| 2k against the density (3.1) with β = 1 and a Gaussian weight gives Clearly, the first term above has a similar structure to that already observed in the GSE and GUE. But the second term, which is a weighted Mellin transform of the error function, is different. It is analytic in the half-plane Re(2k) > −1 and standard properties of Mellin transforms show that it extends as to an analytic function in the entire complex plane except for simple poles when 2k + 1 ∈ {−2, −4, −6, . . .}. Since these simple poles are eliminated on dividing by (k + 1/2), this gives an entire function of k. In fact this analytic continuation can be given in terms of a hypergeometric function: This hypergeometric function reduces to a polynomial whenever k is a positive integer. But its analytic continuation to k ∈ C appears more complicated than in the previously considered cases. Indeed one has the asymptotics (see e.g. [77]): For say n = 4, 6, 8, . . . and so on, this structure persists and follows a similar pattern. As for the GSE, the results of [64][65][66] are again useful here, providing a general formula for the GOE moments: The first two terms in (8.42) have a simple analytic structure, similar to that found in the symplectic case. The term A R k (2n) is the generalization to larger n of the second term in (8.41).

Proposition 8.11
For any positive integer n, the ratio Q R k (2n)/ (k + 1/2) has an analytic continuation to an entire function of k.
Proof It is clear that the first two terms in (8.42) yield a polynomial in k after dividing by (k + 1/2). The third term (8.43) is the Mellin transform of the function φ n (x) = e −x 2 /2 H 2n−1 (x)erf(x/ √ 2) which is analytic in the right-half plane Re(2k) > −3. To extend to the left-half plane Re(2k) ≤ −3 it suffices to notice that φ n (x) has an asymptotic expansion near x = 0 with respect to the sequence {x 2k } k≥1 . Therefore A R 2k (2n) has an meromorphic continuation into the left-half plane except for simple poles when 2k + 1 = −2, −4, −6, . . .. Precisely these poles are eliminated after dividing through by the factor (k + 1/2).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendix A. Orthogonal and symplectic ensembles: duality
The purpose of this section is to prove Proposition 8.8 in the three classical ensembles. The proof is based on explicit results for the eigenvalue density of orthogonal and symplectic ensembles obtained by Adler et al. [2]. To begin with, we introduce the notation and relevant results obtained in [2]. There, the notation e −2V (x) is equivalent to our w 2 (x) as in (3.2). We begin by denoting by p n (x) the unique degree n monic polynomial orthogonal with respect to w 2 (x) and we set h n = I w 2 (x) p n (x) 2 dx. (A.1) An important quantity in the theory is the ratio where f and g are polynomials of minimal degree, with f ≥ 0. For the classical weights, this implies Then we define modified potentials and eigenvalue densitiesρ (1) n (x) with respect to the weight e −V 1 (x) for β = 1 and ρ (4) n (x) with respect to the weight e −2V 4 (x) for β = 4. We have that e −2V 4 (x) =w 4 (x) are precisely the modified weights (8.32). On the other hand, e −V 1 (x) = w 1 (x) and sõ ρ (1) n (x) = ρ (1) n (x). The first result we need is [2,Eq. (4.18)] which writes the density in the n-odd orthogonal ensemble as ρ (1) 2n+1 (x) =ρ Proof For the GOE case the fact thats 2n−1 = 0 follows from symmetry and the formula fors 2n is in [2,Sec. 4]. In the LOE and JOE cases these facts are less obvious, but were derived by Nagao and Forrester in [63, A.2 and A.7] based on evaluations in terms of hypergeometric functions.
Proof of Proposition 8.8 The point of the proof is that (A.5) can be simplified considerably. There are two calculations required in the proof, which in the GOE case were carried out in [2,Eqs 4.19 and 4.13] respectively. The first claim is that the following identity holds in all three cases: This can be verified by direct computation using the above explicit formulae fors 2n and γ n . Then formula (A.5) becomes ρ (1) 2n+1 (x) = ρ (2) 2n (A.17) The second claim is the following identity Differentiating both sides of (A.23) with respect to x reduces the claim to (A.24) Now using standard differential identities for the classical orthogonal polynomials and some routine calculation shows that (A.24) is a consequence of the three term recurrence relation.  W n (x 2 ; a, b, c, d) and In this paper we have considered the less conventional situation when a, b, c, 1 − d > 0. For this range of the parameters, Neretin [71, Section 3.3] found the orthogonality relation 1 2π