Empirical inference on the Majorana mass of the ordinary neutrinos

There is a broad theoretical consensus on the idea that ordinary neutrinos have a Majorana mass, but we have no clear prediction about its value, and direct experimental measurements of this quantity are rather challenging. In this work, we argue that the current cosmological measurements allow us to obtain precise information on the effective Majorana mass, i.e. the electronic-type mass of ordinary neutrinos. We show that the numerical results that we obtain can be accurately reproduced, and hence tested, by a straightforward analytical procedure. We then discuss the stability of the assumptions at the basis of our analysis and the implications of our findings for neutrinoless double beta decay.


I. INTRODUCTION
There is compelling evidence from oscillation phenomena [1][2][3][4][5] that ordinary neutrinos are not massless. However, we do not have empirical knowledge of the absolute values of these masses. Nor do we know whether these are of Dirac type -as for all charged fermions -or of Majorana type.
A Majorana nature for the neutrino mass is considered quite natural from the theoretical point of view [6][7][8][9]. This hypothesis could be experimentally tested via the observation of neutrinoless double beta decay (0νββ), the transition (A, Z) → (A, Z + 2) + 2e − where a nucleus increases its charge by two units, emitting two electrons but no neutrinos [10].
The rate of 0νββ depends quadratically upon a combination of the neutrino masses m i , the mixing matrix U ei and the Majorana phases, which is called the "effective Majorana mass" and, within the three-ordinary-neutrino framework, is defined as A key issue is therefore to understand which is the actual value of m ββ . The theory of fermion masses is not yet sufficiently developed to give definitive answers. For instance, it has been argued that a plausible value for m ββ could be of the order of m·θ n C , where θ C 13 • is the Cabibbo mixing angle, m ≡ ∆m 2 atm 50 meV and n = 1, 2 [11,12]. Still, these indications only offer an "educated guess" pointing * sdelloro@vt.edu † francesco.vissani@lngs.infn.it towards the O(1 meV) scale and might be correct within a factor of a few, which practically means that we do not know m ββ reliably. Principled models, on the other hand, such as SO(10) theories, offer more convincing explanations for light Majorana neutrino masses. However, these motivations have been extensively investigated [13][14][15][16][17][18][19][20], but they remain far from being unique and we lack clear criteria on how to assess them. In this situation, it seems reasonable to fully exploit the information coming from the experimental observations in order to try to quantify m ββ . Neutrino oscillations are a powerful tool in this regard, and the mixing angles and mass splittings are currently measured quite precisely [21][22][23]. In particular, the interpretation of the oscillation phenomena is beginning to prefer the so-called normal hierarchy (NH) over the inverted hierarchy (IH) for the neutrino mass spectrum, now at a level greater than 3σ (see the Appendix). However, neutrino oscillations only measure the differences between the neutrino mass eigenstates m 2 i − m 2 j and not the absolute masses. The most relevant information on absolute masses to date comes from cosmological measurements. It is fair to state that a cautious approach is highly advisable while dealing with the results from cosmological surveys, due to the crucial role of the theoretical assumptions in building the specific model. On the other hand, this cannot provide an excuse in order to ignore a whole category of results. Cosmology is able to probe the value of the sum of the active neutrinos that, consistently with the hypothesis of three light neutrinos, is defined as Σ ≡ m 1 + m 2 + m 3 . Very stringent upper limits Σ < 180 meV at the 2σ C. L. have been obtained by many independent analyses [24][25][26][27][28][29][30][31][32].
In a previous work, we showed the importance of the combination of the results of post-2015 Planck observations and "small scale" cosmology measurements for 0νββ [33]. Here, we demonstrate how to obtain detailed information on m ββ , possibly measuring this important arXiv:1909.05381v2 [hep-ph] 10 Oct 2019  In a frequentist approach, one would assume each value of m ββ lying below the current experimental bound to be equally allowed. In the case of NH, this would imply the impossibility of excluding that this parameter is actually negligibly small. However, this is not necessarily the case. In fact, it is possible to obtain useful information on m ββ independently of that coming from the direct experimental search for 0νββ.
As an illustrative case, let us suppose that one day cosmology will be able to measure a value of Σ, and that its value will turn out to be one of the current 2σ upper limits, i. e. 140 meV [30]. Since the mass splittings are reliably known, the lightest neutrino mass would be set to 38 meV. In turn, still assuming that the neutrinos have a Majorana mass, we could conclude that m ββ lies in the range (13 − 39) meV. In particular, the possibility that m ββ = 0 would be excluded for experimental reasons.
Passing to a more realistic case, today we can dispose of a limit Σ from cosmology, but not a measurement. Still, the range of allowed values of m ββ (for a given Σ) can only be explored by declaring a prior distribution for some relevant parameter -that is, by adopting a Bayesian approach.
A prediction for the possible value of m ββ can by ob-tained by simulating a large series of combinations of neutrino mass parameters, starting from the information available today. We thus developed a Monte Carlo tool that iteratively extracts values of m ββ by providing as input randomly generated values of the mass splittings and mixing angles, and of Σ. The former are picked around the best-fit values of the global analysis of Ref. [21] assuming Gaussian errors. The latter parameter is selected within the limit reported in Ref. [30], in which the authors obtained an almost Gaussian likelihood (see Fig. 10 in the reference) that can be approximated by the expression: Indeed, starting from G(Σ) of Eq. (2), we are able to reproduce the same limit Σ < 140 meV at the 2σ C. L. reported in the reference.
For each event that we simulate, the hierarchy scenario (NH or IH) is initially selected with a preference for NH, according to the result of the analysis of Ref. [21]. A value for the lightest neutrino mass, m lightest , is then calculated starting from Σ and the mass splittings [34]. When m lightest is negative, i. e. it falls outside the physical range, the event is rejected and a new one is generated. The allowed interval for m ββ is fixed once m lightest is known and it is comprised within the extremes [35] The value of m ββ is randomly selected within this range with flat probability. In principle, other solutions are possible. One could adopt the most pessimistic attitude, by setting m ββ = m min ββ or, on the contrary, the most optimistic view, by setting m ββ = m max ββ . These two extreme choices for the prior of m ββ could be of some interest in order to assess the chances of failure or success in the search for 0νββ. However, neither of them is particularly supported by existing theoretical considerations. Also, since our primary goals are to illustrate an analysis methodology and to obtain reasonable expectations on m ββ , we did not consider these options. A valid alternative, that has already been used in the literature, is to take random Majorana phase uniformly distributed in [0; 2π) and to construct m ββ accordingly [36][37][38]. In this latter case, one typically expects larger values for m ββ up to ∼ 10% since there is a larger probability that the two phases will add constructively, rather than destructively (Fig. 3). This situation is more conducive to the experimental searches. However, from a theoretical point of view, it seems easier to motivate a preference for a flat prior on m ββ , since this quantity is more directly connected to the Majorana mass term of the Lagrangian density, and thus we opted for this choice.
The result of the simulation is shown in Fig. 1. In total, one billion events have been generated. The contour lines at different C. L. are drawn by "cutting" the distribution on the (m lightest : m ββ ) plane at fixed ratios of the fraction of surviving events over the total one. The IH region is interested by no 1, 2, 3σ contour lines since this scenario is excluded at more than 3σ by the oscillation studies. The tight limits on Σ push towards smaller values of m lightest and this in turn favors smaller values of m ββ . This appears even more clearly by projecting the plot onto the m ββ axis (Fig. 2). These results will be further discussed in the following sections.

