Anarchy and neutrino physics

The neutrino sector of a seesaw-extended Standard Model is investigated under the anarchy hypothesis. The previously derived probability density functions for neutrino masses and mixings, which characterize the type I-III seesaw ensemble of N × N complex random matrices, are used to extract information on the relevant physical parameters. For N = 2 and N = 3, the distributions of the light neutrino masses, as well as the mixing angles and phases, are obtained using numerical integration methods. A systematic comparison with the much simpler type II seesaw ensemble is also performed to point out the fundamental differences between the two ensembles. It is found that the type I-III seesaw ensemble is better suited to accommodate experimental data. Moreover, the results indicate a strong preference for the mass splitting associated to normal hierarchy. However, since all permutations of the singular values are found to be equally probable for a particular mass splitting, predictions regarding the hierarchy of the mass spectrum remains out of reach in the framework of anarchy.


Introduction
The neutrino sector of the Standard Model (SM) is quite peculiar. Indeed, although the quark and charged lepton mass spectra are quite hierarchical, the neutrino spectrum is simple: all neutrinos are massless. Neutrino oscillations [1][2][3], where neutrinos seemingly change flavor in flight, cannot be accommodated in the SM due to the masslessness of the neutrinos. Neutrino oscillations thus imply massive neutrino eigenstates and the SM must be extended. Moreover, neutrino oscillation experimental data suggest that the neutrino spectrum is not hierarchical, with three massive light neutrinos and a mixing matrix, the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix, exhibiting near-maximal mixing.
The anarchy principle, introduced in [4,5], was put forward to explain the peculiarities of the neutrino sector by postulating a low-energy neutrino mass matrix generated by one of the seesaw mechanisms from randomly-generated high-energy mass matrices with elements distributed following a Gaussian ensemble. In a series of papers [6][7][8][9][10], several numerical analysis of the anarchy principle were performed by randomly generating high-energy mass matrices and computing the corresponding low-energy neutrino mass matrices.
Recently, the low-energy neutrino mass matrix probability density function (pdf) for the anarchy principle was obtained from first principles in [11] following the extensive literature on random matrix theory [12][13][14]. It was shown that the pdfs associated to type I and type III seesaw mechanisms were given by the same complicated integral equation while the pdfs associated to type II seesaw mechanism was simple. A partial investigation of these seesaw ensembles was completed in [11] but an analysis of the physical case of three light neutrinos was not undertaken. This paper closes the gap by studying the implications of the seesaw ensembles for neutrino physics.

JHEP04(2017)131
The low-energy neutrino mass matrix pdf is known to factorize into a singular value pdf and a group variable pdf for all seesaw mechanisms. The singular value pdf corresponds to the pdf for the light neutrino masses while the group variable pdf corresponds to the pdf for the mixing angles and phases of the PMNS matrix. The singular value pdf for type I-III is given in [11] as a multidimensional integral while the singular value pdf for type II is a simple Gaussian-like distribution. The group variable pdf for all types of seesaw mechanisms is the Haar measure, which seems to prefer near-maximal mixing. However, as stressed in [15], the mode of a pdf is not a well-defined quantity (e.g. it is not invariant under change of variables), hence it is preferable to compare pdfs by comparing their probabilities associated to a particular outcome. This probability test, which is hard to perform by randomly generating low-energy neutrino mass matrices, can however be straightforwardly implemented from the analytic pdfs.
Moreover, the factorization of the pdfs for the singular values and the group variables implies that there is no link between the light neutrino masses and the light neutrino mass eigenstates, forbidding an investigation of the preferred mass hierarchy (normal or inverted). From the analytic results for the singular value pdfs, it is however possible to determine which mass splitting, i.e. which ordering ofm 2 med −m 2 min andm 2 max −m 2 med , is favored. This paper is organized as follows: section 2 gives the relevant pdfs for the seesaw ensembles obtained in [11]. In section 3 the pdfs for the complex seesaw ensembles are studied in the N = 2 and N = 3 cases. For both cases, a comparison is made between the analytic results and the numerical results, showing perfect agreement. For the N = 3 case relevant to neutrino physics, a thorough investigation of the implications of the seesaw ensembles is completed. For example, it is shown from the probability test that both type I-III and type II seesaw ensembles prefer the mass splitting associated to normal hierarchy with a neutrino energy scale of O(10 −2 ) eV. Finally, a discussion and a conclusion are presented in section 4.
It is important to note that throughout this paper, the term "analytic results" should be understood as results obtained from the analytic pdfs which are numerically integrated while the term "numerical results" corresponds to results obtained from randomlygenerated mass matrices.

