Moments of discrete orthogonal polynomial ensembles

We obtain factorial moment identities for the Charlier, Meixner and Krawtchouk orthogonal polynomial ensembles. Building on earlier results by Ledoux [Elect. J. Probab. 10, (2005)], we find hypergeometric representations for the factorial moments when the reference measure is Poisson (Charlier ensemble) and geometric (a particular case of the Meixner ensemble). In these cases, if the number of particles is suitably randomized, the factorial moments have a polynomial property, and satisfy three-term recurrence relations and differential equations. In particular, these moment formulae provide simple derivations of the corresponding equilibrium measures. We also briefly outline how these results can be interpreted as Cauchy-type identities for certain Schur measures.


Introduction and definitions
An orthogonal polynomial ensemble is a probability measure on N-tuples x = (x 1 , . . . , x N ) ∈ R N given by (1) dQ where µ is a probability measure on the real line having all moments, ∆(x) is the Vandermonde determinant and the normalization constant Z N is dµ(x j ).
Orthogonal polynomial ensembles arise as joint eigenvalue distributions of certain unitarily invariant models of random matrices [40]. Dyson [12] showed that the measure (1) also have a natural interpretation as a 2d Coulomb gas on the line (also known as log-gas). There are a number of models from statistical physics, probability theory and combinatorics, which are described in terms of orthogonal polynomial ensembles. The reader can consult the survey article by König [33].
The research of PC, FDC and NO'C is supported by ERC Advanced Grant 669306. The research of FDC is partially supported by Gruppo Nazionale di Fisica Matematica GNFM-INdAM.. The measure (1) can be conveniently analysed by using the so-called orthogonal polynomial method pioneered by Mehta [39] in the study of random matrices. Denote by p n (x), n ∈ N, the orthonormal polynomials with respect to µ, i.e. ∫ p m (x)p n (x)dµ(x) = δ mn . Then, a standard calculation shows that (4) dQ where K(x, y) = N−1 n=0 p n (x)p n (y). In fact, all the marginals (4) can be expressed as determinants of the correlation kernel K(x, y). In particular, for every measurable function f , where the normalized one-point function ρ N is the probability measure (6) dρ Therefore the study of the Coulomb gas Q for some reference measure µ amounts to understanding properties of the associated orthogonal polynomials p n (x).
Certainly, one the most studied quantity on orthogonal polynomial ensembles are the so-called moments. Let X = (X 1 , . . . , X N ) be distributed according to the probability measure Q. Then, the kth moment of X is the expectation of the kth power sum (7) E Q N n=1 It is well-known that moments of orthogonal polynomial ensembles admit a 1/N-expansion as N → ∞, whose coefficients satisfy a hierarchy of conditions known as loop equations initially derived in theoretical physics [2,3], and then proved using several analytical techniques [1,14,15,20,25,34,44]. (These expansions are special instances of the asymptotic series and associated topological recursions for more general probability models called β-ensembles; see [5,6,17,18,22] for the state of the art.) When the reference measure µ is the standard Gaussian, Gamma, or Beta distribution, Q is the eigenvalues distribution of the Gaussian (GUE), Laguerre (LUE) or Jacobi (JUE) unitary ensembles of random matrices, with associated Hermite, Laguerre and Jacobi polynomials [40]. These classical ensembles are exactly solvable, in the sense that the coefficients of the 1/N-expansion have an explicit description in terms of enumeration of maps (or factorisations in the symmetric group) on surfaces labelled by their genus [7,8,11,21,26,41]. The 1/N-expansion for the GUE, LUE and JUE is actually convergent; moments (8) of the classical random matrix ensembles are polynomials or rational functions in the variable N (the size of the matrix).
Moreover, the moments (8) of Gaussian, Laguerre and Jacobi unitary ensembles satisfy a three-term recurrences in k (the order of the moment). These remarkable recursions were first discovered by Harer and Zagier [27] for the GUE, and by Haagerup and Thorbjørnsen for the LUE [24] (the extension of the Haagerup-Thorbjørnsen recursion to moments of the inverse LUE was examined later in [9]). Recursions in k for moments of the JUE were obtained by Ledoux [36], and recast as three-term recurrences in [10].
Recent studies by Cunden, Mezzadri, O'Connell and Simm [10] clarified the origin of these three-term recurrences in k for moments of the classical ensembles. Those recursions originate from the second order differential equations (continuous Sturm-Liouville problems) satisfied by the Hermite, Laguerre and Jacobi polynomials. In fact, the Harer-Zagier, Haagerup-Thorbjørnsen and Ledoux recursions can be interpreted as second order difference equation in the variable k (discrete Sturm Liouville problems). It turns out that their solutions are hypergeometric orthogonal polynomials in k belonging to the Askey scheme [31]. The polynomial structure, the orthogonality relation and the hypergeometric representation of the moments as function of k explain several nontrivial symmetries and provide new results. For instance, by duality, moments of classical classical random matrices, if suitably normalized can be view as hypergeometric polynomials in the variable (N − 1) as well; therefore, they also satisfy three-term recurrences in (N − 1), whose nature is different from the general topological recursion. See [10] for the details.
This link between moments of the classical continuous orthogonal polynomial ensembles and hypergeometric orthogonal polynomials suggests that other 'solvable models' in the class of the orthogonal polynomial ensembles might exhibit unexpected representation of the moments in terms of hypergeometric series and/or special orthogonal polynomials.
In this paper we examine the classical discrete orthogonal polynomial (OP) ensembles. Namely, the Charlier, Meixner, and Krawtchouk ensembles, which correspond to the Poisson, negative binomial, and binomial reference measures µ on N, respectively. In these ensembles the reference measure µ is discrete. Hence, the measure Q can be thought of as a discrete Coulomb gas: a system of N particles on the integer lattice N, distributed according to a common distribution µ under the influence of the repelling density ∆(x) 2 . It turns out (as often the case for integer valued random variables) that rather than the moments (8), the factorial moments are more natural to look at in the discrete setting.
The outline of the paper is as follows. In Section 2 we introduce the OP ensembles considered in this paper. Then, we give the factorial moment formulae for the Charlier and Meixner ensembles calculated by Ledoux [37], and we calculate the corresponding formula for the Krawtchouk ensemble, whose statement is omitted in [37].
In Section 3, motivated by recent developments [10] on moments of the classical continuous OP ensembles, we give hypergeometric representations for the factorial moments of the Charlier ensemble, and the special case of Meixner ensemble with γ = 1. From these formulae it is clear that the kth factorial moment (suitably normalized) is a polynomial in k.
The main novel observation reported in the paper is contained in Section 4. It turns out that in the Charlier, and γ = 1 Meixner ensembles, if the number of particles N is randomized according to certain distributions, then the kth factorial moment has an even simpler hypergeometric representation and can be expressed as a Jacobi or Legendre polynomial of degree k. It follows that the randomized factorial moments satisfy threeterm recursions in k (similar to the Harer-Zagier, Haagerup-Thorbjørnsen and Ledoux recurrences in the continuous ensembles), and second order differential equations (Sturm-Liouville problems).
In the last section of this work, as a nice feature of the moment formulae, we outline a simple approach to calculate the equilibrium measures of the Charlier and γ = 1 Meixner ensembles based on comparison with the randomized factorial moments.