III. ANALYTIC TEST OF THE RESULTS
In order to test our Monte Carlo tool, we developed an analytic procedure that allowed us to obtain a distribution of m ββ to be compared to that derived in the previous section (see in particular Fig. 2). We focused on the NH scenario, so that the lightest m i is now m 1 (recall Eq. (1)).
As a starting point, let us consider a priori a flat distribution for each fixed m 1 , namely where the extremes have the known expressions reported in Eqs. (3a)-(3b). As previously stated, we assume a uniform distribution for m ββ within the allowed range. The information on the distribution of m 1 required by Eq.(4) can be extracted from cosmology. The constraints on Σ can be summarized by the likelihood where we can consider Σ as a function of m 1 . We use the expression in Eq. (2) and compute the normalization by integrating G(Σ) down to Σ = 58.5 meV, corresponding to the case m 1 = 0, given the present values of the mass splittings [21]. We thus obtain The distribution of the lightest neutrino mass E(m 1 ) is shown in Fig. 4. The curve has a maximum around m 1 = 10 meV since the Jacobian dΣ/dm 1 , which is a strictly increasing function, disfavors the values close to m 1 = 0, while the likelihood extracted from Ref. [30] does the opposite. It is worth noticing that the propagation of the uncertainties on the mass splitting values introduces a variation in the distribution E(m 1 ). However these are of the order of a few parts per 10,000. Therefore, we can neglect them.
The bidimensional distribution of m ββ and m 1 , which is automatically normalized, is then where The functions describing the extreme values of m 1 given m ββ can be analytically computed, since they correspond to solvable fourth-degree equations. The allowed region in the (m 1 : m ββ ) plane is that shown in Fig. 1, keeping in mind that we are now considering m ββ as the independent variable.
The double integral we need to solve in order to get the probability of a specific m * ββ , i. e. the cumulant distribution L(m * ββ ) ≡ L(m ββ < m * ββ ), can actually be rearranged into a one-dimensional one: .

