Large N (=3) Neutrinos and Random Matrix Theory

The large N limit has been successfully applied to QCD, leading to qualitatively correct results even for N=3. In this work, we propose to treat the number N=3 of Standard Model generations as a large number. Specifically, we apply this idea to the neutrino anarchy scenario and study neutrino physics using Random Matrix Theory, finding new results in both areas. For neutrino physics, we obtain predictions for the masses and mixing angles as a function of the generation number N. The Seesaw mechanism produces a hierarchy of order 1/N^3 between the lightest and heaviest neutrino, and a theta(13) mixing angle of order 1/N, in parametric agreement with experimental data when N goes to 3. For Random Matrix Theory, this motivates the introduction of a new type of ensemble of random matrices, the"Seesaw ensemble."Basic properties of such matrices are studied, including the eigenvalue density and the interpretation as a Coulomb gas system. Besides its mathematical interest, the Seesaw ensemble may be useful in random systems where two hierarchical scales exist.


Introduction
The large N limit is a very useful tool in various theoretical models, both at a qualitative and a quantitative level [1]. Historically, one of the most important examples came from 't Hooft's large N analysis of QCD [2]. Although the number of colors N c = 3 is not very large, many features of mesons and baryons calculated in a 1/N c = 1/3 expansion are in good agreement with experiment [3].
One cannot help but notice that in the Standard Model (SM), the number "3" also appears as the number of fermion generations. Here, we propose to treat the number of families N f = 3 in a large N approximation, and study the properties of fermion masses and mixings as a function of the expansion parameter 1/N f = 1/3.
In this paper, we make the first attempt to use the large N limit to understand some part of the flavor problem (i.e. the origin of masses and mixings of the SM fermions). More specifically, we will apply this idea to the SM neutrinos assuming the mass anarchy scenario of [4,5], where the neutrino couplings are treated as random variables. In light of the recent measurement of θ 13 from Daya Bay [6] and later confirmed by Reno [7], the anarchy hypothesis is consistent with the measured value [8,9].
We will not consider the quark sector here (which requires additional structure), although we will comment briefly on possible large N applications in §6.
Extending the 3 × 3 random matrix for the three active neutrinos to an N × N random matrix (hereafter N = N f denotes the number of generations), we have to analyze the statistical properties of the neutrino spectrum and mixing matrix. This is a classic problem in Random Matrix Theory (RMT); see e.g. [10,11] for reviews. RMT has an incredible range of applications, from medical HIV studies [12] to the string landscape [13,14,15], and now we are adding neutrino physics to this list.
RMT provides analytic results for the correlation functions of eigenvalues of large rank matrices.
For example, the eigenvalues of Hermitian complex matrices obey the well-known Wigner semicircle distribution. Knowing the spectral density for the neutrino masses would already provide valuable information. It will allow us to calculate the expectation value of the ratio of the lightest and heaviest masses, which will in turn determine how hierarchical the mass spectrum is. Although experimentalists have measured all three mixing angles of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix, it is still not clear whether the neutrino mass spectrum is normal or inverted, and hierarchical or quasidegenerate. Our analysis for general N will reveal that the statistical properties of the spectrum depend sensitively on the underlying mass mechanism.
There exist different possibilities for the neutrino mass matrices. Our main focus will be on the Seesaw mechanism, which provides a dynamical explanation for the smallness of neutrino masses. We will demonstrate that the Seesaw mechanism is not only required to explain the hierarchy between the neutrino mass scale and the weak scale, but also preferred by the large N behavior of our RMT analysis and the current experimental results on the neutrino mass-squared differences. Besides providing a useful tool to analyze the spectrum, the connection with RMT will also lead us to define a new type of random matrix ensemble, as we describe in more detail below.
We have attempted to make our paper accessible to both the neutrino physics and RMT communities. We begin in §2 with a short review of the results in RMT that we will need, including the Gaussian unitary ensemble in §2.1 and the Coulomb gas picture in §2.2. Next, we apply RMT techniques to the simplest case of complex Majorana neutrinos in §3. Our main results are presented in §4, where we analyze the large N spectrum of random Seesaw neutrinos, and also define and study the new "Seesaw ensemble." In §5 we apply these results to the phenomenologically relevant N = 3 case and compare with experimental data. Our conclusions and future directions are summarized in §6.