Moment formulae for Charlier, Meixner and Krawtchouk ensembles
We give first the precise definitions of the three classical OP ensembles considered here. These ensembles arise in connection with random partitions [28,29] and conditioned random walks [32,33,42].
(a) The Charlier ensemble corresponds to a Poisson reference measure µ = µ θ with rate θ > 0, given by1 The nth normalized Charlier polynomial c n (x; θ) can be written as [31, 9.14] The Charlier one-point function will be denoted by (11) (b) The Meixner ensemble corresponds to a negative binomial reference measure µ = µ γ q with parameters 0 < q < 1 and γ > 0, The nth normalized Meixner polynomial m n (x; γ, q) is [31, 9.10] We shall denote the Meixner one-point function by (c) The Krawtchouk ensemble corresponds to the binomial distribution µ = µ K p with K ≥ N trials and success probability 0 < p < 1, The nth normalized Krawtchouk polynomial k n (x; p, K) has the hypergeometric representation [31, 9.11] 1In this paper 0 ∈ N.
Note that, unlike in the Charlier and Meixner cases, this is a finite family of (K + 1) polynomials. The binomial distribution (15) corresponds to a negative binomial distribution with negative parameters (12) , and the Krawtchouk polynomials are related to the Meixner polynomials given by (13) in the following way: We denote the one-point function of the Krawtchouk ensemble as In [37], formulae for the factorial moments of the Charlier and Meixner ensembles are calculated. We state the results here.
is given by ). Let ρ γ N,q (x) denote the Meixner one-point function. Then, the factorial moment is given by We can also compute the factorial moments of the Krawtchouk ensemble.
is given by From the general setting of [37], we deduce the following 'integration by parts' formula In order to calculate this integral, we use the forward shift operator [31, Eq. (9.11.6)] and iterate so that Hence we observe that We now make use of the Krawtchouk generating function [31, Eq. (9.11.11)] which, by the generating function (30), is equal to since, after taking out a factor of (1 + t) k−i (1 + s) k−j , we find that the remaining discrete integral can be simplified using the binomial theorem. Thus for any integers n and m, we can calculate the t n s m coefficient by counting the terms in this product with the corresponding coefficients, to get Hence taking n = l − i, m = l − j and relabelling a = r − i, the t l−i s l−j coefficient of (31) gives Combining this with (29) we have a formula for Finally, taking the sum over l from 0 to N − 1, this is which, once we swap the order of the sums, is the desired formula.

