Cornering (quasi) degenerate neutrinos with cosmology

In light of the improved sensitivities of cosmological observations, we examine the status of quasi-degenerate neutrino mass scenarios. Within the simplest extension of the standard cosmological model with massive neutrinos, we find that quasi-degenerate neutrinos are severely constrained by present cosmological data and neutrino oscillation experiments. % % We find that Planck 2018 observations of cosmic microwave background (CMB) anisotropies disfavour quasi-degenerate neutrino masses at $2.4$ Gaussian $\sigma$'s, while adding Baryon acoustic oscillations (BAO) data brings the rejection to 5.9$\sigma$'s. % The highest statistical significance with which one would be able to rule out quasi-degeneracy would arise if the sum of neutrino masses is $\Sigma m_\nu = 60$ \meV (the minimum allowed by neutrino oscillation experiments); % indeed a sensitivity of 15 meV, as expected from a combination of future cosmological probes, would further improve the rejection level up to 17$\sigma$. % We discuss the robustness of these projections with respect to assumptions on the underlying cosmological model, and also compare them with bounds from $\beta$ decay endpoint and neutrinoless double beta decay studies.


I. INTRODUCTION
In analogy to the unification of the Standard Model gauge couplings [1], it has been suggested that perhaps also the masses of the neutrinos may arise from a common seed at high energies. For example, a degenerate neutrino mass spectrum [2] could emerge as a result of some SO(3) family symmetry [3][4][5][6][7] holding at high energies. Alternatively, a simple model with a discrete non-abelian symmetry A4 allows stacking the three lepton families as a triplet [8], leading to quasi-degenerate neutrino masses. In this paper, we remain agnostic as to the underlying theory and consider forms of the mass matrix that could arise from a variety of models.
The goal of this work is to examine the consistency of such quasi-degenerate neutrino mass schemes in light of the improved sensitivity of cosmological observations, as well as improved and upcoming determinations of the neutrino oscillation parameters. We find that degenerate neutrinos are disfavoured by the combination of existing cosmological and oscillation data, being essentially ruled out in the case of inverted neutrino mass ordering, though still allowed within a relatively narrow region of parameters in the case of normal ordering.
Several theoretical scenarios fit within the framework of the models we consider. One simple scenario examined in this paper emerges as a result of the imposition of the A4 family symmetry to the Standard Model [8]. It was originally proposed to provide a symmetry-based realization of the "neutrino unification" scenario suggested in [2]. Within this picture the neutrino mass splittings needed in order to account for neutrino oscillation data [9,10] emerge as calculable radiative corrections [11]. In its simplest form this scenario also predicts maximal atmospheric mixing, θ 23 = π/4, and vanishing reactor mixing angle, θ 13 = 0. However, reactor experiments [12,13] have established that the reactor angle θ 13 is non zero (oscillation data such as those from T2K also indicate a growing hint in favor of leptonic CP violation). This implies that the simplest model needs amendment. Indeed, using the original picture as a "zero-th order" template, the scheme can be easily revamped in order to account for the required value of θ 13 . This makes the CP violating phase δ CP physical, while at the same time providing stringent predictions in the δ CP − θ 23 plane [14]. These will be testable within the upcoming generation of oscillation experiments [15]. Our present dedicated cosmological analysis further constrains the parameters of the mass matrix for this specific scenario, making the degenerate neutrino scenarios strongly disfavoured.
The paper is organized as follows. We first consider a general framework in which the small (solar) mass splitting is neglected, and write the neutrino mass matrix in a simple form that is representative of a wide class of theories. This captures the most relevant features of degenerate schemes. We then extend our analysis to the more general scenario in which both atmospheric and solar mass splittings are taken into account. In all cases, we find that quasidegenerate neutrinos are severely constrained by present cosmological data, at least in the simplest extension of the standard cosmological model with massive neutrinos. We also find that the case of inverted neutrino mass ordering is disfavoured. Future cosmological data would also rule out the surviving parameter regions still allowed for the quasi-degenerate normal ordered neutrino spectrum.