Review of random matrix theory
Let us begin by reviewing some basic results from random matrix theory that will be needed in the following sections. For a pedagogical exposition, we refer the reader to [10,11].

Gaussian unitary ensemble
Random matrix theory is the study of the statistical properties of eigenvalues of very large matrices.
We will first consider N × N Hermitian complex matrices M = A + A † , where each element of A is an independently distributed random variable. The space of such matrices is known as the Gaussian unitary ensemble (GUE).
Writing the probability distribution as P (M )dM , the defining properties of the GUE are that: a) it is invariant under unitary transformations, P (M )dM = P (M ′ )dM ′ , with M ′ = U † M U ; and b) the matrix elements are statistically independent, namely P (M ) = i≤j P ij (M ij ). From this, it follows that [10] where M has been linearly redefined to set the center at zero and to fix the coefficient in the exponential. The matrix M is diagonalized by Here, Λ = diag(m 1 , . . . , m N ) with m i as the eigenvalues of the matrix with the convention m i ≤ m j for i < j and U a special unitary matrix with N (N − 1)/2 independent parameters. Changing variables to Λ and U then gives Here, dU is the invariant Haar measure for U (N ) and C N is a normalization constant.
We see that at large N the distribution of angles and phases contained in U becomes flat in the natural group theory coordinates. This turns out to have important consequences for the neutrino mixing angles, analyzed by [4,5]. On the other hand, the eigenvalues have nontrivial statistical properties. These are characterized by the k-point correlation functions Of particular interest for us will be the one-point function This gives the eigenvalue density, and its inverse 1/ρ(m) determines the mean level spacing around m.
The GUE can be solved exactly in the limit N → ∞. This is done by recognizing that the distribution in Eq. (2.3) can be written as a product of orthogonal Hermite polynomials. The correlators are then of the form where the kernel K N is defined as Here we have introduced the rescaled variables This result implies that at large N the eigenvalue density is given by the Wigner semicircle distribution,

The Coulomb gas picture and generalizations
We will now review the physical interpretation of the GUE as a Coulomb gas, introduced by Dyson [16,17]. This picture will be useful when we study the "Seesaw ensemble" in §4.
The basic idea is that Eq. (2.3) can be interpreted as the thermal partition function for N charged particles with positions (x 1 , . . . , x N ) moving on a fixed line in two dimensions, under the influence of a harmonic potential and electrostatic repulsion. In more detail, the Hamiltonian for this system is Then the partition function Physically, the properties of the eigenvalue distribution appear from a competition between the harmonic potential and the electrostatic repulsion. If the spacing between charged particles decreases with N , at large N the system admits a continuum approximation (2.14) Such an ensemble is equivalent to a system of charged particles with a Hamiltonian The particles interact with each other via the logarithmic Coulomb repulsion, and moreover each of them is subjected to a nontrivial confining potential V (x). The level density ρ(x) of the ensemble can be obtained, as before, by going to a continuum limit at large N and extremizing the analog of Eq. (2.12) with x 2 replaced by V (x). The result is an integral equation for ρ(x), This relation will be applied in §4, where knowledge of the mass distribution generated by the Seesaw mechanism will be used to construct a suitable confining potential V (x).

