A unified fluctuation formula for one-cut $\beta$-ensembles of random matrices

Using a Coulomb gas approach, we compute the generating function of the covariances of power traces for one-cut $\beta$-ensembles of random matrices in the limit of large matrix size. This formula depends only on the support of the spectral density, and is therefore universal for a large class of models. This allows us to derive a closed-form expression for the limiting covariances of an arbitrary one-cut $\beta$-ensemble. As particular cases of the main result we consider the classical $\beta$-Gaussian, $\beta$-Wishart and $\beta$-Jacobi ensembles, for which we derive previously available results as well as new ones within a unified simple framework. We also discuss the connections between the problem of trace fluctuations for the Gaussian Unitary Ensemble and the enumeration of planar maps.


I. INTRODUCTION
Many textbooks on random matrices [2,3,39] begin with the derivation of Wigner's semicircle law for the Gaussian Unitary Ensemble (GUE), 1 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,κ = N −1 TrG κ N , where G N denotes an N × N GUE matrix and κ is an integer, and in applying Wick's formula to evaluate the averages G N,κ . 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 lim N →∞ G N,κ are non-zero only for κ even and given by the Catalan numbers. This computation implies that the empirical density ρ N (x) = N −1 i δ(x − λ i ) of the eigenvalues {λ i } converges in expectation to Wigner's law ρ(x) = (2π) −1 √ 4 − x 2 . The simplicity of this result might leads to believe that the very same approach (known as the method of moments) should be easily extended to higher cumulants of G N,κ . However, it turns out that this problem is highly nontrivial, with very few results available in the literature [8,12,13,19,20,31,32,35]. To date, almost 60 years after Wigner's original paper [48], a self-contained explicit formula for the higher cumulants C v (G N,κ1 , . . . , G N,κv ) for generic v and {κ i } is still lacking 2 . 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 X N whose eigenvalues {λ k } ∈ R have the canonical distribution of a 2D Coulomb gas at inverse temperature β > 0, namely 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 3 includes the classical Gaussian, Wishart and Jacobi ensembles of RMT. For β = 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 β > 0 can be realized from sparse real random matrices, as shown by Dumitriu and Edelman [18] and Killip and Nenciu [28]. Moreover, in view of certain applications to physics (see, e.g., [22]), 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,κ = N −1 TrX κ N for matrices 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 → ∞. More precisely we study the quantities The 1/β-dependence and the N −2 scaling of such covariances are customary in the fluctuations for 2D Coulomb gases [4]. Our main result is a universal formula for the generating function of (2); it is universal in the sense that it depends only on the support 1 For convenience and consistence with the formalism used for the problem of the enumeration of maps (see Appendix A), throughout this paper the GUE is defined by complex Hermitian N × N matrices G N whose diagonal and off-diagonal elements are normal random variables with the same variance, i.e. G ij = 1 √ N x ij , with x ij ∼ N C (0, 1) for i < j , and x ii ∼ N R (0, 1) for all i. 2 Cv(ξ 1 , . . . , ξv) denotes the cumulant of order v of the (not necessarily distinct) random variables ξ 1 , . . . , ξv. For instance C 2 (ξ 1 , ξ 2 ) = ξ 1 ξ 2 − ξ 1 ξ 2 . 3 The Cauchy-Stieltjes transform of the gas distribution has one connected cut on the real line of the complex plane, hence the name one-cut. of the equilibrium density, but not on the potential V (x). This is a direct consequence of 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 (27). 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 β > 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 [6,9,25,26,34]. Some of these results reproduce earlier findings by Tutte [44] for planar enumeration of maps. Johansson [27] 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 [13,31,32,35], 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 [19] and Dumitriu and Paquette [20] 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 [7,8,12].
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 II). In Section III we present our main result. In Section IV we apply our methods to the classical β-ensembles: the Gaussian, Laguerre and Jacobi β-ensembles. Finally, in Section V we discuss some open problems. The paper is complemented by an appendix where we review the connections between the problem of moments fluctuations, enumeration of planar maps and non-crossing pair partitions.