II. PRELIMINARIES
Before moving to the detailed description of the theory setup and the subsequent results, we would like to give some definition and then describe the bottom-line argument of our findings through some general considerations. Although these may fail in capturing the full complexity of the theory model to be discussed later, they serve to highlight the main line of reasoning that eventually will lead us to the results presented in this manuscript.
Strictly speaking, "degenerate neutrino masses" (or, for short, degenerate neutrinos) refers to the case in which the three masses m 1 , m 2 , m 3 are exactly equal, m 1 = m 2 = m 3 , a possibility that is of course excluded by flavour oscillation experiments. In the following we will be concerned about constraining models of quasi-degenerate neutrinos, meaning that the masses are only approximately equal, m 1 m 2 m 3 . This approximate equality amounts to the requirement that the mass differences should be much smaller than the individual masses. In order to make quantitative arguments, like the one that follows, we need to set a criterion as to what "much smaller" means. For definiteness, we will define quasi-degenerate neutrino masses through the requirement that the largest mass difference should be a small fraction ξ (with ξ < 1) of the smallest individual mass.
The squared mass differences ∆m 2 ij ≡ m 2 i − m 2 j are well measured in neutrino oscillation experiments. Global fits [9,10] to oscillation data yield ∆m 2 21 7.6 × 10 −5 eV 2 for the small (solar) squared mass splitting, and |∆m 2 31 | 2.5 × 10 −3 eV 2 for the large (atmospheric) squared mass splitting. Since the sign of |∆m 2 13 | remains unknown, there are two possibilities for the ordering of neutrino masses: the so-called normal (m 1 < m 2 < m 3 ) and inverted (m 3 < m 1 < m 2 ) orderings (NO and IO, respectively). The lowest value of the sum of neutrino masses Σm ν ≡ m 1 + m 2 + m 3 allowed by oscillation measurements is Σm ν > 0.059 eV for normal hierarchy and > 0.10 eV for inverted hierarchy. These minima are found by assuming the lightest neutrino is massless, and using the values reported above for the squared mass differences. However, since we are interested in the case of (quasi) degenerate neutrino masses, the neutrino masses must be substantially higher, so that the difference between the heaviest and lightest neutrino is smaller than any of the neutrino masses.
For the purposes of establishing a quantitative criterion for the definition of "quasi-degenerate", we may ignore the mass difference between m 1 and m 2 , which is so much smaller than the mass difference between m 1 and m 3 . With this approximation m 1 = m 2 , and a straightforward calculation shows that the criterion introduced above reads, for both NO and IO: where m lightest is the mass of the lightest neutrino, i.e. m 1 or m 3 for NO or IO, respectively. Note that neglecting the solar mass splitting is appropriate for the purpose of verifying that the criterion for quasi-degeneracy is satisfied, since the value of the large mass difference |m 3 − m 1,2 | (where one should pick eigenstate 1 or 2 depending on the ordering) mainly depends on |∆m 2 31 |. Moreover, in order to satisfy Eq. (1) with ξ 1, the quantity |∆m 2 31 |/m 2 lightest should be much smaller than unity as well. Expanding the square root in this limit, Eq. (1) becomes Given the measured value of |∆m 2 31 |, for a particular choice for the value of ξ, we can compute the smallest value of m lightest that satisfies the criterion as the value for which the equality in Eq. (1) or (2) holds. Taking ξ = 0.1 yields the condition m lightest > 0.11 eV. In terms of the sum of neutrino masses, a quantity well probed by cosmological observations, this corresponds to Σm ν 0.33 eV.
As a result, an upper bound or a detection from cosmological arguments excluding Σm ν > 0.33 eV at high statistical significance, significantly reduces the parameter region where the neutrino mass spectrum can be degenerate. The latest bounds from the Planck satellite observations of the CMB anisotropies in temperature and polarisation, combined with measurements of the baryon acoustic oscillations (BAO), already corner the degenerate spectrum, providing m ν < 0.12 eV at 95% c.l. 1 . Further improvements are expected from the next generation of cosmological surveys, that will be able to reach a sensitivity of σ( m ν ) ∼ 15 meV. In what follows, we will expand on this basic argument with a more articulated and thorough analysis.

III. THEORY SETUP
In this section, we describe a theoretical setup that might be responsible for quasi-degenerate neutrino mass spectrum. To better illustrate the model, we start in Sec. (III A) from a simple scenario in which the smallest (solar) neutrino mass splitting ∆m 2 12 is neglected, reminiscent of the analysis of Ref. [6]. We then move to Sec. (III B), where we analyse the full-fledged scenario, with both the solar and atmospheric splittings taken into account.

A. Simplest mass matrix
In this section, we begin by neglecting the smallest (solar) mass difference, i.e. we set m 1 = m 2 . We assume that the light neutrino mass matrix M ν , possibly resulting from the seesaw mechanism, has the following form where we use the weak eigenstate basis (ν e , ν µ , ν τ ), m 0 > 0 is the overall neutrino mass scale and |δ| 1 is a small real quantity responsible for the mass difference. This is one of the mixing patterns (specifically, the one called "texture C") appearing in Ref. [6] (see Eq. 7 of that paper), once we identify where m atm and M m atm introduced in Ref. [6] are the scale of atmospheric oscillations, and some higher mass scale, respectively. Notice that, to first order, δ = m atm /M .
A pattern like the one in the mass matrix (3) might emerge, for example, in the presence of a non-Abelian family symmetry, either the small discrete A4 symmetry, or the continuous SO(3) symmetry. These could naturally lead to degenerate neutrinos, with small mass splittings induced as symmetry breaking effects. As an example, Ref. [8] employs an A4 flavour symmetry, with calculable mass differences generated by radiative corrections (arising, for instance, in the context of softly broken supersymmetry). However, for the moment we will be agnostic and just assume that Eq. (3) with δ 1 correctly describes the neutrino mixing pattern. We will, in any case, borrow some useful results from Ref. [8].
Diagonalization of M ν yields the following exact positive eigenvalues: where we have introduced the parameter η ≡ 1 + 2δ, especially in view of the full-fledged scenario that will be discussed in the next section. Note that normal ordering (m 3 > m 1 , m 2 ) is realized for |η| < 1, while inverted ordering (m 3 < m 1 , m 2 ) is realized for |η| > 1. Since the mass eigenstates in Eqs. (5) are independent of the sign of η, in our studies of cosmological bounds we can restrict just to the case η ≥ 0. We can obtain a relation between m 0 and η: where ∆m 2 31 = m 2 3 − m 2 1 is well measured in oscillation experiments. The sum of the individual masses reads m ν = m 1 + m 2 + m 3 = m 0 (1 + 2η), a quantity well probed by cosmological observations. From this one can derive a relation between m ν and η: Note that all results derived so far are exact and do not assume δ 1, i.e., η close to unity. The quasi-degenerate case corresponds to δ 1 i.e. η ∼ 1. Given the measured value of |∆m 2 31 | from oscillations experiments, we have a relation between η and m ν , that allows one to translate cosmological constraints on m ν into constraints on η, and therefore on δ. In other words, using this relation we can use cosmological bounds on neutrino mass to strongly constrain the case of quasi-degenerate neutrinos. The results of this analysis are reported below in Sec. (V A).