Majorana neutrinos and 1/N hierarchies
It is helpful to begin our analysis with the simplest case of neutrinos with random Majorana masses.
Results for Dirac neutrinos are qualitatively similar. We consider N neutrinos ν i with Majorana where C is the charge conjugation operator. M ν is an N × N complex symmetric matrix. For the SM, N = 3 and L m comes from a dimension-five operator (HL) 2 . As in [4,5], M ν is taken to be random. Instead of specializing from the start to N = 3, we will determine the properties of neutrino masses at large N and then take N → 3. At large N the lightest neutrino will be found to have a mass suppressed by 1/N . The results from this section will also help to develop intuition and tools to study in §4 the more interesting and nontrivial random Seesaw mechanism case.
We first analyze the statistical properties of random complex symmetric Majorana mass matrices.
In order to find the probability distribution of eigenvalues, we perform a singular value decomposition The measure becomes [5] As reviewed in §2, statistical independence of the elements (M ν ) ij and invariance under unitary transformations implies that the probability distribution is Gaussian. Therefore, the complex Majorana neutrino masses are distributed according to which is known as the Altland-Zirnbauer CI ensemble [18]. The eigenvalues occur in pairs ±m i , which will be important shortly. (This is also why the exponential contains an extra factor of 1/2 as compared to the GUE.) The statistical distribution of neutrino masses is equivalent to a Coulomb gas system with Hamil- The last term encodes an electrostatic repulsion between the images ±x i . This term dominates near the origin, forcing the density of eigenvalues ρ(x) to vanish as x → 0. For x ≪ N −1/2 , the density is a linear function, while for larger x it is well approximated by the semicircle distribution [18]. In more detail, the approximate distribution of neutrino Majorana masses is  at the origin. However, the more precise Eq. (3.7) shows that the level spacing is parametrically of the same order as that of the semicircle law. Similarly, the median is found to be ∼ 1/2.
In summary, we have We have also checked these results numerically, finding the same type of behavior. The phenomenological predictions of this measure were studied in detail in [4,5].
This type of parametric spectrum is not the one suggested by experimental data (reviewed in §5), at least in the limit where large N results give a good approximation to the N = 3 case. For this reason, we next turn to analyze the statistical properties of masses generated by the Seesaw mechanism.

New hierarchies and ensemble from Seesaw neutrinos
We now focus on the more interesting and nontrivial Seesaw models, which have been used to explain the smallness of the neutrino masses with respect to the electroweak scale. We will find that starting from random masses for the light active and heavy singlet neutrinos, the Seesaw mechanism leads to hierarchies of masses that are parametrically different from those obtained in §3. Furthermore, this mechanism will motivate the definition of a new type of random matrix ensemble, the "Seesaw ensemble", whose interesting mathematical properties we will begin to explore here.

