Global fluctuations for Multiple Orthogonal Polynomial Ensembles

We study the fluctuations of linear statistics with polynomial test functions for Multiple Orthogonal Polynomial Ensembles. Multiple Orthogonal Polynomial Ensembles form an important class of determinantal point processes that include random matrix models such as the GUE with external source, complex Wishart matrices, multi-matrix models and others. Our analysis is based on the recurrence matrix for the multiple orthogonal polynomials, that is constructed out of the nearest neighbor recurrences. If the coefficients for the nearest neighbor recurrences have limits, then we show that the right-limit of this recurrence matrix is a matrix that can be viewed as representation of a Toeplitz operator with respect to a non-standard basis. This will allow us to prove Central Limit Theorems for linear statistics of Multiple Orthogonal Polynomial Ensembles. A particular novelty is the use of the Baker--Campbell--Hausdorff formula to prove that the higher cumulants of the linear statistics converge to zero. We illustrate the main results by discussing Central Limit Theorems for the Gaussian Unitary Ensembles with external source, complex Wishart matrices and Markov processes related to multiple Charlier and multiple Krawtchouk polynomials.


Introduction
Let m ∈ N and a system of Borel measures {µ j } m j=1 on R be given. For each vector k = (k 1 , . . . , k m ) ∈ Z m + , the type II multiple orthogonal polynomial of multi-index k is defined as the monic polynomial p k of smallest possible degree such that p k (x)x l dµ j (x) = 0, ℓ = 0, . . . , k j − 1, j = 1, . . . , m.
We say that the multi-index k is normal if the degree of p k is | k| = k 1 + . . . + k m . A system of measures {µ j (x)} m j=1 is called a perfect system if each k ∈ Z m + is normal. Note that if m = 1 then this definition reduces to that of standard orthogonal polynomials on the real line.
Multiple orthogonal polynomials arise naturally in several different contexts. Originally they were introduced in analytic number theory (for proving irrationality and transcendence). They also play an important role in approximation theory (Hermite-Padé approximations). We refer to [37,42] for surveys on these topics. More recently, they have appeared in random matrix theory and integrable probability, as they integrate the Multiple Orthogonal Polynomial Ensembles [29]. This is an interesting class of point processes that generalize the more classical Orthogonal Polynomial Ensembles [33] and includes many interesting models such as Unitary Ensembles with external source, complex Wishart matrices, products of random matrices, the two matrix model and certain specializations of the Schur processes.

Multiple Orthogonal Polynomial Ensembles
Let us take dµ j (x) = w j (x) dµ(x), j ∈ {1, . . . , m} (1.1) for some Borel measure µ on R and non-negative weight functions w j . Then a probability measure proportional to where n = n 1 + . . . + n m and g j (x) are the functions w 1 (x), xw 1 (x), . . . , x n 1 −1 w 1 (x), · · · , w m (x), xw m (x), . . . , x nm−1 w m (x), (1.3) is called a Multiple Orthogonal Polynomial Ensemble (MOPE). It should be noted that we assume here that (1.2) is non-negative, which is an implicit assumption on the weight functions w j . In the language of [10] we see that (1.2) defines a biorthogonal ensemble, where one of the family of functions consists of polynomials. It is therefore clear that (1.2) is a determinantal point process [26]. The name, MOPE, for these processes was proposed in [29] and is explained by the fact that the correlation kernel can be expressed in terms of multiple orthogonal polynomials. The correlation kernel will however not play a role in the present paper.
A classical example of a MOPE is the Gaussian Unitary Ensemble with external source. This well-studied model was introduced in [16]. In the model, we equip the space of n × n Hermitian matrices with distribution proportional to ∼ e −n Tr(M 2 −AM ) dM, (1.4) where A is a diagonal matrix and dM is the flat Lebesgue measure on the entries. It turns out [7,8] that the eigenvalues of M form a Multiple Orthogonal Polynomial Ensemble where m is the number of distinct values along the diagonal and the n j 's are their multiplicities. The weight functions are the Gaussian w j (x) = e −n(x 2 −a j x) , and µ is the Lebesgue measure on R. The polynomials in this case are called multiple Hermite polynomials. Another example is the complex Wishart Ensemble. Given a diagonal matrix Σ = diag(σ 1 , . . . , σ n ) with σ j > 0, one studies the eigenvalues of an n × n matrix M chosen randomly from ∼ (det M ) α e −n Tr M Σ dM, on the space of positive definite hermitian matrices (for α ∈ N the eigenvalues of M are the square of the singular values of a n × (n + α) where the columns are independent and identically distributed as a complex multivariate normal distribution with covariance matrix Σ −1 ). Similarly to the Gaussian Unitary Ensemble with external source, the eigenvalues values of M form a MOPE, where m is the number of distinct values along the diagonal of Σ and the n j 's are their multiplicities. The multiple polynomials that arise here are called the multiple Laguerre polynomials of the second kind [7,29].
There are many other examples, some of which we address in Section 5. For now, we only briefly mention some important ones and refer to [29] for an extensive survey and a long list of references. In the Hermitian two matrix model, the biorthogonal polynomials can be characterized as multiple orthogonal polynomials [22,30]. MOPE's also appear in the recent activity around the singular values of products of random matrices [32]. The Borodin-Muttalib Ensembles [10,31] with an integer (or, by duality, the reciprocal of an integer) parameter θ lead to multiple orthogonal polynomials with the number of weights θ.
By taking the limit a k → 1 and using l'Hôpital's rule, we see that s λ (a) converges to the Vandermonde determinant in the variables x j = λ j − j.
If, instead, one takes the limits a k → α j such that there are precisely m different values α 1 , . . . , α m for these limits with multiplicities n 1 , . . . , n m , then the Schur polynomial turns into a determinant det g k (λ k ) as in (1.2). This limiting procedure explains why MOPE's to appear in certain specializations of the Schur processes. For example, consider the Schur measure [38] which is the probability measure on λ = (λ 1 , λ 2 , . . . , λ n ) with λ j ≥ λ k for j ≤ k, proportional to ∼ s λ (a)s λ (b).
By taking the limits b k → 1 and a k → α j for m different values α j with multiplicities n k , we see that the Schur measure is then a special case of a (discrete) MOPE of the form (1.2).