B. The full mass matrix
So far we have neglected the mass difference between the states ν 1 and ν 2 characterizing solar neutrino conversions. However, in the realistic case, one needs both mass splitting parameters ∆m 2 21 and ∆m 2 31 to be non-zero, in order to successfully describe the solar and atmospheric neutrino oscillation data. The generalization of the mass matrix in Eq. (3) is given by [8]: This mass matrix reduces to the one discussed in the previous section in the limit in which δ and δ are much smaller than δ. For the moment we take all the parameters to be real. This simplifying approximation amounts to assuming that CP symmetry is conserved in neutrino oscillations, which is sufficient for our purposes 2 . Moreover, without loss of generality we can again take m 0 > 0. The matrix M ν has positive eigenvalues given as One can see that quasi degenerate neutrino models correspond to the case where all three quantities are small: δ, δ , δ 1. For convenience, let us define the following: so that the eigenvalues take the simpler form: Changing sign of either η or η is equivalent to exchanging m 1 and m 2 . Thus we restrict ourselves to the (η > 0, η > 0) quadrant, where m 1 < m 2 . Moreover, exchanging η and η has no effect on the mass eigenvalues, so that in the rest of this section we can further restrict our attention to either of two octants; for definiteness, we take η > η . Quasi-degeneracy corresponds to |η − η | ∼ |η + η | ∼ 1, which, in the first octant, requires η ∼ 1 and η ∼ 0. Cosmological observations bound the sum of the individual masses, thereby placing stringent restrictions on these three parameters where the last equality holds when η > η . As in the two-parameter approximation described in the previous section, one can use the results of neutrino oscillation experiments in order to express m ν as a function of η only. To this aim, we first find the region of the (η, η ) plane that is consistent with current oscillation measurements. The requirement m 2 > m 1 , from the positive There is one region for NO and two disconnected regions for IO. The central dark-shaded region is inconsistent with oscillation data as it corresponds to m1 < m3 < m2. The solid curves are the loci of points that satisfy the constraint from neutrino oscillation measurements on the ratio ∆m 2 12 /|∆m 2 13 |. The whole (η, η ) parameter space is symmetric for reflection around the bisector, i.e. η ↔ η , so that one can focus only on the η > η region (see text for details).
sign of the solar mass splitting, is always satisfied in the first quadrant. As explained above, in the following we consider the η > η region; the results can be easily extended to the rest of the parameter space from the symmetry arguments made above.
It is easy to verify that, when η > η , normal ordering (∆m 2 13 , ∆m 2 23 > 0) is realized for η < 1 − η. Inverted ordering (∆m 2 13 , ∆m 2 23 < 0) is instead realized for η < η − 1. The remaining region of parameter space should be excluded since there m 1 < m 3 < m 2 holds, inconsistent with oscillation experiments. These regions are shown in Fig. (1), where we have used symmetry arguments to reconstruct also the η < η part of the parameter space.
We can further impose that where ∆m 2 1x,obs for x = 2, 3 are the best fit values of the neutrino oscillation global fit analysis [9]. Solving Eq. (14) with Eqs. (12) provides a data-driven relation of the form η = η (η). The relation is one-to-one separately in each of the two regions NO and IO in the first octant. When the two regions are considered together, one finds that, for a given value of η , there are two values of η that satisfy the oscillation constraints, one for each ordering. The curves that satisfy the constraint on ∆m 2 12 /|∆m 2 13 | are shown as solid thick lines in Fig. (1), again after having been extended to the upper octant. Note that the minimum value of η in the first octant that satisfies the oscillation constraint (14) is η 0.088.
We have not yet fully exploited the information coming from neutrino oscillation experiments: given two measured mass differences, we have only required the theory to provide the observed ratio ∆m 2 12 /|∆m 2 13 |. We can use the remaining information to express m 0 in terms of |m 2 3 − m 2 1 | = |∆m 2 31 |. This finally allows us to express the three mass eigenvalues as functions of η only. Thus we can write for the sum of the masses a generalization of Eq. (7): where it is understood that we have used the constraints (14) to express η as a function of η. The corresponding expression in the η > η part of the first quadrant can be obtained from Eq. (15) by the exchange η ↔ η . Constraints on the mass parameters from cosmological data are presented in Sec. (V B).