The seesaw ensembles
This section states without proof the relevant quantities that appear in the seesaw ensembles, which were derived from the anarchy principle applied to the SM extended with the type I-III seesaw mechanism or with the type II seesaw mechanism. The reader interested in the proofs is referred to [11].

Probability density functions
The pdfs for the dimensionless N × N light neutrino mass matrixM ν = M ν /( √ 2Λ ν ) where Λ ν is the (naturally small) light neutrino mass scale, can be expressed in terms JHEP04(2017)131 of the light neutrino mass matrix singular valuesm ν,i and the light neutrino mass matrix group variables U ν with the help of the decompositionM ν = U ν D ν U T ν where D ν = diag(m ν,1 , · · · ,m ν,N ) andm ν,i ≥ 0 for all i. The pdfs are found to factorize into two independent pdfs, one pdf for the singular values and one pdf for the group variables, as in For real (β = 1) and complex (β = 2) matrix elements, the pdfs are given respectively by .

(2.2)
Note that the pdf for the group variables is the normalized Haar measure for all types of seesaw mechanisms. Here the function I β N relevant for the type I-III seesaw mechanism is where the integration is over the full Stiefel manifold V β N ≡ V β N,N if β = 1 but only parts of the full Stiefel manifold if β = 2 (hence the prime, see [11]). The normalization constants and the volume of the Stiefel manifold V β N are .
An important feature of the joint pdfs for the singular values (2.2) is their invariance under permutations of the singular values. Apart from the function I β N , the pdfs are clearly invariant under such a transformation. The function I β N is also invariant under permutations since a permutation only reshuffles the columns and rows of the matrix U , which is integrated over. 1 Hence, the neutrino masses can be reshuffled freely amongst

JHEP04(2017)131
themselves. This invariance, which leads to a complete independence between the light neutrino mass eigenstates and the light neutrino masses, implies that the probability for a specific spectrum of masses is at most 1/N !. This observation has far-reaching consequences in the physical case of neutrino physics.
It is important to note that, at the level of the pdfs, the only difference between the type I-III and the type II seesaw mechanisms lies in the singular value pdfs. Hence, when comparing which ensemble better generates the observed neutrino parameters, only the information on the mass splittings will differ.

A parametrization for unitary matrices
Since the goal of this paper is to compare the implications of the seesaw ensembles with actual neutrino observations, the rest of the paper focuses on complex matrix elements, i.e. β = 2. The analysis of the case with real matrix elements is similar.
With that in mind, it is important to find a convenient parametrization for unitary matrices to proceed with the analysis. Indeed, to write the integral (2.3) more explicitly when β = 2, it is necessary to assume a parametrization for the unitary matrix U . Moreover, the light neutrino group variables pdf is the Haar measure for unitary matrices.
A convenient parametrization, based on [16,17], implies that a N × N unitary matrix U can be written as Here the matrices P j and Σ jk are given by 5) and the N 2 mixing angles θ jk and phases φ jk and ϕ j belong to the following intervals, Moreover, the Haar measure depends only on the mixing angles. For the light neutrino masses, the parametrization (2.4), which gives shows clearly that all dϕ i are not integrated over in (2.3). Indeed, the most interesting quantity in (2.3) is (U † dU ) and it is given by

9)
JHEP04(2017)131 using the wedge product notation [11]. Hence the variables ϕ i are not even part of the integral, implying a N 2 -dimensional integral I β=2 N in the singular value pdf. For the light neutrino mixing matrix, the quantity of interest is simply the Haar measure (2.7). In that case, the corresponding variables ϕ i are the unphysical phases that can be absorbed by a field redefinition. The remaining phases φ ij are the CP-violating Dirac and Majorana phases. All those phases have flat distributions that are uninteresting. The mixing angles on the other hand, have non-trivial distributions. This fact has important consequences for the SM neutrino physics since near-maximal mixings seem highly probable.

