A unified fluctuation formula for one-cut β-ensembles of random matrices

Many textbooks on random matrices [4, 5, 45] begin with the derivation of Wigner's semicircle law for the Gaussian unitary ensemble (GUE)⁵ , which is arguably the best known result in random matrix theory (RMT). One of the most common proofs consists in introducing the moments ${G}_{N,\kappa }={N}^{-1}\mathrm{Tr}{{\mathcal{G}}}_{N}^{\kappa }$, where ${{\mathcal{G}}}_{N}$ denotes an $N\times N$ GUE matrix and $\kappa$ is an integer, and in applying Wick's formula to evaluate the averages $\langle {G}_{N,\kappa }\rangle$. In order to compute these expectation values, certain diagrams are associated to various contributions of Wick's expansion; then, in the large N limit, only certain diagrams (the planar ones) survive, which can be counted explicitly. Eventually, one discovers that the limits ${\mathrm{lim}}_{N\to \infty }\langle {G}_{N,\kappa }\rangle$ are non-zero only for $\kappa$ even and given by the Catalan numbers. This computation implies that the empirical density ${\rho }_{N}(x)={N}^{-1}{\displaystyle \sum }_{i}\delta (x-{\lambda }_{i})$ of the eigenvalues $\{{\lambda }_{i}\}$ converges in expectation to Wigner's law $\rho (x)={(2\pi )}^{-1}\sqrt{4-{x}^{2}}$. The simplicity of this result might lead to believe that the very same approach (known as the method of moments) should be easily extended to higher cumulants of ${G}_{N,\kappa }$. However, it turns out that this problem is highly nontrivial, with very few results available in the literature [11, 15, 16, 25, 26, 38, 39, 42]. To date, almost 60 years after Wigner's original paper [54], a self-contained explicit formula for the higher cumulants ${C}_{v}({G}_{N,{\kappa }_{1}},...,{G}_{N,{\kappa }_{v}})$ for generic v and $\{{\kappa }_{i}\}$ is still lacking⁶ .The situation is even worse for general invariant ensembles, where the method of moments is not applicable.

In this work we consider random self-adjoint matrices ${{\mathcal{X}}}_{N}$ whose eigenvalues $\{{\lambda }_{k}\}\in {\mathbb{R}}$ have the canonical distribution of a 2D Coulomb gas at inverse temperature $\beta \gt 0$, namely

$$\begin{eqnarray}&&{\rm{d}}{{\mathbb{P}}}_{N,\beta }(\{{\lambda }_{k}\})=\displaystyle \frac{1}{{{\mathcal{Z}}}_{N,\beta }}{{\rm{e}}}^{-\beta (-\displaystyle \frac{1}{2}\displaystyle \sum _{i\ne j}\mathrm{ln}| {\lambda }_{i}-{\lambda }_{j}| +N\displaystyle \sum _{i}V({\lambda }_{i}))}\displaystyle \prod _{i=1}^{N}{\rm{d}}{\lambda }_{i}.\end{eqnarray} \tag{ 1 }$$

Here the single particle potential $V(x)$ is assumed to be such that in the thermodynamic (large N) limit the Coulomb gas reaches a stable equilibrium distribution supported on a single bounded interval $[a,b]$ of the real line. This class of one-cut ensembles⁷ includes the classical Gaussian, Wishart and Jacobi ensembles of RMT. For $\beta =1,2$ and 4, (1) is the joint law of the eigenvalues of real symmetric, complex Hermitian or quaternion self-dual invariant matrices, respectively. General ensembles parametrized by non-integer values of $\beta \gt 0$ can be realized from sparse real random matrices, as shown by Dumitriu and Edelman [24] and Killip and Nenciu [36]. Moreover, in view of certain applications to physics (see, e.g., [28]), it is worth studying the general joint probability density function (1) independently of a concrete matrix model realization.

In this paper we consider the moments ${X}_{N,\kappa }={N}^{-1}\mathrm{Tr}{{\mathcal{X}}}_{N}^{\kappa }$ for matrices ${{\mathcal{X}}}_{N}$ with joint probability law of the eigenvalues given by (1), and we derive self-contained formulae for their covariances in the limit $N\to \infty$. More precisely we study the quantities

$$\begin{eqnarray}&&\underset{N\to \infty }{\mathrm{lim}}{N}^{2}\mathrm{Cov}({X}_{N,\kappa },{X}_{N,{\ell }})=\displaystyle \frac{1}{\beta }{\alpha }_{\kappa ,{\ell }}.\end{eqnarray} \tag{ 2 }$$

The $1/\beta$-dependence and the ${N}^{-2}$ scaling of such covariances are customary in the fluctuations for 2D Coulomb gases [6]. The main result of the paper is a universal formula for the generating function of (2); it is universal in the sense that it depends only on the support of the equilibrium density, but not on the potential $V(x)$. This is a direct consequence of the macroscopic universality of the smoothed two-point function in RMT. Therefore, we can always write the covariances (2) in terms of a reference one-cut ensemble—the natural choice is the GUE. This allows us to prove a self-contained and manageable formula for the limiting covariances of an arbitrary one-cut ensemble (28). As particular examples, we derive known as well as new results for the Gaussian (Hermite), Wishart (Laguerre) and Jacobi β-ensembles within a unified simple approach, unveiling the underlying structure.

We will derive our results in the framework of the Coulomb gas analogy. The advantage is twofold. Firstly, the method is insensitive to the specific value of $\beta \gt 0$, so there is no need to produce separate proofs for each symmetry class. Secondly, it easily applies to matrix models whose underlying diagrammatic is cumbersome or unclear (such as the Jacobi ensemble). Moreover, the method is a good candidate for generalizations to higher order cumulants.

Various methods have been employed to compute cumulants of moments of β-ensembles. All of them are combinatorial in nature. There are several works on the Gaussian ensembles scattered throughout the literature of quantum field theory in the limit of large internal symmetry group [8, 12, 33, 34, 41]. Some of these results reproduce earlier findings by Tutte [49, 50] for planar enumeration of maps. Johansson [35] studied the global fluctuations of Hermitian random matrices and was the first one to prove a central limit theorem for β-ensemble with polynomial potentials. More recently, in a series of paper devoted to the second-order freeness [16, 38, 39, 42], the relation between moments of some random matrix ensembles and non-crossing partitions has been employed to compute covariances to leading order in N on complex and real Gaussian and Wishart matrices. Dumitriu and Edelman [25] and Dumitriu and Paquette [26] studied the fluctuations around the equilibrium density in Gaussian, Wishart and Jacobi ensembles at generic Dyson index β. Their approach was based on the realization of these ensembles in terms of sparse matrices. The computation of cumulants of moments for some β-ensembles is possible in principle by solving a recurrence relation associated to the hierarchy of loop equations [10, 11, 15].

The plan of this paper is as follows. We start with some preliminaries and a summary of the main ideas of the Coulomb gas approach (section 2). In section 3 we present our main result. In section 4 we analyse the classical β-ensembles: the Gaussian, Laguerre and Jacobi β-ensembles. Finally, in section 5 we discuss some open problems. The paper is complemented by two appendices. In appendix A we review the connections between the problems of moments fluctuations, enumeration of planar maps and non-crossing pair partitions. Appendix B contains the derivation of a technical result.

Here we briefly summarize (without derivations) some aspects of the 2D Coulomb gas approach, which are more directly relevant to our work. The starting point of the method is the analogy with a gas of Coulomb charges pointed out by Wigner and Dyson [27, 54]. The joint law ${{\mathbb{P}}}_{N,\beta }(\{{\lambda }_{k}\})$ in (1) is the Gibbs–Boltzmann measure of a 2D Coulomb gas confined in a interval $\Lambda \subset {\mathbb{R}}$ at inverse temperature $\beta \gt 0$ (the Dyson index) with a single-particle potential $V(\lambda )$ bounded from below and finite for $\lambda \in \Lambda$ (we adopt the usual physical convention that probabilities vanish in regions of infinite potential). For simplicity we assume $V(x)$ convex with superlogarithmic growth at infinity. It is known [35] that under these assumptions on $V(x)$ the Coulomb gas is supported in the thermodynamic limit on a bounded interval $[a,b]$. These assumptions on $V(x)$ are purely technical, and in fact the correct condition should be the stability of the one-cut configuration under small analytic perturbations of the external potential (this condition goes under the name of 'off-criticality regime').

In the thermodynamic (large N) limit it is convenient to describe the positions $({\lambda }_{1},...,{\lambda }_{N})$ of the gas particles (the eigenvalues of the β-ensemble) in terms of the empirical density of the Coulomb gas

$$\begin{eqnarray}&&{\rho }_{N}(x)=\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}\delta (x-{\lambda }_{i}).\end{eqnarray} \tag{ 3 }$$

The partition function of the gas may be written in terms of the empirical density in the following way. First, notice that for large N the energy of the Coulomb gas admits an integral representation

$$\begin{eqnarray}&&-\displaystyle \frac{1}{2}\displaystyle \sum _{i\ne j}\mathrm{ln}| {\lambda }_{i}-{\lambda }_{j}| +N\displaystyle \sum _{i}V({\lambda }_{i})={N}^{2}{\mathcal{E}}[{\rho }_{N}]+o({N}^{2})\end{eqnarray} \tag{ 4 }$$

in terms of the mean-field energy functional ${\mathcal{E}}$ defined as

$$\begin{eqnarray}&&{\mathcal{E}}[\sigma ]=-\displaystyle \frac{1}{2}{\iint }_{x\ne y}\;\;{\rm{d}}\sigma (x){\rm{d}}\sigma (y)\;\mathrm{ln}| x-y| +\int {\rm{d}}\sigma (x)\;V(x).\end{eqnarray} \tag{ 5 }$$

Using this energy functional we have

$$\begin{eqnarray}&&{{\mathcal{Z}}}_{N,\beta }=\int {\mathcal{D}}[{\rho }_{N}]\;\int \displaystyle \prod _{i=1}^{N}{\rm{d}}{\lambda }_{i}\;{{\rm{e}}}^{-\beta {N}^{2}{\mathcal{E}}[{\rho }_{N}]+o({N}^{2})},\end{eqnarray} \tag{ 6 }$$

where the functional integration $\int {\mathcal{D}}[{\rho }_{N}]$ is restricted to the probability measures compatible with the Coulomb gas configuration $\{{\lambda }_{i}\}$. A spectral linear statistics (or linear statistics for short) is defined by

$$\begin{eqnarray}&&{A}_{N}={N}^{-1}{\displaystyle \sum }_{i=1}^{N}f({\lambda }_{i})=\int {\rm{d}}x\;{\rho }_{N}(x)f(x),\end{eqnarray} \tag{ 7 }$$

where $f(x)$ is a suitably chosen test function. Sometimes (7) is referred to as sum function on the particles positions $\{{\lambda }_{i}\}$. The average of A_N is given by

$$\begin{eqnarray}&&\langle {A}_{N}\rangle =\int {\rm{d}}x\;\langle {\rho }_{N}(x)\rangle f(x).\end{eqnarray} \tag{ 8 }$$

In a similar way one computes the variance of A_N, or in general the covariance between two linear statistics ${A}_{N}^{(1)}={N}^{-1}{\displaystyle \sum }_{i=1}^{N}{f}_{1}({\lambda }_{i})$ and ${A}_{N}^{(2)}={N}^{-1}{\displaystyle \sum }_{i=1}^{N}{f}_{2}({\lambda }_{i})$ as