IV. COSMOLOGICAL DATA
Here we describe in detail the cosmological data samples used to constrain the quasi-degenerate neutrino scenario. We employ the latest results published by the Planck collaboration on the sum of the neutrino masses m ν [16]. These results were obtained assuming a ΛCDM + m ν cosmological model. Such model was tested against the full suite of Planck satellite measurements of the CMB anisotropies in temperature and polarisation, and of the power spectrum of the gravitational lensing potential (a dataset combination labeled as "TTTEEE+lowE+lensing" in the Planck collaboration papers [16]), combined with measurements of the BAO angular scale from 6dF [17], SDSS-MGS [18] and SDSS-DR12 [19] collaborations. We refer to this combination of data as "Planck 2018 + BAO" throughout the manuscript. In some cases, we will also present results for the Planck 2018 dataset alone. For the purposes of this work, we did not run the MCMC analysis. Instead, we have downloaded the MCMC chains provided by the Planck collaboration at the Planck Legacy Archive 3 and reconstructed from them the posterior probability distribution of m ν . Note that this means that we are implicitly assumimng a flat prior on m ν . We make use of the posterior to derive the allowed regions for the parameters in the neutrino mass matrix, using standard statistical techniques.
Future experiments will probe neutrino masses with higher sensitivity. We consider different combinations of future experiments as benchmarks for our projections. CMB observations from the Simons Observatory, combined with the large-scale CMB polarization data from Planck and measurements of the large-scale structure of the Universe, such as those from LSST [20], DESI [21], Euclid [22] are expected to provide a sensitivity on m ν , σ( m ν ) 30 meV. This sensitivity can be improved to σ( m ν ) 20 meV by a cosmic-variance-limited measurement of the reionization optical depth τ , such as that expected from the LiteBIRD satellite [23]. Finally, a sensitivity σ( m ν ) 15 meV is expected from the combination of ultimate CMB experiments, such as CMB-S4 [24], with LiteBIRD and the aforementioned large-scale structure surveys.
To simulate the expected constraints from future cosmology, we interpret the sensitivity on m ν as the square-root of the variance of a Gaussian probability distribution centered in a given fiducial value of m ν . We consider two different fiducial values of m ν , which correspond to two detection scenarios: m ν = 0.06 eV and m ν = 0.1 eV. These correspond to the case in which the "true" value of m ν is the minimal value allowed by measurements of the neutrino mass splittings from oscillation experiments in the normal and inverted ordering, respectively. The highest statistical significance with which one would be able to rule out quasi-degeneracy would be for these cases of minimal allowed neutrino mass. We make use of the Gaussian probability distributions so obtained to derive the expected allowed regions for the parameters in the neutrino mass matrix. The assumption of Gaussianity is expected to provide a good representation of the results from future cosmology, given the expected sensitivity of future surveys.
A final remark concerns the choice of the underlying cosmological model. A well known limitation of the constraints from cosmological probes is the model-dependency issue, i.e. the fact that constraints on cosmological parameters may vary depending on the assumptions on the cosmological model. This happens because there is a certain level of correlation between different cosmological parameters. In other words, the physical effects of one parameter may be compensated by tuning other parameters. Such intrinsic uncertainty of the cosmological analysis can be cured in several ways providing confidence in cosmological results. First of all, one can break the parameter correlation by combining observations of various cosmological probes (CMB and large-scale-structure) that depend differently on the same cosmological parameters. This is the reason why we adopt constraints from combined cosmological probes. Moreover, one can quantify statistically the preference for alternative cosmological models with respect to ΛCDM + m ν . To the best of our knowledge, there is no statistically significant preference reported for extended and/or exotic cosmological models that urges us to consider a different underlying parametrization than the one adopted in this manuscript [25]. Furthermore, the exquisite sensitivity and redundancy of future surveys will help further reduce the impact of model dependency [26]. See e.g. Ref. [27] for a summary concerning the optimal combinations of future cosmological missions. This is the reason why, for the sake of simplicity, we choose to limit our analysis to the ΛCDM + m ν scenario.However, in the conclusions we comment on the impact of considering different cosmological scenarios.