The upper limit of integration is
while we can identify two cases for the inferior one, depending on the value of m * ββ with respect to m min ββ,0 ≡ m min ββ (m 1 = 0), namely If m * ββ → 0, then min m max ββ (m 1 ), m * ββ = m * ββ and the numerator → 0. Instead, if m * ββ → ∞, then min m max ββ (m 1 ), m * ββ = m max ββ (m 1 ), the fraction approaches 1, the inferior integration limit becomes 0, and the integral of E(m 1 ) over all values is indeed 1. In practice, for sufficiently high confidence levels (corresponding to larger values of m * ββ ), we will refer to the latter value of Eq. (13). Therefore, we can rewrite By using Eq. (14), we can compute m ββ for each given confidence level by setting L(m * ββ ) = C. L. and choosing the oscillation parameters at their central (best-fit) values.
The value of m ββ shifts when we include the variation on the mixing angles and mass splittings, since these enter both the integration limits and the integrand of Eq. (14). The individual contribution associated to the uncertainty of each oscillation parameter x i can be analytically computed. For a fixed confidence level, we get: evaluated in the best-fit central region. In practice, the overall effect is numerically small and the impact is negligible for our purposes.
In order check the consistency between the two described procedures, we can compare the resulting differential distributions for m ββ . The distribution from the likelihood of Eq. (9) can directly be compared with that obtained by the Monte Carlo simulation. The result, illustrated in Fig. 2, shows a very good agreement between the two procedures. This is a rather interesting conclusion: on the one hand, we have validated the numerical results. On the other hand, we have learned that the analytical method described here produces accurate results.