Random Seesaw neutrinos
Let us consider the scenario with N active (left-handed) neutrinos, and N singlet (right-handed) neutrinos which have heavy Majorana masses. The case with different numbers N L and N R of leftand right-handed neutrinos also appears to be interesting, offering the possibility of taking the large N R limit while keeping N L fixed. We will comment more on this possibility in §6. The neutrino mass matrix is a 2N × 2N matrix of the form where M R is the symmetric complex Majorana mass matrix for the right-handed neutrinos, and M D is the complex Dirac mass matrix between the left-handed and right-handed neutrinos. For ||M R || ≫ ||M D || the heavy SM singlet neutrinos can be integrated out to obtain the N × N matrix for the active neutrinos, This is the matrix that will be the object of our study.
Before proceeding, it is necessary to point out that the overall scales for M D and M R will not be addressed in the paper. We only assume that ||M D || ≪ ||M R ||, so that the approximation Eq. we will consider an alternative strategy based on an ansatz for the eigenvalue distribution ρ(m ν ) for M ν , which will be verified numerically. Then we will apply the Coulomb gas picture to write down an action for the eigenvalues m i ν , that can be used to derive the correlation functions of the system. The first step is to diagonalize M D and M R separately, rewriting M ν as This is a good approximation away from the edges of the semicircle distributions.
We will find that the distribution of M ν can be described in terms of an action Eq. (2.15), with a confining potential V (m ν ) and a logarithmic Coulomb repulsion. Let us now consider the effects from eigenvalue repulsion, The matrix U rel rotates the relative orientation of M i and m j for the final eigenvalues. Although there are many possible relative orientations, we choose two points for illustration purpose. The first case is to have which has U rel = I N . The second example is with U rel as a permutation matrix that interchanges i ↔ N + 1 − i in M R . In other words, this case corresponds to having the largest suppression for the smallest M ν eigenvalue.
Comparing the repulsion Hamiltonian for these two choices of orientations, we have which shows that the case B has a smaller energy from contribution of the Coulomb repulsion. Concentrating then on the eigenvalue distribution of the case B (we drop the label "B" from now on), we have, for the smallest eigenvalue and the median, This is shown qualitatively in the spectrum in Fig. 2. Based on numerical results, we will now argue that the case B gives a very good approximation to the full eigenvalue distribution.
The spectrum of M ν is studied numerically as follows. Treating all elements of M D and M R as complex numbers, we choose random numbers for both the real and imaginary parts of the elements of M D within (−1, 1) and (−10 10 , 10 10 ) for M R , which incorporates the two hierarchic scales in the Seesaw mechanism. We calculate M ν and diagonalize it. This procedure is repeated several thousand times to calculate the eigenvalue distribution. Let us now focus on the smallest eigenvalue and the median, and compare with the predictions Eq. (4.9). The full distribution will be studied below.
The numerical results for the smallest eigenvalue as a function of N are shown in the blue circles of the left panel of Fig. 3, together with a fitted distribution using a function proportional to N −3 .
One can see that the N −3 is indeed a good description of the large N behavior and also of the N = 3 case. The right panel of Fig. 3 shows the ratio of the median eigenvalue over the largest eigenvalue.
(For simplicity, only odd N numbers are shown in this plot.) The numerical results are well fitted by a N −1 function. Therefore, we find a good agreement between Eq. (4.9) and the numerical eigenvalues; this provides nontrivial evidence that the distribution Eq. (4.7) is a good approximation. Further properties of the eigenvalue distribution will be analyzed shortly. Since we are doing statistical analysis for the expectation values of the eigenvalue ratios, it also important to compute the standard deviation of these observables as a function of N . Defining the large N dependence of σ m i /m j turns out to be the same as that of m i /m j . This is illustrated in Fig. 4, which shows the numerical ratios of the variance over the expectation value. This implies that our large N analysis can not capture detailed order one coefficients in front of the N dependence.
Comparing Fig. 1 and Fig. 2 shows that the Seesaw mechanism generates a larger hierarchical spectrum than the one obtained from random Majorana masses. For a fixed heaviest eigenvalue, random Majorana neutrinos have a lightest eigenvalue suppressed by 1/N , while in the Seesaw case the suppression is by 1/N 3 . The behavior of the median is also parametrically different as a function of N . These statistical properties are then very sensitive to the mass mechanism, providing predictions that distinguish the two scenarios. In §5 we will compare those two spectra to the experimentally observed values, and conclude that the Seesaw scenario is preferred. This effect may also have applications in other systems. For this reason, we now study in more detail the "Seesaw ensemble" by itself, independently of motivations from neutrino physics. As far as we know, this provides a new type of ensemble beyond the ones studied in the literature (see e.g. [11] for a recent review).