The complex seesaw ensembles
In this section the complex seesaw ensemble pdfs (2.2) are analyzed for small values of N . The case N = 1 was studied analytically in [11] for both real and complex matrix elements. Here the case N = 2 with complex matrix elements is investigated and compared to numerical results for 2 ×2 matrices. Then the physical case of N = 3 is analyzed to determine how likely the pdfs are to generate the observed light neutrino masses and mixings.

Consequences of the complex seesaw ensembles
Before discussing specific values of N , it is enlightening to state the implications of the complex seesaw ensemble pdfs (2.2) in general terms.
First, once a decomposition for the light neutrino mass matrixM ν = U ν D ν U T ν is chosen, the implications of near-maximal mixings obtained from the group variable pdf (2.2), i.e the appropriate normalized Haar measure, for some mixing angles seem inescapable. Indeed, using the parameterization (2.4), the Haar measure, given by (2.7), dictates that the most probable value for the mixing angles θ ij is arccot[ 2(j − i) − 1] while all the remaining (unphysical, CP-violating Dirac and Majorana) phases have flat distributions. Thus the preferred value for all mixing angles θ i,i+1 is π/4, which corresponds to maximal mixing, while the preferred value for all mixing angles θ i,i+2 is π/6. It is important however to notice that most probable values are not invariant under change of variables, as pointed out in [15].
Then, the consequences for the light neutrino masses, which are obtained from the singular value pdfs (2.2), are not as sharp. Indeed, although the chosen decomposition is here fixed, the singular value pdfs (2.2) are invariant under permutations of the singular values. Hence, the light neutrino masses and the light neutrino mixing angles and phases are completely independent. In other words, although each singular valuem ν,i has a corresponding singular vector (u ν,i ) j = U ν,ji such thatM ν u * ν,i =m ν,i u ν,i , the probability that a particular neutrino mass spectrum occurs is at most 1/N !. For example, the dimensionless neutrino mass spectrum (m ν,1 , · · · ,m ν,N ) = (µ 1 , · · · , µ N ) is as probable as the spectrum (m ν,1 , · · · ,m ν,N ) = (µ 2 , µ 1 , µ 3 , · · · , µ N ) or any other permutations. It is thus possible to fix an ordering for the singular values, 0 ≤m min ≤ · · · ≤m max , keeping in mind that the relationship between the light neutrino masses and the light neutrino mixing matrix is completely lost. Fixing the ordering implies that the singular value pdfs (2. Therefore, by working with a fixed basis as described above, some mixing angle preferred values correspond to maximal mixing but the ordering of the light neutrino masses for a given spectrum is completely free. It is thus clear that a spectrum exhibiting one of the two hierarchy patterns preferred by the data (normal or inverse) is as probable as the same spectrum but with permuted mass eigenstates. These observations are in the spirit of [15], although with the analytic knowledge of the singular value pdfs (2.2), it is now possible to complete an appropriate statistical test to better check the validity of the anarchy principle. As stressed in [15], the most appropriate statistical test seems to be the probability test which computes the probability that the variables are in a given volume. By choosing the observed values with their error bars for the volume of the light neutrino masses and mixings, the probability that one ensemble leads to the SM is obtained. Clearly, since the variables are continuous, the calculated probability is very small for very precisely-known observed values. One can nevertheless discriminate ensembles, for example the type I-III and type II seesaw ensembles, by comparing their respective probabilities, as shown below.

The case N = 2
As a warm-up exercise, the case N = 2 is studied analytically and compared to randomlygenerated light neutrino mass matrices. Using the parametrization ( Table 1. Location parameters for the marginal singular value pdfs of figure 1.

pdfs (2.2) can be written as
where the subscript ν on the masses was omitted to simplify the equations. The modified Bessel function of the first kind I 0 (z) is generated by the integral over the phase in (2.7). The 4-dimensional integral is thus simplified to a 3-dimensional integral. Figures 1 and 2 show a comparison between the analytic results (3.1) and numerical results for a fixed ordering of the singular values (chosen to be 0 ≤m 1 ≤m 2 , such thatm 1 ≡m min and m 2 ≡m max ). Consequently, the resulting marginal singular value pdfs are obtained by computing the following integrals for both ensembles (i.e. Σ = I-III or II).
Although the analytic behavior of the type I-III singular value pdf (3.1) is hard to see intuitively due to the integral, it is clear that the pdfs (3.1) are correct as seen in figure 1. The behavior of the type I-III singular value pdf at vanishing and large singular values matches the expectations of [11]. The vanishing ofP I-III ν (m 1 ) atm 1 → 0 is not apparent in the histogram due to the bins being too large. The type II singular value pdf is much easier to study. The pdfs for the smallest and the largest singular values can be obtained analytically from (3.1). Again, there is a good match between the analytic results, which are simple exponentials, and the numerical results. In order to make an appropriate comparison between the two ensembles, it becomes essential to fully characterize the previous pdfs. Thus, one is expected to compute their corresponding moments as well as their modes and medians. However, as stated in [11], the only existing moment for the type I-III complex JHEP04(2017)131 Figure 2. Probability density functions for the mixing angle and phases of the complex seesaw ensembles with N = 2. The red curve corresponds to the analytic result while the histogram corresponds to numerical results (with 2.5 × 10 4 dimensionless light neutrino mass matrices generated). The top and bottom rows show the pdfs for the mixing angle θ, the CP-violating phase δ, and the unphysical phases ϕ 1 and ϕ 2 (note the range is halved due to the extra freedom ϕ i → ϕ i ± π [11]).
(β = 2) seesaw ensemble is the first moment (the average singular values). Therefore, a meaningful characterization of the previous pdfs is limited to the first moment, the mode and the median (their location parameters). Their respective values for each distribution are presented in table 1.
From these results, it can be seen that the average singular values coming from the type I-III seesaw ensemble are spread over a wider range than in the type II seesaw ensemble. Moreover, when comparing the mean of a distribution with its respective median, one can quantify the asymmetry of the pdfs presented in figure 1. It turns out that the means are much closer to the medians (and thus the modes) in the type II seesaw ensemble, which leads to more symmetrical pdfs as can already be seen from figure 1.
Moving forward, the group variable pdf (3.1) is easier to analyze. First, all phases have flat distributions as mentioned previously. Moreover, the mixing angle has a non-trivial distribution that prefers near-maximal mixing. From figure 2 the pdf for the mixing angle and the phases agree well with the normalized Haar measure. Since these pdfs were studied extensively in the literature and are easier to analyze, their statistical parameters (mean, median and mode), which are easily obtained from (3.1), are not presented here.
Finally, to provide a global overview of the pdfs with unordered singular values, the density plots of the singular value pdfs for the type I-III and type II seesaw ensembles are shown in figure 3. The symmetry pattern of the pdfs under the exchangem 1 ↔m 2 is easily seen from figure 3. In fact, the plots are constructed in a way that takes advantage of this symmetry to provide a meaningful comparison between analytical and numerical results in both cases. Indeed, by showing only half the data for the analytical (m 1 <m 2 ) and numerical (m 1 >m 2 ) results in each plot, it can be seen that the agreement between the two is once again very good. Moreover, by comparing the distances between modes in each plot, it is possible to determine which ensemble has a stronger repulsion between the singular values. It is found that the modes are 1.7 times farther apart in the type I-III seesaw ensemble, meaning that a stronger repulsion is observed between the singular values for this ensemble. Unfortunately, the origin of this interesting feature is hard to trace back without a completely-integrated analytical expression for P I-III ν (m 1 ,m 2 ). A possible explanation as to why the singular values are closer together in the type II seesaw ensemble is suggested in expression (3.1). Indeed, one can see that the Vandermonde-like contribution |m 2 1 −m 2 2 | in P I-III ν (m 1 ,m 2 ) (responsible for the spreading of the singular values with reference to the symmetry axis), is strongly suppressed by a product of masses to the fifth power, which tends to spread the singular values in a narrow band along the two axes. Eventually, in order to get a better understanding of the type I-III seesaw ensemble, the density plot could be used to guess an analytical form for P I-III ν (m 1 ,m 2 ) by fitting some appropriate functions of the singular values with free parameters to be determined. Such density plots are not really conceivable for the case N = 3 as they would require heavy numerical computation and would be rather hard to illustrate properly. However, the same reasoning and conclusions apply to this case as well.

The case N = 3: SM neutrino physics
With the tools and insights developed in the previous sections, it is now possible to fully analyze the more interesting case of a seesaw-extended SM.

JHEP04(2017)131
Normal Hierarchy Inverted Hierarchy First, recent experimental values of the physical parameters in the neutrino sector are summarized in table 2. The mixing angles, the CP-violating Dirac phase and the squared mass differences ∆m 2 ij = m 2 i − m 2 j are extracted from [2]. 2 The upper bound on the mostly-electronic neutrino m 1 is the 1σ upper bound of [18]. This upper bound on m 1 comes from the study of supernova. It is conservative and quite model-independent. Indeed, it is the weakest bound when compared to the cosmological bound on the sum of the neutrino masses which is somewhat model-dependent or the neutrinoless double βdecay bound and the direct neutrino mass bound which are intertwined with some mixing matrix parameters [19]. These values will be used in the probability test at the end of this section.
For neutrino physics, the most convenient parametrization for the unitary group U(3) is of course the PMNS mixing matrix [1] for the light neutrino U ν ,  phase δ and the two CP-violating Majorana phases α 21 and α 32 are flat. Therefore, any value for the CP-violating phases is equally probable in the complex seesaw ensemble. It is important to note that the unphysical phases are not explicitly included in the PMNS parametrization (3.4).
The full expression for the type I-III singular value pdf is quite long (the explicit expression fills up a few pages) and not enlightening. The form (2.2) with the parametrization (2.4) and N = 3 is sufficient for both type I-III and type II. As discussed in section 3.1, these pdfs are invariant under permutations of the singular valuesm 1 ,m 2 andm 3 . Consequently, the hierarchy of the neutrino mass spectrum of the extended SM cannot be predicted under the anarchy hypothesis. The only claim that can be made is that all hierarchy scenarios (that is to say, every 3! permutations of the three singular values) are equiprobable for a given mass splitting (again, the term mass splitting is understood to be a particular ordering of the two quantitiesm 2 med −m 2 min andm 2 max −m 2 med ). To further study the mass splitting and the resulting marginal pdfs, it is convenient to introduce a particular singular value ordering. From now on, the ordering is chosen to be 0 ≤m 1 ≤m 2 ≤m 3 (such thatm 1 ≡m min ,m 2 ≡m med andm 3 ≡m max ) although it must remain clear that those are not straightforwardly related to the experimental neutrino masses (once again, the choice is completely arbitrary). Thus, one cannot single out any of the two hierarchy scenarios selected by experimental data. Nevertheless, by studying the marginal singular value pdfs of the type I-III and type II seesaw ensembles, some important and interesting results on the regime of low-energy neutrino physics can be obtained.

JHEP04(2017)131
First, the marginal singular value pdfs can be obtained by computing the following integrals,P  figure 4. A first observation that comes to mind when looking at figure 4 is the diversity of the mass spectrum obtained with the type I-III seesaw ensemble compared to the type II seesaw ensemble. This can be traced back to the fact that the pdfsP I-III ν (m i ) are much more complex than the simple Gaussian-like pdfsP II ν (m i ) that arise in the type II seesaw ensemble. For example, there is no internal angular dependence associated to extra variables that needs to be integrated out in the type II pdf. A second observation worth mentioning is the remarkable agreement between analytical and numerical results. For reasons previously mentioned, analytical results derived from the type I-III seesaw ensemble are much more challenging to get than those of the type II seesaw ensemble. Following heavy numerical computation based on adaptive Monte Carlo integration, the 11-dimensional integrals resulting from the marginalization procedure can be obtained for given values ofm 1 ,m 2 orm 3 . The red curves produced in such a way are in very good agreement with their corresponding histograms (generated from a sample of light neutrino mass matrices), which can be viewed as a validation of the Monte Carlo integration method (estimated to be accurate to at least three significant figures) used in this case or a numerical check of the singular value pdfs obtained in [11].
Once the numerical integration method is carefully tested, the next step is to compute the relevant statistical parameters. The results are presented in table 3. Following the same approach as in the N = 2 case, it is found that the average singular values are once again spread over a much wider range in the type I-III seesaw ensemble. Moreover, when compared to the N = 2 case, it can be seen that this range expands significantly more in the type I-III seesaw ensemble as N increases. Next, comparing these values with their respective medians, one can conclude that the pdfs are much more symmetric in the type II seesaw ensemble, as can be expected when looking at figure 4.
Even though determining the hierarchy of the mass spectrum is out of reach in the context of the seesaw ensembles, the previous results can still be used to help identify which of the two possible mass splittings (according to our previous ordering, the mass splittings can be written as ∆m 2 21 =m 2 med −m 2 min and ∆m 2 32 =m 2 max −m 2 med so that the two possibility are either ∆m 2 21 < ∆m 2 32 or ∆m 2 21 > ∆m 2 32 ) is more likely to occur under the anarchy hypothesis.
By studying the distribution of the ratio R, (∆m 2 21 > ∆m 2 32 ) is realized. For the type I-III seesaw ensemble, the probability is 95.8% (4.2%) whereas for the type II seesaw ensemble the probability is 79.0% (21.0%). These probabilities are supported by figure 5. One can thus conclude that the dominant trend for both ensembles is the realization of the mass splitting ∆m 2 21 < ∆m 2 32 reminiscent of the normal hierarchy. Moreover, looking at the very distinct behavior of the two pdfs and comparing the resulting probability, it is possible to state that the type I-III seesaw ensemble is better suited to generate this particular mass splitting. In the context of the type I-III seesaw ensemble, this means that the mass differences are way more likely to coincide with the ones from normal hierarchy, yet the ordering of the masses is still unknown (once again, every 3! permutations are equally probable for this particular splitting).

JHEP04(2017)131
Next, considering the group variable pdf for the mixing angles and phases (3.5), similar conclusions as with the case N = 2 can be drawn. Once again, the phases have flat distributions and two of the three mixing angles prefer near-maximal values as shown in figure 6. Here, the unphysical phases were not considered in the making of figure 6 as they were deemed not interesting for the present discussion. The numerical data coming from a sample of light neutrino mass matrices is once again consistent with the marginal pdfs obtained from the Haar measure. The only non-symmetric pdf for the mixing angles is the one associated to θ 13 . Its mode, located at π/6, is in agreement with the results of section 3.1 obtained with the parametrization (2.4) since the Haar measure is the same for the PMNS matrix (3.4).
To further emphasize the differences between the two ensembles, the probability test is now used to determine how well the complex seesaw ensembles can generate the observed values of the extended SM physical parameters in the neutrino sector. The results of the probability test will then be compared between the type I-III and type II seesaw ensembles to determine which ensemble is more likely to generate the observed values of physical parameters. The experimental values are given in table 2.
With the observed values of table 2, the probabilities (where |detJ| is the Jacobian of the appropriate hierarchy) are P NH m (Λ ν ) = where NH stands for normal hierarchy while IH stands for inverted hierarchy. The first test is achieved by using only the singular value pdfs for both ensembles (see figure 7). For the type I-III seesaw ensemble, the 12-dimensional integrals over the experimental volume defined by the 1σ range are obtained using the same Monte Carlo algorithm. At this point, it is necessary to stress that this test does not require any ordering of the singular values. In fact, all permutations are accounted for in these integrals and there is thus no need for an extra factor of 3! to ensure the normalisation of the pdfs. Since the only free parameter left in the equations is Λ ν , the idea is to plot the probability as a function of Λ ν over a range where the curves reach a maximum. This allows for a simple comparison of their maximum values (by taking appropriate ratios) to determine the likelihood of each ensemble to generate the observed values. It is important to note that beside the fact that the probabilities obtained this way are invariant under a change of basis, the explicit values of the probabilities are not particularly meaningful. In fact, they are bound to shrink further and further as the experimental values get more and more precise. However, the ratios are considered to be relevant quantities since they are subject to only small fluctuations during this process (the order of magnitude should remain the same). The results of the probability test are presented in figure 7.
First, from within the same ensemble, one can compare the probabilities obtained from the normal and inverted hierarchy scenarios. The maximum probability values resulting from a scan over Λ ν reveal that the mass splitting in (3.7) are ∼ 1000 times more likely to originate from the region defined by normal hierarchy (at 1σ) than from the one defined by inverted hierarchy in the type I-III seesaw ensemble. In other words, this means that the type I-III seesaw ensemble is way more likely to generate values for these physical parameters that are contained within the region allowed by the normal hierarchy data set (rather than the inverted hierarchy data set). For the type II seesaw ensemble, the same tendency is observed but with a much smaller ratio between the maximum probability values. Indeed, figure 7 shows that the mass splitting is ∼ 25 times more likely to originate from the region defined by normal hierarchy. It is then possible to conclude that between the two regions scanned in the probability test, both ensembles naturally lead to preferred values for these physical parameters that lie in the region defined by normal hierarchy. Moreover, this preference is strongly accentuated in the type I-III seesaw ensemble.
Second, one can make a comparison between the two ensembles based on the result of figure 7. Since it was shown that one region is actively preferred over the other, it becomes useful to compare the maximum probability values in the case of normal hierarchy for JHEP04(2017)131 both ensembles. This time, the conclusions are not as striking as in the previous case but one can state that the type I-III seesaw ensemble is roughly 2 times better than type II for generating values of these parameters in this particular region. 3 The results obtained from this probability test are thus in agreement with what was found previously by comparing the pdfs of the ratios R and consequently help quantify the underlying trends in both ensembles.
A final point of interest regarding this particular test concerns the energy scale Λ ν . Again from figure 7 one can see that choosing the integration region to be over the accepted experimental values (within the 1σ confidence level) naturally fixes the energy scale of the models. Each scan shows that the maximum probability values are attained for values of Λ ν which are of the same order of magnitude, namely Λ ν ∼ O(10 −2 )eV for both ensembles. Since Λ ν , which corresponds to the light neutrino mass scale, takes the general form Λ ν = v 2 /Λ new for each type of seesaw mechanisms, with v 246 GeV the usual Higgs vacuum expectation value, a quick estimate of the new energy scale Λ new associated to the particle content introduced in the extended SM with type I, type II or type III seesaw mechanism (right-handed neutrinos singlets, Higgs triplet and fermionic triplets respectively) can be made. Indeed, using the previously-mentioned values, one gets Λ new ∼ O(10 15 ) GeV for the new energy scale of the extended SM, which is very close to the energy scale of grand unified theory (GUT). Naturally, Λ new is directly related to the masses of these newlyintroduced particles. However, in order to assess their corresponding mass scales, one needs to specify the order of magnitude of the coupling constants arising from each seesaw scenario. The usual approach is to set the coupling constants to be of O (1) since there is no fundamental principle or symmetry pattern that require particularly small couplings. This in turn suggests that the new particles introduced in the SM are quite heavy since Λ new becomes essentially their corresponding mass scale. In fact, this result is typical of seesaw mechanisms and is often regarded as a prerequisite (when taking the naturalness argument into consideration to avoid seesaw-induced fine-tuning or hierarchy problems) for these mechanisms to give sensible predictions concerning the light neutrino masses. It is however possible to lower the mass scales by simply postulating smaller coupling constants, somewhat disregarding the naturalness argument. Overall, the results of the probability test are therefore consistent with high-energy phenomenology of the seesaw-extended SM.
The second probability test, with results shown in table 4, concerns the mixing parameters of the neutrino sector, namely the mixing angles and phases of table 2. In this case, the analysis is much simpler since there is no free parameter with which to scan a particular region and the pdf (the Haar measure) is also a lot less complicated. Since both ensembles have the same pdf for the mixing angles and phases, a comparison between the two is not possible. However, it is interesting to see how well these ensembles perform when compared to a trivial normalized flat distribution. Here, the idea is simply to test whether there is any improvement when generating parameter values from the Haar measure obtained in the seesaw ensembles as opposed to a less interesting model where there would be 3 It is interesting to note that the type II seesaw ensemble is approximately 18 times more probable than the type I-III seesaw ensemble for the inverted hierarchy.  Table 4. Probability test for the mixing angles and phases of the complex seesaw ensembles with N = 3 and a (trivial) normalized flat distribution. The probability P U and P flat are obtained by integrating the normalized Haar measure and the flat distribution over the experimental volume V exp defined by the data at the 1σ confidence level (see table 2).

JHEP04(2017)131
no information or explicit dependence on the angular part in the pdf. By comparing the probabilities that the generated values lie within V exp in both cases and for the two types of hierarchy, one gets the results of table 4. The first conclusion that can be drawn from these results is that, up to the level of accuracy acknowledged for this test (remember that this analysis remains sensible to the choice of integration volume to some extent), there can be no distinction between normal or inverted hierarchy. Both regions are thus equally probable. However, when comparing P U with P flat , there is indeed improvement as the seesaw ensembles are essentially ∼ 2 times more likely to generate parameters within the allowed experimental region (for both hierarchy scenarios) than the flat distribution.

Large N comparison
To follow up on our previous work [11] regarding the close resemblance of the type I-III singular value pdf at N = 1 and the level density at large N , this section further investigates this connection by adding the comparison with the pdfs at N = 2 and N = 3 for the type I-III seesaw ensemble. The starting point for an appropriate comparison of these quantities is the correlation function 8) for N = 2 and N = 3 respectively. Introducing a convenient rescaling of the variable x → √ Nm ν , the resulting correlation functionŝ withρ(x) = ρ(x)/ √ N , can be compared directly to the large N histogram (N = 60). The resulting 11-dimensional integral for N = 3 is carried out using the previously-mentioned Monte Carlo algorithm.
From figure 8, one can see that the agreement between analytical and numerical results becomes surprisingly good as N reaches 3. However, there is a priori no clue as to why the correlation functions for finite and small N are able to reproduce with great precision the level density at large N since they are independent quantities. This represents the first clear indication that a proper large N analysis would indeed be a good approximation of the physical case N = 3, as was previously suggested (without proof or evidence) in the literature. Thus, there is no doubt that this particular behavior should be investigated further since it motivates the search for an analytical expression (coming from a large N analysis) to better understand the physical case at hand.

Discussion and conclusion
In this work the statistical implications of the seesaw ensembles, following the anarchy principle, for the physical case of three neutrinos were obtained. It is shown that the analytic pdfs computed in [11] are in perfect agreement with the numerical results of randomly-generated light neutrino mass matrices for the complex seesaw ensembles with N = 2 and N = 3. The repulsion between the singular values is stronger in the type I-III seesaw ensemble than in the type II seesaw ensemble, and the strength of the difference between the repulsions of type I-III and type II ensembles increases as N increases.
The loss of correlation between the light neutrino masses and the light neutrino mass eigenstates forbids an investigation of the favored hierarchy pattern (normal or inverted). However, an analysis of the preferred mass splitting, i.e. the preferred ordering ofm 2 med − m 2 min andm 2 max −m 2 med , is completed. The probability test implies that for both seesaw ensembles, the preferred mass splitting is the one associated to normal hierarchy, although any permutation of the mass eigenstates is equally likely. However, a comparison between ensembles shows that the type I-III seesaw ensemble is only twice as likely as the type II seesaw ensemble to generate the neutrino sector experimental data assuming the preferred normal hierarchy.
For all seesaw mechanisms, the preferred neutrino energy scale is of O(10 −2 ) eV, which leads to a scale of new physics similar to the GUT scale when the associated coupling constants are of order one. Smaller coupling constants can partly lower the new physics scale.

JHEP04(2017)131
A comparison of the group variable pdf for all seesaw ensembles (the Haar measure) and a flat distribution shows that the seesaw ensemble is only twice as likely as the flat distribution to lead to the neutrino sector experimental data. One thus concludes that the type I-III seesaw ensemble is marginally favored over the other ensembles, predicting the hierarchy of the neutrino sector to be normal.
Finally, a comparison of the complex type I-III seesaw ensemble level density for N = 3 and large N shows that the properly-normalized N = 3 level density is well approximated by the properly-normalized large N level density. Because of the complexity of the analytic N = 3 singular value pdf and the link between the large N level density and the physical neutrino sector, it would be interesting to obtain an analytical level density at large N . A step in that direction was made in [11] following the usual Coulomb gas technique, but it was shown there that the resulting level density is wrong. The authors hope to return to this question in the near future.