IV. DISCUSSION
The sum of the active neutrino mass coming from cosmology is a key ingredient of the analysis we presented. Therefore, we would like to begin our discussion with some considerations on the robustness of the current measurement of Σ.
To this extent, a first issue to address is the possible existence of additional light neutrinos, beyond those envisaged by the Standard Model. The presence of new neutral fermions was proposed a long time ago [39]. Despite the fact that the values of their masses and mixing with the known neutrinos are a priori unknown, the interpretation of some experiments, when individually taken, would be compatible with the existence of neutrinos with mass of O(1 eV). These neutrinos would give rise to observable oscillation effects. On the other hand, the various analyses that have extensively considered the implications of this hypothesis highlighted the difficulties that arise, and even the inconsistencies, when compared to the available data [40,41]. Therefore, while this could lead to more complicated scenarios for these new neutrinos [42], the "minimal" hypothesis of three ordinary neutrinos considered here is not weakened. 1 Recently, new results on the Hubble constant in the local Universe motivated the idea of rethinking the accepted cosmological model [46]. In principle, this could have an impact on the conclusions of our work. The new observations would be explained by the presence of some "kind" of sterile neutrinos. These neutrinos should have some specific proprieties, namely they should self-interact via strong interactions [47] and should have rather small masses, subject to the stringent limits of cosmology [31]. In other words, while the assumption that there are three active neutrinos in cosmology could be reconsidered, we do not have any evidence that the measurement of Σ should be affected. These considerations can be quantified by an example. Such a type of neutrino, with mass m 0 slightly larger than m 1 , had already been introduced, not only to increase the number of effective neutrinos in cosmology, but also to explain certain features of the observed solar neutrino spectrum [48]. The resulting mixing U 2 e0 ∼ 10 −3 and mass splitting m 2 0 − m 2 1 ∼ 0.2 ∆m 2 [48] would add a contribution U 2 e0 (m 0 − m 1 ) to m ββ of the order of 10 −3 meV at most, and therefore it is completely negligible.
In summary, it seems premature to us to conclude that the the minimal set of assumptions adopted here and the reported results should require significant corrections.
To conclude our discussion, we would like to consider the implications of our findings for 0νββ. Possible predictions for the expected value of m ββ could be of great interest for the experimental community of 0νββ in helping to understand the chances for a positive observation, or to figure out how far from reach a promising target might be.
The current 0νββ experiments have already set limits of the order of (10 25 − 10 26 ) yr on the decay half-life [49][50][51][52][53]. The corresponding upper bounds for m ββ are around 100 meV, with values that actually span a large range, mostly due to the theoretical uncertainties coming from nuclear physics [54], which is a fundamental ingredient in order to extract the information on the neutrino mass. At the same time, the forthcoming generation of experiments is setting very ambitious goals, aiming at sensitivities greater than 10 27 yr. This should allow to explore the parameter space of m ββ of the order of tens of meV [38].
Given the present and near-future experiment sensitivities, we can thus analyze the results shown in Fig. 2. The differential distribution of m ββ is peaked at 4 meV, while the 1σ, 2σ and 3σ intervals extends up to 16,31 and 49 meV, respectively. On the one hand, this means that we should be able to begin to probe the parameter region of interest in the coming years. On the other hand, this means that we are still quite far from fully exploring the "core" of the distribution, i. e. values of m ββ of the order of a few meV. This ultimate investigation would require an extremely challenging multi-tonne experiment [55], but this hypothesis should not be ruled out.
The search for 0νββ remains one of the main ways, if not the only one, to address some major open questions in particle physics: the conservation of the lepton number conservation, and the value of the neutrino mass and its nature!

V. SUMMARY
We developed a Monte Carlo tool to extract the information on the effective Majorana mass starting from the available data on the neutrino masses, coming from the oscillation studies and from cosmology. We found that the distribution of m ββ tends toward low values of the parameter region, with a mode at 4 meV, and a 3σ interval extended up to almost 50 meV. We validated our results with an analytical procedure, whose outcome perfectly matches the numerical one. Finally, we discussed the assumptions at the basis of our analysis and implications of the new information for neutrinoless double beta decay. This work is dedicated to the memory of Simone, a very dear friend and a brilliant scientist, who passed away on August 7th, 2019. He helped in the fulfillment of this article until the very end. We will deeply miss him.

APPENDIX: On the preference for the Normal Hierarchy
In the global analyses, it is customary to present the preference for the hierarchy scenario in terms of the difference ∆χ 2 ≡ χ 2 IH − χ 2 NH > 0 between the two cases. This quantity has a different meaning than in the case of one-dimensional parameters. In order to determine how much is NH preferred by the data, we can estimate the which combines the independent information from cosmology, namely L cosm NH (L cosm IH ≡ 1 − L cosm NH ), and that from neutrino oscillations. The analytical procedure described in the text or, equivalently, the Monte Carlo tool, gives L cosm NH 0.75. In other words, the NH is about 3 times more probable than the IH.