V. RESULTS OF COSMOLOGICAL ANALYSIS
In this section, we report the main findings of our analysis. We use existing and upcoming bounds on the sum of neutrino masses from cosmology to examine the viability of the quasi-degenerate neutrino mass scenario. We continue to study two cases: first, the "simplest mass matrix" presented in Eq. 3, corresponding to negligibly small solar mass splitting and, later, the more general mass matrix presented in Eq. 8.
The basic approach is the following: Taking advantage of relations we have previously found between m ν and parameters η and η in Eqs (7) and (15) above, we will convert cosmological bounds on m ν , into bounds on η and η . As a reminder, quasi-degenerate neutrino models (with small mass differences between species) require η ∼ 1 for all the models we consider as well as η ∼ 0 for the more general mass matrix of Eq. (8). In the latter case, predictions for neutrino masses are unchanged if η and η are exchanged, so that having η ∼ 0 and η ∼ 1 also yields quasi-degenerate neutrino masses. We will show that the combination of oscillations and cosmological bounds is essentially incompatible with such values of η and η .
In addition to studying how well quasi-degenerate neutrino mass can be ruled out from existing data, we also make projections for future data. The highest statistical significance with which one would be able to rule out quasidegeneracy would be for the case of minimal neutrino mass allowed by oscillations data, i.e. Σm ν = 0.06(0.1) eV for NO (IO). Specifically, we examine the bounds on quasi-degeneracy with sensitivity of upcoming experiments ranging from σ (Σm ν ) = 0.030 eV to σ (Σm ν ) = 0.015 eV, assuming the sum of neutrino masses is this minimal allowed value 4 .
A. Results in the simplest mass matrix model Current cosmological data from Planck place limits on the parameters describing the quasi-degenerate neutrino scenario. Within the approximation where the solar neutrino mass splitting is neglected (see Sec. (III A)), the neutrino mass matrix is given in Eq. (3), with the parameters δ and η as defined in Eq. (5). In this simplest case, we have only two parameters, m ν and η, related by Eq. (7); the relation between them depends on ∆m 2 31 , a quantity measured by oscillation experiments, ordering while the dotted orange curve corresponds to the inverted ordering. We recall that in the simplest mass matrix model, η > 1 yields inverted ordering, while η < 1 gives normal ordering.
Cosmological bounds are indicated in Fig. 2 by horizontal bands of various colors. The gray-hatched region is excluded by the current cosmological data (Planck 2018 +BAO). Note that the upper bound on m ν < 0.12 eV quoted in the Planck 2018 parameters paper has been derived assuming a prior m ν > 0. In the present analysis, we are using information from oscillation experiments that require m ν > 0.06 eV, so we should use a prior that reflects this knowledge. Taking this into account yields m ν < 0.15 eV at 95% c.l., which is the value used to produce the gray-hatched exclusion region in Fig. 2. Since quasi-degenerate neutrino masses require δ 1 (η ∼ 1), one can see already by eye that this scenario is ruled out. Roughly, one can see that η 1.8 is required to satisfy the cosmological bounds for inverted ordering, and η 0.7 for normal ordering. This range will be further reduced once data from future cosmological surveys become available. The expected sensitivities of future cosmology are shown as colored horizontal bands for the two fiducial values of m ν introduced in Sec. (IV). These have been chosen as the minimal masses allowed by neutrino oscillations, and correspond, for each ordering, to the strongest rejection for quasi-degeneracy. The viable region of the quasi-degenerate model reduces to the ranges in which the lines overlap with the colored bands.