Universality of global fluctuations
The prime interest of the current paper is to investigate the global behavior of Multiple Orthogonal Polynomial Ensembles. A natural way to study a point process is via its linear statistics. Given a random configuration of points {x j } n j=1 and a smooth function f : R → R, the linear statistic X n (f ) is defined as The first natural question is whether there is a limiting distribution for the random points x j . In other words, does there exists a ρ(x) such that as n → ∞? In many determinantal point processes in random matrix theory and integrable probability this is indeed the case and identifying ρ(x) is an important problem. Naturally, ρ depends on the parameters of the model at hand (in the case of MOPE, these parameters are the weights w j , measure µ, and the indices n j ). It was proved by Hardy in [25] that for MOPE's, ρ(x) is the limiting zero distribution for the multiple orthogonal polynomials p n , and that convergence even holds almost surely (under certain conditions on the recurrence coefficients for the multiple orthogonal polynomials). The characterization of the limiting zero distribution of p n is a classical problem for multiple orthogonal polynomials. This limiting distribution can often be given in terms of the solution of a vector equilibrium problem, which is a good starting point for an asymptotic analysis of the multiple orthogonal polynomials via Riemann-Hilbert techniques in special cases. We refer to [29,39] for general discussions.
In this paper, we will address the phenomenon that the global fluctuations X n (f ) − EX n (f ) are universal. One of the remarkable facts, is that for many random processes in random matrix theory the variance of the linear statistic X n (f ) does not grow with n for functions f that are sufficiently smooth (jump discontinuities are known to contribute with terms that grow logarithmically in n). In fact, it has been proved in a variety of models that for sufficiently smooth f there is limit for some σ 2 f > 0. The models in which this limit holds have the property that the random points accumulate on a single interval. If they accumulate on multiple disjoint intervals, then the variance may have quasi-periodic behavior (see [13] for a discussion and further references). Note that if the points x j were to be independent, then the variance would grow linearly in n. That the variance does not grow is a result of the repulsion between the points that causes rigidity in the location of the points. Moreover, the limiting variance σ 2 f is universal and only depends on a few general properties of the model. Even more is true: for many models, in which the points accumulate on a single interval, it has been established that, as n → ∞, 5) in distribution. This Central Limit Theorem (CLT) together with the universality of σ f , is sometimes referred to as global universality and has been proved in a long list of important models. We will not attempt to give a full list of references, but only mention some. One of the first CLT's in random matrix theory is by Jonsson who proved a CLT for Wishart Matrices [28]. In a seminal work [27], Johansson proved (1.5) for general (continuous) β-ensembles on R with polynomial potentials and certain smoothness conditions on f . The case of discrete β-ensembles was dealt with by Borodin, Gorin and Guionnet [12]. For the special value β = 2, the continuous and discrete log-gases are special cases of Orthogonal Polynomial Ensembles and, consequently, of MOPE's. In [13] [17]. It would be interesting to see how the techniques and results of the present paper would fit into their framework.

Nearest neighbor recurrences
The main tool that we will use is the nearest neighbor recurrence relations for multiple orthogonal polynomials. Similarly to the case of standard orthogonal polynomials on the real line, for each family of multiple orthogonal polynomials (from a perfect system) there exist coefficients {a k,j } k∈N m and {b k,j } k∈N m such that where { e j } m j=1 denotes the standard canonical basis of R m . Equations (1.6) are called the nearest neighbor recurrence relations, see [42,40]. If we subtract two of these relations with ℓ = r and ℓ = s, we arrive at the equality for 1 ≤ r, s ≤ m. Note that for m = 1 these relations are obsolete, so they are a particular feature for multiple orthogonal polynomials.
In general, the map from the system of measures {µ j (x)} m j=1 to its nearest neighbor recurrence coefficients is highly non-trivial. There is an important class of measures though for which the coefficients are known to have a simple limiting behavior. We say that a perfect system {µ j (x)} m j=1 is in the multiple Nevai class if, for every j ∈ {1, . . . , m} and every ν = [0, 1] m with | ν| = 1, there exist a j ( ν) and b j ( ν) such that As was shown in [2, Thm 4.6], multiple Nevai class in particular includes systems for which supp µ j are pairwise disjoint intervals (Angelesco systems) with µ j absolutely-continuous whose a.c. densities dµ j dx are non-vanishing and analytic in a neighborhood of supp µ j . This class is conjectured to be much wider, allegedly containing Angelesco systems with the mere conditions dµ j dx > 0 a.e. on supp µ j . The computation of the constants {a j ( ν)} m j=1 and {b j ( ν)} m j=1 in (1.8) was addressed in [2] and [4]. It should be noted that the behavior (1.8) will typically not hold for measures that have unbounded support. Instead, an alternative formulation of (1.8) is expected to be true for a wide class of measures with unbounded support. Just as in the case of orthogonal polynomials (e.g., the Hermite polynomials), to get reasonable limiting expressions for the multiple orthogonal polynomials and their features, we need an appropriate rescaling. For example, for the Multiple Hermite Polynomials related to the random matrix model (1.4) we have rescaled the Gaussians to be e −n(x−a j ) 2 . This has the consequence that the eigenvalues accumulate on a union of bounded intervals. By formula (2.9) in [7], for | k| ≥ 0. We see that (1.8) does not hold directly, but we do obtain a limit if we let | k| → ∞ and n → ∞ simultaneously.
To accommodate this type of rescaling, we allow the measures to depend on a parameter n and write p if k/| k| → ν and | k|/n → 1 as | k|, n → ∞. The limits (1.8) or (1.10) are easily verified for the generalizations of the classical orthogonal polynomials. We discuss some examples in Section 5.