The Seesaw Ensemble
We define the Seesaw ensemble as the set of N × N matrices of the form This is in the same class as the unitary ensemble that appears in QCD, with no flavors and a nonchiral Dirac spectrum.
Ideally we would want to derive the probability distribution for M ν starting from those of M D and M R . This is, however, a very complicated task. Instead, we now explore the consequences of the ansatz Eq. (4.7) for the spectrum of eigenvalues of M ν , which we found to be in good agreement with numerical results. Our analysis of the Seesaw ensemble will rely heavily on the generalized Coulomb gas picture described in §2.2.
From Eqs. (3.7) and (4.7), it is not difficult to derive the continuous spectral density for the Seesaw ensemble, where the normalization is such that ∞ 0 dxρ(x) = N . For very small and large values of x, the density behaves as (4.14) Both limits reveal interesting properties of the ensemble. Note that the spectral density is divergent when x → 0, but ρ(x) is still integrable. Furthermore, ρ(x) does not have compact support -there is no upper bound on eigenvalues for the Seesaw ensemble. This comes from the inverse of M R in Eq. (4.11), so that it is always possible to have eigenvalues much larger than the light scale m 2 /M . These properties are in sharp contrast with the GUE or AZ-CI ensembles that we discussed before, for which ρ(x) has compact support and is finite as x → 0.
In order to verify our analysis so far, let us compare the spectral density Eq. (4.13) with numerical results. This is done in the left panel of Fig. 5, which shows both the numerically averaged density of eigenvalues for N = 200 and the analytic prediction. One can see a good agreement between those two distributions.
Let us now develop the Coulomb gas model for the Seesaw ensemble. As we reviewed in §2, this representation is given by a system of particles moving along a line in two dimensions, subject to a Coulomb repulsion and a confining potential V (x). We already discussed the effects from repulsion, and now we will determine V (x). Plugging the eigenvalue density into Eq. (2.16) obtains Mapping back to the matrix eigenvalues, this means that the probability distribution for Seesaw matrices is In the limits of small and large values of x, we have Thus, the confining potential V N Seesaw (x) is continuous at x = 0, but its first derivative diverges there. This is much steeper than the harmonic potential for the GUE, explaining why the eigenvalue density is parametrically larger near the origin. 3 The Coulomb gas picture hence provides a physical way of understanding the steeply decreasing behavior of ρ(x) as x increases.
Although we will not study higher correlators in detail, let us mention that given V (x), the correlation functions for the Seesaw ensemble can be obtained as functions of polynomials that are orthogonal with respect to e −V (x) . They are of the form Eq. (2.6) in terms of these new orthogonal polynomials.
Having analyzed the statistical properties of neutrinos for general N , let us now focus more specifically on neutrino physics, comparing our results (extrapolated to N = 3) with the experimental data.
From all neutrino experimental results including solar, atmospheric, reactor and accelerator experiments, the global analysis [21] obtains the following numbers for the neutrino mass squared differences and mixing angles, defined in the PDG [22]: Here, ∆m 2 ij ≡ m 2 i − m 2 j ; "N" means the normal ordering for the neutrino spectrum such that ∆m 2 21 ≪ (∆m 2 32 ≈ ∆m 2 31 > 0); "I" means the inverted ordering for the neutrino spectrum such that ∆m 2 21 ≪ −(∆m 2 31 ≈ ∆m 2 32 < 0). We take their fitted results for treating the reactor fluxes as free parameters and including the results for the reactor experiments with L 100 m. As emphasized in [21], although a zero Dirac CP violation phase δ CP is referred, the significance is less than 1.5σ for the normal ordering and 1.75σ for the inverted ordering. The absolute scale of neutrino masses is currently unknown and bounded from above as m ν ≤ 0.28 eV at the 95% confidence level [23], assuming a flat ΛCDM cosmology.
We should emphasize that our analysis does not explain the overall scale of neutrino masses, but only ratios of masses. Let us concentrate first on the ratio ∆m 2 21 /∆m 2 31 ≈ 0.03. For complex Majorana masses without the Seesaw mechanism (see Fig. 1), the prediction is ∆m 2 21 /∆m 2 31 = O(1) at large N , which does not provide a good fit to the experimental value. On the other hand, the prediction for Seesaw masses (see Fig. 2) is ∆m 2 21 /∆m 2 31 = O(1/N 2 ) = O(0.1) for N = 3, much closer to the experimental result. (Recall that these predictions have order one uncertainties.) Therefore, assuming random masses, the Seesaw mechanism is favored by data.
Focusing then on the Seesaw mechanism, let us make the following observations based on the spectrum in Fig. 2. First, the neutrino spectrum is a normal hierarchical spectrum. The general neutrino mass scale is determined by the heaviest mass, m 3 ≈ ∆m 2 32 ≈ 0.05 eV, which could be tested from the Cosmic Microwave Background by the Planck satellite or the KATRIN experiment [24] in the future. The prediction for the effective Majorana mass | m |, which may be measured from neutrinoless double beta decay experiments, is in the region of (2.3×10 −4 , 5.0×10 −3 ) eV for the normal hierarchical spectrum [22]. Secondly, the lightest neutrino mass is predicted to be m 3 /N 3 = m 3 /3 3 ≈ 0.002 eV, which is unlikely to be tested in the near future unless one can come up a new experimental method to measure the individual neutrino mass more precisely.
Finally, let us discuss the basic properties of the mixing matrix for this case. As mentioned in [5], the basis-independence assumption predicts that the PMNS matrix follows the Haar measure of the Lie group. While in principle the corresponding measure for the Seesaw ensemble can be obtained starting from the Haar measure distributions for M R and M D , for our purpose it is enough to evaluate numerically the expectation values of the elements of the matrix U that diagonalizes M † ν M ν . The left panel of Fig. 6 shows the behavior of the absolute value of the off-diagonal entry U 1N as a function of N , which is proportional to 1/ √ N . This is just what we expect based on the normalization of eigenvectors. For the phases in the rotation matrix, we use the generalized Jarlskog invariants to study the rephasing invariant CP violation [25,26]. In the right panel of Fig. 6, we show the distribution of N 2 |Im[U 11 U * 23 U * 13 U 21 ]| as a function of N . The overall factor N 2 cancels the N dependence of the absolute values of the elements. We can see from this panel that the phase part of this Jarlskog invariant is order 1 and independent of N .
Translating these results into mixing angles and phases, our prediction for sin θ 13 is order 1/N = 1/3, not far from the measured values sin θ 13 ≈ 0.15. For the CP violation phases, we anticipate all the phases in the leptonic mixing matrix in the weak charged current interactions to be order 1. The predictions of our analysis are summarized in Table 1.