B. Cosmological bounds on neutrino mass in the full mass matrix approach
We now move to the full theory setup described in Sec. III B. There, we have shown how we can use neutrino oscillation data to express the three mass eigenvalues, and their sum, as functions of η only, see Eqs. (12) and (15). In the rest of this section we will further assume, unless otherwise stated, that η ≥ η , which implies η ≥ 0.088 once oscillation data are taken into account. From the discussion in Sec. III B, it is clear that one can make this choice without loss of generality. Similarly to what we have done in Sec. (V A) for the simplest model, we show in Fig. (2) the sum of neutrino masses m ν as a function of η. The m ν (η) relation can be used to translate cosmological constraints on m ν into constraints on η. This operation requires some care given the two-valued nature of the m ν − η relation and the multimodality of the posterior. In this section we only report our results; the interested reader is referred to the appendix for technical details on how probabilities and exclusion levels are computed. We only report results in terms of η. Credible intervals for η can be obtained using the η (η) relation built as explained in Sec. III B.
We show the posterior for η from Planck2018 and Planck 2018+BAO in Fig. (3), highlighting the multimodality of the probability distribution. The 95% credible intervals on the parameter η (with η ≥ η ) are reported in Tab. I, for different combinations of current and future cosmological datasets. It is clear that the quasi-degenerate case, corresponding to η 1 is strongly disfavoured by the data. It is useful to quantify the preference of the data for non-degenerate neutrinos. To this purpose, we define the quasi-degenerate scenario as the one in which the large mass splitting is smaller than 10% of the overall mass scale. From Planck 2018+BAO data, we get P deg = 4 × 10 −9 , corresponding to 5.9 Gaussian σ's in favour of nondegenerate neutrinos. See the appendix for more details on how P deg is defined and computed. The preference is relaxed to 2.4 σ's (odds of 64 : 1) if only Planck data are considered. In the last column of Tab. I, we also report values of the Bayes factor B between the quasi-degenerate and non-degenerate scenarios, defined as B ≡ P deg /(1 − P deg ) (see appendix for more details).
Since any given pair (η, η ) allowed by oscillation experiments uniquely corresponds to either NO or IO (see Fig. 1, the posterior for η can be used to assess preference for one ordering or the other. We find that Planck 2018 + BAO prefers normal over inverted ordering with odds 3.3 : 1 (1.2σ's).
Endpoint measurements of the Kurie plot of tritium beta decays, explored at the KATRIN experiment [28], provide an independent probe on the absolute scale of neutrino mass, in terms of the effective mass m e of the electron neutrino. This is complementary to what can be achieved through the cosmological observations considered here. KATRIN currently constrains m e < 1.1 eV at 90% CL [29], and is expected, in case of a nondetection after 5 years of operation, to establish an upper limit m e < 0.2 eV (90% CL). It is instructive to compare the numbers derived above from current cosmological data, to what would be obtained from KATRIN in the latter case. Assuming the KATRIN nominal 90% sensitivity, we approximate the posterior P (Σm ν ) from KATRIN data as a semi-Gaussian peaked in 0 and with standard deviation σ(Σm ν ) = 0.36eV. In this case, odds of 2 : 1 are obtained in favour of nondegenerate masses (log 10 B = −0.3).
Future cosmological data would rule out quasi-degenerate neutrinos at the 17σ level if Σm ν = 0.06 eV, assuming a sensitivity σ(Σm ν ) = 0.015 eV. If, instead, Σm ν = 0.10 eV, a possibility already in mild tension with current data, the same sensitivity would yield an exclusion at the ∼ 14σ level, still basically ruling out the quasi-degenerate hypothesis. In Fig. 4 we compare the sensitivity of future cosmological data to Σm ν to the theoretical expectation Σm ν (η), assuming either Σm ν = 0.06 eV or Σm ν = 0.10 eV.
Finally, in Fig. 5, we show the level at which quasi-degenerate neutrinos can be excluded, as a function of the sensitivity σ(Σm ν ), assuming Σm ν = 0.06 eV.   This value yields the highest statistical significance with which one would be able to rule out quasi-degeneracy.
The neutrinoless nuclear double beta decay (A, Z) → (A, Z +2)+2e − (denoted 0νββ) provides another independent and complementary probe of absolute neutrino mass scale, especially important as it constitutes a unique modelindependent test of the Majorana nature of neutrinos [30].
The effective Majorana mass m ββ characterizing 0νββ decay is given as where m i are the neutrino masses, c 12 and s 13 correspond to the angles measured from oscillations and φ 12 , φ 13 are the Majorana phases. Note that the amplitude is expressed using the original, symmetrical parametrization of the lepton mixing matrix introduced in Ref. [31]. The bounds we have derived above from neutrino oscillation experiments as well as cosmology can be compared also with those that follow from the negative searches for neutrinoless double beta decay [32][33][34][35][36][37][38][39]. Current sensitivity should improve significantly in the future, with good prospects for covering the whole region of parameters associated with the inverted ordering spectrum. The caveat, though, is that all these 0νββ decay bounds rely on nuclear physics calculations of the relevant nuclear matrix elements [40,41], which suffer from non-negligible theoretical uncertainties. For this reason, current bounds on m ββ from 0νββ searches are usually expressed as a range of upper limits. In our study of neutrinoless double beta-decay, we use the full mass matrix of Sec. (III B) with the same η as before in Eqns. (12). We compute the effective Majorana mass m ββ as a function of η, following a similar procedure to the one used for Σm ν . The individual masses are computed as described in Sec. (III B). Note that m ββ depends also on the neutrino mixing angles and on the Majorana phases. As mentioned in the introduction, the scheme we have considered so far, through the mass matrix (8), implies θ 13 = 0, a value now excluded by oscillation measurements [12,13]. However, this can be generalized in order to agree with current oscillation data, as shown in Refs. [14,15], without altering significantly the predictions for neutrino masses. For this reason, we fix the mixing angles to their best-fit values when computing m ββ , as we do for the mass splittings. As far as Majorana phases are concerned, these are treated as free parameters and varied in [0, 2π].
We show in Fig. (6) the effective Majorana mass m ββ as a function of η. Without loss of generality, we restrict our attention to the η > η region of the parameter space, as we did when discussing constraints from cosmological data. Note that, due to the variation of the Majorana phases, a range of theoretical values of m ββ corresponds to a single value of η. We also show the upper limits from KamLAND-Zen, that currently provides the most stringent constraints on m ββ , i.e., m ββ < 0.061 − 0.165 eV (90% confidence level) [33]. We see that the quasi-degenerate region η 1 is disfavoured also by current 0νββ data, at a level depending on the value assumed for the nuclear matrix elements entering the calculation of the decay amplitude. For comparison, we also report, in the same figure, the excluded regions for η derived in Sec. (V) from current cosmological data.

VI. CONCLUSIONS AND DISCUSSION
Degenerate neutrino masses can arise via a high-energy symmetry, that is subsequently broken yielding smaller mass splittings. We have considered a mass matrix of the form in Eq. (8) as a template for this class of models (revamping as in Ref. [14,15] is implicit), and derived constraints on the parameters of the matrix using cosmological data together with information from flavour oscillation experiments. A combination of Planck 2018 and BAO data strongly constrains the model parameters, ruling out quasi-degenerate masses at 5.9 Gaussian σ's (2.4σ's if only Planck 2018 CMB data are included). We define the quasi-degenerate scenario as the one in which the large mass splitting is smaller than 10% of the overall mass scale. Laboratory experiments also allow us to probe the absolute neutrino mass scale and the elements of the mixing matrix. We have also compared the constraining power of cosmological data to that of laboratory experiments, in particular KATRIN for β decay, and KamLAND-ZEN for 0νββ decay.The former only provides a weak preference in favour of nondegenerate neutrino masses. On the other hand, the capability of the latter to provide constraints comparable to those from cosmology is currently hindered by our ignorance of both the matrix elements entering the calculation of the decay amplitude for 0νββ as well as the Majorana phases.
The strongest statistical significance with which one could rule out quasi-degeneracy with upcoming experiments is reached if Σm ν = 0.06 eV (the minimum allowed by oscillation experiments) and assuming a sensitivity σ(Σm ν ) = 0.015 eV. One finds that the exclusion of quasi-degenerate neutrinos from cosmological data would improve to 17 Gaussian σ's. If, instead, Σm ν = 0.10 eV (the minimum allowed for inverted ordering), a possibility already in mild tension with current data, the same sensitivity will yield a 14σ-level exclusion, still strongly disfavouring the quasi-degenerate hypothesis. If the actual sum of neutrino masses is higher than these minimal values, then the quasi-degenerate scenario would still be ruled out, albeit at a slightly lower statistical significance.
We now discuss how robust are our cosmological bounds, commenting briefly on how our conclusions would change by considering a different cosmological model. In most extended models, parameter degeneracies act to degrade constraints on neutrino masses. This is typically the case for models in which the curvature density parameter Ω k , or the equation of state parameter w of dark energy are allowed to vary. For example, using a combination of Planck 2015 temperature and low-polarization data and BAO observations, the 95% constraints on Σm ν degrade from 0.19 eV (ΛCDM+Σm ν ) to 0.30 eV (ΛCDM+Σm ν + Ω k ) or 0.31 eV (ΛCDM+Σm ν + w) [42] (see also Ref. [43] for an analysis using Planck 2018 data and using extended models with up to 6 additional parameters). For such models, the conclusions of this paper -that basically rely on the cosmological upper limit on Σm ν -would be weakened, allowing a larger portion of the parameter space for quasi-degenerate neutrinos with respect to ΛCDM. In particular, from the upper bounds quoted above one has that quasi-degenerate neutrinos would be disfavoured at the 3.2σ level in the minimal extension ΛCDM+Σm ν . This rejection would further weaken to the 2.0σ level if w or Ω k are also allowed to vary. Note that these numbers should not be directly compared with the results presented in the main text, since they use different data combinations. They should, however, illustrate how the level at which quasi-degenerate neutrinos are disfavoured changes when the limits on Σm ν are relaxed. Figure 5 can be used to translate upper bounds on Σm ν obtained in an extended model, into an exclusion level for quasi-degenerate neutrinos. We conclude by stressing that enlarging the parameter space beyond ΛCDM + m ν does not always lead to weaker bounds. As a noticeable example we have models of nonphantom dynamical dark energy (that have w(z) ≥ −1 at all redshifts z), in which the constraints on Σm ν are actually slightly tighter than in ΛCDM [44]. Thus our conclusions would still hold, and get slightly stronger, in such dark energy models.

APPENDIX
In this appendix we discuss in more detail how the constraints on η presented in the main text have been derived. The starting point is the relation between (η, η , m 0 ) and the sum of neutrino masses m ν discussed in Sec. (III B). As explained there we can, without loss of generality, study this relation only in the (η > 0, η > 0) region of parameter space, since changing the sign of either η or η would leave m ν unchanged. Moreover, exchanging η and η also does not change the mass eigenstates, so we can further restrict the analysis to one half of the first quadrant; for definiteness, let us take it to be η ≥ η . We have used the information from neutrino oscillation experiments to write m ν η, η , m 0 as a function of η only, i.e. m ν η, η (η), m 0 (η) ≡ F (η). In the first octant of the (η, η ) plane, this relation is one-valued, so that the value of η uniquely determines Σm ν . The opposite is not true, since all values of Σm ν larger than 0.10 eV can be obtained from two distinct values of η in the first octant. Note that oscillation data constrain η ≥ 0.088 in the first octant, since a smaller value would yield Σm ν < 0.06 eV, which is the minimum value allowed by oscillation experiments. The range η ∈ [0.088, ∞) can be further split in two regions, corresponding to NO (η < 1) and IO (η > 1). The case η = 1 corresponds to exact mass degeneracy, which in turn means Σm ν → ∞, while η = 0.088 and η → ∞ yield the minimum masses allowed in NO (Σm ν = 0.06 eV) and IO (Σm ν = 0.10 eV), respectively. F (η) is always increasing (decreasing) in the NO (IO) region. Then F (η) maps η ∈ [0.088 1] to Σm ν ∈ [0.06, ∞) eV and η ∈ [1, ∞) to Σm ν ∈ [0.10, ∞) eV. Note also that the function F (η) can be inverted over each of the two sub-ranges separately. Referring to Fig. 1, we are looking at the half of the blue curve lying in the lower octant, and at the orange curve.
The probability density for η is thus obtained from the posterior P Σmν (Σm ν ) of the sum of neutrino masses, provided by cosmological data, as up to a proportionality constant that can be obtained a posteriori by requiring that ∞ 0 P (η)dη = 1. Note that the latter requirement amounts to the following normalization for the Σm ν posterior (in the following, it should be understood that Σm ν is measured in eV): This is a consequence of the fact that when η varies from 0.088 to ∞, we are traversing the posterior for Σm ν from 0.06 to ∞ and then again from ∞ to 0.10.
The posterior for η can be used to compute credible intervals for this parameter. In particular, in the main text we quote 95% credible intervals. Such an interval I η is defined as an interval containing 95% of the total probability: and is possibly composed by different disconnected regions. The above requirement does not uniquely defines I η , since there are infinitely many intervals covering 95% of the total probability. Different prescriptions exist for singling out one of these intervals. We choose to quote for η the 95% interval with the property that the probability for Σm ν (η) in every point outside the interval is smaller or equal than the probability inside the interval, i.e.: P Σmν (F (η 1 )) ≥ P Σmν (F (η 2 )) for every η 1 ∈ I η , η 2 / ∈ I η (20) This basically amounts to compute the so-called minimum credible interval for Σm ν and map it to the η parameter space to get I η . Note that this is not the same as computing the minimum credible interval for η, since the reparametrization η → Σm ν = F (η) conserves probability mass but not probability density. In other words, while a 95% credible interval remains a 95% credible interval after reparametrization, the condition that probability outside the interval is lower that the probability inside does not necessarily hold in the new parametrization. We choose to compute the minimum credible interval in Σm ν instead than η because the former is the parameter that is more directly constrained by cosmological observations and that gets a flat prior in our analysis, so we regard it (observation-wise) as a "primary" parameter as opposed to η, that we regard as a derived parameter. As noted above, this procedure for constructing the credible interval can result in an interval composed by several disconnected regions, and, in fact, this is nearly always the case in the case under study. The reason is that, in general, there are two regions of large probability, one for each ordering, corresponding to the values of η that yield values of the mass close to the minimum values allowed by oscillation experiments for NO and IO. Since these regions are separated in η space, corresponding to η → 0.088 (NO) and η → ∞ (IO), the resulting credible interval is the union of two disconnected regions. This explains the intervals quoted in Tab. I. The only case in which only one integral appears is for future experiments with sensitivity σ(Σm ν ) = 0.015 eV and fiducial value Σm ν = 0.06 eV. In that case IO is excluded by the data and only the interval corresponding to NO survives.
The fact that the two regions η < 1 and η > 1 correspond to NO and IO also provides a neat way to quantify the preference for one or the other ordering. In particular, one can compute the probabilities for the two orderings, P NO and P IO , as: We can similarly quantify preference for, or against, degenerate neutrinos from current and future data. In order to do this we need to set a threshold that defines "quasi-degenerate". This is somewhat arbitrary; we choose the criterion |1 − η| < 0.1. Note that the oscillation constraints ensure that η η when this criterion is satisfied. Then, from Eqs (12) it is easy to see that |1 − η| is the ratio between the large mass splitting and the overall mass scale m 0 . We then measure the preference for nondegenerate neutrinos by comparing the probabilities enclosed inside and outside the |1 − η| < 0.1 region. In particular we define the probabilty P deg for quasi-degenerate neutrinos as: P deg ≡ P η (0.9 ≤ η ≤ 1.1) and P non−deg = 1 − P non−deg . Given two mutually excluding scenarios (hypotheses) like NO vs. IO, or quasi-degenerate vs. non-degenerate, the information about the preference for one scenario over the other can be conveyed in different ways. Let us call the two scenarios 1 and 2, with P 1 > P 2 , so that P 1 is favoured by the data. A possibility is to directly quote one of the probabilities P 1 and P 2 . Another possibility is to quote the ratio of P 1 to P 2 in terms of odds, like in "scenario 1 is favoured with odds 3:2", meaning that P 1 /P 2 = 3/2, or P 1 = 0.6 and P 2 = 0.4. One can similarly quote the so-called Bayes factor B ≡ P 1 /P 2 , so that in the example above B = 1.5. Finally, in the main text we also translate these probabilities to an equivalent number x of Gaussian standard deviations. This is defined as the value x such that, considering a normal probability distribution with zero mean and unit variance, the probability in the region −x and +x equals P 1 . This leads to the relation P 1 = 1 − erf(x/ √ 2), where erf is the error function. Note that our use of Gaussian standard deviations should not be meant to imply that we employ frequentist statistics. The statistical analysis presented in this paper is perfomed in the framework of Bayesian statistics.