II. 2D COULOMB GAS AND CUMULANTS OF LINEAR STATISTICS
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 [21,48]. The joint law P N,β ({λ k }) in (1) is the Gibbs-Boltzmann measure of a 2D Coulomb gas confined in a interval Λ ⊂ R at inverse temperature β > 0 (the Dyson index) with a single-particle potential V (λ) bounded from below and finite for λ ∈ Λ (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 [27] 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 (λ 1 , . . . , λ N ) of the gas particles (the eigenvalues of the β-ensemble) in terms of the empirical density of the Coulomb gas, The partition function of the gas may be written in terms of the empirical density where is the mean-field energy functional and the functional integration is restricted to the probability measures σ compatible with the Coulomb gas configuration {λ i }. A spectral linear statistics (or linear statistics for short) is defined by where f (x) is a suitably chosen test function. Sometimes (6) is referred to as sum function on the particles positions {λ i }. The average of A N is given by In a similar way one computes the variance of A N , or in general the covariance between two linear statistics A The same reasoning can be extended to higher order cumulants. In general, the mixed cumulant of v linear statistics (not necessarily distinct) can be recovered by integrating A (j) , also known as v-point connected correlation function: For v = 1 and v = 2 this corresponds to (7) and (8), respectively. In the asymptotic regime as N → ∞, the connected correlation functions of the gas are generated by the free energy − log Z N,β by taking a suitable number of functional derivatives with respect to the external potential. More precisely, for every positive integer v: For instance, the 1-point and 2-point connected correlation functions in the limit as N → ∞ are given by Cov(ρ N (x)ρ N (y)) = 1 βN 2 Under the assumption that the limit lim N →∞ ρ N (x) = ρ(x) exists, from (10) we see that exists too 4 . 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: The quantity ρ(x) is the density of states, or equilibrium density, while K(x, y) is known as smoothed two-point connected correlation kernel (see [5,10] for details). Therefore, the limiting mean and covariances of smooth linear statistics (the lowest order cumulants) are 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 ρ(x) (14) exists, is absolutely continuous and supported on a single bounded interval supp ρ = [a, b], (a < b), and is the unique probability measure that minimizes the energy functional (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, ρ(x) and K(x, y) have the closed form where P denotes Cauchy's principal value. The first result (18) is usually referred to as Tricomi's formula [41]. The beautiful identity (19) is a consequence of Tricomi's formula and (15), and has been independently discovered by several authors [1,5,10,11,23]. It is worth emphasizing that while the one-point density (18) depends explicitly on the potential V (x), the smoothed two-point connected correlation kernel (19) enjoys a very surprising feature: it depends only on the edge points a, b of the spectral density. Hence, it is a universal quantity. This is a key observation that eventually leads to our main result, presented in next section.

III. MAIN RESULT
Let X N,κ = N −1 TrX κ N denote the moments of a generic one-cut β-ensemble of random matrices. For any fixed positive integers κ and , the limits exist and are given by respectively. The identities (22)-(23) are a specialization of (16)-(17)- (19). While the limiting average α κ depends explicitly on the probability measure that defines the ensemble through V (x), and hence ρ(x) (given in (18)), the covariance (1/β)α κ, does not. Indeed, the following theorem is universal, as it depends only on the fact that the density of states is supported on a single bounded interval.
Theorem. Let X N belong to a one-cut β-ensemble (β > 0) and let supp ρ = [a, b] denote the support of the density of states. Then, in the limit N → ∞ the generating function of the covariances (21) is This is the main result of this paper.
Remark. It is worthwhile to emphasize the following properties of F [a,b] (z, ζ).
, hence α κ, = α ,κ , which expresses the fact that the covariance is a symmetric functional.
6. The support [a, b] of the density of states can always be centred at the origin by a shift λ → λ − m, where m = (a + b) /2. This corresponds to a translation in the matrix model X N → 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 [14]. In other words, the covariance structure depends uniquely on the fact that the support is a single bounded interval. The translation X N → X N − mI N induces the change of variables x → x − m and y → y − m into the integral in the r.h.s. of (23), which leads to the relation where m = (a + b)/2, L = (b − a)/2 and α G p,q are the covariances of an equilibrium density supported on [−2, 2] (e.g. of the Gaussian ensemble). The α G p,q are explicitly given in (41). Thus, (27) 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.
In what follows we present the proof of (24). As shown in [15], we can lift (23) to a double integral in the complex plane, where Γ , Γ are two clockwise oriented contours enclosing the cut [a, b] at distance and respectively (see Fig. 1). For κ, ≥ 0, by a residue calculation we get In the above steps the first integration with respect to ζ is performed with z fixed inside the cut. In this way the singularity ζ = z becomes irrelevant. From (30) we have In the last line we have used the observation that if C(z) = κ c κ z κ , then the generating function of κc κ is z∂ z C(z). Hence, the generating function F (z, ζ) = (1/β) κ, ≥0 α κ, z κ ζ is given by Contours of integration and scheme of the prescription used to evaluate (28). First the integral and the limit → 0 + in the ζ-plane with z inside the cut [a, b]. Then the integration in the z-plane.
Carrying out the differentiation and isolating the analytic part, one gets the final expression (24). We also have an alternative derivation of (24), which relies on the chain of loop equations for the v-point resolvent 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 [8,12] 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 [7]). Within this formalism, the generating function (24) can be expressed as is the leading order in 1/N of the 2-point resolvent. Using the notation in [7], we have ω(x, y) = 1 where j a,b (z) = a+b 2 + a−b 4 z + 1 z is the Joukowski transformation.

IV. APPLICATIONS TO THE CLASSICAL β-ENSEMBLES
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 [19] in the framework of the tridiagonal realization of β-ensemble that they introduced [18] (for β = 1, 2 the tridiagonalization of random matrices had already appeared earlier in [42] and [36]). 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 β > 0 (since β appears as just a parameter in the distribution of off-diagonal terms). Later, Dumitriu and Paquette [20] investigated the very same statistics on the Jacobi ensemble. The matrix realization of this ensemble is not tridiagonal [28]; 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 (24). 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 G N are distributed according to where ∆({λ k }) = j<k (λ k − λ j ) is the Vandermonde determinant. Eq. (37) is easily recognized as the canonical measure of a 2D Coulomb gas in a quadratic potential V (x) = x 2 /4. When β = 1, 2, 4 (37) 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 β > 0 the average of the empirical density (3) of the eginvalues {λ i } converges in expectation to Wigner's semicircle law, Therefore, the averages of the moments α G κ = lim N →∞ G N,κ converge to the Catalan numbers The limiting density ρ G is supported on the interval [a, b] = [−2, 2]. Thus, formula (24) reads It is then just a matter of lengthy computations to extract the coefficients α G κ, of this double series and deduce where x is the largest integer smaller than or equal to x. The first values of these covariances are given explicitly in Table I. We mention that (41) can be tracked back to a work by Tutte [44] in combinatorics (see Appendix A).  More generally, using the homogeneity property of F [a,b] (z, ζ) (see Remark III), for any one-cut β ensemble X N whose limiting spectral density is supported on a single symmetric interval [−L, L] we have For L = 1 this expression agrees with the formula derived by Dumitriu and Edelman [19] with a much more complicated combinatorial proof. 5 It is worth emphasizing that all the ensembles with a symmetric support of the equilibrium density share the same Land β-independent correlation matrix (see Remark III), 4 c + c 2 6 c + 3c 2 + c 3 4 c + c 2 4 2c + 5c 2 + 2c 3 12 c + 5c 2 + 5c 3 + c 4 6 c + 3c 2 + c 3 12 c + 5c 2 + 5c 3 + c 4 6 3c + 24c 2 + 46c 3 + 24c 4 + 3c 5   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 W N ≥ 0, characterized by a joint probability density of the eigenvalues For β = 1, 2, 4, this ensemble can be realized as W N = 1 N G N ×M G † N ×M , where G N ×M denotes a Gaussian N × M matrix with real, complex or quaternion entries, respectively. The parameter c ≥ 1 in (44) is the ratio M/N , which is assumed to remain constant as M, N → ∞. We notice that (44) is the canonical measure of a 2D Coulomb gas in the (non-polynomial) external potential V (x) = x/2 − (c − 1) 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 ≥ 1: The limit on average α W κ = lim N →∞ W N,κ of the moments W N,κ = N −1 TrW κ N exists and is described in terms of the Narayana numbers: Formula (27) Few values of α W κ are reported in Table II. The β-Wishart ensemble was also investigated by Dumitru and Edelman [19] 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 [19], Claim 3.15.2) is considerably more involved than (47), and besides, we found that it contains few misprints.
When the support is the interval [0, 2L], the generating function of the α κ, simplifies drastically Extracting the Taylor coefficients of this formula involves lengthy calculations. Eventually, one obtains a self-contained formula for the covariances, Setting 2L = 4 we get the covariance structure α W κ, of the symmetric (c = 1) Wishart ensemble. The correlation coefficients are independent of the length of the interval and are given In Fig. 2 in agreement with an earlier prediction of two of the authors of this paper [14].

C. Jacobi ensemble
As last example we consider the β-Jacobi ensemble. These matrices satisfy the constraint 0 ≤ J N ≤ I and their eigenvalues 0 ≤ λ k ≤ 1 have joint probability density function The density of states of this ensemble is It is compactly supported on the interval [a(γ 1 , γ 2 ), b(γ 1 , γ 2 )], whose edges are By inserting (54) into (24) one obtains the generating function of the limiting covariances of J N,κ = TrJ κ N for any value of γ 1,2 ≥ 0. The α J κ, 's can then be evaluated explicitly via (27) as where α G p,q is given in (41). Furthermore, in the particular case γ 1 = γ 2 = 0, we have [a, b] = [0, 1]; hence, using formula (49) we arrive at Few values are reported in Table III. Dumitriu and Paquette (see [20], Eq. (23)) derived an integral formula for α J κ, that can only be evaluated numerically; one can then check that the numerical integration gives the same values as formula (55). Their proof, however, does not extend to the limiting values γ 1,2 = 0, which are covered by ours.

V. CONCLUSIONS
Using a Coulomb gas approach, we have derived a formula (24) for the generating function of the covariances of the moments in the thermodynamic limit N → ∞ 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 (27) 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 [18,19] with proofs based on sparse matrix realization of β-ensembles. In addition we have proved new formulae for the β-Wishart and the β-Jacobi ensembles.
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 > 2. Concerning the first and second order cumulants, it is natural to look for the next to leading term corrections. In [27], Johansson showed how to compute the subleading terms for the averages at any β > 0 (for V polynomial). We are not aware of similar results for the o(N −2 ) corrections of the covariance structures.
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 correlations kernels are somewhat less universal [16] (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.

Acknowledgments
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.
Appendix A: Gaussian Matrix integrals, Planar diagrams, and 2D Coulomb gas at β = 2 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 [48] and Dyson [21] in the late 1950's 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 1970's by 't Hooft's [40], Koplin et al. [29] and Brézin et al. [9] 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 [6,25,26,34]. (A comprehensive account can be found in the classical references [17,49].) Diagrammatic techniques have been successfully applied to a few specific but important ensembles, namely the Gaussian ensembles, see for instance [45]. The crucial reason is that computing expectation values in this ensembles 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 β = 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. [9]. In order to study the distribution of G N,κ = TrG κ N , for fixed κ, it is customary to consider the function where hereafter the angle brackets denote the average with respect to the GUE measure. From the statistical independence of the Gaussian entries of G N , it is known [9,25] that E N,κ (s) admits a 1/N 2 -expansion where, for all g ≥ 0, E (g) κ (s) is analytic in a neighbourhood of s = 0, In the above expansion, D (g) κ (v) is the number of inequivalent connected diagrams of genus g with v vertices of valence κ. Therefore, in the large N limit only the g = 0 contribution survives, The quantity E For instance, the first two cumulants (average and variance) of G N,κ are given by the number of planar diagrams with one and two κ-valent vertices, respectively, We emphasize the origin of the two main features of these diagrams: the connectedness is intrinsic in the definition of E N,κ (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 (g) κ (s) for higher genera g ≥ 1, that is by considering κ-valent connected diagrams D (g) κ (v) on higher genera surfaces (for instance the first correction comes from legal diagrams D (1) κ (v) on the torus, g = 1). In the one-vertex case (D (0) κ (1) with κ ≥ 1), these planar diagrams are in bijection with non-crossing pair partitions (NCPP) of κ points. A NCPP of κ points can be represented by placing the numbers 1, . . . , κ around a single circle and connecting numbers in pairs without crossings. Hence 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 [46,47], with the notion of free cumulants [33,37]. One can extend the above identification to diagrams with more vertices v ≥ 1 and write D (0) κ (v) ≡ #{Non-crossing pairings of v circles with κ points on each circle such that all circles are connected.} (A8) (See Fig. 3.) The connectedness requirement is clear and the non-crossing condition is the translation of the planarity condition for diagrams. Eqs. (A6) and (A7) thus connect average and variance of a single linear statistics G N,κ of the Gaussian ensemble at β = 2 to the problem of non-crossing pairings on one circle, or on two circles, respectively. 6 In the mathematical literature, these planar diagrams are also known as planar maps with v stars (each star has κ points).
One can generalize the above arguments and consider the joint distribution of (G N,κ1 , . . . , G N,κn ) for n ≥ 2 with distinct κ i 's in general. Using the notation κ = (κ 1 , . . . , κ n ) and s = (s 1 , . . . , s n ), we may introduce the multidimensional analog of (A1) (A9) The large N limit of E N, κ ( s) will be the joint cumulant generating function of (G N,κ1 , . . . , G N,κn ), where D From (A12), one can see that for the GUE ensemble the family {G N,κ } κ≥0 is asymptotically Gaussian (the case κ = 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 [13,31,32], i.e. the limits exist for all κ, ≥ 0, and the higher-order cumulants decay faster to zero for large N . Brézin et al. [9] considered the planar approximation for κ = 3 and κ = 4 (quartic and cubic vertices). They computed exactly the cumulant generating functions E   4 (s) have been computed explicitly in [6]. 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 (0) 3,4 (s 1 , s 2 ) was proposed in [9] but was not worked out explicitly. Later, Harer and Zagier [24] computed the exact finite N value of G N,κ (i.e. the linear term in s of E N,κ (s) in (A1) for all κ).
In Fig. 4 we report few such diagrams. 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,κ (or a fixed n-uple (G N,κ1 , . . . , G N,κn )) in the number of vertices v as in (A4)-(A10), in this paper we have considered the whole family of generating functions of a generic n-uple (G N,κ1 , . . . , G N,κn ) up to quadratic terms. In diagrammatic terms, instead of an expansion of planar diagrams with fixed valence κ in the number of vertices v, we would obtain a series in the valences of vertices κ and of planar diagrams with at most two vertices v ≤ 2.
We mention that the combinatorial problem associated to α G κ and (1/2)α G κ, was solved by Tutte [44] as a particular case of a more general problem on the number of so-called slicings of a band. (This is yet another terminology for planar diagrams or non-crossing pairings of circles.) 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 (0) κ (v 1 , . . . , v n ) of (G N,κ1 , . . . , G N,κn ) (to leading order in N ), with the condition that κ 1 , . . . , κ N should all be even 7 . 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 ϕ κ -interaction for κ odd.