Conclusions and future directions
In this work we proposed to treat the number N = 3 of Standard Model generations as a large number, studying the consequences of this idea for neutrinos. Assuming anarchic neutrino masses, this becomes a problem in Random Matrix Theory. We analyzed the neutrino spectrum and mixing matrices analytically and numerically, and compared with experimental results. The statistical properties of these observables were found to be highly sensitive on the underlying mechanism for neutrino masses.
In particular, the Seesaw mechanism leads to a 1/N 3 hierarchy for the lightest neutrino, and has a spectrum that is qualitatively compatible with experimental data.
Already at this level, the connections between neutrino physics and RMT appear to be interesting.
Based on the Seesaw mechanism, we defined a new type of random matrix ensemble, the "Seesaw ensemble," which, as far as we know, has not been analyzed before in the literature. The properties of this ensemble are qualitatively different from those of Gaussian distributed matrices. We proposed an ansatz for the spectrum that was in good agreement with the numerical results. In the Coulomb gas picture, the confining potential for the Seesaw matrices is much steeper than in the Gaussian case, with a diverging first derivative at the origin.
It will be interesting to determine the correlation functions of the Seesaw ensemble starting from the distributions of the underlying matrices M D and M R in Eq. (4.11). We hope that the analytic and numerical results reported here can shed light on this problem. There are also various possible extensions. First, one could consider directly the ensemble of matrices of the form Eq. (4.1), without restricting from the start to the limit ||M R || ≫ ||M D ||. Furthermore, there is also a natural twoparameter generalization of this ensemble, where M R is a square N R × N R matrix, but M D is an N L × N R matrix. Physically, this corresponds to having N L active and N R singlet neutrinos [27]. One could then study the statistical properties at fixed N L and large N R . More generally, the existence of the new parameter N R /N L can lead to different large N behaviors or multiscale effects. This is reminiscent of QCD with N f fundamental fermions, where taking N f /N c ∼ 1 led to the Veneziano or topological limit [28].
The large N limit can also be applied to neutrino physics with one or more sterile neutrinos or the SM quarks and charged leptons. The dramatically hierarchical pattern of quark masses and angles implies that the couplings in the Lagrangian cannot be randomly distributed -more structure is required to explain the flavor problem. Nevertheless, extra-dimension models such as [29] still allow for random coefficients, and it would be interesting to understand how the dynamics in the quark sector is modified by large N statistical effects.
Finally, given the incredible range of applications of RMT, the Seesaw ensemble may occur in other types of systems. One would need some pattern of the form Eq. (4.11), arising from the interplay of two hierarchical scales.