Main results
Our goal is to show that one can formulate a general CLT for MOPE's. In discussing the limiting behavior X n (f ) as n → ∞ we have to clarify how precisely we want to take the limit. Indeed, for each n we need to choose the indices n j . The choice we make is the following: for each n ∈ Z + we take a vector k n ∈ Z m + such that | k n | = n; k n+1 = k n + e jn for some j n ∈ {1, . . . , m}; there exists a ν ∈ [0, 1] m , | ν| = 1 with k n /n → ν as n → ∞.
One can think of { k n } ∞ n=0 as a path on the lattice Z m + , that connects (0, . . . , 0) to infinity in the direction of ν, see Fig. 1 for an example.
The following two theorems are the main results of the paper. The first one deals with the fixed choice of orthogonality measures and the second allows the measures to vary. Theorem 1.1 (Non-varying systems). Let µ have compact support and assume that the system (1.1) belongs to the multiple Nevai class. Take a path { k n } ∞ n=0 satisfying (1.11) and let a j = a j ( ν) and b j = b j ( ν) be the limits in (1.8). Then, for any polynomial f with real coefficients, we have and γ is a contour around the poles b j with counter-clockwise orientation. Note that the theorem is essentially a result on the universality of the Central Limit Theorem. This result also presents an additional motivation from integrable probability to study which systems of measures are in the multiple Nevai class. However, to have a wide scope of applications, the restriction on the compactness of the support of µ is too severe. Fortunately, it is not necessary.
Let us now allow the weights w (n) j and the measure µ (n) in (1.1) to be n-dependent (where n is the number of points in the process). Then the recurrence coefficients {a } are also n-dependent. Instead of the Nevai condition (1.8) we now only require that for a path (1.11) and for any j ∈ {1, . . . , m}, as n → ∞. Note that these right-limit-type conditions are slightly weaker than (1.10).
Theorem 1.2 (Varying systems). Assume that an n-dependent family of perfect systems satisfies (1.12) along a path with (1.11). Then, for any polynomial f with real coefficients, we have in distribution, where and γ is a contour around the poles b j with counter-clockwise orientation.
Note that Theorem 1.1 is a special case of Theorem 1.2, so it is sufficient to prove Theorem 1.2 which will be the topic of the remainder of this paper.
Our approach for proving Theorem 1.2 is inspired by [13]. In that paper, the authors considered biorthogonal ensembles where the biorthogonal functions satisfy certain recurrence relations. In case that the recurrence coefficients have limits (more precisely, when the recurrence matrix has a right limit along the diagonal that is a banded Toeplitz matrix) then Corollary 2.2 in [13] says that there is a CLT of the type (1.5). The special case m = 1 in Theorem 1.2 is (the varying version of) Theorem 2.5 in [13]. The approach of [13] can also be performed on the unit circle [21]. In further works, it has been extended to prove universality of the fluctuations on the mesoscopic scale for orthogonal polynomial ensembles on the real line [14] and, very recently, on the circle [15]. In [20] a multi-time extension was constructed to prove the global fluctuations for a class of non-colliding processes (such as stationary Dyson's Brownian motion) are governed by the Gaussian Free Field. See also [34] where the results of [13] were revisited and related to several combinatorial identities.
Although the results of [13] also hold for some special MOPE's (in fact, [13] contains examples), it fails to cover important examples such as the ones mentioned previously and which we discuss later in Section 5. The results of [13] assume that there are finite term recurrence relations for the family of polynomials {p kn } ∞ n=0 . By using the nearest neighbor recurrences (1.6) together with the (1.7) we can find the recurrence relations for the family {p kn } ∞ n=0 (cf. Proposition 2.1), but the coefficients in these recurrences depend on the nearest neighbor recurrence coefficients in a complicated way. In general, they will not have limits which is required in order for the main results of [13] to apply. Moreover, the number of terms in this recurrence is not bounded and the recurrence matrix is therefore not banded (nor essentially banded as in [14]). These are the two reasons that our setting falls outside of the scope of [13] and we have to come up with new ideas to deal with them.

Overview of the paper
In Section 2 we construct the recurrence matrix for the family {p kn } ∞ n=0 out of the nearest neighbor recurrences. Moreover, we show that if (1.8) or (1.12) hold, then one can construct a Toeplitz operator with a rational symbol out of the limiting values. The recurrence matrix for the family {p kn } ∞ n=0 and the matrix representation of the Toeplitz operator in a particularly chosen basis share the same asymptotic behavior, as shown in Theorem 2.2. It is important to note that this basis is not the canonical basis and the matrix is therefore not a Toeplitz matrix.
In Section 3 we discuss certain general determinants involving exponentials of one-sided banded matrices. We show that the asymptotic behavior (except for the leading term) of these determinants depend only on a small part of the one-sided banded matrix. We then prove that under certain conditions, the determinant converges to a Gaussian based on the Baker-Campbell-Hausdorff formula. This last part is an important novelty and replaces the use of Ehrhardt's version of the Helton-Howe-Pincus formula in [13,14,20,21]. I The latter is hard to use in our setting since our matrix representation of the Toeplitz operator is not necessarily bounded, leading to non-trivial questions on convergences. By using the Baker-Campbell-Hausdorff formula we avoid such questions and obtain a more direct proof.
In Section 4 we start by recalling that the moment-generating function for the linear statistic is a determinant of the type considered in Section 3. After verifying that the condition of the results in Section 3 are satisfied we then arrive at the proof of Theorem 1.2.
Finally, in Section 5 we present some important applications of Theorem 1.2. We show that Theorem 1.2 gives CLT's for the Gaussian Unitary Ensemble with external source, complex Wishart matrices and certain Markov processes on non-intersecting paths related to multiple Charlier and multiple Krawtchouk polynomials.

Recurrence matrices 2.1 Orthogonal polynomials
Let us remind the reader the basics of orthogonal polynomials on the real line and of Orthogonal Polynomial Ensembles. Let µ be a compactly supported Borel measure on R with infinitely many points in its support. The (monic) orthogonal polynomial p n is the unique monic polynomial of degree n such that One of the important features of orthogonal polynomials on the real line is that they satisfy a three term recurrence.
Here we have set a −1 = 0. We can organize this recurrence into a matrix form The tridiagonal infinite matrix on the right-hand side of (2.1) is called the Jacobi matrix of µ, which we denote by J . A significant part of the literature on orthogonal polynomials is devoted to studying the relationship between (properties of) µ and (properties of) J and the coefficients a n 's and b n 's. An important special class of measures are those that are in the Nevai class. This class is defined as all the measures for which there exists a and b such that a n → a, b n → b, as n → ∞.
The Nevai class of measures is rather large, but it is not easy to explicitly characterize in terms of properties of µ. A famous result of Denisov-Rakhmanov says that if a measure µ has essential support and the absolutely continuous part dµ dx > 0 almost everywhere on this interval, then µ is in the Nevai class (2.2).
We recall that for a(z) = j a j z j the Toeplitz matrix T a with symbol a is defined as Note that the Toeplitz matrix has two standard definitions in the literature, namely the one we have given or its transpose. In our definition we are motivated by (2.12) below.
Observe that if the measure µ is in the Nevai class (2.2), then the matrix J is a compact perturbation of the Toeplitz matrix T s with symbol s(z) = z + b + a/z. In particular, It is important to note that the orthogonality measures for various classical orthogonal polynomials, such as Hermite and Laguerre polynomials, are not in the Nevai class. Indeed, the measures have unbounded support and the recurrence coefficients are unbounded. This is hardly surprising, since it is known that the asymptotic behavior of the Hermite and Laguerre polynomials is best described using a rescaling (so that their zeros accumulate on compact intervals), in other words, varying orthogonality. That is, while taking the limit n → ∞ as in (2.3), we introduce the n-dependence in the measure µ (n) , the recurrence coefficients a n , and the Jacobi matrix J (n) . If done appropriately, the following generalization for (2.3) is expected to be universal: as n → ∞, for s(a) = z + b + a/z and some a > 0. Very roughly speaking, the limit (2.4) can be expected to hold if the zeros of p If they accumulate on several intervals, the limit is not true. This has been proved for varying measures of the type e −nV (x) dx where V is real analytic function with sufficient growth at infinity [18].

Recurrence matrices for multiple orthogonal polynomials
We now return to multiple orthogonal polynomials: let {µ j } m j=1 be a perfect system of measures. Then the type II multiple orthogonal polynomials satisfy the nearest neighbor recurrence relations (1.6) and the relations (1.7). Now let us choose any path with (1.11). Note that for each n ∈ Z + we have that xp kn (x) is a polynomial of degree n + 1. Therefore there exist coefficients J n,j such that These coefficients can be organized into a matrix J = (J n,j ) ∞ n,j=0 (we take J n,n+1 = 1 for all n and J n,j = 0 if j ≥ n + 2), and the above relation can be written as The matrix J has a (lower) Hessenberg form. The lower triangular part of J can be computed from the nearest neighbor recurrence (1.6) and (1.7), as we show in the next proposition.
We see that J is more complicated than the Jacobi matrix J for orthogonal polynomials. It is natural to ask whether (1.10) implies that there exists a Laurent polynomial s(z) = p j=−q s j z j such that as n → ∞, similarly to (2.4). If (2.7) holds, then Theorem 1.2 follows from the results in [13] and we are done. Although it is true for some special cases (and we refer to [13] for examples), in general it is too much to ask for. First of all, note that J depends on the specific path {k n } and, in general, J will not be a banded matrix but only one-sided banded (Hessenberg). The matrix J will happen to be banded in very special cases such as the step-line path k n = k n−1 + e n mod m which results in the so-called step-line recurrence relation.
Secondly, J depends on the nearest neighbor recurrences in a complicated fashion and there is no reason to expect that the entries of J behave in a simple way, even if (1.12) holds. In the examples in Section 5 we will see cases where one can obtain perturbations of block Toeplitz matrices, which is likely as good as it gets.
Our first main result of this paper is that (2.7) does hold if T s is replaced by the matrix representation of a Toeplitz operator in a non-standard basis, which we discuss in the next subsection.

Limiting behavior of recurrence matrices
Denote the space of rational functions on C by R and the (sub-)space of all polynomials by P. We denote the operator that projects any rational function to its polynomial part by P : R → P and denote the embedding of P into R by P * : P → R. Note that where γ is a simple counter-clockwise oriented contour that goes around all poles of g and around the pole at z, and that with this integral representation it is easily verified that P p = p for any polynomial p and that P (1/p) = 0 if the degree of p is 1 or greater. Then for any rational function r we define the multiplication operator M r : R → R by M r g = rg for all g ∈ R. The Toeplitz operator τ r : P → P, with symbol r ∈ R, is now defined as the operator τ r = P M r P * .
We will be specifically interested in the Toeplitz operator τ c with the symbol with m ∈ N, b 1 , . . . , b m , ∈ R and a 1 , . . . a m ∈ R. Thus, for any polynomial p, Note that τ c acts on such functions in a rather simple way. In fact, we have the following relation for ℓ = 1, . . . , m and k ∈ Z m + . By subtracting any pair of equations we also find π k+ e j = π k+ e ℓ for ℓ, j = 1, . . . , m. It is illuminating to compare these two identities with (1.6) and (1.7) for multiple orthogonal polynomials.
The collection of all functions π k can not be a basis since it has too many elements. We need to choose one polynomial of degree n for each n ∈ Z + , to obtain a basis. We do this using the path { k n } ∞ n=0 with (1.11). In analogy with (2.5), we define T c = (T c ) ∞ j,k=0 as follows  Theorem 2.2. Assume that an n-dependent family of perfect systems satisfies (1.12) along a path with (1.11). Let T c be as in (2.11) and J (n) be the lower Hessenberg matrix corresponding to multiplication by x in the basis {p . Then, as n → ∞, where c is (2.8) with a ℓ = a ℓ ( ν) and b ℓ = b ℓ ( ν).
Proof. Repeating the proof of Proposition 2.1 with (2.9)/(2.10) instead of (1.6)/(1.7), we arrive at τ c π kn = π k n+1 + b jn π kn + m ℓ=1 a ℓ π k n−1 (x) Using this and (2.6) (with n-dependence), for s ≥ r + 2, we get a ℓBn+s,ℓBn+s−1,ℓ · · ·B n+r+2,ℓ Now the condition (1.12) implies that the last expression goes to 0. For s = r and s = r + 1 the proof is similar, while for s < r there is nothing to prove due to the Hessenberg structure of J (n) and T c .
Note that the operator τ c is a Toeplitz operator, but T c is not a Toeplitz matrix. To get a Toeplitz matrix, we change the basis in the equality (2.11). We obtain  where the matrix T c is the Toeplitz matrix (see (2.1)) with the symbol c(z) (2.8). The matrices T c and T c are related as indicated in the following lemma.
Lemma 2.3. For any r ∈ R we have where S is the triangular matrix where γ 0 is simple counter-clockwise oriented contour that goes around the origin.
Proof. By the Cauchy residue theorem, The lemma follows from the fact that τ r is linear on the space of polynomials P.

One-sided band matrices
In this section, we recall and extend parts of [13] that are necessary for our proofs. The main results are Theorems 3.2 and 3.3. The most important differences are: (i) the analysis of [13] is restricted to band matrices, where as here we only require the matrices to be banded from one side; (ii) we will use the Baker-Campbell-Hausdorff formula to establish that the fluctuations are Gaussian in a more direct way than in [13].

Determinants and traces of one-sided banded matrices
Consider the space B of one-sided banded matrices B = (B i,j ) ∞ i,j=1 . More precisely, B ∈ B if and only if there exists b > 0 (depending on B) such that B ij = 0 for j > i + b. It is not hard to see that the space of such one-sided band matrices is closed under addition and matrix multiplication. Therefore, if B ∈ B then for any polynomial p the matrix p(B) is well-defined and belongs to B. In particular, For t ∈ R and B ∈ B, we will be interested in det (P n exp r (tB)P n + Q n ) , where Q n and P n are the complementary projection matrices: P n is the projection (x 0 , x 1 , x 2 , . . .) T → (x 0 , x 1 , . . . , x n−1 , 0, 0, . . .) T , and Q n = I − P n .
Proof. This is a standard identity that can be proved by writing log det (P n exp r (tB)P n + Q n ) = log det(I + P n (exp r (tB) − I)P n ) = Tr log(I + P n (exp r (tB) − I)P n ).
Expanding the logarithm and then the exponential near t = 0, we get the statement. See, for example, [20,Lemma 4.1] (with N = 1).
The following theorem is an extension of the results from [13] for (twosided) banded matrices.
Proof. We first recall the statement of Lemma 4.1 in [13]. In equality t = log(1 + (e t − 1)) we expand the logarithm first and then the exponential, and then observing that all coefficients of t m in the expansion vanish, we get m j=1 (−1) j+1 j for m ≥ 2. As a consequence, it follows that for any infinite matrix B = (B ij ) ∞ i,j=1 , and m ≥ 2, 4) which is Lemma 4.1 in [13].
Let us now fix m ≥ 2.
If we can show that for n sufficiently large, then, by (3.1) and (3.5), we obtain (3.2) and this finishes the proof. To prove (3.5), we use the representation (3.4), and analyze each summand on the right hand side. By definition of the trace, we have where I j,n = {(r 1 , . . . , r j−1 ) : r i > n for some i = 1, . . . , j − 1} .
Since B is one-sided banded, it follows that if r i > r i−1 + l i b, with r 0 = r j = s. Since l i ≤ m, it follows that (3.7) is 0 if r i > r i−1 + bm. Combined with the fact that s = r 0 ≤ n and that there is an i such that r i > n, it follows that the only possible contributions in (3.6) can come from the summands where |r i − n| < jmb ≤ m 2 b for all i (including r 0 = r j = s). Finally, for a chosen i = 0, . . . , j, let us write Note that each of these summands is zero unless t k ≤ t k−1 + b (we adopt t 0 = r i−1 and t l i = r i ). Combining it with |r i − n| < m 2 b, this leads to m (B) only depends on the entries B ik of B for which |i − n|, |k − n| < 2m 2 b. This proves (3.5).
The following is an important improvement over [13,21]. Then, for m ≥ 3, C (n) m (B) → 0 as n → ∞. The proof will make use of the Baker-Campbell-Hausdorff formula. We postpone the proof for a moment and proceed with an intermezzo on several well-known identities around the Baker-Campbell-Hausdorff formula and some of their consequences.

Baker-Campbell-Hausdorff formula
The Baker-Campbell-Hausdorff formula gives an expression of Z in exp Z = exp X exp Y in terms of nested commutators of X and Y .
Before we state the form that we will need, we start with a more straightforward expansion stated in the next lemma. (3.8) Remark 3.5. The series on the right-hand side (3.8) should be considered as formal series. As noted in the previous subsection, the space B forms an algebra, so each of the coefficients on the right-hand side of (3.8) is well defined. Therefore the statement is to be interpreted in the following standard way: if, on the left-hand side, the exponential is replaced by its truncated series exp r and similarly for the logarithm, then each coefficient on the left-hand side matches the corresponding coefficient on the right-hand side for sufficiently large r.
Proof. This follows by expanding first the logarithm and then the exponentials in the following sequence of equalities: The sum can be further reorganized to which gives the statement.
This immediately implies the following.
Corollary 3.6. Let A 1 and A 2 be infinite matrices with only finitely many non-zero entries. Then Proof. The trace of the left-hand side of (3.8) is linear in t for finite matrices. Therefore the trace of the coefficients of t m for m ≥ 2 on the right-hand side must vanish and this proves the statement.
The Baker-Campbell-Hausdorff formula now states the far from obvious results that each coefficient of t m on the right-hand side of (3.8) can be written as a nested commutator as stated in the next proposition. The form that we have chosen here is due to Dynkin [23].
Let X 1 , X 2 , . . . , X j be in B. We define the nested commutator by For example, We can also take "powers" within nested commutator: For example, [X 1 , X (−1) j+1 j (3.9) As before, we point out that the series on the right-hand side of (3.9) is in general not convergent and has to be interpreted as a formal series. The point here is that the coefficients on the right-hand side of (3.9) match the ones on the right-hand side of (3.8). The proof of this fact is purely algebraic, see, for example, [9,23].

Proof of Theorem 3.3
Now we return to the proof of Theorem 3.3.
We will need the following two lemmas.
Lemma 3.9. Let B ± be as in Theorem 3.3. Then, for ℓ ∈ N, where A is a matrix such that AQ s = 0 for some s ∈ N.
Proof. Choose any σ ∈ {+, −} ℓ and let ℓ + and ℓ − be the number of + and − entries in σ, respectively. By successively changing the order of B − 's and B + 's, it is not hard to see that where A σ is a linear combination of terms with each term being a product of B − 's, B + 's and [B − , B + ] in a certain order. By assumption we have that [B − , B + ] has only finitely many non-trivial columns. Clearly, this structure is preserved if we multiply this commutator by an arbitrary matrix from the left. This also holds true if we multiply the commutator from the right by a lower triangular matrix. Finally, if we multiply from the right by an upper triangular matrix that is banded, then the number of non-trivial columns increases with the bandwidth b. Still, there are only finitely many non-trivial columns. Then there exists s σ such that AQ sσ = 0. Hence the statement follows by taking s = max σ∈{+,−} ℓ s σ . for any matrix Asuch that AQ s = 0 for some s ∈ N and for n sufficiently large.
Proof. Let us start with Tr P n [B + , A]P n = 0. Then for n sufficiently large, On the other hand, Therefore, Tr P n [B + , A]P n = 0 for n sufficiently large. The proof of Tr P n [B − , A]P n = 0 is even easier. Since A = AP n for n ≥ s and P n B − = P n B − P n we have P n [B − , A]P n = P n B − P n AP n − P n AP n B − P n . Each term is a product of two finite square matrices, and for such products the trace is cyclic.
Corollary 3.11. Let B ± be as in Theorem 3.3. Let m ≥ 3, j ∈ N and set u 1 + v 1 + · · · + u j + v j = m where u i + v i ≥ 1 for i = 1, . . . , j. Then for n sufficiently large.
Proof. If v j > 1, the nested commutator is trivially zero. If v j = 1, then where r = u 1 + v 1 + · · · + u j + v j , and where each X i = B − or B + , and X r−1 = B − and X r = B + . Let Then Y r−1 = [B − , B + ], and by assumption it has a finite number of nontrivial columns. We assume that this property holds true for Y i , and prove the same for Y i−1 . By definition, Since Y i has only a finite number of non-trivial columns by assumption, and X i−1 is either lower triangular or upper triangular with finite bandwidth, it follows that both X i−1 Y i and Y i X i−1 have only a finite number of non-trivial columns. The same property therefore holds true for Y i−1 and by induction for Y 2 . By Lemma 3.10 we obtain TrP n Y 1 P n = 0, which proves the corollary for v j = 1.
If v j = 0 and u j > 1, then the nested commutator is trivially zero. If v j = 0 and u j = 1, the arguments above apply, with the right-most part of As argued in Lemma 3.2, the one-side banded structure shows that this only depends on s and r j that are within a certain distance of n. In particular, one has s > n − m 2 b and r i > n − m 2 b for all i = 1, 2, . . . , j − 1, and thus by Lemma 3.9, it follows that, for sufficiently large n, where v j = l j − u j . Note that we used that contributions from various matrices A (see lemma 3.9) disappear since s > n − m 2 b and r i > n − m 2 b.

Therefore
Tr P n B l 1 P n . . . B l j P n − P n B m P n l 1 ! · · · l j ! = Tr By the triangular structure, we have P n B − = P n B − P n and B + P n = P n B + P n , and hence we can write for n sufficiently large. Now, by Corollary 3.6 we can ignore the contribution of the first term to C for m ≥ 2. By (3.8) and the Baker-Campbell-Hausdorff theorem we can rewrite this as nested commutators. Then by applying Corollary 3.11 we see that C In this section we prove our main result, Theorem 1.2. The starting point is that the moment-generating function for the linear statistic can be written as a determinant involving the exponential of a one-sided banded matrix. We then employ the results from the previous section to prove Theorem 1.2.

Cumulants
We start by recalling the definition of the cumulants for X n (f ) = n j=1 f (x j ). For r = 1, 2, . . . , consider where M (n) j (f ) = E (X n (f )) j are the moments (we chose to work with exp r since the moment-generating function is not necessarily well-defined. An alternative is to work with formal expressions). The cumulants C m (X n (f )) are defined by for λ in a neighbourhood of 0. It is easily verified that C (n) m (X n (f )) defined in such a manner is indeed independent of r > m.
The first step is to give a formula for the cumulants C m (X n (f )) in terms of f (J). We start with the following observation. f (x j )     = det P n exp r (λf (J))P n , with J as defined in (2.5), with respect to the measures w 1 (x)dµ(x), . . . , w m (x)dµ(x), and k i = (k i,1 , . . . , k i,m ) is any path so that k n = (n 1 − 1, . . . , n m − 1).

Remark 4.2.
We will allow J, µ and w to depend on n. To avoid cumbersome notation we will supress the n-dependence and just write J, µ and w from now on.
Proof. This lemma can be found in [20,Lemma 4.1] in a more general setup. For clarity and completeness we include a proof here.
Since the system is perfect, it follows that where j i is given as in (1.11) by the condition k i+1 − k i = e j i (if this were not the case p k i+1 would not be unique, and the system would not be perfect).
Thus there is a biorthogonal system satisfying the following.
where q i,j is a polynomial of degree k i,j − 1 for j = 1, . . . , m and q i,j i is monic, satisfying for all j = i. Then by standard row and column operations, the measure (1.2) is given by for a suitable normalizing constant Z n . By Andreief's identity, it follows that . (4.2) Taking ℓ = 0, and noting that in this case the right hand side of (4.2) is the determinant of a diagonal matrix, it follows that Z n in (4.1) is given by Since the system is perfect, it follows that κ j = 0 for all j = 1, 2, . . . . Taking the expectation of (4.1) and using (4.2), it follows that .
It is easily verified that for ν = 1, 2, 3, . . . , and since exp r (f (x)) is a sum of monomials, the lemma is proven.
Together with Lemma 3.1 we thus find the following expression for the cumulants.
and P n is the projection onto the first n dimensions of ℓ 2 (N).

Comparison with the right limit
Our next step is to show that the results of Section 3 imply the following lemma. for any polynomial f , with and T f •c is given by the matrix defined in (2.11) but with τ c replaced by τ f •c .
Before we come to the proof of this lemma, we first mention an easier consequence of Lemma 3.2.
for any polynomial f .
Lemma 4.6. Let f be a polynomial. Then, for any polynomial p, is a polynomial of degree less than or equal to deg f − 1.
Proof. We recall that M c is the multiplication operator on R with multiplier and since f is a polynomial, it suffices to show that is a polynomial of degree less than or equal to j − 1, for any polynomial p. It is clear that this holds for j = 1. We assume that it is true for j − 1 (j ≥ 2), and prove it for j. By assumption where q j−2 is some polynomial of degree less than or equal to j − 2. The first term P M c M j−1 c − P * P M j−1 c p is a constant (indeed, for any rational function g we have that P M c (g − P * P g) is a constant) and P M c q j−2 has degree less than or equal to j − 1, and thus we have proven the lemma.
Proof of Lemma 4.4. By Lemma 4.6 we see that the matrices T f •c and f (T c ) only differ possibly in the first deg f − 1 columns. Hence, as n → ∞, for any fixed ℓ, m. Thus, the statement follows from Theorem 3.2 and (4.3).

Computing the cumulants for the right-limit
We now compute the limiting behavior of C for some given polynomial f . We start with the following lemma.
Lemma 4.7. Split r = r + + r − where r + is the polynomial part of r, i.e. r + = P r. Then 2. T r + is upper triangular and banded (more precisely, (T r + ) j,j+k = 0 for k < 0 or k > deg r + ), 3. T r − is strictly lower triangular, Proof. 1. This is trivial.
2. By construction of the basis it follows that r + (z)π kn can be expressed in term of π km with only n ≤ m ≤ n + deg r + . Therefore, T r + is upper triangular and banded.
3. The polynomial part of r − p for any polynomial p is a polynomial of strictly lower degree.
4. Note that where γ is a simple counter-clockwise oriented contour that goes around all poles of r − p. Now, r + (w)−r + (z) w−z is a polynomial of degree deg r + − 1 in z. After going to the basis {π kn } n we see therefore that only the first deg r + − 1 columns of [T r − , T r + ] contain non-zero entries. This proves the statement with s ≥ deg r + − 1.
By Theorem 3.3 and Lemma 4.7 we readily find the following corollary. It remains to compute the limiting behavior of C (n) 2 (T r ). For this, we will change basis and work with C  Proof. We will use the decomposition r = r + + r − with r + = P r. We wish to rely on the connection between T r and T r from Lemma 2.3. First observe that C (n) 2 (T r ) = − Tr P n T r P n T r P n + Tr P n (T 2 r )P n = Tr P n T r Q n T r P n .
Because of the triangular structure, P n T r − Q n = 0, Q n T r + P n = 0, and thus we can write C (n) 2 (T r ) = Tr P n T r + Q n T r − P n . Now, recalling Lemma 2.3, we change basis from {π kn } n to the canonical basis {z n } n . This change of basis shows that T r ± and T r ± are related by conjugation with a lower triangular matrix S, Thus, we can write It remains to argue that we can drop the S ±1 everywhere on the right-hand side. This is is in fact not obvious since P n and S ±1 do not commute. By the triangular structure we have P n S ±1 = P n S ±1 P n and S ±1 Q n = Q n S ±1 Q n .
This implies that P n SP n S −1 P n = P n and Q n SQ n S −1 Q n = Q n . (4.5) This allows us to rewrite (4.4) to C (n) 2 (T r ) = Tr P n SP n T r + Q n S −1 Q n ST r − S −1 P n . Now insert I = P n + Q n after the T r − at the right-hand side to obtain C (n) 2 (T r ) = Tr P n SP n T r + Q n S −1 Q n ST r − P n S −1 P n + Tr P n SP n T r + Q n S −1 Q n ST r − Q n S −1 P n .
Since Tr AB = Tr BA if A, B are of finite rank and by (4.5), we can simplify the first term at the right-hand side to By inserting I = P n + Q n before T r − in the first term at the right-hand side we find In the first term at the right-hand side we use (4.5). For the third term, we note that since T r − is lower triangular we have T r − Q n = Q n T r − Q n and combining this with (4.5) gives C (n) 2 (T r ) = Tr P n T r + Q n T r − P n + Tr P n T r + Q n S −1 Q n SP n T r − P n + Tr P n SP n T r + Q n T r − Q n S −1 P n .
Using the cyclicity of the trace in the second and third term we find C (n) 2 (T r ) = Tr P n T r + Q n T r − P n + Tr P n T r − P n T r + Q n S −1 Q n SP n + Tr P n T r + Q n T r − Q n S −1 P n SP n .
By the triangular structure we have P n T r − P n = P n T r − and Q n T r − Q n = T r − Q n , so we can drop the Q n and P n between the matrices T r ± and get C (n) 2 (T r ) = Tr P n T r + Q n T r − P n + Tr P n T r − T r + Q n S −1 Q n SP n + Tr P n T r + T r − Q n S −1 P n SP n .
Since Q n P n = O we have Q n S −1 Q n SP n = −Q n S −1 P n SP n and thus C (n) 2 (T r ) = Tr P n T r + Q n T r − P n + Tr P n [T r + , T r − ]Q n S −1 P n SP n .
Now, by the same argumentation as in Lemma 4.7, we have [T r + , T r − ]Q n = O for n large enough, and thus, which proves the statement.
The benefit of working with C  Proof. This follows directly by rewriting (4.6) min(j, n)r j r −j , and then taking the limit as n → ∞.

Proof of Theorem 1.2
Proof. To prove that X n (f )−EX n (f ) converge to a Gaussian in distribution, it is sufficient (1) to show that C m (X n (f )) → 0 as n → ∞ for m ≥ 3 and (2) to compute C 2 (X n (f )). Corollary 4.3 shows that these cumulants are given by C m (T f •c )) as n → ∞. Corollary 4.8 tells us that the higher cumulants indeed tend to zero as n → ∞. The second cumulant is computed using Lemmas 4.9 and 4.10.

Applications
In this section we will illustrate the main results with some applications. We will prove CLT's for the MOPEs related to the multiple Hermite, Laguerre, Charlier and Krawtchouk polynomials and also discuss them in the context of random matrix theory and integrable probability. The recurrence coefficients of the multiple Hermite polynomials of type II, denoted by p k for k = (k 1 , . . . , k n ) are then given by (1.9). We readily verify that the conditions of Theorem 1.1 are satisfied. To the best of our knowledge, this is the first CLT for this model. Although one may argue that the point process is determinantal and that there are explicit double integral formulas for the kernel [7], which means that proving a CLT using steepest descent techniques on these integrals should be possible, we emphasize that the many technical details in that approach increase with m. Our approach makes a very direct verification possible, without all the cumbersome technical details. We end this example with a comment on the fact that the results of [13] are not sufficient for proving the above CLT. As mentioned before, the main assumption for the CLT in [13] was that the right limit of the recurrence matrix J in (2.5) is constant along the diagonals. We claim that this is only possible for the multiple Hermite polynomials in case m = 1. This is most easily understood by an example: Consider the external source model with m = 2, the values a 1 = a 2 on the diagonal and n = (n/2, n/2) with n even. Then take the family of polynomials {p k j } n j=0 where k j = (⌊(j +1)/2⌋, ⌊j/2⌋) and ⌊x⌋ stands for the largest integer less than x. From the recurrences (1.9) and (1.7), we then see that Therefore there is a 2-periodic structure along the diagonals in the recurrence matrix J, and every right limit will have this 2-periodicity. In fact, the right limit is a particular block Toeplitz matrix. The proof of the CLT in [13] only apply to cases where the right limit is a scalar Toeplitz matrix and the extension to block Toeplitz is in general false. The point is that the block Toeplitz matrix in this example is of a very special type.

Wishart ensembles and Multiple Laguerre polynomials
In the next example, we consider the measure (det M ) α e −n Tr M Σ dM, on the space of positive definite matrices, where α > 0 and Σ is a diagonal matrix with strictly positive entries. We will study the case where Σ has precisely m different values σ 1 , . . . , σ m on the diagonal with multiplicities n 1 , . . . , n m . It is known that the eigenvalues of M form a multiple orthogonal polynomial ensemble, with weights now given by The multiple orthogonal polynomials are called multiple Laguerre polynomials of the second kind [41] and are denoted by L α, σ k where σ = (σ 1 , . . . , σ m ). The nearest neighbor recurrences now are given by and [41, 3.6 We see directly that also here the conditions of Theorem 1.2 are satisfied. Note that we can even let α depend linearly on n which will change the parameters in the CLT in an obvious way. To the best of our knowledge, the CLT for Wishart ensembles has not appeared in the literature before.

Discrete multiple orthogonal polynomials
The examples above are well-known in the literature. It is lesser known that discrete multiple orthogonal polynomials also appear in integrable probability. We will discuss two families of examples related to multiple Charlier and multiple Krawtchouk polynomials, based on Markov processes for noncolliding particles.

Markov processes on Weyl chambers
Let P t (x − ξ) be the transition kernel (for the probability to jump from ξ to x after time t) for a single particle Markov process (time may be discrete or continuous) on Z. We will mainly be interested in the cases of the Poisson process P t (x) = e −tλ (tλ) x x! , x = 0, 1, 2, 3 . . . , 0, x < 0, t > 0. (5.1) and the simple random walk 0, There is a standard construction for defining non-colliding processes with n particles, starting from these Markov processes, using Doob's h-transform [33]. The key to this construction is the fact that k → α k is a harmonic function for the transiton kernel P t (x − ξ). That is, for some constant c t,α . On the Weyl Chamber we can then, using Andreief's identity, define a Markov process by taking the transition function For now the coefficients α j are arbitrary, but all with distinct values, in (0, 1]. The process depends on the values of α j . One can construct certain multiple orthogonal polynomial ensembles in the following way.
In case of (5.1) and (5.2) this will allow us to simplify the factor det P t (x i − ξ j ) to a product of Vandermonde determinant and certain weight functions. The arguments of this part depend however on the precise form of P t and we leave that to the concrete examples.
• Fix m and for ℓ = 1, . . . , m take (distinct) values γ ℓ and integers n ℓ such that n 1 + . . . + n m = n. Then we take the limit in which n ℓ of the α j 's converge to γ ℓ for ℓ = 1, . . . , m. By taking this limit, using L'Hôpital's rule, we find where the g j 's are given by (1.3) with w j (x) = γ x j . Then, as we will see, in case of (5.1) and (5.2) the process on the Weyl Chamber at time t forms a multiple orthogonal polynomial ensemble.
Before we come to the more explicit constructions in case of (5.1) and (5.2), we first mention an interpretation of the parameter γ j . This is best understood using the construction in [11]. In that paper, the authors construct a more general evolution on a system of interlacing particles on a two dimensional lattice Z × N. On the first level (or horizontal section) (Z, 1) there is a single particle moving along Z, on the second level (Z, 2) there are two particles moving along Z, and on the n'th level (Z, n) there are n particles moving along Z. Within each level the particles interact by being conditioned not to take the same coordinate (non-intersecting), and between levels the particles interact by the condition that the particles on the n-th and n + 1-th level are interlacing.. When time evolves, particles randomly hop in horizontal direction, according to certain rules that ensure that the interlacing condition is satisfied at all times. In the construction, particles on the same level have the same jump rate or speed. If one considers the marginal distribution of the particles on the n-th level, then the evolution of the particles is exactly as described above in (5.3). One may now think of α j as the speed of the particles on the j-th horizontal section. Multiple orthogonal polynomials thus arise in the construction of [11] when there are only m different jump rates to be distributed among these sections and n ℓ is the number of sections (below or including the n-th section) that have speed γ ℓ .

Multiple Charlier Ensemble
Consider P t (u) as in (5.1) and let us start the process at consecutive integers ξ i = i − 1, for i = 1, . . . , n. The heart of the matter is the fact that where (a) j stands for the Pochhamer symbol (a) j = a(a + 1)(a + 2) · · · (a + j − 1). Thus, up to the factor (tλ) x /x!, we see that P t (i − 1 − x) is a polynomial of degree i − 1 in x. Therefore (after using some standard rules for determinants), we see that det (P t (x j − i + 1)) n i,j=1 = c n (λ, t) , for some c n (λ, t) independent of the x j 's. Concluding, we find that the process (5.3) at time t is indeed a MOPE on the non-negative integers, with n = (n 1 , . . . , n m ), and w j (x) = γ x j , j = 1, . . . , m, dµ(x) = (λt) x x! .
The multiple orthogonal polynomials are called the Multiple Charlier polynomials, denoted by C k (x). Observe that in the case m = 1, these are indeed the classical Charlier polynomials. The nearest neighbor recurrence coefficients for C k have been computed explicitly in [40]. By taking the parameter a j in [40] to be a j = λtγ j we find xC k (x) = C k+ e ℓ (x) + b k,ℓ C k+ e ℓ (x) + m j=1 a k,j C k− e j (x).
To get a reasonable limit theorem, we need to rescale the time parameter and space parameters. The explanation of this is that after long time the particles will spread out over a large interval. It turns out the correct rescaling is t = nτ and x = nξ. By settingC k (ξ) = n −| k| C k (nξ) we obtain ξC k (ξ) =C k+ e ℓ (ξ) +b k,ℓC k+ e ℓ (ξ) + m j=1ã k,jC k− e j (ξ).
Now one readily verifies that the conditions of Theorem 1.2 are satisfied. The multiple Charlier ensemble with m = 2 has appeared before (although not explicitly mentioned). Indeed, in [19], the dynamics on interlacing particles systems of [11] with two speeds was studied. By taking the marginal density at the horizontal level with n particles, we obtain exactly the process described here. The different coefficients γ 1 and γ 2 represent two different speeds in the evolution of the particle system. The main result of [19] was the computation of the long time behavior of global fluctuation for the two dimensional system, which turns out to be described by the Gaussian Free Field. Observe that the proof of the fluctuations of the present paper avoids many of the technical details that one needs to overcome in the approach of [19]. In fact, in an effort to keep the number of technical details to a minimum, [19] deals with only two possible speeds, whereas here we can allow m different values without much extra effort.

Multiple Krawtchouk Ensemble
Let us now consider the case (5.2). Now note that for x = 0, . . . , t + n − 1. This implies that for x j = 0, . . . , t + n − 1, where c m,i,t does not depend on x j and q i is a polynomial of degree n − 1. In fact, all q 1 , . . . , q n can be shown to be linearly independent and thus, after some standard rules for determinants, we find det (P t (x j + 1 − i)) n i,j=1 = c n (p, t) n j=1 (p/(1 − p)) x j x j !(t + n − 1 − x j )! det x i−1 j n i,j=1 , where c n (p, t) is a constant independent of the x j 's. Concluding, we find that the process (5.3) at time t is indeed a MOPE on {0, 1, . . . , t + n − 1}, with n = (n 1 , . . . , n m ), and w j (x) = γ x j , j = 1, . . . , m, dµ(x) = (p/(1 − p)) x x!(t + n − 1 − x)! .
These are the weights for the Multiple Krawtchouk Polynomials K p,t+n−1 k that were studied in [5,24,42]. Naturally, when m = 1, these reduce to the standard Krawtchouk polynomials. The nearest neighbor recurrence relation reads (see [