R
1. The reader can check that the moments (25) of the Krawtchouk ensemble can be obtained formally from the Meixner case (23) by the formal substitution q → −p/(1 − p) and γ → −K:

Hypergeometric moment formulae
Propositions 1, 2 and 3 are the starting points of our analysis on the factorial moments of the classical discrete OP ensembles as functions of k. Specifically we look for properties that mirror the nice structure recently found on the classical continuous OP ensembles [10]. A 'nice property' in the continuous setting is that the moments have a hypergeometric representation. We found similar hypergeometric series for the factorial moments of the classical discrete OP ensembles when the reference measure is Poisson µ θ (x) = θ x e −θ /x!, θ > 0 (Charlier ensemble) and geometric µ 1 q (x) = q x (1 − q) (Meixner ensemble with parameter γ = 1).
In the following we will use the hypergeometric notation If one of the parameters a 1 , . . . , a p is a nonpositive integer, then the series terminates.
Recall that a series ∞ i=0 t i is hypergeometric if the ratio of consecutive terms t i+1 t i is a rational function of i, and from (34) we see that We will make use of the following identity, which can be proved by induction on N: T 1 (C ). The factorial moment of the Charlier ensemble can be written as In particular, θ −k M θ (k, N) can be extended to the entire complex plane to a polynomial in k of degree 2(N − 1). As a function of N − 1, it is a polynomial of degree k.

P
. Firstly, by (36) (with k = i) and the identity 1 If we denote we have t 0 = θ k , and for all 0 ≤ i ≤ k, So, by (35), the factorial moment M θ (k, N) = i ≥0 t i has the required hypergeometric representation.
The first few moments as functions of k are Consider now the negative binomial distribution with the special parameter γ = 1, so that the reference measure is simply a geometric distribution dµ 1 q (x) = q x (1 − q). The corresponding measure Q is the Meixner ensemble with parameter γ = 1.
The factorial moment of the Meixner ensemble with γ = 1 can be written as In particular, M 1 q (k, N) can be extended to the entire complex plane to (q/(1 − q)) k (N+1) k If we denote then we have Using (35) we conclude the proof.
The first few moments as functions of k are

Randomized factorial moments
In the classical continuous OP ensembles, the moments are themselves hypergeometric orthogonal polynomials in the variable k [10]. This is not the case in the discrete setting. However, if we suitably randomize the number of particles N, then the factorial moments simplify dramatically and can be expressed as Jacobi polynomials.
Recall that the Jacobi polynomials with shape parameters (α, β), are orthogonal with respect to the Beta distribution (x − 1) α (x + 1) β on the interval x ∈ [−1, 1]. They have the hypergeometric representation [31, 9.8] (Here we do not require them to be have squared norm 1.) When α = β = 0, this is a special case known as the Legendre polynomials, orthogonal with respect to the uniform measure on x ∈ [−1, 1], and we use the standard notation [31, 9.8.3] P n (x) = P (0,0) n (x).
For a given reference measure µ, the OP ensemble Q can be thought of as a random point configuration where the number of points (or particles) is N. We now show that the formulae of the factorial moments in the Charlier ensemble and in the Meixner ensemble with γ = 1 simplify if we randomize N.