$$\begin{eqnarray}&&\mathrm{Cov}({A}_{N}^{(1)},{A}_{N}^{(2)})=\int {\rm{d}}x{\rm{d}}y\;\mathrm{Cov}({\rho }_{N}(x){\rho }_{N}(y)){f}_{1}(x){f}_{2}(y).\end{eqnarray} \tag{ 9 }$$

The same reasoning can be extended to higher order cumulants. In general, the mixed cumulant ${C}_{v}({A}_{N}^{(1)},{A}_{N}^{(2)},...,{A}_{N}^{(v)})$ of v linear statistics (not necessarily distinct) can be recovered by integrating ${A}_{N}^{(j)}={N}^{-1}{\displaystyle \sum }_{i=1}^{N}{f}_{j}({\lambda }_{i})$ against the cumulant of the gas density at v points ${C}_{v}({\rho }_{N}({x}_{1}),{\rho }_{N}({x}_{2}),...,{\rho }_{N}({x}_{v}))$, also known as v-point connected correlation function:

$$\begin{eqnarray}&&{C}_{v}({A}_{N}^{(1)},{A}_{N}^{(2)},...,{A}_{N}^{(v)})=\int (\displaystyle \prod _{j=1}^{v}{f}_{j}({x}_{j}){\rm{d}}{x}_{j}){C}_{v}({\rho }_{N}({x}_{1}),{\rho }_{N}({x}_{2}),...,{\rho }_{N}({x}_{v})).\end{eqnarray} \tag{ 10 }$$

For v = 1 and v = 2 this corresponds to (8) and (9), respectively. In the asymptotic regime as $N\to \infty$, the connected correlation functions of the gas are generated by the free energy $-\mathrm{log}{{\mathcal{Z}}}_{N,\beta }$ by taking a suitable number of functional derivatives with respect to the external potential. More precisely, for every positive integer v:

$$\begin{eqnarray}&&{C}_{v}({\rho }_{N}({x}_{1}),{\rho }_{N}({x}_{2}),...,{\rho }_{N}({x}_{v}))={(-\displaystyle \frac{1}{\beta {N}^{2}})}^{v}[\displaystyle \prod _{i=1}^{v}\displaystyle \frac{\delta }{\delta V({x}_{i})}]\mathrm{ln}{{\mathcal{Z}}}_{N,\beta }[1+o(1)].\end{eqnarray} \tag{ 11 }$$

For instance, the one-point and two-point connected correlation functions in the limit as $N\to \infty$ are given by

$$\begin{eqnarray}&&\langle {\rho }_{N}(x)\rangle =-\displaystyle \frac{1}{\beta {N}^{2}}\displaystyle \frac{\delta \mathrm{ln}{{\mathcal{Z}}}_{N,\beta }}{\delta V(x)}[1+o(1)]\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\mathrm{Cov}({\rho }_{N}(x){\rho }_{N}(y))={(\displaystyle \frac{1}{\beta {N}^{2}})}^{2}\displaystyle \frac{{\delta }^{2}\mathrm{ln}{{\mathcal{Z}}}_{N,\beta }}{\delta V(x)\delta V(y)}[1+o(1)].\end{eqnarray} \tag{ 13 }$$

Under the assumption that the limit ${\mathrm{lim}}_{N\to \infty }\langle {\rho }_{N}(x)\rangle =\rho (x)$ exists, from (11) we see that

$$\begin{eqnarray}&&\underset{N\to \infty }{\mathrm{lim}}{N}^{2(v-1)}{C}_{v}({\rho }_{N}({x}_{1}),{\rho }_{N}({x}_{2}),...,{\rho }_{N}({x}_{v}))={(-\displaystyle \frac{1}{\beta })}^{v-1}[\displaystyle \prod _{i=1}^{v-1}\displaystyle \frac{\delta }{\delta V({x}_{i})}]\underset{N\to \infty }{\mathrm{lim}}\langle {\rho }_{N}({x}_{v})\rangle \end{eqnarray} \tag{ 14 }$$

exists too⁸ . Hence, in principle a complete answer to the general problem of computing spectral linear statistics in the large N limit could be achieved if all functional derivatives with respect to the external confining potential were known. At present the complete solution to this problem is out of reach. Fortunately, in some circumstances we have a partial but explicit answer. In the following we focus to first and second order effects:

$$\begin{eqnarray}&&\rho (x)=\underset{N\to \infty }{\mathrm{lim}}\langle {\rho }_{N}(x)\rangle \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\mathcal{K}}(x,y)=-\underset{N\to \infty }{\mathrm{lim}}{N}^{2}\mathrm{Cov}({\rho }_{N}(x){\rho }_{N}(y))=\displaystyle \frac{1}{\beta }\displaystyle \frac{\delta \rho (x)}{\delta V(y)}.\end{eqnarray} \tag{ 16 }$$

The quantity $\rho (x)$ is the density of states, or equilibrium density, while ${\mathcal{K}}(x,y)$ is known as smoothed two-point connected correlation kernel (see [7, 13] for details). Therefore, the limiting mean and covariances of smooth linear statistics (the lowest order cumulants) are

$$\begin{eqnarray}&&\underset{N\to \infty }{\mathrm{lim}}\langle {A}_{N}\rangle =\int {\rm{d}}x\;\rho (x)f(x)\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\underset{N\to \infty }{\mathrm{lim}}{N}^{2}\mathrm{Cov}({A}_{N}^{(1)},{A}_{N}^{(2)})=-\int {\rm{d}}x{\rm{d}}y\;{\mathcal{K}}(x,y){f}_{1}(x){f}_{2}(x).\end{eqnarray} \tag{ 18 }$$

From now on, we will restrict ourselves to one-cut ensembles, i.e. we will always assume that the potential $V(x)$ is convex and with superlogarithmic growth at infinity. Under these hypotheses, the density of states $\rho (x)$ (15) exists, is absolutely continuous and supported on a single bounded interval $\mathrm{supp}\;\rho =[a,b]$, ($a\lt b$), and is the unique probability measure that minimizes the mean-field energy functional ${\mathcal{E}}$ in (5). We notice that β does not appear in the minimization problem, and hence the density of states (and in particular its support) does not depend on β. Remarkably, under the condition that the support is a single interval, $\rho (x)$ and ${\mathcal{K}}(x,y)$ have the closed expressions

$$\begin{eqnarray}&&\rho (x)=\displaystyle \frac{1}{\pi }\displaystyle \frac{{{\bf{1}}}_{a\lt x\lt b}}{\sqrt{(x-a)(b-x)}}[1+\displaystyle \frac{1}{\pi }\;P{\displaystyle \int }_{a}^{b}{\rm{d}}y\displaystyle \frac{V\prime (y)\sqrt{(y-a)(b-y)}}{y-x}],\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{\mathcal{K}}(x,y)=\displaystyle \frac{1}{\beta {\pi }^{2}}\displaystyle \frac{{{\bf{1}}}_{a\lt x,y\lt b}}{\sqrt{(y-a)(b-y)}}\displaystyle \frac{{\rm{d}}\;}{{\rm{d}}x}\displaystyle \frac{\sqrt{(x-a)(b-x)}}{x-y},\end{eqnarray} \tag{ 20 }$$

where P denotes Cauchy's principal value.

The first result (19) is usually referred to as Tricomi's formula [47]. The beautiful identity (20) is a consequence of Tricomi's formula and (16), and has been independently discovered by several authors [3, 7, 13, 14, 29]. It is worth emphasizing that while the one-point density (19) depends explicitly on the potential $V(x)$, the smoothed two-point connected correlation kernel (20) enjoys a very surprising feature: it depends only on the edge points $a,b$ of the spectral density. The dependence of ${\mathcal{K}}(x,y)$ on the details of the model only through the edges of the spectral density is often referred to as to long-range (or macroscopic) universality. Conversely, the microscopic universality emerges on a scale ${\mathcal{O}}(1/N)$ of the interior points of the limiting support of $\rho (x)$ and is governed by the celebrated sine-kernel [13]. The universality of ${\mathcal{K}}(x,y)$ is a key observation that eventually leads to our main result, presented in next section.

Let ${X}_{N,\kappa }={N}^{-1}\mathrm{Tr}{{\mathcal{X}}}_{N}^{\kappa }$ denote the moments of a generic one-cut β-ensemble of random matrices. For any fixed positive integers $\kappa$ and ${\ell }$, the limits

$$\begin{eqnarray}&&{\alpha }_{\kappa }=\underset{N\to \infty }{\mathrm{lim}}\langle {X}_{N,\kappa }\rangle \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{\beta }{\alpha }_{\kappa ,{\ell }}=\underset{N\to \infty }{\mathrm{lim}}{N}^{2}\mathrm{Cov}({X}_{N,\kappa },{X}_{N,{\ell }})\end{eqnarray} \tag{ 22 }$$

exist and are given by

$$\begin{eqnarray}&&{\alpha }_{\kappa }={\displaystyle \int }_{a}^{b}{\rm{d}}x\;\rho (x){x}^{\kappa }\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{\alpha }_{\kappa ,{\ell }}=\displaystyle \frac{1}{{\pi }^{2}}\;P{\iint }_{a}^{b}{\rm{d}}x{\rm{d}}y\;\sqrt{\displaystyle \frac{(x-a)(b-x)}{(y-a)(b-y)}}\displaystyle \frac{\kappa }{y-x}{y}^{{\ell }}{x}^{\kappa -1},\end{eqnarray} \tag{ 24 }$$

respectively. The identities (23)–(24) are a specialization of (17)–(20). While the limiting average ${\alpha }_{\kappa }$ depends explicitly on the probability measure that defines the ensemble through $V(x)$, and hence $\rho (x)$ (given in (19)), the covariance $(1/\beta ){\alpha }_{\kappa ,{\ell }}$ does not. Indeed, the following result is universal, as it depends only on the fact that the density of states is supported on a single bounded interval.

Theorem. Let ${{\mathcal{X}}}_{N}$ belong to a one-cut β-ensemble ($\beta \gt 0$) and let $\mathrm{supp}\;\rho =[a,b]$ denote the support of the density of states. Then, in the limit $N\to \infty$ the generating function of the covariances (22) is

$$\begin{eqnarray}&&{F}_{[a,b]}(z,\zeta )=\displaystyle \frac{1}{\beta }\displaystyle \sum _{\kappa ,{\ell }=0}^{\infty }{\alpha }_{\kappa ,{\ell }}{z}^{\kappa }{\zeta }^{{\ell }}=\displaystyle \frac{1}{\beta }\displaystyle \frac{z\zeta }{{(z-\zeta )}^{2}}[\displaystyle \frac{2{abz}\zeta -(a+b)(z+\zeta )+2}{2\sqrt{(1-{az})(1-{bz})(1-a\zeta )(1-b\zeta )}}-1].\end{eqnarray} \tag{ 25 }$$

This is the main result of this paper.

Remark. It is worthwhile to emphasize the following properties of ${F}_{[a,b]}(z,\zeta )$.