T 3. Define the Poissonized Charlier factorial moment as
Then, M θ (k; t) can be written in terms of a Jacobi polynomial of degree k in the variable (1 + t)(1 − t) −1 , with parameters α = 1 and β = 0: In the proof we shall make use of the following generating function from the theory of hypergeometric functions.
P . We first rewrite the hypergeometric functions and exponential on the left hand side as series, i.e.
Then, we change the variables in the double sum by defining m = k + j and eliminating j, so we have which is exactly the form of the hypergeometric function required. P T 3. We use Lemma 1 with a = b = −k and c = 2, so that by definition of the Poisson measure µ t θ . Hence by the factorial moment formula of Theorem 1, we have To show the second equality (44) we use the hypergeometric representation (41) of the Jacobi polynomial.
The first few Poissonized moments are Note that (44) are polynomials in t. In particular, they are regular at t = 1, corresponding to the Poissonized Charlier ensemble with parameter θ. C 1. When we set t = 1 (i.e. Poissonize with parameter θ), we have k is the kth Catalan number.
P . The second equality is clear. To show the first equality, we observe that Then, by the Chu-Vandermonde identity, we have the binomial term as required.
The Jacobi polynomials satisfy a three-term recurrence [31,Eq. (9.8.4)] (in the degree) and a second order differential equation [31,Eq. (9.8.6)] (in the argument). In the special case (α, β) = (1, 0) they are These, combined with Theorem 3, show that the Poissonized factorial moments M θ (k, t) of the Charlier ensemble satisfy a three-term recurrences in k and a second order differential equation in t.

C
2. The Poissonized moments M θ (k, t) satisfy the three-term recurrence and the second order differential equation The recursion (49) is an immediate consequence of (47). For the proof of the differential equation (50), note that from (48) we have A straightforward manipulation concludes the proof.

Negative binomialized Meixner moments.
For the Meixner ensemble with γ = 1 it turns out that negative binomializing the number of particles provides tractable formulae. More precisely, consider the OP ensemble Q with geometric reference measure µ 1 q , 0 < q < 1. Suppose that N − 1 has 2-negative binomial distribution µ 2 t q , t > 0, i.e. (52) Then, M 1,2 q (k; t) can be written in terms of a a Legendre polynomial in the variable (1 + t)(1 − t) −1 : P . We first follow a similar method to Exton's lemma to find the hypergeometric representation in (54), and then compare with (41) to show that this is a Legendre polynomial.
By (39) and the definition of µ 2 t q , we have where in the third line we relabelled l = m − i. Simplifying the terms we get the first equality in (54). By (41), and the second equality in (54) follows.
C 3. When we set t = 1, i.e. negative binomialize with parameter q, we have P . An immediate consequence of As in the Charlier case, the negative binomialized factorial moments of the γ = 1 Meixner ensemble satisfy three term recurrences in k and differential equations in the parameter t. 4. The negative binomialized factorial moments M 1,2 q (k; t) satisfy the three term recurrences and the second order differential equation In place of the Legendre polynomials, we could also use the Legendre rational functions R k (x), which are defined by and form an orthogonal system of rational functions on the positive half-line with respect to the weight ω(x) = (1 + x) −2 . See [23].

Interpretations in terms of Schur measures.
It is well known that the Charlier and Meixner ensembles can be interpreted in terms of Schur measures on integer partitions. For a particle configuration x 1 > · · · > x N on the nonnegative integers, we can associate an integer partition λ via λ i = x i + i − N. Under this identification, the Charlier ensemble with N particles and parameter θ corresponds to the probability distribution on integer partitions given by and, assuming γ is a positive integer, the Meixner ensemble with N particles and parameters γ, q corresponds to the probability distribution For λ an integer partition with at most N parts, define Then, in the notations of Theorems 1 and 2, we have Thus the statements of Theorems 1-4 can be interpreted as Cauchy-type identities.

Equilibrium measures
Let X = (X 1 , . . . , X N ) be distributed according to an OP ensemble Q. The mean spectral measure is the probability measure ρ N = E Q 1 N N n=1 δ Xn N , and it is a simple rescaling of the one-point function It is well-known that in the limit N → ∞, under some quite general conditions, the mean spectral measure of OP ensembles converges to a limiting probability measure known as equilibrium measure (it is the unique minimizer of a certain energy functional). For the classical OP ensembles (continuous and discrete), the equilibrium measure is known explicitly in terms of elementary functions [13,[28][29][30]35,43]. In the context of the corresponding Schur measures, these equilibrium measures describe limit shapes of random partitions and have close connections to asymptotic representation theory and free probability [4].
In [36,37], Ledoux gave a precise description of the equilibrium measures of the classical ensembles as the distribution of adapted mixtures of an arcsine random variable with an independent uniform random variable.
A nice feature of Theorems 3 and 4 is that the randomized factorial moments formulae provide an easy way to compute the limit N → ∞ of the one-point function of the Charlier and γ = 1 Meixner ensembles.
In the following, ξ will be a random variable with the arcsine distribution on (−1, 1), and U a uniform random variable on (0, 1), independent of ξ,

Charlier ensemble.
The equilibrium measure can be described as follows.
T 5 (L [37]). Consider the mean spectral measure of the Charlier ensemble ρ θ N N , and suppose that θ N ∼ hN as N → ∞. Then, ρ θ N N converges to the distribution of We now explain how the Poissonized Charlier moment formulae can be used to prove Theorem 5. When we Poissonize the Charlier ensemble, we take N − 1 a Poisson random variable with rate tθ, see (42). Therefore, EN = 1 + θt and Var(N) = θt, and we see that N/θ → t in probability, as θ → ∞. In other words, θ ∼ N t , which is the assumption needed for Ledoux's theorem above, with h = 1 t . L 2. The kth moment of the random variable U + 2 √ Uhξ + h is . First we substitute sin x = s and then use a multinomial expansion to get For the integral in x we use , and the double sum simplifies to 2 F 1 −k, −k 2 ; 1 h which completes the proof.
Comparing the moments (68) of the equilibrium measure with Theorem 3 we readily get the following result.

C
5. For all θ, h > 0, the kth Poissonized factorial moments are related to the moments of (67) in the following way: After rescaling, the factorial moments converge to the moments as θ → ∞. Therefore, the moments of the Poissonized Charlier mean spectral measure converge to the moments of (U + 2 √ Uhξ + h). Clearly, the distribution (U + 2 √ Uhξ + h) has compact support, so it is determined by its the moments. Hence, from Corollary 5 we deduce an alternative proof of Theorem 5.
Here we consider the special case γ = 1, that is c = 0 (in fact for any fixed γ this is the case), when the limiting random variable takes the simpler form of As in the Charlier case, we now calculate the kth moment of this random variable.
The kth moment is the product where the first integral is (1 − q) −k (k + 1) −1 , and the second can be transformed by a substitution and binomial expansion to We can once again use (69) to see that k! (k − 2 j)!( j!) 2 (1 + q) −2j q j , and the sum is equal to q 1+q k 2 F 1 −k, −k 1 ; 1 q .
In the Meixner case, the connection between factorial moments and moments of the equilibrium measure is not as straightforward as in the Charlier ensemble. Nevertheless, there is a natural way to obtain a match by 2-negative binomialized (N − 1), weighting the factorial moments, and then taking the limit t → 1/q. We state this as a corollary of Theorem 2. (l + 2) j (l + k + 2) j (l + 1)(tq) l .
The factorial terms cancel and we get This is exactly the hypergeometric representation of Lemma 3.
We have been unable to find a hypergeometric representation for this. If it were possible, it would suggest how (if possible at all) to randomize N in a way to obtain tractable results and recover the equilibrium measure for general γ ∼ cN.
We could calculate the moments of (75) as but, as in the general Meixner case, we have not been able to properly identify a useful hypergeometric representation.