• (1)  
  ${F}_{[a,b]}(z,\zeta )={F}_{[a,b]}(\zeta ,z)$, hence ${\alpha }_{\kappa ,{\ell }}={\alpha }_{{\ell },\kappa }$, which expresses the fact that the covariance is a symmetric functional.
• (2)  
  One can verify that $\underset{z,\zeta \to 0}{\mathrm{lim}}{F}_{[a,b]}(z,\zeta )=0$. Therefore ${\alpha }_{\kappa ,{\ell }}=0$ whenever $\kappa$ or ${\ell }$ are zero. This is a consequence of the fact that ${X}_{N,0}=1$ is a constant. Moreover, ${F}_{[a,b]}(z,\zeta )$ is continuous in $(z,\zeta )$ and analytic in each variable (keeping the other variable fixed) in a neighbourhood of the origin that does not intersect the cuts of the square roots in the denominator. If a function of complex variables is continuous and separately analytic in each variable, then it is also jointly analytic (Hartogs' theorem [32]); therefore, ${F}_{[a,b]}(z,\zeta )$ is analytic in a neighbourhood of $(z,\zeta )=(0,0)$.
• (3)  
  The limiting covariances are proportionals to $1/\beta$, which is true whenever $V(x)$ is independent of β (see, e.g. [7, 28]).
• (4)  
  The generating function (25) satisfies the homogeneity identity ${F}_{[{ta},{tb}]}(z,\zeta )={F}_{[a,b]}({tz},t\zeta )$ for all real t. Therefore, under the map $(a,b)\mapsto ({ta},{tb})$ the limiting covariances transform as ${\alpha }_{\kappa ,{\ell }}\mapsto {t}^{\kappa +{\ell }}{\alpha }_{\kappa ,{\ell }}$. This implies that the correlation coefficients ${r}_{\kappa ,{\ell }}={\alpha }_{\kappa ,{\ell }}/\sqrt{{\alpha }_{\kappa ,\kappa }{\alpha }_{{\ell },{\ell }}}$ are invariant in absolute value. More precisely

  $$\begin{eqnarray}&&{r}_{\kappa ,{\ell }}\mapsto {(\mathrm{sign}(t))}^{\kappa +{\ell }}{r}_{\kappa ,{\ell }}\quad {\rm{as}}\quad (a,b)\mapsto ({ta},{tb}).\end{eqnarray} \tag{ 26 }$$
• (5)  
  When the support of the equilibrium density is centred at the origin, i.e. $\mathrm{supp}\;\rho =[-b,b]$ with $b\gt 0$, then the generating function simplifies considerably:

  $$\begin{eqnarray}&&{F}_{\pm b}(z,\zeta )=\displaystyle \frac{1}{\beta }\displaystyle \frac{z\zeta }{{(z-\zeta )}^{2}}[\displaystyle \frac{1-{b}^{2}z\zeta }{\sqrt{(1-{b}^{2}{z}^{2})(1-{b}^{2}{\zeta }^{2})}}-1].\end{eqnarray} \tag{ 27 }$$

  Furthermore, it satisfies ${F}_{\pm b}(z,\zeta )={F}_{\pm b}(-z,-\zeta )$; as a consequence ${\alpha }_{\kappa ,{\ell }}=0$ whenever $\kappa \ne {\ell }\mathrm{mod}2$.
• (6)  
  The support $[a,b]$ of the density of states can always be centred at the origin by a shift $\lambda \mapsto \lambda -m$, where $m=(a+b)/2$. This corresponds to a translation in the matrix model ${{\mathcal{X}}}_{N}\mapsto {{\mathcal{X}}}_{N}-{{mI}}_{N}$ where I_N is the identity matrix. It is known that the 2D Coulomb gas interaction is invariant under conformal transformations [17]. In other words, the covariance structure depends uniquely on the fact that the support is a single bounded interval. The translation ${{\mathcal{X}}}_{N}\mapsto {{\mathcal{X}}}_{N}-{{mI}}_{N}$ induces the change of variables $x\mapsto x-m$ and $y\mapsto y-m$ into the integral in the rhs of (24), which leads to the relation

  $$\begin{eqnarray}&&{\alpha }_{\kappa ,{\ell }}={m}^{\kappa +{\ell }}\displaystyle \sum _{p=0}^{\kappa }\displaystyle \sum _{q=0}^{{\ell }}\displaystyle (\displaystyle \genfrac{}{}{0em}{}{\kappa }{p})\displaystyle (\displaystyle \genfrac{}{}{0em}{}{{\ell }}{q}){(\displaystyle \frac{L}{2m})}^{p+q}{\alpha }_{p,q}^{{\mathcal{G}}},\end{eqnarray} \tag{ 28 }$$

  where $m=(a+b)/2$, $L=(b-a)/2$ and ${\alpha }_{p,q}^{{\mathcal{G}}}$ are the covariances of an equilibrium density supported on $[-2,2]$ (e.g. of the Gaussian ensemble). The ${\alpha }_{p,q}^{{\mathcal{G}}}$ are explicitly given in (42). Thus, (28) gives a tool to evaluate the limiting covariances for an arbitrary one-cut β-ensemble. This is somehow analogous to the usual process of centring the variables in the theory of probability.
• (7)  
  If ${\alpha }_{\kappa ,{\ell }}\in {\mathbb{N}}$ for all $\kappa ,{\ell }$ then $\mathrm{max}(| a| ,| b| )\gt 1$. This is a necessary condition to have positive integers ${\alpha }_{\kappa ,{\ell }}$.

In what follows we present the proof of (25). As in [19], the real integral (24) can be lifted to a double integral in the complex plane

$$\begin{eqnarray}&&{\alpha }_{\kappa ,{\ell }}=\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{\epsilon \prime \to {0}^{+}}{\mathrm{lim}}\displaystyle \frac{1}{{\pi }^{2}}\;{\displaystyle \int }_{{\Gamma }_{\epsilon }}{\rm{d}}z{\displaystyle \int }_{{\Gamma }_{{\epsilon }^{\prime }}}{\rm{d}}\zeta \;\sqrt{\displaystyle \frac{(z-a)(z-b)}{(\zeta -a)(\zeta -b)}}\displaystyle \frac{\kappa }{\zeta -z}{\zeta }^{{\ell }}{z}^{\kappa -1},\end{eqnarray} \tag{ 29 }$$

where ${\Gamma }_{\epsilon }$, ${\Gamma }_{{\epsilon }^{\prime }}$ are two clockwise oriented contours enclosing the cut $[a,b]$ at distance and $\epsilon \prime$ respectively (see figure 1). The first integration with respect to ζ is performed with z fixed inside the cut. In this way the singularity $\zeta =z$ becomes irrelevant. One verifies that the result of the first integration (a function of z) is analytic in a neighbourhood of the cut and hence integrable on the contour ${\Gamma }_{\epsilon }$ enclosing the cut. With this prescription in mind, for $\kappa ,{\ell }\geqslant 0$, we get using residues

$$\begin{eqnarray}&&{\alpha }_{\kappa ,{\ell }}=\kappa \times \mathrm{Res}(\sqrt{\displaystyle \frac{(z-a)(z-b)}{(\zeta -a)(\zeta -b)}}\displaystyle \frac{1}{z-\zeta }{\zeta }^{{\ell }}{z}^{\kappa -1};\zeta =\infty ,z=\infty )\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&=\;\kappa \times \mathrm{Res}(\sqrt{\displaystyle \frac{(1-{az})(1-{bz})}{(1-a\zeta )(1-b\zeta )}}\displaystyle \frac{\zeta }{z-\zeta }\displaystyle \frac{1}{{\zeta }^{{\ell }+1}}\displaystyle \frac{1}{{z}^{\kappa +1}};\zeta =0,z=0).\end{eqnarray} \tag{ 31 }$$

Zoom In Zoom Out Reset image size

From (31) we have

$$\begin{eqnarray}&&{\alpha }_{\kappa ,{\ell }}=\kappa \times {\rm{coefficient}}\;{\rm{of}}\;{\rm{}}{z}^{\kappa }{\zeta }^{{\ell }}{\rm{in}}\;{\rm{the}}\;{\rm{expansion}}\;{\rm{of}}\\ &&\quad (\sqrt{\displaystyle \frac{(1-{az})(1-{bz})}{(1-a\zeta )(1-b\zeta )}}\displaystyle \frac{\zeta }{z-\zeta }){\rm{about}}\;{\rm{}}(z,\zeta )=(0,0)\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&={\rm{coefficient}}\;{\rm{of}}\;{\rm{}}{z}^{\kappa }{\zeta }^{{\ell }}{\rm{in}}\;{\rm{the}}\;{\rm{expansion}}\;{\rm{of}}\;{\rm{}}z\;{\partial }_{z}\\ &&\quad (\sqrt{\displaystyle \frac{(1-{az})(1-{bz})}{(1-a\zeta )(1-b\zeta )}}\displaystyle \frac{\zeta }{z-\zeta }){\rm{}}\;{\rm{about}}\;{\rm{}}(z,\zeta )=(0,0).\end{eqnarray} \tag{ 33 }$$

In the last line we have used the observation that if $C(z)={\displaystyle \sum }_{\kappa }{c}_{\kappa }{z}^{\kappa }$, then the generating function of $\kappa {c}_{\kappa }$ is $z{\partial }_{z}C(z)$. Hence, the generating function $F(z,\zeta )=(1/\beta ){\displaystyle \sum }_{\kappa ,{\ell }\geqslant 0}{\alpha }_{\kappa ,{\ell }}{z}^{\kappa }{\zeta }^{{\ell }}$ is given by

$$\begin{eqnarray}&&{F}_{[a,b]}(z,\zeta )=\displaystyle \frac{1}{\beta }\times {\rm{analytic}}\;{\rm{part}}\;{\rm{at}}(z,\zeta )=(0,0)\\ &&\;{\rm{of}}\;{\rm{the}}\;{\rm{function}}\;{\rm{}}z\;{\partial }_{z}(\sqrt{\displaystyle \frac{(1-{az})(1-{bz})}{(1-a\zeta )(1-b\zeta )}}\displaystyle \frac{\zeta }{z-\zeta }).\end{eqnarray} \tag{ 34 }$$

Carrying out the differentiation and isolating the analytic part, one gets the final expression (25).

We also have an alternative derivation of (25), which relies on the chain of loop equations for the v-point resolvent

$$\begin{eqnarray}&&W({z}_{1},...,{z}_{v})={C}_{v}(\displaystyle \frac{1}{N}\mathrm{Tr}\displaystyle \frac{1}{{z}_{1}-{{\mathcal{X}}}_{N}},...,\displaystyle \frac{1}{N}\mathrm{Tr}\displaystyle \frac{1}{{z}_{v}-{{\mathcal{X}}}_{N}}).\end{eqnarray} \tag{ 35 }$$

When the equilibrium density has square root singularities at the edges of the support (soft edges) these quantities admit a power expansion in $1/N$ (we refer to [11, 15] for a more detailed exposition) and, under certain hypotheses, one can compute the v-point resolvents recursively order by order by solving a hierarchy of loop equations (for an algorithmic approach see [10]). Within this formalism, the generating function (25) can be expressed as ${F}_{[a,b]}(z,\zeta )=z\zeta {W}_{2}^{(0)}(1/z,1/\zeta )$ where ${W}_{2}^{(0)}$ is the leading order in $1/N$ of the two-point resolvent. Using the notation in [10], we have

$$\begin{eqnarray}&&{F}_{[a,b]}(z,\zeta )=\displaystyle \frac{2}{\beta z\zeta }\;\omega ({j}_{a,b}^{-1}(1/z),{j}_{a,b}^{-1}(1/\zeta )),\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&\omega (x,y)=\displaystyle \frac{1}{{({j}_{a,b}^{\prime }(x)-{j}_{a,b}^{\prime }(y))}^{2}}\displaystyle \frac{1}{{({xy}-1)}^{2}},\end{eqnarray} \tag{ 37 }$$

where ${j}_{a,b}$ is the Joukowski transformation defined as ${j}_{a,b}(z)=\frac{a+b}{2}+\frac{a-b}{4}(z+\frac{1}{z})$ (in (36) and (37) ${j}_{a,b}^{-1}$ and ${j}_{a,b}^{\prime }$ of denote the inverse map and the derivative of ${j}_{a,b}$, respectively).

In this section we focus on the classical β-ensembles: the Gaussian, Wishart and Jacobi β-ensembles. Fluctuations of traces in the Gaussian and Wishart β-ensembles were considered by Dumitriu and Edelman [25] in the framework of the tridiagonal realization of β-ensemble that they introduced [24] (for $\beta =1,2$ the tridiagonalization of random matrices had already appeared earlier in [48] and [43]). Using a combinatorial technique simplified by the sparseness of their matrix models, they managed to obtain formulae for the averages and covariances of the moments (among other results). Their formulae are insensitive to the particular value of $\beta \gt 0$ (since β appears as just a parameter in the distribution of off-diagonal terms). Later, Dumitriu and Paquette [26] investigated the very same statistics on the Jacobi ensemble. The matrix realization of this ensemble is not tridiagonal [36]; this turns out to be a serious obstruction to obtain explicit formulae.

Here we show that our approach allows to treat these cases on the same footing: neither the specific matrix model nor the concrete matrix representation do matter. We will recover previous results, and produce new self-contained formulae for the limiting covariances of the Wishart and Jacobi β-ensemble as particular cases of the general formula (25).

4.1. Gaussian ensemble and symmetric intervals $[a,b]=[-L,L]$

The density of states of the β-Gaussian ensemble is symmetric with respect to the origin. It is, therefore, a natural place to start. The eigenvalues of β-Gaussian matrices ${{\mathcal{G}}}_{N}$ are distributed according to

$$\begin{eqnarray}&&{\rm{d}}{{\mathbb{P}}}_{N,\beta }(\{{\lambda }_{k}\})\propto {| \Delta (\{{\lambda }_{k}\})| }^{\beta }\displaystyle \prod _{k=1}^{N}{{\rm{e}}}^{-N\beta {\lambda }_{k}^{2}/4}{\rm{d}}{\lambda }_{k},\end{eqnarray} \tag{ 38 }$$

where $\Delta (\{{\lambda }_{k}\})={\displaystyle \prod }_{j\lt k}({\lambda }_{k}-{\lambda }_{j})$ is the Vandermonde determinant. Equation (38) is easily recognized as the canonical measure of a 2D Coulomb gas in a quadratic potential $V(x)={x}^{2}/4$. When $\beta =1,2,4$ (38) is the joint probability density function of the eigenvalues of real symmetric (GOE), complex Hermitian (GUE) and quaternion self-dual (GSE) gaussian random matrices, respectively.

For any $\beta \gt 0$ the average of the empirical density (3) of the eigenvalues $\{{\lambda }_{i}\}$ converges in expectation to Wigner's semicircle law

$$\begin{eqnarray}&&{\rho }_{{\mathcal{G}}}(x)=\displaystyle \frac{1}{2\pi }\sqrt{4-{x}^{2}}\;{{\bf{1}}}_{| x| \lt 2}.\end{eqnarray} \tag{ 39 }$$

Therefore, the averages of the moments ${\alpha }_{\kappa }^{{\mathcal{G}}}={\mathrm{lim}}_{N\to \infty }\langle {G}_{N,\kappa }\rangle$ converge to the Catalan numbers

$$\begin{eqnarray}{\alpha }_{\kappa }^{{\mathcal{G}}}=\int {\rm{d}}x{\rho }_{{\mathcal{G}}}(x){x}^{\kappa }=\{\begin{array}{cc}\displaystyle \frac{2}{\kappa +2}\displaystyle (\displaystyle \genfrac{}{}{0em}{}{\kappa }{\displaystyle \frac{\kappa }{2}}) & {\rm{}}\;{\rm{for}}\;{\rm{}}\;\kappa =0{\rm{}}\;{\rm{even}}\;{\rm{}},\\ 0 & {\rm{}}\;{\rm{for}}\;{\rm{}}\;\kappa =0{\rm{}}\;{\rm{odd.}}\;{\rm{}}\end{array}\end{eqnarray} \tag{ 40 }$$

The limiting density ${\rho }_{{\mathcal{G}}}$ is supported on the interval $[a,b]=[-2,2]$. Thus, formula (25) reads

$$\begin{eqnarray}&&{F}_{\pm 2}(z,\zeta )=\displaystyle \frac{1}{\beta }\displaystyle \frac{z\zeta }{{(z-\zeta )}^{2}}[\displaystyle \frac{1-4z\zeta }{\sqrt{(1-4{z}^{2})(1-4{\zeta }^{2})}}-1].\end{eqnarray} \tag{ 41 }$$

The coefficients ${\alpha }_{\kappa ,{\ell }}^{{\mathcal{G}}}$ of this double series can be first extracted by taking kth and ${\ell }\mathrm{th}$ order derivatives of (41) symbolically, and recognizing an easy pattern. This leads eventually to the well-supported ansatz

$$\begin{eqnarray}{\alpha }_{\kappa ,{\ell }}^{{\mathcal{G}}}=\underset{N\to \infty }{\mathrm{lim}}\beta {N}^{2}\mathrm{Cov}({G}_{N,\kappa },{G}_{N,{\ell }})=\{\begin{array}{cc}\displaystyle \frac{4\kappa {\ell }}{\kappa +{\ell }}\displaystyle (\displaystyle \genfrac{}{}{0em}{}{\kappa -1}{\lfloor \displaystyle \frac{\kappa }{2}\rfloor })\displaystyle (\displaystyle \genfrac{}{}{0em}{}{{\ell }-1}{\lfloor \displaystyle \frac{{\ell }}{2}\rfloor })\quad & {\rm{}}\;{\rm{if}}\;{\rm{}}\;\kappa ={\ell }\;\;\mathrm{mod}\;\;2,\\ 0 & {\rm{}}\;{\rm{if}}\;{\rm{}}\;\kappa \ne {\ell }\;\;\mathrm{mod}2,\end{array}\end{eqnarray} \tag{ 42 }$$

where $\lfloor x\rfloor$ is the largest integer smaller than or equal to x. Then, it remains to prove that the guessed form (42) is indeed correct by resumming the series ${\displaystyle \sum }_{\kappa ,{\ell }}{\alpha }_{\kappa ,{\ell }}^{{\mathcal{G}}}{z}^{\kappa }{\zeta }^{{\ell }}$—this is outlined in appendix B.

The first values of these covariances are given explicitly in table 1. We mention that (42) can be tracked back to a work by Tutte [50] in combinatorics (see appendix A).

Table 1.  A block of limiting covariance structure ($1\leqslant \kappa ,{\ell }\leqslant 8$) for the β-Gaussian ensemble from (42).

$$\begin{eqnarray*}\underset{N\to \infty }{\mathrm{lim}}{N}^{2}\mathrm{Cov}({G}_{N,\kappa },{G}_{N,{\ell }})=\displaystyle \frac{1}{\beta }[\begin{array}{cccccccc}2 & 0 & 6 & 0 & 20 & 0 & 70 & 0\\ 0 & 4 & 0 & 16 & 0 & 60 & 0 & 224\\ 6 & 0 & 24 & 0 & 90 & 0 & 336 & 0\\ 0 & 16 & 0 & 72 & 0 & 288 & 0 & 1120\\ 20 & 0 & 90 & 0 & 360 & 0 & 1400 & 0\\ 0 & 60 & 0 & 288 & 0 & 1200 & 0 & 4800\\ 70 & 0 & 336 & 0 & 1400 & 0 & 5600 & 0\\ 0 & 224 & 0 & 1120 & 0 & 4800 & 0 & 19600\end{array}]\end{eqnarray*}$$

More generally, using the homogeneity property of ${F}_{[a,b]}(z,\zeta )$ (see remark III), for any one-cut β ensemble ${{\mathcal{X}}}_{N}$ whose limiting spectral density is supported on a single symmetric interval $[-L,L]$ we have

$$\begin{eqnarray}&&\underset{N\to \infty }{\mathrm{lim}}{N}^{2}\mathrm{Cov}({X}_{N,\kappa },{X}_{N,{\ell }})\\ \quad =\displaystyle \frac{1}{\beta }\{\begin{array}{cc}{(\displaystyle \frac{L}{2})}^{\kappa +{\ell }}\displaystyle \frac{4\kappa {\ell }}{\kappa +{\ell }}\displaystyle (\displaystyle \genfrac{}{}{0em}{}{\kappa -1}{\lfloor \displaystyle \frac{\kappa }{2}\rfloor })\displaystyle (\displaystyle \genfrac{}{}{0em}{}{{\ell }-1}{\lfloor \displaystyle \frac{{\ell }}{2}\rfloor }) & {\rm{}}\;{\rm{if}}\;{\rm{}}\;\kappa ={\ell }\;\;\mathrm{mod}2\\ 0\quad & {\rm{}}\;{\rm{if}}\;{\rm{}}\;\kappa \ne {\ell }\;\;\mathrm{mod}2,\end{array},\quad {\rm{if}}\;[a,b]=[-L,L].\end{eqnarray} \tag{ 43 }$$

For L = 1 this expression agrees with the formula derived by Dumitriu and Edelman [25] with a much more complicated combinatorial proof ⁹ .

It is worth emphasizing that all the ensembles with a symmetric support of the equilibrium density share the same L- and β-independent correlation matrix (see remark III),

$$\begin{eqnarray}{r}_{\kappa ,{\ell }}=\displaystyle \frac{{\alpha }_{\kappa ,{\ell }}}{\sqrt{{\alpha }_{\kappa ,\kappa }{\alpha }_{{\ell },{\ell }}}}=\{\begin{array}{cc}2\displaystyle \frac{\sqrt{\kappa {\ell }}}{\kappa +{\ell }}\quad & {\rm{}}\;{\rm{if}}\;{\rm{}}\;\kappa ={\ell }\;\;\mathrm{mod}2,\\ 0 & {\rm{if}}\;\kappa \ne {\ell }\;\mathrm{mod}2,\end{array}\qquad {\rm{if}}\;{\rm{}}[a,b]=[-L,L].\end{eqnarray} \tag{ 44 }$$

4.2. Wishart ensemble and intervals $[a,b]=[0,2L]$

The same direct method that we used for the Gaussian β-ensemble can be applied to the β-Wishart (Laguerre) ensemble of positive semi-definite random matrices ${{\mathcal{W}}}_{N}\geqslant 0$, characterized by a joint probability density of the eigenvalues

$$\begin{eqnarray}&&{\rm{d}}{{\mathbb{P}}}_{N,\beta }(\{{\lambda }_{k}\})\propto {| \Delta (\{{\lambda }_{k}\})| }^{\beta }\displaystyle \prod _{k=1}^{N}{\lambda }_{k}^{\beta (c-1)N}{{\rm{e}}}^{-\displaystyle \frac{\beta }{2}N{\lambda }_{k}}\theta ({\lambda }_{k})\;{\rm{d}}{\lambda }_{k},\qquad (c\geqslant 1).\end{eqnarray} \tag{ 45 }$$

For $\beta =1,2,4$, this ensemble can be realized as ${{\mathcal{W}}}_{N}=\frac{1}{N}{{\mathcal{G}}}_{N\times M}{{\mathcal{G}}}_{N\times M}^{\dagger }$, where ${{\mathcal{G}}}_{N\times M}$ denotes a Gaussian $N\times M$ matrix with i.i.d. real, complex or quaternion entries, respectively. The parameter $c\geqslant 1$ in (45) is the ratio $M/N$, which is assumed to remain constant as $M,N\to \infty$. We notice that (45) is the canonical measure of a 2D Coulomb gas in the (non-polynomial) external potential $V(x)=x/2-(c-1)\mathrm{ln}x$ for nonnegative x and infinite on the negative half-line.

The average of the empirical spectral density of this ensemble converges to the Marčenko–Pastur law with parameter $c\geqslant 1$:

$$\begin{eqnarray}&&{\rho }_{{\mathcal{W}}}(x)=\displaystyle \frac{1}{2\pi x}\sqrt{(x-a)(b-x)}\;{{\bf{1}}}_{x\in (a,b)},\qquad {\rm{with}}\;{\rm{}}\quad [a,b]=[{(1-\sqrt{c})}^{2},{(1+\sqrt{c})}^{2}].\end{eqnarray} \tag{ 46 }$$

The limit on average ${\alpha }_{\kappa }^{{\mathcal{W}}}={\mathrm{lim}}_{N\to \infty }\langle {W}_{N,\kappa }\rangle$ of the moments ${W}_{N,\kappa }={N}^{-1}\mathrm{Tr}{{\mathcal{W}}}_{N}^{\kappa }$ exists and is described in terms of the Narayana numbers:

$$\begin{eqnarray}&&{\alpha }_{\kappa }^{{\mathcal{W}}}=\displaystyle \sum _{p=1}^{\kappa }{c}^{p-1}\mathrm{Nar}(\kappa ,p),\qquad \mathrm{Nar}(\kappa ,p)=\displaystyle \frac{1}{p}\displaystyle (\displaystyle \genfrac{}{}{0em}{}{\kappa }{p-1})\displaystyle (\displaystyle \genfrac{}{}{0em}{}{\kappa -1}{p-1}).\end{eqnarray} \tag{ 47 }$$

Formula (28) combined with (42) gives

$$\begin{eqnarray}{\alpha }_{\kappa ,{\ell }}^{{\mathcal{W}}} & = & {\mathrm{lim}}_{N\to \infty }\beta {N}^{2}\mathrm{Cov}({W}_{N,\kappa },{W}_{N,{\ell }})\\ & = & 4{(1+c)}^{\kappa +{\ell }}\displaystyle \sum _{\begin{array}{c}0\leqslant p\leqslant \kappa \\ 0\leqslant q\leqslant {\ell }\\ p=q\mathrm{mod}2\end{array}}{(\displaystyle \frac{\sqrt{c}}{1+c})}^{p+q}\displaystyle \frac{{pq}}{p+q}\displaystyle (\displaystyle \genfrac{}{}{0em}{}{\kappa }{p})\displaystyle (\displaystyle \genfrac{}{}{0em}{}{{\ell }}{q})\displaystyle (\displaystyle \genfrac{}{}{0em}{}{p-1}{\lfloor \displaystyle \frac{p}{2}\rfloor })\displaystyle (\displaystyle \genfrac{}{}{0em}{}{q-1}{\lfloor \displaystyle \frac{q}{2}\rfloor }).\end{eqnarray} \tag{ 48 }$$

Few values of ${\alpha }_{\kappa ,{\ell }}^{{\mathcal{W}}}$ are reported in table 2.

Table 2.  Block of the limiting covariance structure ($1\leqslant \kappa ,{\ell }\leqslant 3$) of the β-Wishart ensemble from (48).

$$\begin{eqnarray*}&&\underset{N\to \infty }{\mathrm{lim}}{N}^{2}\mathrm{Cov}({W}_{N,\kappa },{W}_{N,{\ell }})\\ =\displaystyle \frac{1}{\beta }[\begin{array}{ccc}2c & 4(c+{c}^{2}) & 6(c+3{c}^{2}+{c}^{3})\\ 4(c+{c}^{2}) & 4(2c+5{c}^{2}+2{c}^{3}) & 12(c+5{c}^{2}+5{c}^{3}+{c}^{4})\\ 6(c+3{c}^{2}+{c}^{3}) & 12(c+5{c}^{2}+5{c}^{3}+{c}^{4}) & 6(3c+24{c}^{2}+46{c}^{3}+24{c}^{4}+3{c}^{5})\end{array}]\end{eqnarray*}$$

The β-Wishart ensemble was also investigated by Dumitru and Edelman [25] within the framework of the tridiagonal realization of the ensemble. Their proof requires an ad hoc combinatorial analysis, which is different from the β-Gaussian ensemble. Their formula (see [25], claim 3.15.2) is considerably more involved than (48), and besides, we found that it contains few misprints.

When the support is the interval $[0,2L]$, the generating function of the ${\alpha }_{\kappa ,{\ell }}$ simplifies drastically

$$\begin{eqnarray}&&{F}_{[\mathrm{0,2}L]}(z,\zeta )=\displaystyle \frac{1}{\beta }\displaystyle \frac{z\zeta }{{(z-\zeta )}^{2}}[\displaystyle \frac{1-L(z+\zeta )}{\sqrt{(1-2{Lz})(1-2L\zeta )}}-1].\end{eqnarray} \tag{ 49 }$$

Extracting the Taylor coefficients of this formula involves lengthy calculations. However, for c = 1 the moments of Wishart matrices are related in a simple way to the moments of Gaussian matrices and therefore one obtains a self-contained formula for the covariances

$$\begin{eqnarray}&&\underset{N\to \infty }{\mathrm{lim}}{N}^{2}\mathrm{Cov}({X}_{N,\kappa },{X}_{N,{\ell }})=\displaystyle \frac{1}{\beta }{(\displaystyle \frac{L}{2})}^{\kappa +{\ell }}\displaystyle \frac{4\kappa {\ell }}{\kappa +{\ell }}\displaystyle (\displaystyle \genfrac{}{}{0em}{}{2\kappa -1}{\kappa })\displaystyle (\displaystyle \genfrac{}{}{0em}{}{2{\ell }-1}{{\ell }}),\qquad {\rm{if}}\;{\rm{}}[a,b]=[0,2L].\end{eqnarray} \tag{ 50 }$$

Setting $2L=4$ we get the covariance structure ${\alpha }_{\kappa ,{\ell }}^{{\mathcal{W}}}$ of the symmetric ($c=1$) Wishart ensemble.

The correlation coefficients are independent of the length of the interval and are given

$$\begin{eqnarray}&&{r}_{\kappa ,{\ell }}=\displaystyle \frac{{\alpha }_{\kappa ,{\ell }}}{\sqrt{{\alpha }_{\kappa ,\kappa }{\alpha }_{{\ell },{\ell }}}}=2\displaystyle \frac{\sqrt{\kappa {\ell }}}{\kappa +{\ell }},\qquad {\rm{if}}\;{\rm{}}[a,b]=[0,2L].\end{eqnarray} \tag{ 51 }$$

In figure 2 we report Monte Carlo simulations of correlation coefficients for ensembles of type $[0,2L]$. From (50) we can derive the large $\kappa ,{\ell }$ behaviour of the covariances

$$\begin{eqnarray}&&\underset{N\to \infty }{\mathrm{lim}}{N}^{2}\mathrm{Cov}({X}_{N,\kappa },{X}_{N,{\ell }})\sim \displaystyle \frac{{(2L)}^{\kappa +{\ell }}}{\beta \pi }\displaystyle \frac{\sqrt{\kappa {\ell }}}{\kappa +{\ell }}\qquad {\rm{if}}\;{\rm{}}[a,b]=[0,2L],\end{eqnarray} \tag{ 52 }$$

in agreement with an earlier prediction of two of the authors of this paper [17].

Zoom In Zoom Out Reset image size

4.3. Jacobi ensemble

As last example we consider the β-Jacobi ensemble. These matrices satisfy the constraint $0\leqslant {{\mathcal{J}}}_{N}\leqslant I$ and their eigenvalues $0\leqslant {\lambda }_{k}\leqslant 1$ have joint probability density function

$$\begin{eqnarray}&&{\rm{d}}{{\mathbb{P}}}_{N,\beta }(\{{\lambda }_{k}\})\propto {| \Delta (\{{\lambda }_{k}\})| }^{\beta }\displaystyle \prod _{k=1}^{N}{\lambda }_{k}^{\displaystyle \frac{\beta }{2}N{\gamma }_{1}}{(1-{\lambda }_{k})}^{\displaystyle \frac{\beta }{2}N{\gamma }_{2}}\theta ({\lambda }_{k}(1-{\lambda }_{k}))\;{\rm{d}}{\lambda }_{k},\qquad ({\gamma }_{\mathrm{1,2}}\geqslant 0).\end{eqnarray} \tag{ 53 }$$

The density of states of this ensemble is

$$\begin{eqnarray}&&{\rho }_{{\mathcal{J}}}(x)=\displaystyle \frac{{\gamma }_{1}+{\gamma }_{2}+2}{2\pi x(1-x)}\sqrt{(x-a)(b-x)}\;{{\bf{1}}}_{x\in (a,b)}.\end{eqnarray} \tag{ 54 }$$

It is compactly supported on the interval $[a({\gamma }_{1},{\gamma }_{2}),b({\gamma }_{1},{\gamma }_{2})]$, whose edges are

$$\begin{eqnarray}a({\gamma }_{1},{\gamma }_{2}) & = & {[\displaystyle \frac{\sqrt{{\gamma }_{2}+1}-\sqrt{({\gamma }_{1}+1)({\gamma }_{1}+{\gamma }_{2}+1)}}{{\gamma }_{1}+{\gamma }_{2}+2}]}^{2},\\ b({\gamma }_{1},{\gamma }_{2}) & = & {[\displaystyle \frac{\sqrt{{\gamma }_{2}+1}+\sqrt{({\gamma }_{1}+1)({\gamma }_{1}+{\gamma }_{2}+1)}}{{\gamma }_{1}+{\gamma }_{2}+2}]}^{2}.\end{eqnarray} \tag{ 55 }$$

By inserting (55) into (25) one obtains the generating function of the limiting covariances of ${J}_{N,\kappa }=\mathrm{Tr}{{\mathcal{J}}}_{N}^{\kappa }$ for any value of ${\gamma }_{\mathrm{1,2}}\geqslant 0$. The ${\alpha }_{\kappa ,{\ell }}^{{\mathcal{J}}}$'s can then be evaluated explicitly via (28) as

$$\begin{eqnarray}{\alpha }_{\kappa ,{\ell }}^{{\mathcal{J}}} & = & {[\displaystyle \frac{{\gamma }_{1}^{2}+{\gamma }_{1}{\gamma }_{2}+2({\gamma }_{1}+{\gamma }_{2}+1)}{{({\gamma }_{1}+{\gamma }_{2}+2)}^{2}}]}^{\kappa +{\ell }}\\ & & \times \displaystyle \sum _{\begin{array}{c}0\leqslant p\leqslant \kappa \\ 0\leqslant q\leqslant {\ell }\end{array}}{[\displaystyle \frac{\sqrt{({\gamma }_{1}+1)({\gamma }_{2}+1)({\gamma }_{1}+{\gamma }_{2}+1)}}{{\gamma }_{1}^{2}+{\gamma }_{1}{\gamma }_{2}+2({\gamma }_{1}+{\gamma }_{2}+1)}]}^{p+q}\displaystyle (\displaystyle \genfrac{}{}{0em}{}{\kappa }{p})\displaystyle (\displaystyle \genfrac{}{}{0em}{}{{\ell }}{q}){\alpha }_{p,q}^{{\mathcal{G}}},\end{eqnarray} \tag{ 56 }$$

where ${\alpha }_{p,q}^{{\mathcal{G}}}$ is given in (42). Furthermore, in the particular case ${\gamma }_{1}={\gamma }_{2}=0$, we have $[a,b]=[0,1]$; hence, using formula (50) we arrive at

$$\begin{eqnarray}{\alpha }_{\kappa ,{\ell }}^{{\mathcal{J}}} & = & \underset{N\to \infty }{\mathrm{lim}}\beta {N}^{2}\mathrm{Cov}({J}_{N,\kappa },{J}_{N,{\ell }})\\ & = & {4}^{1-\kappa -{\ell }}\displaystyle \frac{\kappa {\ell }}{\kappa +{\ell }}\displaystyle (\displaystyle \genfrac{}{}{0em}{}{2\kappa -1}{\kappa })\displaystyle (\displaystyle \genfrac{}{}{0em}{}{2{\ell }-1}{{\ell }})\quad {\rm{if}}\;{\rm{}}{\gamma }_{\mathrm{1,2}}=0\quad ([a,b]=[0,1]).\end{eqnarray} \tag{ 57 }$$

Few values are reported in table 3.

Table 3.  A block of the limiting covariance structure ($1\leqslant \kappa ,{\ell }\leqslant 5$) of the Jacobi ensemble for ${\gamma }_{\mathrm{1,2}}=0$ from (57).

$$\begin{eqnarray*}\underset{N\to \infty }{\mathrm{lim}}{N}^{2}\mathrm{Cov}({J}_{N,\kappa },{J}_{N,{\ell }})=\displaystyle \frac{1}{\beta }[\begin{array}{ccccc}\displaystyle \frac{1}{8} & \displaystyle \frac{1}{8} & \displaystyle \frac{15}{128} & \displaystyle \frac{7}{64} & \displaystyle \frac{105}{1024}\\ \displaystyle \frac{1}{8} & \displaystyle \frac{9}{64} & \displaystyle \frac{9}{64} & \displaystyle \frac{35}{256} & \displaystyle \frac{135}{1024}\\ \displaystyle \frac{15}{128} & \displaystyle \frac{9}{64} & \displaystyle \frac{75}{512} & \displaystyle \frac{75}{512} & \displaystyle \frac{4725}{32\;768}\\ \displaystyle \frac{7}{64} & \displaystyle \frac{35}{256} & \displaystyle \frac{75}{512} & \displaystyle \frac{1225}{8192} & \displaystyle \frac{1225}{8192}\\ \displaystyle \frac{105}{1024} & \displaystyle \frac{135}{1024} & \displaystyle \frac{4725}{32\;768} & \displaystyle \frac{1225}{8192} & \displaystyle \frac{19\;845}{131072}\end{array}]\end{eqnarray*}$$

Dumitriu and Paquette (see [26], equation (23)) derived an integral formula for ${\alpha }_{\kappa ,{\ell }}^{{\mathcal{J}}}$ that can be evaluated numerically; one can then check that the numerical integration gives the same values as formula (56). Their proof, however, does not extend to the limiting values ${\gamma }_{\mathrm{1,2}}=0$, which are covered by ours.

Using a Coulomb gas approach, we have derived a formula (25) for the generating function of the covariances of the moments in the thermodynamic limit $N\to \infty$ for one-cut β-ensembles of random matrices. This result is universal, as it depends only on the support of the limiting spectral density, but not on the potential $V(x)$ that defines the ensemble, nor on the inverse temperature β; this is a direct consequence of the universality of the smoothed two-point spectral correlation function. Given a β-ensemble whose equilibrium density has support $[a,b]$, the covariances of the moments can always be expressed in terms of a set of centred covariances, i.e. corresponding to a limiting density with support centred at the origin. This allows us to write an explicit self-contained formula (28) for the limiting covariances of any one-cut β-ensemble. Since the proof is based on the Dyson Coulomb gas analogy it is independent of the particular matrix realization of the ensemble and as such is simpler and more transparent than methods previously available.

The classical β-Gaussian, β-Wishart and β-Jacobi ensembles can be treated in a unified way. In this way we recover results previously established [24, 25] with proofs based on sparse matrix realization of β-ensembles. In addition we have proved new formulae for the β-Wishart and the β-Jacobi ensembles. We mention that the issues of second order cumulants and joint behavior of spectral linear statistics have recently attracted attention in the context of the theory of quantum transport in mesoscopic cavities [19, 20, 30] as well as in practical spectral inference [18, 23].

We conclude by mentioning several related open problems. The challenge posed by higher cumulants is still open and the promising features of the 2D Coulomb gas approach suggest that the same line of reasoning could be pursued to investigate the universality of the v-point correlation kernel for $v\gt 2$. Concerning the first and second order cumulants, it is natural to look for the next to leading term corrections. In [35], Johansson showed how to compute the subleading terms for the averages at any $\beta \gt 0$ (for V polynomial). We are not aware of similar results for the $o({N}^{-2})$ corrections of the covariance structures. Coming back to leading order effects, it would be interesting to obtain an analogue of (25) for multi-cut β-ensembles where the universality of the two-point smoothed kernel is more delicate [1, 2, 9]. As a last remark, we point out that the Coulomb gas approach is tailored to invariant ensembles. Wigner matrices (other than Gaussian) escape from the unitary invariance and this is a serious problem when discussing the v-point correlation kernels. For Wigner matrices it is known that, while the density of states is universal (the Wigner law), the higher-order correlation kernels are somewhat less universal [21] (they depend on more details of the entries distribution). The solution to these problems may well uncover yet another layer of 'universality' features of 'classical' RMT, a field from which many lessons are surely still to be learnt.

FDC and FM are grateful to G Borot for enlightening discussions on the loop equations technique and for suggestion of a few relevant references. This work was partially supported by EPSRC Grant number EP/L010305/1. FDC acknowledges partial support of Gruppo Nazionale di Fisica Matematica GNFM-INdAM. No empirical or experimental data were created during this study. PV acknowledgs support from EPSRC Center for Doctoral Training in Cross-Disciplinary Approaches to Non-Equlilibrium Systems (CANES).

Enumeration of maps, i.e. graphs drawn on certain surfaces, random matrices and 2D Coulomb gases are inter-related topics in mathematical physics. The electrostatic analogy suggested by Wigner [54] and Dyson [27] in the late 1950s has proved extremely useful in many branches of RMT. On the other hand, extensive elaborations on the connections between random matrices and the enumeration of maps did not begin until the works in the 1970s by 't Hooft's [46], Koplin et al [37] and Brézin et al [12] on the planar approximation in various field theories. Once this connection was brought to the surface, matrix models and field theoretical techniques have been exploited to make progress on unsolved graphical enumeration problems [8, 33, 34, 41]. (A comprehensive account can be found in the classical [22, 55].)

Diagrammatic techniques have been successfully applied to a few specific but important ensembles, namely the Gaussian ensembles, see for instance [51]. The crucial reason is that computing expectation values in this ensemble is tantamount to applying Wick's formula in various settings; furthermore, when pairing the Gaussian entries, a nice geometrical structure emerges. These patterns have been rediscovered many times by different communities and mastering the relevant vast literature can be a daunting task for the uninitiated. The aim of this appendix is twofold. Firstly, we summarize the general ideas connecting matrix integrals and planar enumeration of maps, and at the same time we try to clarify different technical terminologies from various communities. Secondly, we collect results available in the mathematical physics literature and compare them with our findings. Throughout this Appendix, we will concentrate on the GUE ensemble and on the diagrammatic techniques associated to the computation of its moments. For GUE, the joint law of the eigenvalues (1) is the canonical distribution of a 2D Coulomb gas confined in a harmonic potential $V(x)={x}^{2}/4$ at inverse temperature $\beta =2$.

The idea that the large N asymptotic expansion of a matrix integral around a saddle point can be represented by diagrams identifiable with certain maps was pioneered by Brézin et al [12]. In order to study the distribution of ${G}_{N,\kappa }=\mathrm{Tr}{{\mathcal{G}}}_{N}^{\kappa }$, for fixed $\kappa$, it is customary to consider the function

$$\begin{eqnarray}&&\mathrm{exp}\{-{N}^{2}{E}_{N,\kappa }(s)\}=\langle \mathrm{exp}\{-{N}^{2}s\;{G}_{N,\kappa }\}\rangle ,\end{eqnarray} \tag{ A.1 }$$

where hereafter the angle brackets denote the average with respect to the GUE measure. From the statistical independence of the Gaussian entries of ${{\mathcal{G}}}_{N}$, it is known [12, 33] that ${E}_{N,\kappa }(s)$ admits a $1/{N}^{2}$-expansion

$$\begin{eqnarray}&&{E}_{N,\kappa }(s)=\displaystyle \sum _{g\geqslant 0}\displaystyle \frac{1}{{N}^{2g}}{E}_{\kappa }^{(g)}(s),\end{eqnarray} \tag{ A.2 }$$

where, for all $g\geqslant 0$, ${E}_{\kappa }^{(g)}(s)$ is analytic in a neighbourhood of s = 0,

$$\begin{eqnarray}&&{E}_{\kappa }^{(g)}(s)=-\displaystyle \sum _{v\geqslant 1}\displaystyle \frac{{(-s)}^{v}}{v!}{D}_{\kappa }^{(g)}(v).\end{eqnarray} \tag{ A.3 }$$

In the above expansion, ${D}_{\kappa }^{(g)}(v)$ is the number of inequivalent connected diagrams of genus g with v vertices of valence $\kappa$.

Therefore, in the large N limit only the g = 0 contribution survives

$$\begin{eqnarray}&&\underset{N\to \infty }{\mathrm{lim}}{E}_{N,\kappa }(s)={E}_{\kappa }^{(0)}(s)=-\displaystyle \sum _{v\geqslant 1}\displaystyle \frac{{(-s)}^{v}}{v!}{D}_{\kappa }^{(0)}(v).\end{eqnarray} \tag{ A.4 }$$

The quantity ${E}_{\kappa }^{(0)}(s)$ is often referred to as cumulant generating function of the κth moment ${G}_{N,\kappa }$. The coefficients in this expansion count the connected diagrams of a ${\varphi }^{\kappa }$-theory in which all non planar diagrams are omitted: ${D}_{\kappa }^{(0)}(v)$ is the number of connected diagrams with v vertices all of valence $\kappa$ that can be drawn on a plane (genus g = 0)¹⁰ . The probabilistic meaning of these numbers is clear from (A.1): if we denote by ${C}_{v}({G}_{N,\kappa },...,{G}_{N,\kappa })$ the vth cumulant of ${G}_{N,\kappa }$, then by differentiating (A.1) we arrive at

$$\begin{eqnarray}&&\underset{N\to \infty }{\mathrm{lim}}{N}^{2(v-1)}{C}_{v}({\underbrace{{G}_{N,\kappa },...,{G}_{N,\kappa }}_{v\;{\rm{times}}}})={D}_{\kappa }^{(0)}(v).\end{eqnarray} \tag{ A.5 }$$

For instance, the first two cumulants (average and variance) of ${G}_{N,\kappa }$ are given by the number of planar diagrams with one and two $\kappa$-valent vertices, respectively

$$\begin{eqnarray}&&\langle {G}_{N,\kappa }\rangle ={D}_{\kappa }^{(0)}(1)+o(1),\end{eqnarray} \tag{ A.6 }$$

$$\begin{eqnarray}&&{N}^{2}\mathrm{Var}({G}_{N,\kappa })={D}_{\kappa }^{(0)}(2)+o(1).\end{eqnarray} \tag{ A.7 }$$

We emphasize the origin of the two main features of these diagrams: the connectedness is intrinsic in the definition of ${E}_{N,\kappa }(s)$ (the logarithm of an exponential generating function); the planar approximation g = 0 emerges in the large N limit. Correction terms are encoded in ${E}_{\kappa }^{(g)}(s)$ for higher genera $g\geqslant 1$, that is by considering $\kappa$-valent connected diagrams ${D}_{\kappa }^{(g)}(v)$ on higher genera surfaces (for instance the first correction comes from legal diagrams ${D}_{\kappa }^{(1)}(v)$ on the torus, g = 1). In the one-vertex case (${D}_{\kappa }^{(0)}(1)$ with $\kappa \geqslant 1$), these planar diagrams are in bijection with non-crossing pair partitions (NCPP) of $\kappa$ points. A NCPP of $\kappa$ points can be represented by placing the numbers $1,...,\kappa$ around a single circle and connecting numbers in pairs without crossings. Hence

$$\begin{eqnarray*}&&{D}_{\kappa }^{(0)}(1)\equiv \#\{{\rm{non-crossing}}\;{\rm{pairings}}\;{\rm{of}}\;\kappa \;{\rm{points}}\;{\rm{on}}\;{\rm{a}}\;{\rm{circles}}\}.\end{eqnarray*}$$

The theory of non-crossing partitions is an important topic in modern combinatorics and their role in random matrices has been elucidated by free probability [52, 53], with the notion of free cumulants [40, 44].

One can extend the above identification to diagrams with more vertices $v\geqslant 1$ and write

$$\begin{eqnarray}&&{D}_{\kappa }^{(0)}(v)\\ &&\equiv \#\{{\rm{Non-crossing}}\;{\rm{pairings}}\;{\rm{of}}\;v\;{\rm{circles}}\;{\rm{with}}\;\kappa \;\\ &&{\rm{points}}\;{\rm{on}}\;{\rm{each}}\;{\rm{circle}}\;{\rm{such}}\;{\rm{that}}\;{\rm{all}}\;{\rm{circles}}\;{\rm{are}}\;{\rm{connected}}.\}\end{eqnarray} \tag{ A.8 }$$

(See figure A1 .) The connectedness requirement is clear and the non-crossing condition is the translation of the planarity condition for diagrams. Equations (A.6) and (A.7) thus connect average and variance of a single linear statistics ${G}_{N,\kappa }$ of the Gaussian ensemble at $\beta =2$ to the problem of non-crossing pairings on one circle, or on two circles, respectively.

Zoom In Zoom Out Reset image size

One can generalize the above arguments and consider the joint distribution of $({G}_{N,{\kappa }_{1}},...,{G}_{N,{\kappa }_{n}})$ for $n\geqslant 2$ with distinct ${\kappa }_{i}$'s in general. Using the notation $\vec{\kappa }=({\kappa }_{1},...,{\kappa }_{n})$ and $\vec{s}=({s}_{1},...,{s}_{n})$, we may introduce the multidimensional analog of (A.1)

$$\begin{eqnarray}&&\mathrm{exp}\{-{N}^{2}{E}_{N,\vec{\kappa }}(\vec{s})\}=\langle \mathrm{exp}\{-{N}^{2}\displaystyle \sum _{i=1}^{n}{s}_{i}{G}_{N,{\kappa }_{i}}\}\rangle .\end{eqnarray} \tag{ A.9 }$$

The large N limit of ${E}_{N,\vec{\kappa }}(\vec{s})$ will be the joint cumulant generating function of $({G}_{N,{\kappa }_{1}},...,{G}_{N,{\kappa }_{n}})$,

$$\begin{eqnarray}&&{E}_{\vec{\kappa }}^{(0)}(\vec{s})=-\displaystyle \sum _{v\geqslant 1}\displaystyle \sum _{\stackrel{{v}_{1},...,{v}_{n}\geqslant 0}{{v}_{1}+\cdots +{v}_{n}=v}}[\displaystyle \prod _{i=1}^{n}\displaystyle \frac{{(-{s}_{i})}^{{v}_{i}}}{{v}_{i}!}]{D}_{\vec{\kappa }}^{(0)}({v}_{1},...,{v}_{n}),\end{eqnarray} \tag{ A.10 }$$

where ${D}_{\vec{\kappa }}^{(0)}({v}_{1},...,{v}_{n})$ is now the number of connected planar diagrams with v₁ '${\kappa }_{1}$-valent' vertices, v₂ '${\kappa }_{2}$-valent' vertices, and so on. The relation with the theory of non-crossing partitions is now

$$\begin{eqnarray}{D}_{\vec{\kappa }}^{(0)}({v}_{1},...,{v}_{n}) & \equiv & \#\{{\rm{Non-crossing}}\;{\rm{pairings}}\;{\rm{of}}\;{\rm{}}v={v}_{1}+\cdots +{v}_{n}\\ & & {\rm{circles}}\;{\rm{with}}\;{\kappa }_{1}\;{\rm{points}}\;{\rm{on}}\;{\rm{}}{v}_{1}{\rm{circles,}}\\ & & {\kappa }_{2}\;{\rm{points}}\;{\rm{on}}\;{v}_{2}\;{circles},\ldots ,\\ & & {\kappa }_{n}{\rm{}}\;{\rm{points}}\;{\rm{on}}\;{v}_{n}{\rm{circles}}\;{\rm{such}}\;{\rm{that}}\;{\rm{all}}\;{\rm{circles}}\;{\rm{are}}\;{\rm{connected.}}\}\end{eqnarray} \tag{ A.11 }$$

Clearly, no such pairing is possible if the total number of points ${\displaystyle \sum }_{i=1}^{n}{v}_{i}{\kappa }_{i}$ is odd. These numbers provide the leading order of the joint mixed cumulants $({G}_{N,{\kappa }_{1}},...,{G}_{N,{\kappa }_{n}})$ of order $v={v}_{1}+\cdots +{v}_{n}$,

$$\begin{eqnarray}&&\underset{N\to \infty }{\mathrm{lim}}{N}^{2(v-1)}{C}_{v}(\underset{{v}_{1}\;{\rm{times}}}{\underbrace{{G}_{N,{\kappa }_{1}},...,{G}_{N,{\kappa }_{1}}}},...,\underset{{v}_{n}\;{\rm{times}}}{\underbrace{{G}_{N,{\kappa }_{n}},...,{G}_{N,{\kappa }_{n}}}})={D}_{\vec{\kappa }}^{(0)}({v}_{1},...,{v}_{n}).\end{eqnarray} \tag{ A.12 }$$

From (A.12), one can see that for the GUE ensemble the family ${\{{G}_{N,\kappa }\}}_{\kappa \geqslant 0}$ is asymptotically Gaussian (the case $\kappa =0$ is degenerate) as N goes to infinity. In the language of free probability, this means that the GUE ensemble has a second order limit distribution [16, 38, 39], i.e. the limits

$$\begin{eqnarray}&&{\alpha }_{\kappa }^{{\mathcal{G}}}=\underset{N\to \infty }{\mathrm{lim}}\langle {G}_{N,\kappa }\rangle \quad {\rm{and}}\quad (1/2){\alpha }_{\kappa ,{\ell }}^{{\mathcal{G}}}=\underset{N\to \infty }{\mathrm{lim}}{N}^{2}[\langle {G}_{N,\kappa }{G}_{N,{\ell }}\rangle -\langle {G}_{N,\kappa }\rangle \langle {G}_{N,{\ell }}\rangle ]\end{eqnarray} \tag{ A.13 }$$

exist for all $\kappa ,{\ell }\geqslant 0$, and the higher-order cumulants decay faster to zero for large N.

Brézin et al [12] considered the planar approximation for $\kappa =3$ and $\kappa =4$ (quartic and cubic vertices). They computed exactly the cumulant generating functions ${E}_{3}^{(0)}(s)$ and ${E}_{4}^{(0)}(s)$ thus providing the full family of cumulant of $\mathrm{Tr}{{\mathcal{G}}}^{3}$ and $\mathrm{Tr}{{\mathcal{G}}}^{4}$ to leading order in N. For the quartic interaction, the higher genus corrections ${E}_{4}^{(1)}(s)$ and ${E}_{4}^{(2)}(s)$ have been computed explicitly in [8]. The problem of mixed interaction (diagrams whose vertices have mixed valence) has been less explored even in the planar regime g = 0—the case of the combined cubic and quartic interactions, i.e. the computation of ${E}_{3,4}^{(0)}({s}_{1},{s}_{2})$ was proposed in [12] but was not worked out explicitly. Later, Harer and Zagier [31] computed the exact finite N value of $\langle {G}_{N,\kappa }\rangle$ (i.e. the linear term in s of ${E}_{N,\kappa }(s)$ in (A.1) for all $\kappa$).

Our results (when restricted to the GUE ensemble) provide a complete answer to the following problem: for any finite family of power traces $({G}_{N,{\kappa }_{1}},...,{G}_{N,{\kappa }_{n}})$ compute the joint cumulant generating function ${E}_{\vec{\kappa }}^{(0)}(\vec{s})$ in (A.10) up to quadratic terms

$$\begin{eqnarray}&&{E}_{\vec{\kappa }}^{(0)}(\vec{s})=-\displaystyle \sum _{1\leqslant v\leqslant 2}\displaystyle \sum _{\stackrel{{v}_{1},...,{v}_{n}\geqslant 0}{{v}_{1}+\cdots +{v}_{n}=v}}[\displaystyle \prod _{i=1}^{n}\displaystyle \frac{{(-{s}_{i})}^{{v}_{i}}}{{v}_{i}!}]{D}_{\vec{\kappa }}^{(0)}({v}_{1},...,{v}_{n})+\cdots \end{eqnarray} \tag{ A.14 }$$

$$\begin{eqnarray}&&=\displaystyle \sum _{i=1}^{n}{s}_{i}{D}_{{\kappa }_{i}}^{(0)}(1)-\displaystyle \frac{1}{2}\displaystyle \sum _{i,j=1}^{n}{s}_{i}{s}_{j}{D}_{{\kappa }_{i},{\kappa }_{j}}^{(0)}(1,1)+\cdots \end{eqnarray} \tag{ A.15 }$$

This corresponds to considering the contribution from all types of planar diagrams (every kind of combined interactions) with at most two vertices (of arbitrary valence). Here ${D}_{{\kappa }_{i}}^{(0)}(1)$ and ${D}_{{\kappa }_{i},{\kappa }_{j}}^{(0)}(1,1)$ give the leading order behaviour of the averages and covariances of the family $({G}_{N,{\kappa }_{1}},...,{G}_{N,{\kappa }_{n}})$. The results of section 4.1, restricted to $\beta =2$, give

$$\begin{eqnarray}&&\underset{N\to \infty }{\mathrm{lim}}\langle {G}_{N,\kappa }\rangle ={D}_{\kappa }^{(0)}(1)={\alpha }_{\kappa }^{{\mathcal{G}}},\end{eqnarray} \tag{ A.16 }$$

$$\begin{eqnarray}&&\underset{N\to \infty }{\mathrm{lim}}{N}^{2}\mathrm{Cov}({G}_{N,\kappa },{G}_{N,{\ell }})={D}_{\kappa ,{\ell }}^{(0)}(1,1)=(1/2){\alpha }_{\kappa ,{\ell }}^{{\mathcal{G}}}.\end{eqnarray} \tag{ A.17 }$$

The explicit formulae were given in (40) and (42), respectively. Therefore, from the above discussion we have

$$\begin{eqnarray}&&{\alpha }_{\kappa }^{{\mathcal{G}}}=\#\{{\rm{Non-crossing}}\;{\rm{pairings}}\;{\rm{of}}\;\kappa \;{\rm{points}}\;{\rm{on}}\;{\rm{the}}\;{\rm{circle}}\},\end{eqnarray} \tag{ A.18 }$$

$$\begin{eqnarray}(1/2){\alpha }_{\kappa ,{\ell }}^{{\mathcal{G}}} & = & \#\{{\rm{Non-crossing}}\;{\rm{pairings}}\;{\rm{of}}\;{\rm{two}}\;{\rm{circles}}\;{\rm{with}}\;\kappa \;{\rm{points}}\;{\rm{on}}\;{\rm{the}}\;{\rm{1st}}\;{\rm{circle}}\\ & & {\rm{and}}\;{\ell }\;{\rm{points}}\;{\rm{on}}\;{\rm{the}}\;{\rm{2nd}}\;{\rm{circle}}\;{\rm{such}}\;{\rm{that}}\;{\rm{the}}\;{\rm{two}}\;{\rm{circles}}\;{\rm{are}}\;{\rm{connected}}\}.\end{eqnarray} \tag{ A.19 }$$

In figure A2 we report few such diagrams.

Zoom In Zoom Out Reset image size

The difference of our method compared to previous works on the planar approximation in field theories is the following. Instead of expanding the cumulant generating function of a fixed ${G}_{N,\kappa }$ (or a fixed n-uple $({G}_{N,{\kappa }_{1}},...,{G}_{N,{\kappa }_{n}})$) in the number of vertices v as in (A.4)–(A.10), in this paper we have considered the whole family of generating functions of a generic n-uple $({G}_{N,{\kappa }_{1}},...,{G}_{N,{\kappa }_{n}})$ up to quadratic terms In diagrammatic terms, instead of an expansion of planar diagrams with fixed valence $\kappa$ in the number of vertices v, we would obtain a series in the valences of vertices $\kappa$ and ${\ell }$ of planar diagrams with at most two vertices $v\leqslant 2$.

We mention that the combinatorial problem associated to ${\alpha }_{\kappa }^{{\mathcal{G}}}$ and $(1/2){\alpha }_{\kappa ,{\ell }}^{{\mathcal{G}}}$ was solved by Tutte [50] as a particular case of a more general problem on the number of so-called slicings of a band. We note that Tutte's formula enumerates general slicings, with a strong restriction on the parity of the vertices. In the language of RMT, Tutte provided a formula for the generic mixed cumulant ${D}_{\vec{\kappa }}^{(0)}({v}_{1},...,{v}_{n})$ of $({G}_{N,{\kappa }_{1}},...,{G}_{N,{\kappa }_{n}})$ (to leading order in N), with the condition that ${\kappa }_{1},...,{\kappa }_{N}$ should all be even ¹¹ . We find interesting to note that this obstruction (that might appear unjustified in the combinatorial problem) can be explained in field-theory with the instability of any ${\varphi }^{\kappa }$-interaction for $\kappa$ odd.

It is convenient to cast ${\alpha }_{\kappa ,{\ell }}^{{\mathcal{G}}}$ in (42) for $\kappa ,{\ell }\gt 0$ in term of Euler Gamma functions

$$\begin{eqnarray}{\alpha }_{\kappa ,{\ell }}^{{\mathcal{G}}}=\{\begin{array}{cc}\displaystyle \frac{16\Gamma (\kappa )\Gamma ({\ell })}{(\kappa +{\ell })\Gamma {(\displaystyle \frac{\kappa }{2})}^{2}\Gamma {(\displaystyle \frac{{\ell }}{2})}^{2}} & {\rm{}}\;{\rm{for}}\;{\rm{}}{\ell }{\rm{}}\;{\rm{and}}\;{\rm{}}\kappa {\rm{}}\;{\rm{even}}\;{\rm{}},\\ \displaystyle \frac{\kappa {\ell }{4}^{\kappa +{\ell }-1}\Gamma {(\displaystyle \frac{\kappa }{2})}^{2}\Gamma {(\displaystyle \frac{{\ell }}{2})}^{2}}{{\pi }^{2}(\kappa +{\ell })\Gamma (\kappa )\Gamma ({\ell })} & {\rm{}}\;{\rm{for}}\;{\rm{}}{\ell }{\rm{}}\;{\rm{and}}\;{\rm{}}\kappa {\rm{}}\;{\rm{odd}}\;{\rm{}}.\end{array}\end{eqnarray} \tag{ B.1 }$$

Setting $\kappa =2r$ and ${\ell }=2s$ (in the even–even case) and $\kappa =2r+1$ and ${\ell }=2s+1$ (in the odd–odd case), we can evaluate the double sum in (25) as

$$\begin{eqnarray}\displaystyle \sum _{\kappa ,{\ell }=1}^{\infty }{\alpha }_{\kappa ,{\ell }}^{{\mathcal{G}}}{z}^{\kappa }{\zeta }^{{\ell }} & = & \displaystyle \sum _{r,s=1}^{\infty }\displaystyle \frac{8\Gamma (2r)\Gamma (2s)}{(r+s)\Gamma {(r)}^{2}\Gamma {(s)}^{2}}{z}^{2r}{\zeta }^{2s}\\ & & +\displaystyle \sum _{r,s=0}^{\infty }\displaystyle \frac{(2r+1)(2s+1){4}^{2r+2s+1}\Gamma {(\displaystyle \frac{2r+1}{2})}^{2}\Gamma {(\displaystyle \frac{2s+1}{2})}^{2}}{{\pi }^{2}(2r+2s+2)\Gamma (2r+1)\Gamma (2s+1)}{z}^{2r+1}{\zeta }^{2s+1}.\end{eqnarray} \tag{ B.2 }$$

The two series on the rhs can be computed in closed form using the following identities

$$\begin{eqnarray}&&\displaystyle \frac{1}{r+s}={\displaystyle \int }_{0}^{\infty }{\rm{d}}t\;{{\rm{e}}}^{-t(r+s)},\end{eqnarray} \tag{ B.3 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{r=1}^{\infty }\displaystyle \frac{\Gamma (2r)}{\Gamma {(r)}^{2}}{(z{{\rm{e}}}^{-t/2})}^{2r}=\displaystyle \frac{{{\rm{e}}}^{-t}{z}^{2}}{{({{\rm{e}}}^{-t}({{\rm{e}}}^{t}-4{z}^{2}))}^{3/2}},\end{eqnarray} \tag{ B.4 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{r=0}^{\infty }\displaystyle \frac{(2r+1)\Gamma {(\displaystyle \frac{2r+1}{2})}^{2}}{\Gamma (2r+1)}{(4z{{\rm{e}}}^{-t})}^{2r}=\displaystyle \frac{\pi {{\rm{e}}}^{2t}}{({{\rm{e}}}^{2t}-4{z}^{2})\sqrt{1-4{{\rm{e}}}^{-2t}{z}^{2}}}.\end{eqnarray} \tag{ B.5 }$$

Therefore we can rewrite the double sum in integral form

$$\begin{eqnarray}\displaystyle \sum _{\kappa ,{\ell }=1}^{\infty }{\alpha }_{\kappa ,{\ell }}^{{\mathcal{G}}}{z}^{\kappa }{\zeta }^{{\ell }} & = & 8{\displaystyle \int }_{0}^{\infty }{\rm{d}}t\displaystyle \frac{{{\rm{e}}}^{-t}{z}^{2}}{{({{\rm{e}}}^{-t}({{\rm{e}}}^{t}-4{z}^{2}))}^{3/2}}\displaystyle \frac{{{\rm{e}}}^{-t}{\zeta }^{2}}{{({{\rm{e}}}^{-t}({{\rm{e}}}^{t}-4{\zeta }^{2}))}^{3/2}}\\ & & +\displaystyle \frac{4}{{\pi }^{2}}z\zeta {\displaystyle \int }_{0}^{\infty }{\rm{d}}t\ {{\rm{e}}}^{-2t}\displaystyle \frac{\pi {{\rm{e}}}^{2t}}{({{\rm{e}}}^{2t}-4{z}^{2})\sqrt{1-4{{\rm{e}}}^{-2t}{z}^{2}}}\displaystyle \frac{\pi {{\rm{e}}}^{2t}}{({{\rm{e}}}^{2t}-4{\zeta }^{2})\sqrt{1-4{{\rm{e}}}^{-2t}{\zeta }^{2}}}.\end{eqnarray} \tag{ B.6 }$$

The integrals can be carried out in closed form, and after simplifications we obtain

$$\begin{eqnarray}&&\displaystyle \sum _{\kappa ,{\ell }=1}^{\infty }{\alpha }_{\kappa ,{\ell }}^{{\mathcal{G}}}{z}^{\kappa }{\zeta }^{{\ell }}=-\displaystyle \frac{\zeta z(\sqrt{1-4{\zeta }^{2}}\sqrt{1-4{z}^{2}}+4\zeta z-1)}{\sqrt{1-4{\zeta }^{2}}\sqrt{1-4{z}^{2}}{(z-\zeta )}^{2}},\qquad {\rm{}}\;{\rm{for}}\;{\rm{}}-1/2\lt z,\zeta \lt 1/2,\end{eqnarray} \tag{ B.7 }$$

which (after multiplying by $1/\beta$) coincides with ${F}_{\pm 2}(z,\zeta )$ as given in (41).