Neutrino Mass Ordering from Oscillations and Beyond: 2018 Status and Future Prospects

The ordering of the neutrino masses is a crucial input for a deep understanding of flavor physics, and its determination may provide the key to establish the relationship among the lepton masses and mixings and their analogous properties in the quark sector. The extraction of the neutrino mass ordering is a data-driven field expected to evolve very rapidly in the next decade. In this review, we both analyze the present status and describe the physics of subsequent prospects. Firstly, the different current available tools to measure the neutrino mass ordering are described. Namely, reactor, long-baseline (accelerator and atmospheric) neutrino beams, laboratory searches for beta and neutrinoless double beta decays and observations of the cosmic background radiation and the large scale structure of the universe are carefully reviewed. Secondly, the results from an up-to-date comprehensive global fit are reported: the Bayesian analysis to the 2018 publicly available oscillation and cosmological data sets provides \emph{strong} evidence for the normal neutrino mass ordering versus the inverted scenario, with a significance of 3.5 standard deviations. This preference for the normal neutrino mass ordering is mostly due to neutrino oscillation measurements. Finally, we shall also emphasize the future perspectives for unveiling the neutrino mass ordering. In this regard, apart from describing the expectations from the aforementioned probes, we also focus on those arising from alternative and novel methods, as 21~cm cosmology, core-collapse supernova neutrinos and the direct detection of relic neutrinos.

The Royal Swedish Academy of Sciences decided to award the 2015 Nobel Prize in Physics to Takaaki Kajita and Arthur B. McDonald "for the discovery of neutrino oscillations, which shows that neutrinos have mass. [...] New discoveries about the deepest neutrino secrets are expected to change our current understanding of the history, structure and future fate of the Universe", see [1][2][3][4][5][6][7] for essential publications. These discoveries robustly established that neutrinos are massive particles. However, neutrinos are massless particles in the Standard Model (SM) of particle physics: in the absence of any direct indication for their mass available at the time, they were introduced as fermions for which no gauge invariant renormalizable mass term can be constructed. As a consequence, in the SM there is neither mixing nor CP violation in the lepton sector. Therefore, neutrino oscillations and masses imply the first known departure from the SM of particle physics.
Despite the good precision that neutrino experiments have reached in the recent years, still many neutrino properties remain unknown. Among them, the neutrino character, Dirac versus Majorana, the existence of CP violation in the leptonic sector, the absolute scale of neutrino masses, and the type of the neutrino mass spectrum. Future laboratory, accelerator and reactor, astrophysical and cosmological probes will address all these open questions, that may further reinforce the evidence for physics beyond the SM. The main focus of this review is, however, the last of the aforementioned unknowns. We will discuss what we know and how we could improve our current knowledge of the neutrino mass ordering.
Neutrino oscillation physics is only sensitive to the squared mass differences (∆m 2 ij = m 2 i − m 2 j ). Current oscillation data can be remarkably well-fitted in terms of two squared mass differences, dubbed as the solar mass splitting (∆m 2 21 7.6 × 10 −5 eV 2 ) and the atmospheric mass splitting (|∆m 2 31 | 2.5 × 10 −3 eV 2 ) [8,9]. Thanks to matter effects in the Sun, we know that ∆m 2 21 > 0 1 . Since the atmospheric mass splitting ∆m 2 31 is essentially measured only via neutrino oscillations in vacuum, which exclusively depend on its absolute value, its sign is unknown at the moment. As a consequence, we have two possibilities for the ordering of neutrino masses: normal ordering (NO, ∆m 2 31 > 0) or inverted ordering (IO, ∆m 2 31 < 0). The situation for the mass ordering has changed a lot in the last few months. The 2017 analyses dealing with global oscillation neutrino data have only shown a mild preference for the normal ordering. Namely, the authors of Ref. [10], by means of a frequentist analysis, found χ 2 IO − χ 2 NO = 3.6 from all the oscillation data considered in their analyses. Very similar results were reported in the first version of [8] 2 , where a value of χ 2 IO − χ 2 NO = 4.3 was quoted 3 . Furthermore, in Ref. [13], the authors verified that the use of a Bayesian approach and the introduction of cosmological or neutrinoless double beta decay data did not alter the main result, which was a weak-to-moderate evidence for the normal neutrino mass ordering according to the Jeffreys' scale (see Table II). The most recent global fit to neutrino oscillation data, however, reported a strengthened preference for normal ordering that is mainly due to the new data from the Super-Kamiokande [14], T2K [15] and NOνA [16] experiments. The inclusion of these new data in both the analyses of Ref. [17] and the 2018 update of Refs. [8,9] increases the preference for normal ordering, which now lies mildly above the 3σ level. In this review we will comment these new results (see section II) and use them to perform an updated global analysis, following the method of Ref. [13] (see section V).
The two possible hierarchical 4 neutrino mass scenarios are shown in Figure 1, inspired by Ref. [18], which provides a graphical representation of the neutrino flavor content of each of the neutrino mass eigenstates given the current preferred values of the oscillation parameters [8], see section II. At present, even if the current preferred value of δ CP for both normal and inverted mass orderings lies close to 3π/2 [8], the precise value of the CP violating phase in the leptonic sector remains unknown. Consequently, in Figure 1, we have varied δ CP within its entire range, ranging from 0 to 2π.
Given the two known mass splittings that oscillation experiments provide us, we are sure that at least two neutrinos have a mass above ∆m 2 21 8 meV and that at least one of these two neutrinos has a mass larger than |∆m 2 31 | 50 meV. For the same reason, we also know that there exists a lower bound on the sum of the three active neutrino masses ( m ν = m 1 + m 2 + m 3 ): m NO ν = m 1 + m 2 1 + ∆m 2 21 + m 2 1 + ∆m 2 31 , m IO ν = m 3 + m 2 3 + |∆m 2 31 | + m 2 3 + |∆m 2 31 | + ∆m 2 21 , 1 Note that the observation of matter effects in the Sun constrains the product ∆m 2 21 cos 2θ 12 to be positive. Therefore, depending on the convention chosen to describe solar neutrino oscillations, matter effects either fix the sign of the solar mass splitting ∆m 2 21 or the octant of the solar angle θ 12 , with ∆m 2 21 positive by definition. 2 See the "July 2017" version in [9]. 3 A somewhat milder preference in favor of normal mass ordering was obtained in the corresponding version of the analysis in Refs. [11,12]. 4 A clarification about the use of "hierarchy" and "ordering" is mandatory. One talks about "hierarchy" when referring to the absolute scales of neutrino masses, in the sense that neutrino masses can be distinguished and ranked from lower to higher. This does not include the possibility that the lightest neutrino mass is much larger than the mass splittings obtained by neutrino oscillation measurements, since in this case the neutrino masses are degenerate. On the other hand, the mass "ordering" is basically defined by the sign of ∆m 2 31 , or by the fact that the lightest neutrino is the most (least) coupled to the electron neutrino flavor in the normal (inverted) case.  mν as a function of the mass of the lightest neutrino, m1 (m3) for the normal (inverted) ordering, in red (blue) respectively. The (indistinguishable) width of the lines represents the present 3σ uncertainties in the neutrino mass splittings from the global fit to neutrino oscillation data [8]. The horizontal bands illustrate two distinct 95% Confidence Level (CL) limits on mν from cosmology, see the text for details.
where the lightest neutrino mass eigenstate corresponds to m 1 (m 3 ) in the normal (inverted) ordering. Using the best-fit values for the neutrino mass splittings in Table I one finds that m ν 0.06 eV in normal ordering, while m ν 0.10 eV in inverted ordering. Figure 2 illustrates the values of m ν as a function of the lightest neutrino mass for the two possible ordering schemes. We also show the two representative bounds on the sum of the neutrino masses from cosmology (discussed later in section IV) which is currently providing the strongest limits on m ν thanks to the fact that neutrinos affect both the evolution of the cosmological background and perturbation quantities (see e.g. the excellent detailed reviews of Refs. [19][20][21][22][23]).
The state-of-knowledge of cosmological observations [24] points to a flat Universe whose mass-energy density includes 5% of ordinary matter (baryons), 22% non-baryonic dark matter, and that is dominated by the dark energy, identified as the motor for the accelerated expansion. This is the so-called ΛCDM Universe, which fits extremely well the Cosmic Microwave Background (CMB) fluctuations, distant Supernovae Ia and galaxy clustering data.
Using the known neutrino oscillation parameters and the standard cosmological evolution, it is possible to compute the thermalization and the decoupling of neutrinos in the early universe (see e.g. [25,26]). While neutrinos decoupled as ultra-relativistic particles, currently at least two out of the three neutrino mass eigenstates are non-relativistic. Neutrinos constitute the first and only known form of dark matter so far. Indeed, neutrinos behave as hot dark matter particles, possessing large thermal velocities, clustering only at scales below their free streaming scale, modifying the evolution of matter overdensities and suppressing structure formation at small scales. The CMB is also affected by the presence of massive neutrinos, as these particles may turn non-relativistic around the decoupling period. However, the strong degeneracy between the Hubble constant and the total neutrino mass requires additional constraints (from Baryon Acoustic Oscillations, Supernovae Ia luminosity distance data and/or direct measurements of the Hubble constant) to be added in the global analyses. In this regard, CMB lensing is also helpful and improves the CMB temperature and polarization constraints, as the presence of massive neutrinos modify the matter distribution along the line of sight through their free streaming nature, reducing clustering and, consequently, CMB lensing. The most constraining cosmological upper bounds to date on m ν can be obtained combining CMB with different large scale structure observations and range from m ν < 0.12 eV to m ν < 0.15 eV at 95% CL [23,[27][28][29][30][31][32], as illustrated in Figure 2.
If the massive neutrino spectrum does not lie in the degenerate region, the three distinct neutrino masses affect the cosmological observables in a different way. For instance, the transition to the non-relativistic period takes place at different cosmic times, and the associated free-streaming scale is different for each of the neutrino mass eigenstates. However, the effect on the power spectrum is very small (permille level) and therefore an extraction of the neutrino mass hierarchy via singling out each of the massive neutrino states seems a very futuristic challenge. This will be possibly attainable only via huge effective volume surveys, as those tracing the 21 cm spin-flip transition in neutral hydrogen, see sections VI D and VI E. On the other hand, should the cosmological measurements of m ν be strong enough to rule out the m ν parameter space corresponding to the inverted ordering (i.e. strong enough to establish in a very significant way that m ν < 0.1 eV), we would know that the neutrino mass ordering must be normal. A word of caution is needed here when dealing with Bayesian analyses, usually performed when dealing with cosmological data: a detection of the neutrino mass ordering could be driven by volume effects in the marginalization, and therefore the prior choice can make a huge difference, if data are not powerful enough [33].
Another way to probe the neutrino mass ordering, apart from direct determinations of the sign of the atmospheric mass splitting ∆m 2 31 in neutrino oscillation experiments and, indirectly, from cosmological bounds on the sum of the neutrino masses, is neutrinoless double β decay [34][35][36][37]. This process is a spontaneous nuclear transition in which the charge of two isobaric nuclei would change by two units with the simultaneous emission of two electrons and without the emission of neutrinos. This process is only possible if the neutrino is a Majorana particle and an experimental signal of the existence of this process would constitute evidence of the putative Majorana neutrino character. The non-observation of the process provides bounds on the so-called effective Majorana mass m ββ , which is a combination of the (Majorana) neutrino masses weighted by the leptonic flavor mixing effects (see section III). Figure 3 illustrates the (Bayesian) 95.5% and 99.7% credible intervals for m ββ as a function of the lightest neutrino mass in the case of three neutrino mixing, considering a logarithmic prior on the lightest neutrino mass. The picture differs from the plot that is usually shown, which features an open band towards increasingly smaller values of m ββ for m lightest 5 meV, due to cancellations which depend on the values of the Majorana phases α i (see section III). In the Bayesian sense of credible intervals, the values of α i which produce such a suppression of m ββ represent an extremely small fraction of the parameter space, which is therefore not relevant when computing the 95.5% and 99.7% credible intervals. In other words, given our knowledge of the neutrino mixing parameters, having m ββ 2 × 10 −4 eV would require some amount of fine tuning in the Majorana phases. This figure is in perfect agreement with the results shown in Figure 1 of Ref. [38], which shows that most of the allowed parameter space is not concentrated at small m ββ if one considers a linear prior on the lightest neutrino mass. We also show the most competitive current limits, as those from KamLAND-Zen (m ββ < 61 − 165 meV at 90% CL) [39], GERDA Phase II (m ββ < 120 − 260 meV at 90% CL) [40] and CUORE (m ββ < 110 − 520 meV at 90% CL) [41]. Please note that a detection of the effective Majorana mass will not be sufficient to determine the mass ordering if the lightest neutrino mass is above ∼ 40 meV: in this case, indeed, the normal and the inverted ordering become indistinguishable from the point of view of neutrinoless double beta decay. Similarly to the case of the cosmological bounds on the neutrino mass m ν , in which only constraining m ν to be below 0.1 eV could be used to disfavor the inverted mass ordering, only a limit on m ββ below ∼ 10 meV could be used to rule out the inverted ordering scheme, and only assuming that neutrinos are Majorana particles.
Since neutrino oscillation measurements, cosmological observations and neutrinoless double beta decay experiments are cornering the inverted mass ordering region, it makes sense to combine their present results. Indeed, plenty of works have been recently devoted to test whether a preference for one mass ordering over the other exists, given current oscillation, neutrinoless double beta decay and cosmological data. A number of studies on the subject [10,[43][44][45][46]    , taking into account the current uncertainties on the neutrino mixing parameters (angles and phases), when three neutrinos are considered. The horizontal regions refer to the most conservative (obtained with the most disfavorable value for the nuclear matrix element of the process) upper bounds from KamLAND-Zen [39], GERDA Phase II [40] and CUORE [41]. The vertical band in the lower panel indicates the strongest limit reported by Planck [42], using the Planck TT,TE,EE + SimLow + lensing data combination.
found that the preference for the normal versus the inverted mass scenario is rather mild with current data, regardless the frequentist versus Bayesian approach. In the latter case, however, the results may be subject-dependent, as a consequence of different possible choices of priors and parameterizations when describing the theoretical model, for example in the case of sampling over the three individual neutrino mass states. Therefore, one must be careful when playing with different priors, as recently shown in Ref. [13]. The current status of the preference of normal versus inverted ordering will be further investigated carefully throughout this review. Furthermore, as it will be carefully detailed in section V, the Bayesian global fit to the 2018 publicly available oscillation and cosmological data points to a strong preference (3.5 standard deviations) for the normal neutrino mass ordering versus the inverted one.
To summarize and conclude this introductory part, we resume that the current available methods to determine the neutrino mass ordering can be grouped as: a) neutrino oscillation facilities; b) neutrinoless double beta decay experiments, with the caveat that the results will only apply in case neutrinos are Majorana fermions; c) CMB and large scale structure surveys.
For each of these three categories we will review the current status and also analyse the future prospects, with a particular focus on the existing experiments which will be improved in the future and on new facilities which aim at determining the neutrino mass ordering in the next ten-to-twenty years. In the second part of this review we will also focus on possible novel methods that in the future will enable us to determine the neutrino mass ordering, as for example future cosmological observations of the 21 cm line, the detection of neutrinos emitted by core-collapse supernovae, measurements of the electron spectrum of β-decaying nuclei and the direct detection of relic neutrinos. We shall exploit the complementarity of both cosmology and particle physics approaches, profiting from the highly multidisciplinary character of the topic. We dedicate sections II, III and IV to explain the extraction of the neutrino mass ordering via neutrino oscillations, β and neutrinoless double β decays and cosmological observations, which will be combined in section V where we present the analysis of current data related to these three data sets. Future perspectives are described throughout section VI and its subsections, while the final remarks will be outlined in section VII.

II. NEUTRINO OSCILLATIONS
Our current knowledge on the neutrino mass ordering comes mainly from the analysis of the available neutrino oscillation data. The sensitivity to the neutrino mass spectrum at oscillation experiments is mostly due to the presence of matter effects in the neutrino propagation. Therefore, one can expect that this sensitivity will increase with the size of matter effects, being larger for atmospheric neutrino experiments, where a fraction of neutrinos travel through the Earth. For long-baseline accelerator experiments, matter effects will increase with the baseline, while these effects will be negligible at short-baseline experiments.
When neutrinos travel through the Earth, the effective matter potential due to the electron (anti)neutrino chargedcurrent elastic scatterings with the electrons in the medium will modify the three-flavor mixing processes. The effect will strongly depend on the neutrino mass ordering: in the normal (inverted) mass ordering scenario, the neutrino flavor transition probabilities will get enhanced (suppressed). In the case of antineutrino propagation, instead, the flavor transition probabilities will get suppressed (enhanced) in the normal (inverted) mass ordering scenario. This is the Wolfenstein effect [47], later expanded by Mikheev and Smirnov [48,49], and commonly named as the Mikheev-Smirnov-Wolfenstein (MSW) effect (see e.g. Ref. [50] for a detailed description of neutrino oscillations in matter).
Matter effects in long-baseline accelerator or atmospheric neutrino oscillation experiments depend on the size of the effective mixing angle θ 13 in matter, which leads the transitions ν e ↔ ν µ,τ governed by the atmospheric mass-squared difference ∆ 31 = ∆m 2 31 /2E. Within the simple two-flavor mixing framework, the effective θ 13 angle in matter reads as sin 2 2θ m 13 = sin 2 2θ 13 where the minus (plus) sign refers to neutrinos (antineutrinos) and N e is the electron number density in the Earth interior. The neutrino mass ordering fixes the sign of ∆ 31 , that is positive (negative) for normal (inverted) ordering: notice that, in the presence of matter effects, the neutrino (antineutrino) oscillation probability P (ν µ → ν e ) [P (ν µ → ν e )] gets enhanced if the ordering is normal (inverted). Exploiting the different matter effects for neutrinos and antineutrinos provides therefore the ideal tool to unravel the mass ordering.
Matter effects are expected to be particularly relevant when the following resonance condition is satisfied: The precise location of the resonance will depend on both the neutrino path and the neutrino energy. For instance, for ∆m 2 31 ∼ 2.5 × 10 −3 eV 2 and distances of several thousand kilometers, as it is the case of atmospheric neutrinos, the resonance effect is expected to happen for neutrino energies ∼ 10 GeV.
In the case of muon disappearance experiments, in the ∼ GeV energy range relevant for long-baseline and atmospheric neutrino beams, the P µµ survival probabilities are suppressed (enhanced) due to matter effects if the ordering is normal (inverted). If the matter density is constant, the P µµ survival probability at terrestrial baselines is given by P µµ = 1 − cos 2 θ m 13 sin 2 2θ 23 × sin 2 1.27 − sin 2 θ m 13 sin 2 2θ 23 × sin 2 1.27 where A = 2 √ 2G F N e E, θ m 13 is that of Equation (2) and The dependence of the survival probability P µµ on the neutrino energy E and the cosine of the zenith angle cos θ z , related to the distance the atmospheric neutrinos travel inside the Earth before being detected at the experiments, is shown in Figure 4 for normal (left panel) and inverted (right panel) ordering. There, we can see that reconstructing the oscillation pattern at different distances and energies allows to determine the neutrino mass ordering (see also section VI A). Until very recently, oscillation experiments were not showing a particular preference for any of the mass orderings, not even when combined in a global analysis, see for instance Ref. [51]. Lately, however, the most recent data releases from some of the experiments have become more sensitive to the ordering of the neutrino mass spectrum. In particular, the long-baseline experiments T2K and NOνA on their own obtain a slight preference in favor of normal mass ordering, with ∆χ 2 ≈ 4 each [15,16]. Note that these results have been obtained imposing a prior on the mixing angle θ 13 , according to its most recent determination at reactor experiments. Relaxing the prior on the reactor angle results in a milder preference for normal over inverted mass ordering. The latest atmospheric neutrino results from Super-Kamiokande also show some sensitivity to the neutrino mass ordering. In this case, the collaboration obtains  a preference for normal ordering with ∆χ 2 ≈ 3.5, without any prior on the reactor angle. Constraining the value of θ 13 , the preference for normal mass ordering increases up to ∆χ 2 ≈ 4.5 [14]. The full sensitivity to the ordering of the neutrino mass spectrum from oscillations is obtained after combining the data samples described above with all the available experimental results in a global fit [8]. This type of analysis exploits the complementarity among the different results as well as the correlations among the oscillation parameters to obtain improved sensitivities on them. In the global analysis to neutrino oscillations, the parameters sin 2 θ 12 and ∆m 2 21 are rather well measured by the solar experiments [52][53][54][55][56][57][58][59][60][61] and the long-baseline reactor experiment KamLAND [62]. The short-baseline reactor neutrino experiments Daya Bay [63], RENO [64] and Double Chooz [65] are the most efficient ones in measuring the reactor angle θ 13 and also measure very well the atmospheric mass splitting, ∆m 2 31 . This mass splitting is also measured, together with the atmospheric angle θ 23 , by the atmospheric experiments IceCube-DeepCore [66], ANTARES [67] and Super-Kamiokande [14], where the latter shows some sensitivity to θ 13 and δ CP , too. The long-baseline accelerator experiments are also sensitive to these four parameters through their appearance and disappearance neutrino channels. Apart from the already mentioned T2K [15] and NOνA [16], the global fit also includes the previous experiments K2K [68] and MINOS [69].
The result of the global analysis is summarized in Table I and Figure 5. Before discussing the sensitivity to the neutrino mass ordering, we shall briefly discuss some other features of this global fit. Notice first that now the best-fit value for the atmospheric mixing angle θ 23 lies in the second octant, although values in the first octant are still allowed with ∆χ 2 = 1.6 (3.2) for normal (inverted) ordering. Therefore, the octant problem remains unsolved so far. Note also that, for the first time, the CP violating phase δ CP is determined with rather good accuracy. The best-fit values for this parameter lie close to maximal CP violation, being δ CP = 1.32π for normal ordering and δ CP = 1.56π for inverted ordering. As can be seen from the ∆χ 2 profile in Figure 5, values around δ CP ≈ 0.5π are now highly disfavored by data. Indeed, only around 50% of the parameter space remains allowed at the 3σ level, roughly the interval [0.9π, 1 for normal and [1.1π, 1.9π] for inverted ordering. In the case of normal ordering, CP conservation remains allowed at 2σ, while it is slightly more disfavored for inverted ordering. For the remaining oscillation parameters, one clearly sees that neutrino oscillations are entering the precision era, with relative uncertainties on their determination of 5% or below. For a more detailed discussion about these parameters we refer the reader to Refs. [8,9]. Concerning the neutrino mass ordering, we obtain a global preference of 3.4σ (∆χ 2 = 11.7) in favor of normal ordering. This result emerges from the combination of all the neutrino oscillation experiments, as we explain in the following. Starting with long-baseline data alone, the inverted mass ordering is disfavored with ∆χ 2 = 2.0, when no prior is considered on the value of θ 13 . However, as explained above, the separate analysis of the latest T2K and NOνA data independently report a ∆χ 2 ≈ 4 among the two possible mass orderings when a prior on the reactor angle is imposed. This comes from the mismatch between the value of θ 13 preferred by short-baseline reactor and long-baseline accelerator experiments, which is more important for inverted ordering. Besides that, the combination of T2K and reactor data results in an additional tension relative to the preferred value of the atmospheric mass splitting ∆m 2 31 , which is again larger for the inverted mass ordering. This further discrepancy results in a preference for normal ordering with ∆χ 2 = 5.3 for the combination of "T2K plus reactors" and ∆χ 2 = 3.7 for the combination of "NOνA plus reactors". From the combined analysis of all long-baseline accelerator and short-baseline reactor data one obtains a ∆χ 2 = 7.5 between normal and inverted ordering, which corresponds to a preference of 2.7σ in favor of normal mass ordering. By adding the atmospheric data to the neutrino oscillations fit, we finally obtain ∆χ 2 = 11.7 5  The most reliable method to determine the absolute neutrino masses in a completely model independent way is to measure the spectrum of β-decay near the endpoint of the electron spectrum. The reason for this is related to the fact that, if neutrinos are massive, part of the released energy must go into the neutrino mass and the electron spectrum endpoint shifts to lower energies. When there are more than one massive neutrino, each of the separate mass eigenstates contributes to the suppression of the electron energy spectrum and it becomes possible to study the pattern of the neutrino masses. Nowadays none of the β-decay experiments can reach the energy resolution required to be able to determine the mass hierarchy 6 , but we will explain in the following how, in principle, future experiments may aim at such result.
The best way to depict the effects of the separate mass eigenstates is to compute the Kurie function for β-decay. The complete expression can be written as (see e.g. Ref. [71]): where Q β is the Q-value of the considered β-decay, T is the electron kinetic energy, Θ is the Heaviside step function and |U ei | 2 is the mixing matrix element that defines the mixing between the electron neutrino flavor and the i-th mass eigenstate with a mass m i . The standard scenario features N = 3, but the formula is valid also if a larger number of neutrinos exists (i.e. if there are sterile neutrinos, as explained for example in Ref. [72]). The Kurie function of Equation (6) is depicted in Figure 6, where we show in red (blue) the result obtained using a massless lightest neutrino and the current best-fit mixing angles and mass splittings for normal (inverted) ordering, as described in the previous section. As a reference, we also plot K(T ) for a case with massless neutrinos only (in black). Should we consider higher values for the lightest neutrino mass, the detection of the mass ordering would be increasingly more difficult, since the separation of the mass eigenstates would decrease, eventually becoming negligible in the degenerate case. For this reason we will only discuss the case of a massless lightest neutrino from the perspective of the β-decay experiments.
Given the unitarity of the mixing matrix ( N i=1 |U ei | 2 = 1), the normalization of the Kurie function is the same at Q β −T m i . Since we are interested in the small differences which appear near the endpoint, the plot only focuses on the very end of the electron spectrum and the common normalization is not visible for the inverted ordering case. In the considered range, however, the effect of the different correspondence between the mass eigenstates and the mixing matrix elements introduces a difference which in principle would allow to determine the mass ordering through the observation of the β-decay spectrum. The observation of the kinks in the electron spectrum is very challenging, especially in the case of normal ordering, for which even the more pronounced kink (at Q β − T ∆m 2 21 8 meV) is barely visible in Figure 6. In the case of inverted ordering, since the mass difference between the two lightest mass eigenstates is the largest possible one ( ∆m 2 31 50 meV), and the lightest neutrino is the one with the smallest mixing with the electron neutrino, the amplitude of the kink is much larger. As a consequence, an experiment with enough energy resolution to measure the spectrum in the relevant energy range can directly probe the mass ordering observing the presence of a kink. Note that this measurement can be obtained even without a detection of the lightest neutrino mass. As we show in Figure 6, however, it is crucial to have a non-zero observation of the electron spectrum between Q β and Q β − ∆m 2 31 , otherwise one could confuse the inverted ordering spectrum with a normal ordering spectrum obtained with a larger lightest neutrino mass m lightest ∆m 2 31 (green curve). Another consideration is due. One could think to probe the neutrino mass ordering just using the fact that the expected number of events is smaller in the inverted ordering than in the normal ordering case. As we discussed above, this could be possible only if some independent experiment would be able to determine the mass of the lightest neutrino, in order to break the possible degeneracy between m lightest and the mass ordering depicted by the blue and green curves in Figure 6, otherwise the conditions required to observe the electron spectrum between Q β and Q β − ∆m 2 31 would be probably sufficient to guarantee a direct observation or exclusion of the kink. The best way to determine the neutrino mass ordering, however, may be to use an estimator which compares the binned spectra in the normal and inverted ordering cases, as proposed for example in Ref. [73] in the context of reactor neutrino experiments. The authors of the study, indeed, find that a dedicated estimator can enhance the detection significance with respect to a standard χ 2 comparison.
To conclude, today the status of β-decay experiments is far from the level of determining the mass ordering, since the energy resolution achieved in past and current measurements is not sufficient to guarantee a precise probe of the interesting part of the spectrum. KATRIN, for example, aims at a sensitivity of 0.2 eV on the effective electron neutrino mass [74,75], only sufficient to probe the fully degenerate region of the neutrino mass spectrum.

B. Mass ordering from neutrinoless double beta decay
In the second part of this section we shall discuss instead the perspectives from the neutrinoless double beta decay (see e.g. the reviews [34,35]), a process allowed only if neutrinos are Majorana particles [76], since it requires the lepton number to be violated by two units. Neutrinoless double beta decay experiments therefore aim at measuring the life time T 0ν 1/2 of the decay, which can be written as: where m e is the electron mass, G N 0ν is the phase space factor, M 0ν N is the nuclear matrix element (NME) of the neutrinoless double beta decay process, N indicates the chemical element which is adopted to build the experiment and m ββ is the effective Majorana mass, see below. In case no events are observed, a lower bound on T 0ν 1/2 can be derived. Recent constraints on the neutrinoless double beta decay half-life come from the EXO-200 [77], KamLAND-Zen [39], CUORE [41], Majorana [78], CUPID-0 [79], Gerda [40] and NEMO-3 [80] where N is the number of neutrino mass eigenstates, each with its mass m k , α k are the Majorana phases (one of which can be rotated away, so that there are N − 1 independent phases), and U ek represents the mixing between the electron flavor neutrino and the k-th mass eigenstate. Notice that the conversion between the half-life of the process and the effective Majorana mass depends on the NMEs (see e.g. [81,82]), which are typically difficult to compute. Several methods can be employed and there is no full agreement between the results obtained with the different methods. As a consequence, the quoted limits on T 0ν 1/2 can be translated into limits on m ββ which depend on the NMEs. If the most conservative values for the NMEs are considered, none of the current constraints reaches the level required to start constraining the inverted ordering in the framework of three neutrinos, see Figure 3.
If we compute m ββ as a function of the lightest neutrino mass with the current preferred values of the mixing parameters and in the context of three neutrinos, we discover that the value of m ββ depends on the mass ordering only for m lightest 40 meV, see Figure 3. For this reason, neutrinoless double beta experiments can aim to distinguish the mass ordering only for the smallest values of the lightest neutrino mass. Please note that this also means that if the lightest neutrino has a mass above ∼ 40 meV, perfectly allowed by all the present constraints on the neutrino mass scale, the two mass orderings will never be distinguished in the context of neutrinoless double beta decay experiments.
When going to smaller m lightest , the situation changes, as m ββ becomes independent of m lightest . In the region m lightest 10 meV, a difference between the two mass orderings appears, since the effective Majorana mass is constrained by the mass splittings to be larger than ∼ 10 meV for inverted ordering, while it must be below ∼ 7 meV for normal ordering. This means that experiments which can test the region m ββ < 10 meV can rule out the inverted scenario. Note that a positive detection of T 0ν 1/2 in the range that corresponds to m ββ 10 meV, on the other hand, would not give sufficient information to determine the mass ordering without an independent determination of m lightest . To resume, in the context of three neutrino mixing, neutrinoless double beta decay experiments alone will be able to determine the neutrino mass ordering only ruling out the inverted scheme, that is to say if the ordering is normal and m lightest 10 meV.
In any case, we should remember that if no neutrinoless double beta decay candidate event will ever be observed we will not have determined the mass ordering univocally: Dirac neutrinos escape the constraints from this kind of process, so that it would be still perfectly allowed to have an inverted ordering scheme and no Majorana fermions in the neutrino sector. Due to the presence of the Majorana phases in Equation (8), unfortunately, there is a small window for m lightest in normal ordering that can correspond to almost vanishing values of m ββ , which will possibly never be observable. As we show in Figure 3, however, the region of parameter space where this happens has a very small volume if one considers the phases to vary between 0 and 2π, so that the credible region for m ββ in a Bayesian context shows that it is rather unlikely to have m ββ 2 × 10 −4 eV, as a significant amount of fine tuning would be needed in the (completely unknown) Majorana phases. Our statement, which arises from assuming a logarithmic prior on m lightest , is in perfect agreement with the results of Ref. [38], where a linear prior on m lightest is assumed.
Please note that the situation depicted in Figure 3 is only valid if there are only three neutrinos. If, as the current DANSS [83] and NEOS [84] experiments may suggest, a sterile neutrino with a mass around 1 eV exists (see e.g. [85][86][87][88]), the situation would be significantly different. The allowed bands for m ββ as a function of the lightest neutrino mass when a light sterile neutrino is introduced are reported for example in Ref. [89] (see also Ref. [72]). In this three active plus one sterile neutrino case (3+1), the contribution of the fourth neutrino mass eigenstate (mainly mixed with the sterile flavor) must be added in Equation (8), with the consequence that the allowed bands are located at higher m ββ . In Figure 7, adapted from Ref. [90], we reproduce the dependence of the effective Majorana mass on the lightest neutrino mass when one assumes the 3+1 neutrino scenario, compared with the standard three neutrino case. As we can see, with the introduction of an extra sterile neutrino state, m ββ is significantly increased for the normal ordering case, reaching the level of the inverted ordering bands, which are less shifted towards higher values of m ββ . Furthermore, in the 3+1 scenario, also in the inverted ordering case it is possible to have accidental cancellations due to the three independent Majorana phases in Equation (8) (see the detailed discussion of Ref. [89]), so that a non-detection of the neutrinoless double beta decay process would never be sufficient to rule out the inverted ordering. The opposite situation may occur in case the lightest neutrino mass will be independently constrained to be below ∼ 10 meV while m ββ 10 meV: in this case, however, we would rule out normal ordering. Consequently, if a light sterile neutrino exists, neutrinoless double beta experiments will never be able to determine the mass ordering if the mass ordering is normal, while some possibility remains if the ordering of the three active neutrino masses is inverted, provided that the lightest neutrino is very light and the Majorana phases are tuned enough. The KamLAND-Zen, Gerda and CUORE experiments, using three different materials, may very soon start probing the inverted ordering region in the case of 3+1 neutrino mixing for all the possible values of the NMEs, see Figure 7, where the current KamLAND-Zen [39] constraints are reported.
To conclude and summarize the current status: neutrinoless double beta decay cannot yet provide constraints on the neutrino mass ordering. Depending on the lightest neutrino mass and on the existence of a fourth (sterile) neutrino, it would be possible that not even far-future experiments could be able to reach this goal. Adapted from Ref. [90]. The green band represents the 90% CL bounds from KamLAND-Zen [39], given the uncertainty on the NME.

IV. RESULTS FROM COSMOLOGY
Massive neutrinos affect the cosmological observables in different ways, that we shall summarize in what follows. For a comprehensive review of the effects of neutrino masses in cosmology, we refer the reader to the recent work presented in [23].
A very important epoch when discussing the impact of massive neutrinos in the cosmological expansion history and in the perturbation evolution is the redshift at which neutrinos become non-relativistic. This redshift is given by with m i referring to the mass of each massive neutrino eigenstate. Current bounds on neutrino masses imply that at least two out of the three massive eigenstates became non-relativistic in the matter dominated period of the universe. As stated in the introductory section, and as we shall further illustrate along this section, cosmological measurements are currently unable to extract individually the masses of the neutrino eigenstates and the ordering of their mass spectrum and, therefore, concerning current cosmological data, all the limits on the neutrino mass ordering will come from the sensitivity to the total neutrino mass m ν . Consequently, in what follows, we shall mainly concentrate on the effects on the cosmological observables of m ν , providing additional insights on the sensitivity to the ordering of the individual mass eigenstates whenever relevant.

A. CMB
There are several imprints of neutrino masses on the CMB temperature fluctuations pattern once neutrinos become non-relativistic: a shift in the matter-radiation equality redshift or a change in the amount of non-relativistic energy density at late times, both induced by the evolution of the neutrino background, that will, respectively, affect the angular location of the acoustic peaks and the slope of the CMB tail, through the Late Integrated Sachs Wolfe (ISW) effect. The former will mostly modify Θ s , i.e. the angular position of the CMB peaks, which is given by the ratio of the sound horizon and the angular diameter distance, both evaluated at the recombination epoch. Massive neutrinos enhance the Hubble expansion rate, with a consequent reduction of the angular diameter distance and an increase of Θ s , which would correspond to a shift of the peaks towards larger (smaller) angular scales (multipoles). The latter, the Late ISW effect, is related to the fact that the gravitational potentials are constant in a matter-dominated universe. The inclusion of massive neutrinos will delay the dark energy dominated period and consequently reduce the time variation of the gravitational potential at late times, suppressing the photon temperature anisotropies in the multipole region 2 < < 50. A very similar effect occurs at early times through the so-called Early ISW effect, which governs the height of the first CMB peak. Light active neutrino species, indeed, reduce the time variations of the gravitational potential also around the recombination period, due to the different evolution of these potentials in radiation/matter dominated epochs, leaving a signature on the CMB photon fluctuations when they become non-relativistic. Massive neutrinos will therefore decrease the temperature anisotropies by ∆C /C ∼ (m ν,i /0.1 eV) % in the multipole range 20 < < 500 [20]. Unfortunately, the Late ISW effect affects the CMB spectrum in a region where cosmic variance does not allow for very accurate measurements. From what regards the other two effects, i.e. the shift in the location of the acoustic peaks and the Early ISW effect, they can both easily be compensated varying other parameters which govern the expansion of the universe. For example, within the minimal ΛCDM framework, the total amount of matter in the universe and the Hubble constant H 0 can be tuned in order to compensate the effects of massive neutrinos. Therefore, CMB primary anisotropies alone can not provide very tight bounds on the neutrino masses, due to the strong parameter degeneracies. This automatically implies that CMB measurements alone are unable to extract any information concerning the neutrino mass ordering, as shown in Figure 8, obtained by means of the publicly available Boltzmann solver Cosmic Linear Anisotropy Solving System (CLASS) [91][92][93][94]. In the figure we can notice that the difference between normal and inverted neutrino mass orderings, for m ν = 0.12 eV 7 is almost negligible. Moreover, the largest differences appear in the multipole range where cosmic variance dominates.
Among the secondary CMB anisotropies, i.e. those generated along the photon line of sight and not produced at recombination, there are two effects that can notably improve the sensitivity to the total neutrino mass m ν from CMB observations. One of them is CMB lensing, that is, a distortion of the photon paths because of the presence of matter inhomogeneities along the line of sight. Due to such distortion, the CMB acoustic oscillation features will be smeared out in a scale-dependent way, mostly due to matter overdensities at z 5. By measuring the nongaussianities of CMB polarization and temperature maps it is possible to extract the power spectrum of the lensing potential. This, in turn, contains very useful information on the integrated matter distribution along the line of sight. Since massive neutrinos behave differently from a pure cold dark matter component, characterized by zero velocities, the small-scale structure suppression induced by the non-negligible neutrino dispersion velocities will decrease the CMB lensing signal expected in the absence of neutrino masses [95][96][97][98][99][100], leaving unchanged the power spectrum of the lensing potential at large scales, and suppressing it at small scales. Furthermore, since CMB lensing involves high redshifts, non-linearities do not enter in the calculation of the matter density field. Therefore, CMB lensing enhances the capabilities to bound the neutrino masses using CMB data. In the future, this technique may even surpass weak lensing capabilities, based on statistical analyses of the ellipticity of remote galaxies, see below and section VI D. Indeed, nowadays, measurements from the Planck satellite constrain the neutrino masses dominantly through CMB gravitational lensing. As stated in Ref. [101], increasing the neutrino mass implies an increase on the expansion rate at redshifts z ≥ 1, corresponding to a suppression of clustering at scales below the size of the horizon at the non-relativistic transition. This effect leads to a decrease in CMB lensing that, at multipoles = 1000, is ∼ 10% for m ν = 0.66 eV.
On the other hand, we have the reionization process in the late universe, when the first generation of galaxies emitted ultraviolet (UV) photons that ionized the neutral hydrogen, leading to the end of the so-called dark ages. Reionization increases the number density of free electrons n e which can scatter the CMB with a probability given by a quantity named reionization optical depth, τ , which can be computed as an integral over the line of sight of n e . The consequence of an increase of τ on the CMB temperature fluctuations is the suppression of the acoustic peaks by a factor exp(−2τ ) at scales smaller than the Hubble horizon at the reionization epoch. Even if from the point of view of CMB temperature anisotropies this effect is highly degenerate with a change in the amplitude of the primordial power spectrum A s , which governs the overall amplitude of the CMB spectra, reionization induces linear polarization on the CMB spectrum, leading to a "reionization bump" in the polarization spectra at large scales, which otherwise would vanish. Even if the reionization signal is rather weak, as it amounts to no more than ∼ 10% of the primary polarization signal [102], very accurate measurements of the reionization optical depth τ sharpen considerably the CMB neutrino mass bounds [28], as they alleviate the degeneracy between A s and τ and consequently the existing one between m ν and A s .

B. Large scale structure of the universe
The largest effect of neutrino masses on the cosmological observables is imprinted in the matter power spectrum [103,104]. Neutrinos are hot dark matter particles and, therefore, due to the pressure gradient, at a given redshift z, the non-relativistic neutrino overdensities can only cluster at scales for which the wavenumber of perturbations is below the neutrino free streaming scale k fs (i.e. at scales k < k fs ), with being Ω m the ratio of the total matter energy density over the critical density at redshift zero. The free-streaming nature of the neutrino will be directly translated into a suppression of the growth of matter fluctuations at small scales. One could then conclude that extracting the neutrino relic masses and their ordering is a straightforward task, once that measurements of the matter power spectrum at the relevant scales are available at a different number of redshifts. The former statement is incorrect, not only because it does not consider the existence of degeneracies with the remaining cosmological parameters, but also because a number of subtleties must be taken into account, as we shall explain in what follows. The decrease of the matter power spectrum due to the total neutrino mass m ν is, in principle, currently measurable. Nonetheless, when fixing m ν , the total mass could be splitted differently among the three neutrino mass eigenstates (i.e. m 1 , m 2 and m 3 ), modifying slightly the relativistic to non-relativistic transition. This will affect both the background evolution and the perturbation observables [105]: the different free-streaming scales associated to each of the three neutrino mass eigenstates will be imprinted in the matter power spectrum. Figure 9 shows the ratios of the matter power spectrum for normal over degenerate, inverted over degenerate, and inverted over normal neutrino mass spectra for a total neutrino mass of 0.12 eV. We illustrate such ratios at different redshifts. Notice that the differences among the possible neutrino mass schemes are tiny, saturating at the 0.06% level at k > 0.2h Mpc −1 . Therefore, only very futuristic means of measuring the matter power spectrum could be directly sensitive to the neutrino mass ordering, and, eventually, be able to isolate each of the free-streaming scales associated to each individual neutrino mass eigenstate. We shall comment on these future probes in section VI E.
Since, at present, matter power spectrum data constrain exclusively m ν , it is only via these bounds, combined with CMB or other external data sets, that nowadays a limit on the neutrino mass ordering can be obtained, see section IV C. Nevertheless, and as aforementioned, there are a number of problems which may interfere with a proper understanding of the scale-dependent neutrino mass suppression of clustering. The first of them is due to the fact that observations measure the spatial distributions (or their Fourier transforms, the power spectra) of galaxies, clusters, or quasars, e.g. of given tracers, mapping the large scale structure of the Universe at a number of redshifts, by measuring the growth of fluctuations at different scales. However, the matter distribution is not directly measured, i.e. it needs to be inferred from the tracers observed. A simple model of structure formation suggests that at large scales and, therefore, when the perturbation evolution is still in the linear regime, the galaxy power spectrum is related to the matter one by a constant named b, the light-to-mass bias [106]. The galaxy bias can be determined either separately by independent methods or to be considered as an additional free parameter to be measured together with the neutrino mass m ν . This approach has been followed in many studies in the past [28][29][30]. However, when dealing with neutrino masses, the relationship between the tracers and the underlying matter field may be more complicated, as neutrinos themselves may induce scale-dependent features in the bias [107][108][109] due to their free-streaming nature (see also the recent work of Ref. [110] for a new method to extract a scale-dependent bias, based on the cross-correlation of CMB lensing maps and galaxy samples).
Another additional complication when extracting the neutrino mass from clustering observations arises from to the presence of non-linearities at scales k k 0 The effect of neutrino masses is very-well understood on linear scales, i.e. scales below k NL at z = 0 (or located at slightly larger values of k but at higher redshifts). Massive neutrinos induce a suppression in the linear matter power spectrum below their free streaming scale ∆P/P ∝ −8f ν , with f ν the fraction of matter in the form of massive neutrinos [104]. Accurate descriptions of the matter power spectrum in the non-linear regime are therefore mandatory in order to be sensitive to the full neutrino mass signature. This is particularly important in the case of galaxy surveys, in which the information depends on the number of independent modes available, and where going to smaller scales (i.e. larger values of k) has a profound impact on the sensitivity to neutrino masses. Several approaches have been followed in the literature to account for the effect of massive neutrinos in the non-linear regime, most of them relying on N-body cosmological simulations, which have been upgraded to include the effects of neutrino clustering in the evolution of the cosmological structures. Methods range from perturbative attempts [111][112][113][114][115][116][117][118] to the fully non-linear inclusion [118][119][120][121][122][123] of neutrinos as an extra set of particles. A conservative alternative consists on using exclusively power spectrum measurements within the linear regime (i.e., k < 0.1 h Mpc −1 ). Some of the cosmological constraints have also been obtained using the mildly non-linear regime (k < 0.2 h Mpc −1 ) by means of the so-called Halofit formalism [124,125]. The Halofit prescription models the non-linear matter power spectrum, and it has been calibrated against a wide range of CDM simulations. It has also been extended for massive neutrino cosmologies [126]. Other predictions for the non-linear matter power spectrum include the Coyote emulator [127], which is based on a set of high-accuracy N-body simulations.
However, there is also another avenue to use large scale structure information, the geometrical approach, which exploits the so-called Baryon Acoustic Oscillations (BAO) rather than the measurements of the broad-band shape of the galaxy power spectrum. The BAO signal appears as a peak in the two-point mass correlation function corresponding to the distance a sound wave can travel in the photon-baryon fluid from very early in the universe until the drag epoch, when the baryon optical depth equals one. The BAO signature provides a standard ruler to measure the distance to various redshifts, and it can be measured either along the line of sight, in which the radial distance is inversely proportional to the Hubble expansion rate H(z), or across the line of sight, in which case the angular distance is proportional to an integral of H(z), the angular diameter distance d A (z). To use the BAO method, one must, therefore, extract the acoustic scale from the clustering of some tracer of the baryon distribution (galaxies, quasars). This is typically done statistically using the two-point correlation function of the spatial distribution of tracers, or from its Fourier transform, the power spectrum. From these functions, it is possible to measure two different quantities corresponding to the oscillations parallel and perpendicular to the line of sight, that is r s H(z) and d A (z)/r s , with r s the sound horizon at the drag epoch. Many of the BAO analyses to date have used spherically averaged clustering statistics, measuring an effective distance D V ≡ (zd A (z) 2 /H(z)) 1 3 , which is the volume-averaged distance. Some of the most recent BAO extractions by the Sloan Digital Sky Survey III (SDSS-III) [128] Baryon Oscillation Spectroscopic Survey (BOSS) [129] have achieved, by measuring the clustering of 1.2 million galaxies with redshifts 0.2 < z < 0.75, 1.8% precision on the radial BAO distance and 1.1% precision on the transverse distance in the z < 0.75 redshift region [130][131][132]. These results improve former determinations from previous data releases of BOSS and SDSS [133][134][135][136][137] or other galaxy surveys [138][139][140][141][142], see also the recent works of Ref. [143] for a 2.6% measurement of D V at 2.8σ significance with the extended Baryon Oscillation Spectroscopic Survey (eBOSS) from SDSS-IV [144]. The Dark Energy Survey (DES) has also achieved a 4.4% accuracy on the measurement of d A (z)/r s at z = 0.81 [145].
Galaxy clustering measurements can also be exploited to constrain, at a number of redshifts, the product of the linear growth rate f × σ 8 8 , by means of the so-called redshift space distortions, caused by galaxy peculiar velocities, see the recent analyses of Refs. [146][147][148].
Apart from the spatial distribution of galaxies, there are also other ways of mapping the large scale structure of the universe at different cosmic times. The Lyman-α forest power spectrum from distant quasars plays a major role for constraining the neutrino masses, as it is sensitive to smaller scales, where the effect of neutrino masses is more pronounced. We refer the reader to the seminal works of Refs. [149][150][151][152][153][154][155]. In addition, since the redshifts at which Lyman-α forest probes are sensitive to are higher than those corresponding to galaxy surveys, a fixed scale k will be closer to the linear regime in the Lyman-α case. An additional benefit of going to higher redshifts is that uncertainties related to the evolution of the dark energy fluid will be sub-dominant, as dark energy effects are expected to be more prominent at very low redshifts. However, modeling the neutrino mass effect in the Lyman-α forest power spectrum is highly non-trivial as it may strongly rely on hydrodynamical simulations [155]. These numerical calculations try to properly account for the late time non-linear evolution of the intergalactic medium (IGM), including reionization processes [155,156]. The BAO signature can also be measured in the flux correlation function of the Lyman-α forest of quasars, first detected at a mean redshift z = 2.3 in Ref. [157], see also Refs. [158][159][160][161][162] and Ref. [163], in which joint constraints from the BAO signature from galaxies and quasars have been presented.
Galaxy clusters provide yet another test which allows us to trace the clustering of matter perturbations and, therefore, to test the suppression due to the presence of a non-zero m ν . Galaxy clusters are, by far, the largest virialised objects in the universe, providing a measurement of the so-called cluster number count function dN/dz. This function gives the number of clusters of a certain mass M within a redshift interval (bin) z + δz and, for a given survey: The quantity f sky = ∆Ω/4π refers to the fraction of sky covered by the survey and the unit volume is given by While the redshift is relatively easy to measure, the main uncertainty of this method comes from the cluster mass estimates, determined through four main available methods: X-rays, velocity dispersion, Sunyaev-Zeldovich (SZ) effect 9 and weak lensing. The overall error in the cluster mass determination is usually around ∆M/M ∼ 10%. Furthermore, in order to relate the cluster number count function to the underlying cosmological parameters, one needs as an input a mass function dn(z, M )/dM describing the abundance of virialised objects at a given redshift, usually obtained by means of N -body simulations [164,165]. This mass function depends on both the matter massenergy density and on the standard deviation (computed in linear perturbation theory) of the density perturbations: where P (k) is the matter power spectrum, W (kR) is the top-hat window function, R is the comoving fluctuation size, related to the cluster mass M as R = (3M/4πρ m ) 1/3 , and There are still some degeneracies in the cosmological parameters probed by cluster surveys, whose results are reported by means of a relationship between the matter clustering amplitude σ 8 (obtained from Equation (13)), and the matter mass-energy density Ω m parameters. More concretely, cluster catalogues measure the so-called cluster normalization condition, σ 8 Ω γ m , where γ ∼ 0.4 [166,167]. Current cluster catalogs include X ray clusters (see e.g. [168,169] and references therein), the optically detected SDSS photometric redMaPPer cluster catalog [170][171][172] and the Planck SZ galaxy clusters (PSZ2) [173], which contains more than a thousand confirmed clusters. Other SZ cluster catalogs are those detected from the Atacama Cosmology Telescope (ACT) [168] and from the South Pole Telescope (SPT) [174].
Last but not least, weak lensing surveys are also an additional probe of the large scale structure effects of massive neutrinos [175][176][177][178][179][180][181][182][183]. Light rays from distant galaxies are bent by the matter density perturbations between the source galaxies and the observer, thereby inducing distortions in the observed images of the source galaxies, see the reviews [184,185]. Commonly, the deformations in the source galaxies are rather weak and to extract the lensing signature one needs a correlation among different galaxy images, the so-called shear-correlation functions. By measuring the angular correlation of these distortions, one can probe the clustering statistics of the intervening matter density field along the line of sight, without relying strongly on bias assumptions, setting therefore independent constraints on the neutrino masses. Weak lensing surveys usually report their cosmological constraints in terms of the clustering amplitude σ 8 and the current matter energy density Ω m . More specifically, they make use of the combination S ≡ σ 8 √ Ω m as an accurate description of the amplitude of structure growth in the universe. The most recent weak lensing cosmological analyses profiting of weak lensing data from DES and from the Kilo Degree Survey (KiDS), consisting of ∼ 450 deg 2 of imaging data, are presented in Refs. [186] and [187][188][189], respectively.

C. Cosmological bounds on neutrino masses and their ordering
In the following, we shall review the current cosmological bounds on neutrino masses and on their ordering, firstly in the standard ΛCDM scenario and then when considering extended cosmological models.

Constraints within the ΛCDM universe
Focusing on bounds exclusively from the Planck collaboration, making use of their CMB temperature anisotropies measurements in the multipole range 2500 (Planck TT) and of their low-multipole (up to = 29) polarization data, lowP, a bound of m ν < 0.72 eV at 95% CL [24] is reported. When high-multipole (i.e. small scale, > 30) polarization measurements are included in the analyses (Planck TT,TE,EE + lowP), the quoted constraint is m ν < 0.49 eV at 95% CL. As the Planck TT,TE,EE data combination may still have some systematics due to temperature-to-polarization leakage [24], the bounds including these measurements provide the less conservative approach when exploiting CMB data. In 2016, the Planck collaboration presented a series of new results based on a new analysis, in which the modeling and removal of unexplained systematics in the large angular polarization data were accounted for [42]. The value of the optical depth τ found in these refined analyses (using the SimLow likelihood) was smaller than that quoted in previous analyses [24]: while the lowP data was providing τ = 0.067 ± 0.022, the SimLow likelihood results in τ = 0.055 ± 0.009. The most important consequence of this lower value of τ on the CMB bounds on m ν is related to the degeneracy between the amplitude of the primordial power spectrum, A s , and τ , as already introduced in section IV A: a lower value of τ will imply a lower value of A s , thus implying a lower overall normalization of the spectrum, leading therefore to tighter constraints on neutrino masses. The 95% CL limits of m ν < 0.72 eV and m ν < 0.49 eV, respectively from the Planck TT + lowP and Planck TT,TE,EE + lowP analyses, are updated to m ν < 0.59 eV and m ν < 0.34 eV when using Planck TT + SimLow and Planck TT,TE,EE + SimLow, respectively. These constraints are clearly located away from the region in which a preference for a given mass ordering (normal versus inverted) may show up. Indeed, the CMB data alone were used by the authors of Ref. [44] which, by means of a novel approach to quantify the neutrino mass ordering, have shown that the odds favoring normal ordering versus inverted ordering are 1 : 1 and 9 : 8 in the case of the Planck TT + lowP and Planck TT,TE,EE + lowP data combinations, respectively. These results point to an inconclusive strength of evidence, see Table II. Based on a full Bayesian comparison analysis, Ref. [13] has shown, using Planck TT,TE,EE + lowP measurements together with global neutrino oscillation data, that the Bayes factor for such a combination is log(B NO,IO ) 2.5 for almost all the possible parameterizations and prior choices considered. This value of the Bayes factor, which only points to weak preference for normal ordering, is entirely due to neutrino oscillation data, in agreement with the results of Ref. [46]. Therefore, Planck temperature and polarization measurements alone can not further improve our current knowledge of the neutrino mass ordering from global oscillation data.
The CMB limits on neutrino masses can also include the lensing likelihood, which leads to m ν < 0.59 eV at 95% CL for the case of Planck TT,TE,EE + lowP + lensing measurements [24]. Notice that the bound with the lensing likelihood is less tight than that obtained without the lensing potential extraction ( m ν < 0.49 eV at 95% CL from Planck TT,TE,EE + lowP). The reason is related to the fact that, while the Planck CMB power spectra favor a larger lensing amplitude, the lensing potential reconstructions prefer a lower one. Since increasing the neutrino masses reduces the lensing amplitude, the one dimensional posterior distribution of m ν arising from the combination of CMB temperature, polarization and lensing data sets shifts the neutrino mass constraints away from zero, so that less posterior volume is found near zero than when constraining m ν only with CMB temperature and polarization data.
A significant strengthening on the aforementioned limits can be obtained by means of additional data sets, which help enormously in breaking the degeneracies which are allowed when only CMB data are considered. Among them, the one existing between m ν and the Hubble constant H 0 (see e.g [190]). Large scale structure data from galaxy clustering are of great help in breaking degeneracies. When exploited in the geometrical (BAO) form, the Planck collaboration quotes 95% CL limits of m ν < 0.17 eV from the combination Planck TT,TE,EE + SimLow + lensing + BAO data [42] 10 . Concerning the neutrino mass ordering, the addition of BAO measurements to CMB Planck measurements leads to odds for the normal versus the inverted ordering of 4 : 3 and of 3 : 2, in the case of the Planck TT + lowP + BAO and Planck TT,TE,EE + lowP + BAO respectively, suggesting only very mild evidence for the normal ordering case [44]. These results confirmed the previous findings obtained in Ref. [43]. The authors of Ref. [28] reported odds for the normal versus the inverted ordering of 2.4 : 1 from the combination of Planck TT,TE,EE + BAO plus the SimLow prior on the reionization optical depth, i.e. τ = 0.05 ± 0.009. Notice that if data are not informative enough, the choice of prior on m lightest will make a difference in the odds ratio 11 .
Another possible avenue to exploit galaxy clustering data is to use the information contained in the full-shape of the galaxy power spectrum (see e.g. Refs. [28-30, 154, 192-210]). Notice however that using BAO is currently a more robust method, as the effects of the galaxy bias and non-linearities are not as severe as in the shape approach. In the minimal ΛCDM scenario, the BAO geometrical approach can supersede the neutrino mass constraints obtained from the shape one, see e.g. Refs. [190,204]. Indeed, a dedicated analysis has been devoted in Ref. [28] to compare the constraining power of these two different approaches to large scale structure data with the SDSS-III BOSS measurements. The conclusions are that, even if the latest measurements of the galaxy power spectrum map a large volume of our universe, the geometric approach is still more powerful, at least within the minimal ΛCDM + m ν cosmology. The better performance of BAO measurements is partly due to the upper cutoff applied in the scale k of the power spectrum when dealing with shape analyses (mandatory to avoid non-linearities), and partly due to the fact that two additional nuisance parameters are further required to relate the galaxy power spectrum to the matter power one 12 . As an example, the 95% CL bound of m ν < 0.118 eV obtained with Planck TT,TE,EE + BAO plus SimLow is degraded to m ν < 0.177 eV when replacing part of the BAO data [more concretely, the high redshift BOSS CMASS Data Release 11 (DR11) sample by the full-shape power spectrum measurements from the BOSS CMASS Data Release 12 (DR12)].
An alternative tracer to map out the large scale structure in our universe and improve the CMB-only bounds on the sum of the three active neutrinos is the Lyman-α forest, leading to neutrino mass bounds which turn out to be among the most constraining ones. By means of the one-dimensional Lyman-α forest power spectrum extracted by Ref. [211] and combining these measurements with Planck TT,TE,EE + lowP + BAO, the authors of Ref.
[31] find a 95% CL upper limit of m ν < 0.12 eV. It is also remarkable the fact that, even without the addition of CMB data, the combination of the Lyman-α forest power spectrum of [211], together with those from the XQ-100 quasars at z 3.5 − 4.5 and the high-resolution HIRES/MIKE spectrographs at z = 4.2 and z = 4.6 [212], is already able to provide a limit of m ν < 0.8 eV [213], showing clearly the enormous potential of small-scale probes to extract the neutrino masses.
The degeneracies among m ν and the other cosmological parameters that appear when considering CMB data only can also be strongly alleviated by the addition of Supernova Ia luminosity distance data and/or local measurements of the Hubble parameter 13 . Concerning Supernovae Ia data, the most complete photometric redshift calibrated sample joins the SuperNova Legacy Survey (SNLS) and SDSS supernova catalogs. This Joint Light-Curve Analysis (JLA) catalogue [216][217][218] has been used by the Planck collaboration and by other analyses to improve the constraints on m ν , being its impact particularly crucial in non-minimal cosmologies [27], as we shall explain towards the end of this section. Concerning the value of H 0 , as there exists a strong anti-correlation between the Hubble constant and m ν when considering CMB measurements, larger mean values of H 0 will lead to tighter constraints on the neutrino mass and consequently on the inverted mass ordering. When performing combined analyses of CMB and  [220]], which was in mild (2.5σ) tension with the value of the Hubble parameter derived from 2013 Planck CMB data [101]. This reanalysis [219] considers the original Cepheid data of Ref. [220] and uses a new geometric maser distance estimate to the active galaxy NGC 4258 [221], which is used as a distance anchor to find a value of the Hubble constant H 0 = (70.6 ± 3.3) km s −1 Mpc −1 14 . The limit on the sum of the three active neutrino masses reported by the Planck collaboration using this value of H 0 is m ν < 0.23 eV at 95% CL, when combined with Planck TT + lowP + lensing + BAO + SNIa data. Other estimates of the Hubble constant, however, exist. The 2.4% determination of Ref. [222] profits from new, near-infrared observations of Cepheid variables, and it provides the value H 0 = (73.02 ± 1.79) km s −1 Mpc −1 [222]. As the former mean H 0 value is higher than the one considered by the Planck collaboration, it will lead to tighter limits on m ν . Indeed, the work of Ref. [28] quotes the 95% CL bounds of m ν < 0.196 eV and m ν < 0.132 eV when combining with external data sets using the priors H 0 = (70.6 ± 3.3) km s −1 Mpc −1 and H 0 = (73.02 ± 1.79) km s −1 Mpc −1 , respectively. Focusing on the less conservative choice H 0 = (73.02 ± 1.79) km s −1 Mpc −1 , odds for the normal versus the inverted neutrino mass ordering of 3.3 : 1 were found for both the Planck TT,TE,EE + BAO + SimLow + H 0 and the Planck TT,TE,EE + BAO + SimLow + H 0 + Planck SZ data sets [28]. The 95% CL upper bounds on the neutrino mass for these two combinations are m ν < 0.094 eV and m ν < 0.093 eV, respectively. These results indicate, once again, very mild evidence for the normal mass ordering, even within these more aggressive and less conservative scenarios, in which the very tight limit on m ν is mostly due to the tension between CMB and direct measurements of the Hubble constant H 0 , together with the strong degeneracy between m ν and H 0 . Using these results, we stress that having an upper bound m ν 0.1 eV at 95% CL is not equivalent to having a 95% CL preference for normal ordering: the probabilities for normal ordering and inverted ordering, as computed from the odds 3.3 : 1, are approximately 77% and 23% (see also section V A).
In general, the combination of data sets that are inconsistent is potentially dangerous. Apart from the constraining effect on the neutrino mass limits when considering a particular prior on the Hubble constant H 0 , there have been also other related cases in which the neutrino masses were a tool to accommodate tensions among different data sets. For instance, in the case of galaxy cluster counts, a larger neutrino mass could in principle fit both CMB and low-redshift universe constraints on the power spectrum normalization σ 8 [196]. The effect of combining CMB and BAO observations with clusters and/or shear data is presented in Ref. [223], where it is shown that the inclusion of either cluster or shear measurements in the Planck + BAO joint analysis indicates a preference for m ν > 0 at more than 2σ. However, the authors clearly state that these results can not be interpreted as a claim for a cosmological detection of the neutrino mass, but rather as a remedy to palliate the existing tension between clusters/shear data and Planck/BAO observations. Finally, weak lensing constraints from the Dark Energy Survey Year 1 results [186] (DES Y1), have also recently provided bounds on the sum of the total neutrino mass. Based on 1321 deg 2 imaging data, DES Y1 analyses exploit the galaxy correlation function (from 650.000 luminous red galaxies divided into five redshift bins) and the shear correlation function (from twenty-six million source galaxies from four different redshift bins) as well as the galaxyshear cross-correlation. The 95% CL upper bound reported on the neutrino mass after combining their measurements with Planck TT,TE,EE + lowP + BAO +JLA is m ν < 0.29 eV, ∼ 20% higher than without DES measurements. The reason for this higher value of m ν is that the clustering amplitude in the case of DES Y1 is mildly below the one preferred by Planck measurements. Since larger values of the neutrino mass will decrease the value of the clustering amplitude, the upper limit on the total neutrino mass is loosened by ∼ 20% after the DES results are also considered.

Extensions to the minimal ΛCDM universe
So far we have discussed the neutrino mass and neutrino mass ordering sensitivities within the minimal ΛCDM universe. However, these limits will change when additional parameters are introduced in the analyses.
The first and most obvious scenario one can consider is to test the stability of the neutrino mass limits when new physics is added in the neutrino sector. As already mentioned in section III B, short baseline neutrino experiments indicate that a light sterile neutrino at the eV scale may exist. These extra sterile species will contribute to the effective number of relativistic degrees of freedom, N eff , defined by where ρ rad (ρ γ ) is the total radiation (CMB photons) energy density. In the standard picture N eff = 3.046 [25,26]. This number accounts for the three active neutrino contribution and considers effects related to non-instantaneous neutrino decoupling and QED finite temperature corrections to the plasma evolution 15 . Variations in N eff , apart from the light sterile neutrino, may be related to the existence of additional relativistic particles, as thermally-produced axions (see below). Analyses in which both the active neutrino masses and the number of additional massless or massive species are varied simultaneously have been extensively carried out in the literature [190,[224][225][226][227][228][229][230][231], showing that the bounds on the active neutrino mass are relaxed when additional sterile species are added to the fermion content of the SM of particle physics. The constraints on the total neutrino mass m ν are less stringent than in the standard three neutrino case due to the large degeneracy between m ν and N eff , which arises from the fact that a number of massless or sub-eV sterile neutrino species contributing to the radiation content of the universe will shift both the matter-radiation equality era and the location of the CMB acoustic peaks. This effect could be compensated by enlarging the matter content of the universe, implying therefore that larger values for the neutrino masses could be allowed. Consequently, a priori, the constraints on m ν when N eff is also a free parameter in the analyses are not very competitive. Fortunately, CMB measurements from the Planck collaboration help enormously in sharpening the measurement of N eff , especially when considering polarization measurements at small scales: including data at high multipoles, one obtains ∆N eff < 1 at more than 4σ significance from Planck CMB observations alone. Indeed, the limit on the sum of the three active neutrino species considering also additional radiation neutrino species (i.e. massless sterile neutrino species) is m ν < 0.178 eV at 95% CL from Planck TT,TE,EE + lowP + BAO data, very similar to the bound m ν < 0.168 eV at 95% CL arising from the very same dataset within the minimal ΛCDM scenario with three active massive neutrinos. Another possible way of relaxing (or even avoiding) the cosmological neutrino mass limits is via the addition of non-standard interactions in the active neutrino sector [232][233][234][235][236][237][238][239][240][241][242][243][244][245].
Furthermore, additional relics different from sterile neutrinos, as thermal axions [246][247][248][249], contributing to both N eff at early times and to the hot dark matter component in the late-time universe, suppress small-scale structure formation and show effects very similar to those induced by the (active) three massive neutrino species. Therefore, the cosmological bounds on the three active neutrino masses are modified in scenarios with thermal axions, see Refs. [250][251][252][253][254][255][256][257][258], as these two species have to share the allowed amount of dark matter. Nonetheless, there are non-negligible differences among neutrinos and thermal axions: (a) axions are colder than neutrinos, as they decouple earlier; (b) since the axion is a scalar particle, an axion mass larger than the neutrino one is required in order to make identical contributions to the current mass-energy density of the universe; (c) in the case of axions, the contribution to N eff is related to their mass, while for neutrinos this is usually not true. Consequently, the bounds on the axion mass are always less constraining than for the neutrino, and m ν is slightly more constrained in scenarios in which thermal axions are also present. For instance, Ref. [32] quotes m ν < 0.175 eV at 95% CL from the Planck TT,TE,EE + lowP + BAO data combinations when considering only neutrinos, while the analyses in Ref. [257], including massive axions, find m ν < 0.159 eV and m a < 0.763 eV, both at 95% CL, for the very same data combination. There are also other ways of relaxing the cosmological neutrino mass bounds, related either to the early or the late-time accelerating periods in the universe. In the former case one can play with inflationary processes. There have been a number of studies devoted to explore their degeneracies with the neutrino sector, see the recent works of Refs. [201,210,231,[258][259][260][261]. The authors of Ref. [258] have considered a non-standard and parametric form for the primordial power spectrum, parameterized with the PCHIP (piecewise cubic Hermite interpolating polynomial) formalism with twelve nodes between k 1 = 5 × 10 −6 Mpc −1 and k 2 = 10 Mpc −1 and derived the neutrino mass constraints within this more general scenario. When only Planck TT + lowP measurements were considered, the 95% CL mass bound of m ν < 0.75 eV obtained with the usual power-law description of the primordial power spectrum was relaxed to m ν < 2.16 eV. This large value is explained in terms of the strong degeneracy between m ν and the PCHIP nodes corresponding to the wave-numbers where the contribution of the Early ISW effect is located, in such a way that the effect induced by a non-zero neutrino mass is easily compensated by an increase of the primordial power spectrum at these scales only. BAO information improves considerably the limits in the PCHIP prescription, but it is the addition of high-polarization data what further constrains the effect. The 95% CL upper limit in the PCHIP scenario from the Planck TT,TE,EE + lowP + BAO data combination is m ν < 0.218 eV, quite close to the bound found when the usual power-law description is applied ( m ν < 0.175 eV). Reference [261] deals instead with the robustness of the constraints on the scalar spectral index n s under several neutrino physics scenarios. The authors have explored the shifts induced in the inflationary parameters for different choices of the neutrino mass ordering, comparing the approximate massive neutrino case (one massive eigenstate plus two massless species when the total mass is close to the minimum allowed value by oscillation data, and three degenerate massive neutrinos otherwise) versus the exact case (normal or inverted mass orderings). While the mass-ordering assumptions are not very significant when m ν is fixed to its minimum value, there is a shift in n s when m ν is a free parameter, inherited from the strong degeneracies in the m ν , H 0 and Ω m h 2 parameter space. Fortunately, BAO measurements revert the m ν -n s anti-correlation present with CMB data only, and the shift in the spectral index turns out to be negligible.
The other possibility is to play with the late-time acceleration period and study how the neutrino mass bounds change. The current accelerated expansion of the universe, explained in terms of a cosmological constant in the minimal ΛCDM scenario, may be due to a dynamical dark energy fluid with a constant equation of state w = −1 or a time-dependent w(z) [262,263], or to quintessence models, based on the existence of a cosmic scalar field [264][265][266][267][268][269], which provide a dynamical alternative to the cosmological constant scenario with w = −1. It is naturally expected that the neutrino mass bounds will increase when enlarging the parameter space. Indeed, when the dark energy equation of state is allowed to vary within the phantom region w < −1, there is a very well-known degeneracy between the dark energy equation of state w and the sum of the three active neutrino masses, as first noticed in Ref. [270] (see also Refs. [27,231,259,[271][272][273]) 16 . It has been pointed out that for very high neutrino masses only dark energy models lying within the phantom region will be allowed. The reason for that is the following: a larger m ν can be compensated by a larger Ω m , which in turn can be compensated by a smaller equation of state of the dark energy component, i.e. w < −1. Interestingly, the recent work of Ref. [27] shows that the cosmological bounds on m ν become more restrictive in the case of a dynamical dark energy component with w(z) ≥ −1. Following the usual dynamical dark energy description, whose redshift dependence is described by the standard Chevallier-Polarski-Linder (CPL) parametrization [262,263], the authors of [27] have shown that the combination of Planck TT,TE,EE + BAO + JLA plus the SimLow prior on the reionization optical depth provides, at 95% CL, m ν < 0.11 eV in the CPL case when restricting w(z) ≥ −1 (within the physical, non-phantom region), while m ν < 0.13 eV in the ΛCDM case. When w(z) is also allowed to be in the phantom region (w(z) < −1) within the CPL parameterization, the resulting 95% CL constraint on the three active neutrino masses is m ν < 0.37 eV. These results have a direct impact on the cosmological preference for a given neutrino mass ordering. Following Refs. [28,43], it is found that the normal ordering is mildly preferred over the inverted one, with posterior odds 3 : 1 for the data combination quoted above when w(z) ≥ −1. On the contrary, if there is no such a restriction and w(z) can also take values in the phantom region, the odds are 1 : 1. The odds in the non-phantom dynamical dark energy case show a mild preference for normal ordering. Therefore, if neutrino oscillation experiments or neutrinoless double beta decay searches find that the neutrino mass ordering is the inverted one, if the current cosmic acceleration is due to a dynamical dark energy component, one would require this component to be phantom.
As a final point in this section, we would like to note that also in scenarios in which the current accelerated expansion is explained by means of modifications of gravity at ultra-large length scales, the cosmological limits on neutrino masses will differ from those in the standard ΛCDM model, see e.g. [279][280][281][282][283][284][285][286][287].

V. GLOBAL 2018 DATA ANALYSIS
In this section we shall combine the available measurements that allow us to constrain the neutrino mass ordering, updating the results presented in Ref. [13].

A. Bayesian model comparison
Before performing the analysis, we will briefly summarize the method we will adopt to compare the two possible orderings. We will follow a Bayesian approach to model comparison, which makes use of the Bayesian evidence Z. This quantity, which is also known as the marginal likelihood, is defined as the integral over the entire parameter space Ω M of the prior π(θ) ≡ p(θ|M) times the likelihood L(θ) ≡ p(d|θ, M), where θ is the set of parameters that describe the model M and d represents the available data: The posterior probability of the model M can be written in terms of its prior probability π(M) times the Bayesian evidence Z M : where the proportionality constant depends only on the data. In our case we will be interested in comparing normal ordering (NO) and inverted ordering (IO), which can be considered as two different competing models M 1 ≡ NO and M 2 ≡ IO. The ratio of the posterior probabilities of the two models can be written as having defined the Bayes factor as Assuming the same prior probabilities for normal and inverted ordering, the Bayes factor is what determines the odds in favor of one of the competing models. In particular we will indicate the results in terms of its natural logarithm ln B NO,IO , which will be positive when data will prefer normal ordering and negative otherwise. Quantitatively, the preference is given in terms of posterior odds, which are always |B NO,IO | : 1 in favor of the preferred model. The strength of the preference can be also translated into an empirical scale, which in our case is summarized in the third column of Table II. Let us briefly discuss the correspondence of the quoted levels that classify the strength of the preference in favor of one of the competing models. In the case of the neutrino mass ordering, we have only two possibilities (normal or  inverted), so that p(NO|d) + p(IO|d) = π(NO) + π(IO) = 1. If we assign the same prior probability to the two cases, π(NO) = π(IO) = 1/2, it is easy to compute the posterior probability for each of the two cases, which will be having used Equations (18) and (19). The confidence levels for the rejection of the disfavored (e.g. inverted) mass ordering will then be x = 100 × (1 − |B NO,IO | −1 ) %. For example, a Bayes factor B NO,IO = 10 corresponds to a rejection of the inverted ordering at 90% CL. If, instead, we want to reproduce the probability levels P = erf(N/ √ 2) that are usually associated to the classical N σ levels for a Gaussian measurement, being erf the error function and considering, for example, N ∈ (1, 2, 3, 4, 5), the corresponding Bayes factors B can be computed to be B = P/(1 − P ), which gives us ln B N σ 0.77, 3, 5.9, 9.7, 14.37. Therefore, our "strong", "very strong" and "decisive" levels roughly correspond to the > 3σ, > 4σ and > 5σ probabilities, as indicated in the fourth column of Table II.

B. Parameterization and data
Our two competing models are described by the same number of parameters, listed with their priors in Table III: the three neutrino mixing angles (sin 2 θ 12 , sin 2 θ 13 , sin 2 θ 23 ), the CP violating phase δ CP and the parameters associated with neutrino masses, neutrinoless double beta decay (0νββ) and cosmology, as we shall describe now.
We consider in our analysis the parameterization that uses the two mass splittings (∆m 2 21 and ∆m 2 31 ) and the lightest neutrino mass m lightest with logarithmic priors. This parameterization, strongly motivated by the physical observables, was shown to provide the optimal strategy to successfully explore the neutrino parameter space, see Ref. [13] 17 . Within the other possible choice, that is, within the parametrization that uses the three neutrino masses as free parameters, most of the parameter space at high neutrino masses is useless for the data fit. Therefore, this second possibility is penalized by the Occam's razor and we shall not explore it here.
The neutrino mixing parameters are constrained using the same data we described in section II. The complete oscillation data set is indicated with the label "OSC" in the following.
For the cosmological part, we will describe the universe using the ΛCDM model and its six parameters: the baryon and cold dark matter densities, Ω b h 2 and Ω c h 2 ; the optical depth to reionization, τ ; the angular scale of the acoustic peaks through Θ s and the amplitude log(10 10 A s ) and tilt n s of the power spectrum of initial curvature perturbations. In addition, we add the effect of the three massive neutrinos computing the evolution of the cosmological observables assuming three independent mass eigenstates, which, in terms of the parameters involved in our analyses, read as m 1 = m lightest m 1 = m 2 lightest + |∆m 2 31 | , m 2 = m 2 lightest + ∆m 2 21 m 2 = m 2 lightest + |∆m 2 31 | + ∆m 2 21 and m 3 = m 2 lightest + ∆m 2 31 (m 3 = m lightest ) for normal (inverted) neutrino mass orderings. When considering cosmological data, we will focus on the Planck measurements of the CMB spectrum and on the most recent results from BAO observations. For the former we consider the 2015 Planck release [24,290] of the highlikelihood [291], together with a prior on τ as obtained in the 2016 intermediate results [42]. For the purposes of our analyses, this will be sufficient to mimic the final Planck release which is expected within the next few months. Complementary to the CMB, we include in our calculations the final constraints from the SDSS BOSS experiment, the DR12 release, in the form denoted as "final consensus" [292], which provides constraints from observing 1.2 million massive galaxies in three separate bands at effective redshifts 0.38, 0.51 and 0.61, plus results from the 6DF survey at z = 0.106 [142] and from the SDSS DR7 MGS survey at z = 0.15 [136]. The combined dataset including the mentioned CMB and BAO data will be denoted as "Cosmo".
In addition, we shall impose a prior on the Hubble parameter as obtained in the recent results from Ref. [222]: H 0 = (73.24 ± 1.74) km s −1 Mpc −1 . We will denote the data combinations including this prior with the label "H 0 ".
Finally, concerning neutrinoless double beta decay, we vary the two Majorana phases in the entire available range (0-2π) and the NMEs according to the range allowed by recent theoretical calculations. We revised the NME ranges adopted in [13], which were the ones suggested in [293]. Here we use these new ranges: [3.3 -5.7] for 76 Ge and [1.5 -3.7] for 136 Xe, following the 1σ range proposed in [81].
We use 0νββ data from the 136 Xe experiments KamLAND-Zen [39] and EXO-200 [77] and from the 76 Ge experiment Gerda, for which we use the results in Ref. [294], since the latest publication [40] does not contain enough information that allows us to parameterize a likelihood function. The most stringent bounds, anyways, still come from KamLAND-Zen, so that not including the new Gerda results does not affect significantly our results. For the very same reason we do not include the results of CUORE [41], for which the uncertainty on the NME of 130 Te is very large and the constraints corresponding to most of the values of M 0ν 130 Te are much looser than the ones from KamLAND-Zen, and of CUPID-0 [79], which establishes a much less stringent limit on the 82 Se half-life. The complete neutrinoless double beta set of data will be denoted as "0νββ".
All the previously listed data are coded as likelihood terms in a full Bayesian analysis. We compute the cosmological quantities using the Boltzmann solver CAMB [295], the likelihoods using the interface provided by CosmoMC [296], modified in order to take into account the oscillation and neutrinoless double beta decay data, while the calculation of the Bayesian evidence is committed to PolyChord [297,298].

C. Constraints on the mass orderings
The main results are depicted in Figure 10. The first data point corresponds to the Bayesian evidence from oscillation data only. Notice that the Bayes factor [ln(B NO,IO ) = 6.5 ± 0.2 for concreteness] indicates strong evidence for the normal mass ordering from oscillation data only. This Bayes factor is translated into a ∼ 3.2σ evidence favoring normal mass ordering. This result was expected in light of the results presented in section II, arising from the frequentist joint analysis. There it was reported a ∆χ 2 = 11.7 in favor of the normal mass ordering from the combination of all long baseline, reactor and atmospheric data, which corresponds, roughly, to ∼ 3.4σ. Adding information from neutrinoless double beta decay searches does not affect the Bayesian analysis, as shown by the second data point in Figure 10, and as expected from previous work [13].
Once CMB and BAO measurements are also added in the Bayesian analysis, ln(B NO,IO ) = 7.4 ± 0.3 is obtained (see the third point in Figure 10), improving the significance of the preference for normal ordering from ∼ 3.2σ to ∼ 3.4σ. Notice that, even if the preference for the normal neutrino mass ordering is mostly driven by oscillation data, the information provided by cosmological observations is more powerful than that in the analysis carried out in Ref. [13], as the Bayesian analyses here also include BAO measurements, together with CMB data. Indeed, from the two Bayes factors obtained considering oscillation data only [ln(B NO,IO ) = 6.5 ± 0.2] and oscillation plus cosmological measurements [ln(B NO,IO ) = 7.4 ± 0.3], it is straightforward to infer the probability odds for normal ordering arising exclusively from cosmology. By doing so, one obtains odds of 2.7 : 1 for the normal ordering against the inverted one, in perfect agreement with the analyses of Ref. [28], where odds of 2.4 : 1 with cosmological data only were reported when considering the very same data sets adopted here (albeit the odds were derived with an alternative method).
Finally, the addition of the prior on the Hubble constant raises the evidence for the normal ordering to ln(B NO,IO ) = 7.7 ± 0.3 (i.e. ∼ 3.5σ). This improvement is expected, as previously explained in section IV, since a prior on the Hubble constant breaks the degeneracy between m ν and H 0 and, therefore, sharpens the neutrino mass bounds from cosmology. By performing a similar exercise to the one previously quoted, one finds that the odds for normal versus inverted ordering from cosmology data only are 3.3 : 1 for the combination of CMB, BAO plus the H 0 prior, again in excellent agreement with the results obtained in Ref. [28].

A. Prospects from oscillations
As we have seen in section II, the combination of all current neutrino experiments leads to a preference for normal ordering of 3.4σ, within the context of the latest frequentists global data analyses. The Bayesian analysis described in the previous section confirms these results, as we have reported a 3.2σ evidence for normal mass ordering. In principle, one may expect to achieve further sensitivity on the neutrino mass ordering from more precise data by the current long-baseline and atmospheric neutrino experiments. However, indications above the 3σ level from a single experiment seem very unlikely. This will be one of the main goals of the next-generation neutrino oscillation experiments, including new long-baseline, reactor, and atmospheric neutrino detectors. These upcoming experiments will be able to measure the neutrino mass ordering with astonishing precision. In this section we briefly discuss some of the proposed projects and their physics potential.

Long-baseline experiments
The Deep Underground Neutrino Experiment (DUNE) [299][300][301][302] will be a new long-baseline accelerator experiment, with a small near detector and a huge far detector with a fiducial mass of 40 kton located 1300 km away from the neutrino source at Fermilab. With its powerful 1.1 MW beam, it will be exposed to around 15 × 10 20 POTs (protons on target) per year, which will lead to a huge number of events and therefore to high precision measurements of the neutrino oscillation parameters. As explained in section II, the presence of matter affects differently the neutrino appearance probabilities for normal and inverted mass orderings. DUNE, with the longest baseline ever for an accelerator neutrino experiment, will be able to measure the neutrino mass ordering with a significance above 5σ for any set of the oscillation parameters (θ 23 , δ CP ) after 7 years of data taking. Note that this sensitivity could be further increased by using an improved energy reconstruction method, as shown in Ref. [303]. On the other hand, the sensitivities could also be biased by the potential presence of new physics beyond the SM, such as non-standard neutrino interactions [304][305][306][307], deviations from unitarity [308][309][310] or the presence of light-sterile neutrinos [311]. Indeed, besides providing very precise information about the neutrino oscillation mechanism, the DUNE experiment will also be very useful to test different models for neutrino masses and mixings [312][313][314][315][316] as well as to check for various effects of new physics such as the ones mentioned above, neutrino decay scenarios [317,318] or even CPT invariance [319,320].
There are also plans to build a larger version of the Super-Kamiokande detector, Hyper-Kamiokande [321], that will be very similar to its predecessor but with a fiducial mass of 560 kton, 25 times larger than Super-Kamiokande. The Hyper-Kamiokande detector will be a requirement for the upgrade of T2K, the T2HK (Tokai-to-Hyper-Kamiokande) experiment [322]. The very massive detector together with the upgraded neutrino beam from J-PARC will guarantee a huge number of neutrino events and therefore larger statistics. As a consequence, T2HK will be able to determine the neutrino mass ordering after few years of running time with very high significance, as well as to explore new physics scenarios, see for instance Refs. [321,323,324]. In combination with atmospheric data from Hyper-Kamiokande, a 3σ rejection of the wrong mass ordering would be expected after five years of data taking. A third project has been proposed as an extension of T2HK to Korea, the T2HKK (Tokai-to-Hyper-Kamiokande-and-Korea) experiment [325]. This proposal includes a second far detector facility for the J-PARC neutrino beam, located at 1000-1300 km from the source. The longer path traveled within the Earth by the neutrinos detected in T2HKK will result in an enhanced sensitivity to the neutrino mass ordering if compared to T2HK alone.
The synergies and complementarities among the three long-baseline proposals above, DUNE, T2HK and T2HKK, have been discussed in Ref. [326]. It is found that the combination of their experimental results may significantly mitigate the limitations of a given experiment, improving the precision in both the determination of the mass ordering and the measurement of CP violation.
Note that, although here we have focused on the long-baseline side of DUNE and Hyper-Kamiokande, they are actually designed as multi-purpose experiments, with a rich physics program aiming to study the neutrino oscillations with accelerator, atmospheric and solar neutrinos as well as to detect neutrinos from astrophysical sources and proton decay.

Atmospheric experiments
In atmospheric neutrino experiments, the sensitivity to the mass ordering comes from the matter effects that distort the pattern of neutrino oscillations inside Earth, see Equation (4). Based on the oscillatory pattern that depends on the reconstructed neutrino energy and zenith angle, an ideal experiment would observe a given number of events in each energy and zenith angle bin as shown in Figure 11. Comparing the observed two-dimensional histograms with the theoretical ones for normal (left panel) or inverted ordering (right panel) allows to determine the true mass ordering that is realized in nature. In the following we list some of the future projects with this aim.
The Oscillation Research with Cosmics in the Abyss (ORCA) experiment [327] will be a large neutrino telescope placed deep inside the Mediterranean sea. It will detect the Cherenkov light emitted by the muons and electrons created by the interactions of atmospheric neutrinos in the sea and that propagate into water. Unlike its precursor, ANTARES, with 12 lines and a separation of 70 meters between neighbouring optical modules, ORCA will have 60 lines with modules separated by 9 meters. Due to the matter effects on the propagation of atmospheric neutrinos, the ORCA experiment will be able to measure the neutrino mass ordering with very good precision. In particular, a 3σ determination of the mass ordering can be expected after only three years of data taking, with even higher significance for the case in which nature has chosen normal ordering and the upper octant for the atmospheric mixing angle. Several studies have been performed in order to analyze the sensitivity of ORCA to the standard oscillation parameters [328,329]. Its potential to determine the Earth matter density through neutrino oscillation tomography [330] or to test new physics scenarios [331,332] have also been extensively discussed. PINGU (Precision IceCube Next Generation Upgrade) [333] is a planned upgrade of the IceCube DeepCore detector, an ice-Cherenkov neutrino telescope which uses the antarctic ice as a detection medium. The IceCube design aims at the detection of very high energy neutrinos, with an energy threshold above the relevant energy range for neutrino oscillations. However, the denser instrumented region DeepCore allows IceCube to decrease its energy threshold down to E th = 6.3 GeV. A further improvement with an even denser zone, PINGU, could lower E th to only a few GeV. With this very low-energy threshold, one of the main purposes of PINGU is the determination of the neutrino mass ordering, with expected sensitivities similar to the ORCA experiment 18 . Besides that, PINGU is expected to have the best sensitivity to ν τ appearance and to determine accurately the octant of the atmospheric mixing angle. The PINGU capabilities to detect high-energy supernova neutrinos [335], and to investigate scenarios beyond the Standard Model, such as non-standard interactions [336] or dark matter self-interactions [337,338] have been also analyzed in the literature.
The India-based Neutrino Observatory (INO) is a very ambitious project, aiming to detect atmospheric neutrinos with a 50 kton magnetized iron calorimeter (ICAL) [339]. The most outstanding feature of the INO experiment will be its capability to distinguish neutrinos from antineutrinos in an event by event basis. As a result, the identification of the matter effects in the neutrino and antineutrino propagation will be much cleaner in comparison with the sea water/ice Cherenkov detectors. Indeed, one of the main scientific goals of INO will be the determination of the neutrino mass ordering [340]. According to the Physics White Paper of the ICAL (INO) Collaboration [339], after 10 years run, INO will be able to identify the correct neutrino mass ordering with a significance larger than 3σ. As the experiments discussed above, the atmospheric neutrino results from INO can also be used to test the presence of new physics beyond the SM, such as CPT-or Lorentz violation [341], sterile neutrinos [342,343], dark matter related studies [344,345], non-standard neutrino interactions [346] or decaying neutrinos [347].

Medium-baseline reactor experiments
We have focused so far on extracting the neutrino mass ordering from matter effects in the neutrino propagation through the Earth interior. An alternative technique is that provided by medium-baseline reactor neutrino experiments. For baselines of the order of 50 km, the survival probability for reactor antineutrinos exhibits a pattern that may allow the discrimination between normal and inverted mass orderings. Indeed, for such distances, the electron antineutrino survival probability is given by the following expression: P νe→νe = 1 − cos 4 θ 13 sin 2 2θ 12 sin 2 ∆ 21 − sin 2 2θ 13 sin 2 ∆ 31 + sin 2 θ 12 sin 2 ∆ 21 cos 2∆ 31 ∓ sin 2 θ 12 2 sin 2∆ 21 sin 2|∆ 31 | , 18 The effect of statistic and systematic uncertainties on the PINGU sensitivity to the mass ordering has been presented in Ref. [334]. where ∆ ij = ∆m 2 ij L 4E and the minus (plus) sign in the last term corresponds to normal (inverted) mass ordering. This probability contains a main oscillatory term with a frequency given by the solar neutrino mass splitting ∆m 2 21 , plus an additional term whose frequency depends on the sign of the atmospheric splitting ∆m 2 31 , i.e. on the neutrino mass ordering. The effect of the ordering over the neutrino survival probability in a medium-baseline reactor experiment is illustrated in Figure 12. There, we depict in black the oscillatory term corresponding to the solar splitting frequency. The red (blue) line corresponds to the full neutrino survival probability for normal (inverted) mass ordering. Note that this plot was obtained using the best-fit values from Table I for each ordering.
The Jiangmen Underground Neutrino Observatory (JUNO) [348] is a 20 kton multi-purpose underground liquid scintillator detector. The site of the experiment, located 53 km away from the Yangjiang and Taishan nuclear power plants in China, was chosen to optimize its sensitivity to the neutrino mass ordering, one of its main physics goals. Like any other reactor neutrino experiment, JUNO will be sensitive to the disappearance of electron antineutrinos, with about 10 5 events expected after six years of run time. From this high-statistics data sample, JUNO will try to reconstruct with extremely good precision the neutrino oscillation spectrum and to discriminate the different highfrequency behavior for normal and inverted mass ordering, as illustrated in Equation (22) and Figure 12. For a projected energy resolution of 3% at 1 MeV, JUNO will be able to establish the neutrino mass ordering at the level of 3-4σ in 6 years. Apart from the mass ordering, JUNO will also provide precision measurements of the solar oscillation parameters, θ 12 and ∆m 2 21 , with an accuracy of around 1%. In this sense, JUNO might help to solve the observed disagreement between the mass splitting measured at solar experiments and at the reactor experiment KamLAND. If the discrepancy persists after new measurements by JUNO and future solar results by Super-Kamiokande, it could be considered as an indication of new physics [304]. Moreover, JUNO will be sensitive to different types of new physics scenarios beyond the SM, as studied in Refs. [349][350][351][352][353][354][355][356][357].
In parallel to JUNO, there is a proposal to extend the already existing experiment RENO with a third mediumbaseline detector located at a distance of 47 km. This new project is known as RENO-50 [358], given its location, at approximately 50 km from the Hanbit power plant, in South Korea. The detector would consist of a 18 kton ultralow-radioactive liquid scintillator instrumented with 15000 high quantum efficiency photomultiplier tubes. Using the same technique described above, RENO-50 will be able to determine the neutrino mass ordering as well as the solar oscillation parameters with extremely good precision. Conceived as multi-purpose detectors, JUNO and RENO-50 will have a wide physics program, including not only the observation of reactor and solar neutrinos, but also neutrinos from supernova bursts, the diffuse supernova neutrino background, atmospheric neutrinos and geoneutrinos.

B. Prospects from beta-decay experiments
As already mentioned in section III, the determination of the mass ordering through the observation of the energy spectrum near the endpoint of β-decay or similar will be extremely challenging, because an impressive energy resolution is required to distinguish the kink due to the second and third mass eigenstates in the spectrum. We list here the main projects that aim at detecting the neutrino mass in the future and comment their perspectives for the mass ordering determination.
The first experiment we will comment on is KATRIN, which has recently started operations and aims at a detection of the effective electron antineutrino mass with a sensitivity of 0.2 eV [74,75]. The first results from KATRIN are expected in early 2019, but the final target statistics will be reached after 3 yr of data taking. Thanks to the detailed study of the detector systematics which has been carried out, it is possible that the final mass determination will reach a better sensitivity than the nominal one of 0.2 eV, eventually reaching something closer to 0.1 eV [359]. Even with the more optimistic sensitivity, however, it will be impossible for KATRIN to determine the mass ordering.
Other tritium experiments exploiting different technologies include the Project-8 [360][361][362] experiment, which will use the Cyclotron Radiation Emission Spectroscopy (CRES) in order to determine the mass of the electron antineutrino. The technique consists in measuring the frequency of cyclotron radiation emitted by the electrons released during tritium decay and spiralling into a magnetic field. The frequency can then be related with the electron energy and consequently the energy spectrum can be determined. At the moment, Project-8 is in the calibration phase (phase-II) [363] for a small prototype which will not have enough sensitivity to be competitive in the determination of the neutrino mass. Next phases include a large volume system using molecular tritium (phase-III), starting in 2020, which will be competitive in determining the neutrino mass and will serve as an intermediate step before moving to phase-IV, which will use atomic tritium, required in order to avoid uncertainties related to the existence of excited molecular tritium states. Project-8 in its atomic tritium phase is expected to reach the sensitivity mν e 40 meV with an exposure of 10 − 100 m 3 yr, sufficient to probe the values of mν e allowed in the context of inverted ordering [362], so that in case of no observation we will know that the ordering of neutrino masses must be normal.
Another interesting class of the experiments includes the HOLMES [364,365] and ECHo [366] experiments, which both aim at the determination of the electron neutrino mass through observations of the endpoint of the electron capture decay of 163 Ho, which practically proceeds through the measurement of de-excitation transitions of the Dy atoms, which are produced in the process 163 Ho + e − → 163 Dy * + ν e [367]. As for the tritium β-decay, also the endpoint of the 163 Ho electron capture spectrum depends on the value of the neutrino masses and, in principle, it would be possible to determine the mass ordering in this way. Besides the experimental and theoretical problems that the HOLMES and ECHo collaborations must face, however, it seems that the current technology is not yet at the level of precision required for the mass ordering determination. The HOLMES demonstrator, currently running, should reach a sensitivity of m νe 10 eV by the end of 2018, while the full-scale experiment, possibly starting in 2019, has a target sensitivity of m νe 1 eV [368]. ECHo, on the other hand, is running a first phase (ECHo-1k) which has also a target of m νe 10 eV in 1 yr, while the full scale ECHo-100k will reach m νe 1.5 eV in 3 yr of data taking, expected to start in 2019 [368]. Both results are impressive when compared with the current upper limit on the electron neutrino mass using the same isotope, which is 225 eV [369], more than two orders of magnitude larger.
Finally, to conclude this subsection we want to mention that the PTOLEMY proposal [370,371], aiming at the detection of the relic neutrino background and recently approved by the Scientific Committee of the Laboratori Nazionali del Gran Sasso (LNGS), will be able to study and possibly determine the mass ordering through the observation of the β-spectrum of tritium decay. PTOLEMY will be discussed later in section VI G.

C. Prospects from neutrinoless double beta decay
We list here the future perspectives for neutrinoless double beta decay experiments in terms of sensitivity to the half-life for the processes of interest. As we already commented in section III B, the conversion between the half-life T 0ν 1/2 and the effective Majorana mass m ββ depends on the NME and the phase space factor of the process of interest, see Equation (7). In order to exclude the inverted ordering allowed range for m ββ (in case there is no sterile neutrino), one would need to constrain m ββ 10 meV, which corresponds to T 0ν 1/2 1 × 10 28 yr, with some dependence on the material (phase space and NME). This means that none of the current generation experiments will be able to reach the required sensitivity, and we will have to wait for next-generation upgrades and new projects. Many of the information listed in the following has been taken from Refs. [38,372].

Current generation experiments
Let us firstly address the current generation of experiments, which at most will be able to start exploring the three-neutrino inverted mass ordering regime or to probe the upper range for m ββ allowed within the 3+1 neutrino scenario. The experiments will be listed in alphabetical order.
AMoRE [373] is an experiment devoted to determine the life-time of 100 Mo. After a first pilot run, the current status (AMoRE-I) is to test the technology with a 100 Mo mass of 5-6 kg, in order to demonstrate the scalability before moving to the full scale (AMoRE-II) detector, which will use 200 kg of material and is expected to start around 2020, with a final target sensitivity of T 0ν 1/2 5 × 10 26 yr.
CUORE [374][375][376], already mentioned in section III B, works with 130 Te and is already taking data with the full scale detector, which will have as ultimate sensitivity T 0ν 1/2 9 × 10 25 yr after 5 yr of data taking [377]. The KamLAND-Zen experiment [39,378], after the previous successful data taking period, is now upgrading the detector for a new observation run with approximately 750 kg of 136 Xe and a new balloon inside the KamLAND detector. The target sensitivity for the upcoming phase is around T 0ν 1/2 5 × 10 26 yr, a factor of five larger than the current limit [39].
A smaller experiment is NEXT [379], which is running background studies in the Canfranc laboratories in Spain. NEXT will use high pressure 136 Xe TPCs, which will allow an impressive tracking of the emitted particles through scintillation and electroluminescence. A prototype with 10 kg of natural Xenon will start data taking this year to demonstrate that the expected background control and particle tracking have been achieved. In 2019 NEXT is expected to start a new phase with 100 kg of 136 Xe, which will reach T 0ν 1/2 1 × 10 26 yr with 5 yr of data. A similar project is called Panda-X-III [380], which is based in the Jinping underground laboratories in China. Panda-X-III will run the first phase using 200kg of 136 Xe to reach T 0ν 1/2 1 × 10 26 yr in 3 yr. Going to a different concept, SNO+ [381] will feature a detector of 760 ton of ultra-pure liquid scintillator. SNO+ will be a multipurpose detector, as it will be capable of studying reactor, solar, supernova and geoneutrinos, and also to probe proton decay [382]. After the background studies will be completed, a 0.5% loading will be performed, inserting 130 Te in the detector to measure double beta decay processes. The target sensitivity after 5 yr is T 0ν 1/2 2 × 10 26 yr. Future plans for the SNO+ experiment include the further 130 Te loading to 1%, or even more, of the detector mass, with the advantage that increasing the 130 Te amount will not influence the backgrounds but only the signal. The final target for this second phase is to reach T 0ν 1/2 1 × 10 27 yr, thus starting to cover the inverted ordering allowed range. Let us finally comment the SuperNEMO experiment [383,384], which uses 82 Se for its study. SuperNEMO is particularly interesting because it will be able to perform a full topological reconstruction of the events, which is extremely important in case of detection because it opens the possibility to directly test the mechanism that underlies neutrinoless double beta decay and, in principle, to determine the lepton-number violating process. A first demonstrator of about 7 kg is expected to start in 2018 and to reach T 0ν 1/2 6 × 10 24 yr with 2.5 yr of data. The subsequent plans include an extension with a ∼ 100 kg scale detector with 20 modules, which will be able to probe T 0ν 1/2 up to 1 × 10 26 yr, and the possibility to use the 150 Nd isotope, for two reasons: to have a more favorable phase space when converting T 0ν 1/2 to m ββ and to get rid of the Rn background which affects the 82 Se measurements [372].
As a summary, some of the current generation experiments will be able to probe the inverted ordering range of m ββ within the standard three neutrino framework and assuming an exchange of light Majorana neutrinos. However, none of them will be able to rule out completely the inverted mass ordering, because of the uncertainty related to the NMEs.

Next generation experiments
The situation will be different for the following generation of experiments, which are mostly the natural evolution of current experiments to the ton-scale of decaying material. With the increased amount of material, a larger statistics will be achieved and stronger bounds, of the order of T 0ν 1/2 1 × 10 28 yr, will be feasible. We briefly discuss here the main current proposals for the next 10-20 years. The time schedules for these projects will be necessarily vague, as they will depend on the results of the present ones.
Let us start with CUPID (CUORE Upgrade with Particle ID) [79,385], which will be the evolution of the previously discussed CUORE experiment. The goal of CUPID is to use particle tracking in order to have a better discrimination of background and ultimately allow a background-free experiment: the target is < 0.1 counts/(ton yr) [377]. A first demonstrator, named CUPID-0 [79], is already running with about 5 kg of 82 Se, and already obtained the strongestto-date constraint on the life-time on this isotope. In order to reach the target sensitivity T 0ν 1/2 1 × 10 27 yr, however, further improvement in the crystals quality and radio-purity is required. A full development plan for CUPID is currently under discussion.
The successor of KamLAND-Zen, KamLAND2-Zen [372,378,386] will benefit the upgrades of KamLAND into KamLAND2, including the improved light collection and better energy resolution guaranteed by the new photomultipliers, together with an increased amount of 136 Xe, to reach at least 1 ton of material. These upgrades will be performed after the completion of KamLAND-Zen 800, expected to start this year. The target sensitivity after 5 yr will be T 0ν 1/2 1 × 10 28 yr 19 , sufficient to give a complete coverage of the inverted ordering range for m ββ . Future studies will also test the possibility to accommodate scintillating crystals inside the detector and run a multi-isotope experiment.
Back to 76 Ge-based experiments, the efforts of the Gerda and Majorana collaborations will join to work on the LEGEND (Large Enriched Germanium Experiment for Neutrinoless Double beta decay) experiment. Learning from both its precursors, LEGEND will need further background rejection and will be built in different phases. The first module, LEGEND-200, made of 200 kg of Germanium and expected to start in 2021, will be built on top of the existing Gerda infrastructures and will have a target sensitivity T 0ν 1/2 1 × 10 27 yr in 5 yr. The full scale detector, LEGEND-1000, consisting in several modules summing up to a total of 1 ton of material, will have as an ultimate goal T 0ν 1/2 1 × 10 28 yr in 10 yr [387], giving a full coverage of the inverted mass ordering region. Even larger in size, nEXO [388,389] will replace the EXO-200 experiment after its completion, expected this year. The new detector will use 5 ton of Xenon in order to reach T 0ν 1/2 1 × 10 27 yr with just 1 yr of data and T 0ν 1/2 1 × 10 28 yr with the full statistics, after 10 yr.
After the completion of the upcoming phase, NEXT-100 will be possibly upgraded into NEXT 2.0, which will need a 1.5 ton of Xenon to obtain the statistics for achieving T 0ν 1/2 1 × 10 27 yr after 5 yr of running [38,372]. In the same way, the Panda-X-III collaboration is also planning a 1 ton scale phase II with a target of T 0ν 1/2 1 × 10 27 yr [380].
The last comment regards another interesting possibility related to the SNO+ experiment. The THEIA proposal [390] is a concept study for a gigantic detector of something around 30-100 kton of target material which will use waterbased liquid scintillator. Such target allows to track both Cherenkov and delayed scintillation light, thus enabling high light yield and low-threshold detection with attenuation close to that of pure water. The result is that such a detector would be able to achieve excellent background rejection thanks to directionality, event topology, and particle ID, with very large statistics. Loading of metallic ions which can undergo neutrinoless double beta decay would enable to use the THEIA detector for studying the Dirac/Majorana nature of neutrinos. Given the size of the detector, a 0.5% loading will allow to store several tons of decaying material, which naturally result in huge statistics when compared with current experiments. A 3% loading with natural (not enriched) Tellurium will be sufficient to reach T 0ν 1/2 1 × 10 28 yr in 10 yr [372,391].

D. Prospects from cosmology
There are a number of studies in the literature focused on forecasting the expected sensitivity from both future CMB and large scale structure surveys to the total neutrino mass m ν [99,100,[392][393][394][395][396][397][398][399][400][401]. Awaiting for very futuristic measurements which may allow for the extraction of each of the individual masses associated to the neutrino mass eigenstates (see section IV), the extraction of the neutrino mass ordering strongly relies on the error achieved on m ν for a chosen fiducial value of the neutrino mass. A complete, updated and useful summary is provided in Table II of Ref. [23], which shows the expected sensitivity [σ( m ν )] from different future cosmological probes, assuming the fiducial value m ν = 0.06 eV. Nevertheless, the authors of Ref. [44] considered different fiducial values for the total neutrino mass and computed the odds for the normal versus the inverted ordering for possible combinations of future cosmological probes including the current information from oscillation experiments. We shall comment on these results towards the end of this section.

CMB prospects
Two main missions are expected to lead the next decade generation of CMB experiments, albeit a number of other experiments are in progress between now and then. The latter list includes ground-based observatories as the ACT (Atacama Cosmology Telescope) [402], the SPT-3G (South Pole Telescope-3G) [403], the Simons Array [404], 19 The collaboration usually quotes the sensitivity in terms of m ββ 20 meV, we converted it into a value of T 0ν 1/2 assuming the most pessimistic value for the NME.
CLASS [405], BICEP3 [406] and the Simons Observatory 20 . The two main missions are expected to be the CMB-Stage IV project [408] and CORE (Cosmic Origin Explorer) [409]. The former, the CMB-Stage IV project [408], expected to be the definitive ground-based CMB experiment, aims at 250000 detectors operating for four years, covering a 40% fraction of the sky. Depending on the beam size and on the effective noise temperature, CMB Stage IV could reach sensitivities of σ( m ν ) = 0.073 − 0.11 eV, assuming m ν = 0.058 eV as the fiducial model and an external prior on the reionization optical depth of τ = 0.06 ± 0.01, see Ref. [408] for the precise configuration details. The latter, CORE, a medium-size space mission proposed to the European Space Agency (ESA) [409], is expected to have an one order of magnitude larger number of frequency channels and a twice better angular resolution than Planck. With these improved capabilities, CORE could achieve a sensitivity of σ( m ν ) = 0.044 eV [23,399], for a fiducial total neutrino mass of 0.06 eV. As it is evident from these estimates, future CMB experiments alone will not be able to determine the neutrino mass ordering.

Large scale structure prospects
From the large scale structure perspective, in analogy to the future CMB probes, there are also two main surveys, DESI (Dark Energy Spectroscopic Instrument) [410,411], a ground-based telescope which will improve the SDSS-III and IV legacies (BOSS [129] and eBOSS galaxy surveys [144]), and the Euclid space mission [398]. The baseline design of DESI assumes that it will run over five years, covering 14000 deg 2 of the sky targeting four different tracers: Bright, Luminous Red and Emission Line Galaxies plus quasars in the redshift interval (0.05 < z < 1.85), and a Lyman-α survey in the 1.9 < z < 4 redshift interval. The expected error in m ν from DESI and Planck data is 0.02 eV. This number corresponds, approximately, to a 2σ determination of the neutrino mass ordering in case neutrinos have the minimal mass within the normal ordering scenario [411]. The authors of Ref. [397] have also explored a number of possible combinations of DESI with other surveys. Namely, combining DESI measurements with the final results from DES, an error of 0.017 eV in m ν could be achieved. Their most constraining result, σ( m ν ) = 0.011 eV, however, arises from an extension of the DESI survey, together with data from Euclid and LSST (Large Synoptic Survey Telescope) [412,413] (see below). In case this small error is achieved, the neutrino mass ordering can be determined with a high accuracy, again assuming a massless lightest neutrino and normal ordering. Other analyses have also reduced the nominal σ( m ν ) = 0.02 eV expected from the DESI survey replacing the Planck CMB information with that expected from the future CMB Stage IV [408] or CORE [399] probes.
Euclid, an ESA mission expected to be launched early in the upcoming decade, mapping ∼ 15000 deg 2 of the sky, has also been shown to provide excellent capabilities to test the neutrino properties [398]. Euclid will focus on both galaxy clustering and weak lensing measurements, which, combined with Planck CMB data, will provide errors on the sum of the neutrino masses of σ( m ν ) = 0.04 eV [392] and σ( m ν ) = 0.05 eV [178], respectively, albeit exploiting the mildly non-linear regime could highly reduce these errors [414]. While these errors are large to extract useful information concerning the neutrino mass ordering, the weak gravitational lensing abilities from Euclid have also been considered to extract the neutrino mass ordering when it lies far enough from the degenerate region, see e.g. Ref. [398]. The addition of future CMB measurements, as those from CORE, could notably improve the expected Euclid sensitivity. The authors of [400] have shown that CMB measurements from CORE, combined with full shape measurements of the galaxy power spectrum and weak lensing data from Euclid, could reach σ( m ν ) = 0.014 eV. This result clearly states the complementarity of cosmic shear and galaxy clustering probes, crucial to test the neutrino mass ordering. Further improved measurements of the reionization optical depth τ could strengthen this bound and consequently the sensitivity to the ordering of the neutrino masses [400,401,415], see the following section. Other future large scale structure surveys are the aforementioned LSST and WFIRST [416,417], that will lead as well to accurate measurements of the total neutrino mass. Their combination with e.g. Euclid could provide an error of a few meV on the total neutrino mass, σ( m ν ) 0.008 eV [418].
The above neutrino mass (neutrino mass ordering) projected errors (sensitivities), even if strongly constraining, are highly dependent on the fiducial value of m ν , in the sense that the majority of the forecasts (a) are usually carried out assuming the minimal neutrino mass allowed within the normal ordering scheme, i.e. m ν 0.06 eV 21 ; (b) the quoted sensitivities in the neutrino mass ordering are computed via an extrapolation of the error on the sum of neutrino masses rather than from proper Bayesian comparison tools. The authors of Ref. [44] found that a future CMB CORE-like satellite mission, even combined with a 1% measurement of the Hubble constant H 0 and with the future DESI survey [397,411] can not extract the ordering if nature has chosen a value for the neutrino masses of m ν = 0.1 eV. Odds for the normal versus the inverted ordering of 1 : 1 were reported [44]. When considering the minimum allowed value for the total neutrino mass set by neutrino oscillation experiments, i.e. m ν = 0.06 eV, they quote odds of 3 : 2 (9 : 1) for the case in which CORE and the prior on H 0 without (with) DESI measurements are considered 22 . Therefore, the next generation of CMB and large scale structure surveys will be sensitive to the mass ordering only if it is normal and the lightest neutrino mass is close to zero. The significance of such a measurement will crucially depend on how far m ν lies from its minimum allowed value from oscillation probes.

E. Prospects from 21 cm surveys
Cosmological measurements of the redshifted 21 cm hydrogen line provide a unique test of the Epoch of Reionization (EoR) and the "dark ages", the period before the first stars formed. The 21 cm line is due to spin-flip transitions in neutral hydrogen between the more energetic triplet state and the ground singlet state, and its intensity depends on the ratio of the populations of these two neutral hydrogen hyperfine levels. At a given observed frequency ν, the 21 cm signal can be measured in emission or in absorption against the CMB. The so-called differential brightness temperature δT b therefore refers to the contrast between the temperature of the hydrogen clouds and that of the CMB, which, for small frequencies and up to first order in perturbation theory, reads as [421][422][423][424] where x HI is the fraction of neutral hydrogen, δ b is the baryon overdensity, Ω b h 2 and Ω m h 2 the present baryon and matter contributions to the mass-energy budget of the Universe, H(z) the Hubble parameter and ∂v r /∂r the comoving peculiar velocity gradient along the line of sight. Therefore, 21 cm cosmology aims to trace the baryon overdensities via transitions in neutral hydrogen.
There are a number of current and future experimental setups devoted to detect the 21 cm global signal averaged over all directions in the sky, as EDGES (Experiment to Detect the Reionization Step) [425], the future LEDA (Large Aperture Experiment to Detect the Dark Ages) [426] or DARE (Moon space observatory Dark Ages Radio Experiment) [427]. The EDGES experiment has quoted the observation of an absorption profile located at a frequency of 78 ± 1 MHz, corresponding to a redshift of z ∼ 17, with an amplitude of about a factor of two larger than the maximum expected in the canonical ΛCDM framework [428]. This recent claim has led to a number of studies aiming either to explain the effect or to constrain some non-standard scenarios .
Fluctuations in the redshifted 21 cm signal can be used to compute the power spectrum of the differential brightness temperature. This is the major goal of experiments as GMRT (Giant Metrewave Radio Telescope) [456,457], LOFAR (LOw Frequency ARray) [458], MWA (Murchison Widefield Array) [459] and PAPER (Precision Array for Probing the Epoch of Reionization) [460][461][462], targeting statistical power-spectrum measurements of the 21 cm signal employing large radio interferometers. Even if current experiments have not yet detected the 21 cm cosmological signature, the PAPER collaboration has recently improved the previous upper limits at z = 8.4 [461]. Next decade, high-redshift 21 cm experiments include the SKA (Square Kilometre Array) [463] and HERA (Hydrogen Epoch of Reionization Array) [464]. A three-dimensional map of the 21 cm signal could also be obtained by means of the so-called intensity mapping technique, which measures the collective emission from neutral hydrogen in dense clumps, targeting large regions without resolving individual galaxies in the post-reionization era (z 3) [465][466][467][468]. The experimental efforts for this technique include the GBT-HIM project, with the GBT (Green Bank Telescope) [469], CHIME (Canadian Hydrogen Intensity Mapping Experiment) [470], the Tianlai project [471] and SKA-mid frequency [472], see e.g. [473].
Despite the fact that the primary task of future 21 cm experiments is to improve our current knowledge of the reionization history, they provide as well an additional tool for fundamental cosmology [415,[474][475][476][477][478][479][480][481][482][483][484][485], complementary to CMB missions and galaxy surveys. Indeed, 21 cm cosmological observations will play a very important role concerning neutrino physics. As previously stated, there are two types of experiments. First of all, we will have observations focused on the pre-reionization and EoR periods, that can probe very large volumes (where the non-linear scale is small). Remember that the largest signal from relic neutrino masses and their ordering appears at scales which, at the redshifts attainable at galaxy clustering surveys, lie within the mildly non-linear regime. Therefore one needs to rely on either N-body simulations or on analytical approximations for the matter power spectrum to simulate the massive neutrino signature. EoR 21 cm experiments will achieve the scales required to observe the neutrino signature within the linear regime, avoiding the simulation problems described in section IV B. In this regard, these probes may widely surpass the constraints on neutrino masses expected from even very large galaxy surveys [394,480,482,[486][487][488][489][490]. Furthermore, the neutrino constraints will be largely independent of the uncertainties in the dark energy fluid, which, as we have seen in section IV C 2, have instead a non-negligible impact in lower redshift, galaxy survey measurements. This is a byproduct of using the 21 cm line to trace the matter overdensities: at redshifts z 2, the universe starts to be dominated by the dark energy fluid and the growth of matter perturbations is modified depending on the dark energy equation of state w(z), whose precise time-evolution remains unknown. Consequently, for a given perturbation in the matter fluid, a suppression in its structure growth could be either due to the presence of massive neutrinos or to an evolving dark energy fluid. Focusing at higher redshifts, the neutrino mass constraints from 21 cm probes will be largely independent of the uncertainties in the dark energy fluid properties.
Expectations from MWA, SKA and FFTT (Fast Fourier Transform Telescope) [486] were considered in Ref. [482]. Focusing on 4000 hours of observations of two areas in the sky in a range of z = 6.8 − 8.2 (divided into three redshift bins) and a value of k max = 2 Mpc −1 , the reported errors on m ν are 0.19 (0.027), 0.056 (0.017), 0.007 (0.003) for MWA, SKA and FFTT, respectively, in their middle (optimistic) scenarios 23 , when combined with Planck measurements. These forecasts were performed for a fiducial Ω ν h 2 = 0.0875, which corresponds to a quite high value for the neutrino mass, lying in the fully degenerate neutrino mass spectrum.
The authors of Ref. [488] devoted a dedicated analysis to establish the potential for extracting the neutrino mass ordering combining the FFTT capabilities with future CMB polarization measurements. Based exclusively on the induced effect of the neutrino mass ordering on the cosmic expansion rate, a robust 90% CL neutrino mass ordering extraction was reported if m ν < 0.1 eV, regardless the underlying true ordering (i.e. normal or inverted). In Ref. [490], the authors propose to combine ground-based CMB polarization observations, SKA Phase 2 and BAO measurements from DESI. With these data sets, a 2σ extraction of the neutrino mass ordering seems feasible, unless the neutrino spectrum is degenerate. Notice that these results arise from the signature induced by the neutrino mass ordering in the cosmic expansion rate, as the minimum cutoff of the wavenumber in the 21 cm observations is k min = 0.06h Mpc −1 , while the wavenumber corresponding to the neutrino free-streaming scale is k min 0.02h Mpc −1 for a 0.05 eV massive neutrino.
More futuristic 21 cm experiments, as FFTT, may open the possibility of going beyond measurements of the total neutrino mass m ν and measure the individual neutrino masses, revealing the uniqueness of such experiments for constraining the neutrino properties. As shown in Figure 9 in section IV, the differences in the power spectra for the two possible mass orderings are tiny. Therefore, exquisite precision measurements are required to identify such signatures. Galaxy surveys, already discussed in the previous section, are limited by two facts. The first one is related to non-linearities, which will not allow for a measurement of the power spectrum at scales k > 0.2h Mpc −1 at small redshifts, see section IV B. Since the non-linear scale at z = 8 is k 3h Mpc −1 , both SKA and FFTT can measure the entire linear region and be more sensitive to the scale-dependent suppression, which is different in the two neutrino mass orderings. The second one is related to the fact that a galaxy survey requires a large number density of tracers to ensure a good sensitivity at small scales, while for 21 cm surveys, tracing the ubiquitous permeating hydrogen, a high-density antennae distribution will already warrant excellent small-scale sensitivities. One drawback of 21 cm probes are foregrounds, which should be kept under control.
The authors of Ref. [487] have studied the perspectives for extracting the individual neutrino masses with SKA and FFTT, finding that FFTT could be able to distinguish all the three neutrino masses from zero at the 3σ level, due to its enormous effective volume (see Figure 3 of Ref. [487]). Extracting the neutrino mass ordering directly from the individual masses, however, was shown to be a very difficult achievement. Our calculations show that, for the total neutrino mass we use here as a reference, m ν = 0.12 eV, the differences among the lightest (l), medium (m) and heaviest (h) neutrino mass eigenstates between the normal and inverted orderings are (|∆m l |, |∆m m |, |∆m h |) = (0.015, 0.0209, 0.0059) eV, which, especially for the case of |∆m h | = 0.0059 eV, are tiny and very difficult to resolve, even with very futuristic 21 cm measurements. While increasing the exposure of FFTT may improve its capabilities for this purpose (the error in the most optimistic FFTT scenario of Ref. [482] is 0.003 eV), it seems an extremely challenging task. Figure 13 depicts the differences in the values of the three neutrino masses as a function of the total neutrino mass between inverted and normal orderings. We show with a dashed vertical line our representative case m ν = 0.12 eV (the present most constraining 95% CL upper limit) and another one for m ν = 0.34 eV (the most recent 95% CL bound from the Planck collaboration after the removal of systematics in their polarization data at high angular scales [42]). Notice that, as expected, the differences between the values of the three neutrino masses decrease with the total neutrino mass. In this regard, the lower the neutrino mass, the easier it could be to single out the three neutrino mass eigenstates, because they are more separated. However, an extraction of the mass ordering in the non-degenerate region via the values of the individual neutrino masses seems very difficult. Indeed, Figure 14 illustrates the values of the individual neutrino masses for the heaviest, medium and lightest states for the normal and inverted orderings as a function of the total neutrino mass. The bands, from top to bottom panels, depict the errors  Figure 14 shows the results if we assume the futuristic value of σ(m i ) = 0.005, expected to be achieved by FFTT. In this case, a measurement of the three neutrino masses will be achieved. Furthermore, in this (very optimistic) situation, the error bars will be, in principle, sufficiently small to detect the presence (or the lack) of two massive neutrino states with masses in the 0.02-0.03 eV range, required if the ordering is normal to explain m ν 0.1 eV, which would strongly confirm the normal (or inverted) neutrino mass ordering. If σ(m i ) = 0.005, the detection of the mass ordering will still be possible even if m ν 0.1 eV, since the error on m ν will allow to exclude the inverted ordering with great accuracy. As already mentioned, another possibility is the so-called 21 cm intensity mapping, which will focus on low redshifts z 3 and will measure, with low angular resolution, the integrated 21 cm flux emitted from unresolved sources observing large patches of the sky. The lack of high angular resolution will result in a less precise measurement of non-linear scales. On the other hand, low angular resolution will imply a much faster survey. Future planned intensity mapping surveys are developed within the Phase 1 of the SKA experiment, which will include a wide and deep survey at low redshifts (z 3, the SKA1-MID array) and a narrow and deep survey at higher redshift (3 z 6, the SKA1-LOW array), and within the Phase 2 of SKA (SKA2). Since, in some sense, these intensity mapping probes will be complementary to future planned optical surveys, as DESI or Euclid, it makes sense to combine their expected results. The intensity mapping technique, as galaxy clustering, is also affected by bias uncertainties and non-linearities at small scales. Several studies have been carried out in the literature to unravel the perspectives of the intensity mapping technique in unveiling the neutrino properties. Some of them include the combination of the expectations from future large scale structure and intensity mapping surveys [394,400,401,467,483,491]. Notice that all these studies rely on different assumptions on the cosmological parameters, on the foregrounds and on the systematic uncertainties, therefore we can not do comparisons among them. Instead, we quote the most recent findings and the impact for an eventual future detection of the neutrino mass ordering.
The authors of Ref. [491] found that, by combining SKA1-LOW with Planck measurements, the 95% CL error on m ν could be ∼ 0.089 eV. It is remarkable that such a combination could potentially rule out the inverted ordering scenario, assuming that normal ordering is realized in nature. These authors also find that, under identical assumptions in the forecasted analyses, their combination of intensity mapping surveys (SKA1-LOW and MID) should be regarded as competitive with future spectroscopic surveys concerning neutrino mass properties. The authors of Ref. [400] showed that constraints of the future CORE CMB mission and galaxy redshift/weak lensing large scale structure surveys (as Euclid) on the neutrino mass can be improved if a prior on the reionization optical depth from 21 cm probes as HERA or SKA is also included. A prior of σ(τ ) = 0.001 will reduce the freedom in the amplitude of the primordial power spectrum A s , as CMB measurements mostly constrain the combination A s exp(−2τ ), see section IV A. Therefore, the direct correlation between m ν and A s , both modifying the amplitude of the matter power spectrum (although the change induced by m ν is, obviously, scale dependent), is largely affected by the presence of a precise determination of τ . The 1σ sensitivity they find for the combination of CORE, Euclid plus the prior on the optical depth from future 21 cm observations is σ( m ν ) = 0.012 eV 24 .
Nevertheless, as carefully detailed above, even if these tiny errors on m ν will be reached and extrapolated to an error on the individual neutrino mass eigenstates, the possibility of extracting the neutrino mass ordering via singling out the neutrino mass eigenstates with cosmological observables remains unfeasible, unless very visionary scenarios, as FFTT under the most optimistic assumptions, are envisaged.

F. Prospects from core-collapse supernova
Neutrinos from core-collapse supernovae offer an independent and complementary way to test neutrino physics. The existence of these neutrinos was robustly confirmed by the detection of twenty-five events from Supernova 1987A in the Large Magellanic Cloud [492][493][494], located at ∼ 50 kpc from our Milky Way galaxy. Such a detection allowed to set very compelling bounds on a number of neutrino properties [495,496]. Even if laboratory experiments have surpassed some of these limits, the eventual detection of supernovae neutrinos will still provide precious information about the details of the explosion process (see e.g. [497][498][499] and references therein), and also of neutrino mixing effects in dense media, see also Ref. [500].
Neutrino production in core-collapse supernovae occurs in a number of different stages. The first one is the infall, in which electron neutrinos are produced, confined, as a result of the process e − + p → n + ν e . When electrons are converted, the outwards pressure they generate disappears and the gravity forces are no more balanced: the core will start to collapse until its density reaches that of matter inside atomic nuclei, i.e. nuclear densities. Once these densities are reached, matter becomes incompressible, and a hydrodynamic shock is formed. As this shock wave propagates outwards, it heats up the nuclei and disintegrates them, releasing neutrinos. This initial neutrino release is commonly known as neutronization burst, and it is mainly composed of ν e and may last for a few tens of milliseconds. After the neutronization burst, the remnant proto-neutron star may evolve into a neutron star or collapse to a black hole, depending on the mass of the progenitor star. During this phase of explosion and accretion, which lasts for one to two seconds, the ν e contribution is still the dominant one, albeit there is also a contribution from other (anti)neutrino flavors, in particularν e . The neutrinos produced in the cooling stage give the main contribution to the total flux, as it is in this phase when the supernova releases its energy via all-flavor neutrino-antineutrino pair production, reaching its final cold state. This process lasts for about tens of seconds. The differences in the mean temperature of the neutrino fluxes of ν e ,ν e and ν x (ν x ) are due to the different medium opacity of each species. The larger the opacity, the lower the temperature that the (anti)neutrino will have at decoupling. The neutrino fluxes read as [497] φ(E ν ) = N 0 (α + 1) (α+1) E ν Γ(α + 1) where N 0 is the total number of emitted neutrinos, and both α and the mean energy E ν are flavor dependent. The supernova neutrino energy spectra peaks around the 10 − 20 MeV region. The most popular process for supernova neutrino detection is inverse beta decay on protons (ν e + p → n + e + ). Other possibilities include elastic scattering on electrons (ν +e − → ν +e − ), whose kinematics may provide information on the supernova location. Supernova neutrinos can also interact with nuclei via charged current or neutral current interactions, giving rise to charged leptons and/or excited nuclei which may provide flavor tagging. A very important process on argon nuclei is ν e + 40 Ar → e − + 40 K * , which allows for electron neutrino tagging. In practice, water Cherenkov and scintillator detectors are mostly sensitive to electron antineutrinos via inverse beta decay, while the liquid argon technique mainly detects electron neutrinos. While other flavors may also be detected, the two processes above are the dominant ones. Large detector volumes (dozens of kilotons) are required to detect neutrinos from core-collapse supernovae located at ∼ O(10) kpc. A convenient way to scale the total number of supernova neutrino events in a detector of given effective mass is [501,502] (25) 24 More recently, this very important synergy between Euclid and future 21 cm surveys, concretely with the intensity mapping survey SKA1, has been further assessed in Ref. [401].
In the expression above, E B is the gravitational binding energy of the collapsing star and D OS the distance between the observer and the supernova. Assuming sensitivity to all reactions, the reference rate is N 0 = O(10 4 ) for the Super-Kamiokande water Cherenkov detector with 32 kton and 5 MeV energy detection threshold. References [497,503] give an estimate of the number of neutrino events for a number of ongoing and future facilities, based on different detection techniques: water Cherenkov (including also those with long string photosensors in ice, as Icecube and PINGU), liquid argon time projection chambers, and liquid scintillators. Upcoming neutrino detectors, already described in section VI A and crucial for oscillation physics measurements, such as the JUNO liquid scintillator [348], the liquid argon DUNE [299][300][301][302] and the water Cherenkov Hyper-Kamiokande [322] can lead to a number of 6000, 3000 and 75000 supernova neutrino events respectively, assuming that the explosion occurs at 10 kpc from our position. Flavor transitions inside a supernova have been carefully reviewed in Refs. [497,498] (see also [504][505][506][507][508]). Here we summarize the most relevant results. As we have seen in section II, when neutrinos propagate through matter their mixing effects undergo the so-called MSW mechanism, feeling a matter potential which is proportional to the electron number density N e . If the supernova matter density has a profile which varies slowly, the neutrino matter eigenstates will propagate adiabatically and their final flavor composition will depend on the neutrino mass ordering, which will establish whether or not resonant transitions associated to each neutrino mass squared difference (solar and atmospheric) take place 25 . In the normal ordering case, the neutrino fluxes will have a significantly transformed spectrum, while the electron antineutrino one will only be partially transformed (F final νe = F initial νx and F final νe = cos 2 θ 12 F initial νe + sin 2 θ 12 F initial νx ). In the inverted ordering case, the effects on the electron neutrino and antineutrino fluxes will be approximately the opposite ones (F final νe = sin 2 θ 12 F initial νe + cos 2 θ 12 F initial νx and F final νe = F initial νx ). Once neutrinos exit from supernovae, they can still undergo flavor transitions if they traverse the Earth. Their final flavor composition at the detector location will again depend on the neutrino mass ordering, as matter effects in Earth depend on it, see e.g. Ref. [497] and references therein.
Furthermore, collective effects from neutrino self-interactions, due to ν e +ν e → ν x +ν x flavor processes, can lead to departures from the above summarized three-flavor oscillation picture [498,[509][510][511][512][513][514]. The effective potential, proportional to the difference between the electron antineutrino and the muon/tau antineutrino fluxes, and inversely proportional to the supernova radius, should dominate over the standard matter one, leading to spectral swaps or splits [515][516][517][518]. In the early stages, these self-interacting effects are sub-leading for mass ordering signatures, albeit we shall comment on possible non-thermal features in the neutrino or antineutrino spectra which depend on the mass ordering [519].
In the following, we shall summarize the most relevant available methods to extract the neutrino mass ordering using the mentioned fluxes. For a recent and thorough review of the mass ordering signatures from supernovae neutrinos, we refer the reader to Ref. [497]. The electron neutrinos produced in the neutronization burst undergo the MSW effect being fully (only partially) transformed, i.e. F final νe = F initial νx (F final νe = sin 2 θ 12 F initial νe + cos 2 θ 12 F initial νx ) if the mass ordering is normal (inverted), respectively. Therefore, detectors with good ν e tagging, such as liquid argon or water Cherenkov ones, will detect a neutronization burst only in the inverted neutrino mass ordering case. Concerning the accretion phase, and once electron antineutrinos are also produced, as they are almost unchanged in the MSW resonance, the largest signature is expected to occur for the normal ordering case for the three type of aforementioned detector types (liquid argon, water Cherenkov and scintillator), although the Icecube detector, with its excellent capabilities to reconstruct the time dependence of the signal, could also distinguish between the normal and inverted mass orderings [520]. On the other hand, collective effects, which lead to spectral swaps in the electron (anti)neutrino spectra, show very sharp features at fixed energy values which depend, among other factors, on the neutrino mass ordering. However, these signatures are not as robust as the ones existing in the neutronization and accretion phases. Finally, a very significant imprint of the neutrino mass ordering on the supernovae neutrino fluxes is that due to their propagation through the Earth interior, where the standard MSW effect will induce a few percent-level oscillatory pattern in the 10 − 60 MeV energy range, in the electron (anti)neutrino spectra in case of (normal) inverted mass ordering. The detection of these wiggles requires however excellent energy resolution.

G. Prospects from relic neutrino direct detection
In the early Universe, neutrinos decoupled from the cosmic plasma during the cool down, in a process similar to the one leading to the formation of the CMB but at an earlier time, when the universe was seconds to minutes old. These neutrinos have been free-streaming for such a long time that they have decohered and are currently propagating as mass eigenstates. The decoupling of neutrinos occurred just before e ± annihilated and reheated photons, leading to the following ratio between the photon (T γ ) and neutrino (T ν ) temperatures, see Equation (15): Today, the temperature of the neutrino background is T 0 ν 1.6 × 10 −4 eV. Their mean energy is E ν 3 T ν 5 × 10 −4 eV, much smaller than the minimal mass of the second-to-lightest neutrino as required by flavor oscillations, so that at least two out of three neutrinos are non-relativistic today. The cosmic neutrino background (CνB) is the only known source of non-relativistic neutrinos and it has never been detected directly.
Apart from the imprints that relic neutrinos leave in the CMB (see section IV A), which allow to have an indirect probe of their existence through the determination of N eff , the direct detection of the CνB would offer a good opportunity to test neutrino masses and their ordering. Capturing relic neutrinos is not only rewarding from the point of view of what we can learn about neutrino properties, but also because it would be a further confirmation of the standard Big Bang cosmological model. Different ideas on how to achieve such a detection have been proposed [521][522][523][524][525][526][527][528][529][530][531][532], ranging from absorption dips in the ultra-high-energy (UHE) neutrino fluxes due to their annihilation with relic neutrinos at the Z boson resonance, to forces generated by coherent scattering of the relic bath on a pendulum and measured by laser interferometers. Most of these proposed methods are impractical from the experimental point of view. The one exploiting UHE neutrinos [522,523,525,526] has two problems, one related with the fact that it is difficult to think about a source that produces such UHE neutrinos, of energies and another one regarding the difficulties of detecting a large enough sample of UHE neutrinos in order to resolve the dips. The method based on interferometers [532] is even more complicated to address. At interferometers, current sensitivities to accelerations are of the order of a 10 −16 cm/s 2 , with an optimistic estimation of a 3 · 10 −18 cm/s 2 [532] for the incoming generation. However, expected accelerations due to relic neutrino interactions are of the order of 10 −27 − 10 −33 cm/s 2 [524,532], many orders of magnitude below the sensitivity of the nextgeneration interferometers.
The most promising approach to detect relic neutrinos is to use neutrino capture in a β-decaying nucleus A where the signal for a positive detection is a peak located about 2m ν above the true β-decay endpoint (see below). In particular, tritium is considered as the best candidate since it has a high neutrino capture cross section, low Q-value and it is long-lived [529,[533][534][535][536]. The proposal for an experiment chasing this purpose was made in [529]. Currently, efforts are put for such experiment, the PonTecorvo Observatory for Light Early-Universe Massive-Neutrino Yield (PTOLEMY) [370,371], to be built. The experiment has recently been approved by the Scientific Committee of the Italian National Laboratories of Gran Sasso and, in the following months, the existing prototypes for various components are expected to be moved from Princeton, where the R&D has been performed up to now, to Gran Sasso. The idea is to implant the tritium source on graphene layers, to avoid the problems related to a gaseous source, then collect and measure the energy of the emitted electrons using a combination of MAC-E filter, radio-frequency tracking and micro-calorimetry to obtain a determination of the β-decay and neutrino capture spectrum of tritium with an energy resolution of the order ∆ 0.05 − 0.1 eV. The total expected event rate from relic neutrino capture for a PTOLEMY-like experiment, assuming the estimated tritium mass of 100 g, is where n 0 (ν h R,L ) is the averaged number density of relic neutrinos with right (R) or left (L) helicity, N T = M T /m( 3 H) is the approximated number of tritium atoms in the source,σ 3.834 × 10 −45 cm 2 [536], and f c (m i ) is a massdependent overdensity factor that accounts for the clustering of relic neutrinos under the gravitational attraction of the matter potential (mostly from the dark matter halo) of our galaxy. This last factor was originally computed in Refs. [537,538] and later updated in Ref. [539] (see also Ref. [540]), where smaller masses were considered for the neutrinos, and the treatment of the matter potential of the Milky Way was improved. The values of f c (m ν ) range from 1.1 − 1.2 for a neutrino with m ν = 60 meV to 1.7 − 2.9 for m ν = 150 meV [539].
For unclustered neutrinos (i.e. f c = 1) and 100 g of tritium, the expected number of events per year is [536] Γ  where the upperscripts D and M stand for the possible Dirac and Majorana neutrino character. If neutrinos are Majorana particles, the expected number of events is doubled with respect to the Dirac case. The reason is related to the fact that, during the transition from ultra-relativistic to non-relativistic particles, helicity is conserved, but not chirality. The population of relic neutrinos is then composed by left-and right-helical neutrinos in the Majorana case, and only left-helical neutrinos in the Dirac case. Since the neutrino capture can only occur for left-chiral electron neutrinos, the fact that in the Majorana case the right-handed neutrinos can have a left-chiral component leads to a doubled number of possible interactions. While this means that in principle it is possible to distinguish the Dirac or Majorana neutrino nature with a precise determination of the event rate, there are two problems. First of all, even without assuming new physics, the factor of two coming from the neutrino nature is degenerate with the clustering factor, see Equation (29), so that a precise calculation of f c is required to determine if neutrinos are Dirac or Majorana particles through the direct detection of relic neutrinos [539]. Moreover, non-standard interactions can increase the event rate in the Dirac case by a factor larger than two, canceling the difference with Majorana neutrinos in some scenarios [541]. Let us come back to the PTOLEMY proposal. Instead of considering the total event rate, for this kind of experiment it is much better to study the energy spectrum, as the direct detection of relic neutrinos can only be possible if one can distinguish the signal events due to neutrino capture from the background events due to the β-decay of tritium. A crucial issue for such an experiment, actually more important than the event rate, is therefore the energy resolution. In order to distinguish the peak due to the captured relic neutrinos from the β-decay background, a full-width half maximum (FWHM) energy resolution ∆ 0.7m ν is needed [536]. If neutrinos are non-degenerate in mass, the neutrino capture signal has a peak for each of the separate neutrino mass eigenstates. The full expression of the energy spectrum of neutrino capture, given an energy resolution σ = ∆/ √ 8 ln 2, can be written as: where E end is the energy of the β-decay endpoint, E end = E end,0 − m lightest , being E end,0 the endpoint energy when m lightest = 0. If the energy resolution is good enough, the three peaks coming from the three neutrino mass eigenstates could be resolved, each of them with an expected number of events modulated by |U ei | 2 . This might lead to a positive detection of the neutrino mass ordering, since the electron-flavor component of ν 1 is larger than the one of ν 2 and ν 3 , and therefore the furthest peak from the β-decay endpoint (again if neutrinos are non-degenerate) is enhanced if the ordering of neutrino masses is inverted. This can be seen in the three panels of Figure 15, which also show the effect of changing the mass ordering on the β-decay spectrum. Dashed lines represent the spectrum which would be determined by an experiment capable of measuring the β spectrum with zero energy uncertainty, while solid lines represent the shape of the spectrum that one would observe in a real experiment. We plot in red (blue) the spectrum obtained using normal (inverted) ordering, a FWHM resolution ∆ = 10 meV (top), ∆ = 20 meV (middle), ∆ = 50 meV (bottom) and a lightest neutrino mass m lightest = 10 meV. As we can see from the figure, the kink commented in section III is clearly visible when one observes the huge number of events that come from the 100 g of decaying tritium with a sufficient energy resolution. While for distinguishing the relic neutrino events from the β-decay background and for having a direct detection of the CνB the energy resolution is a crucial requirement, in principle even a worse energy resolution may allow to determine the neutrino mass scale and the mass ordering, thanks to the fact that we expect less events near the endpoint when the ordering is inverted. A direct observation of the amplitude of all the CνB peaks, however, would give a much cleaner signal, because the peak corresponding to the heaviest neutrino would be always higher in the inverted ordering case, independently of any other factor.
In summary, the CνB capture event rate in a PTOLEMY-like experiment (Equation (31)), even within SM physics and without considering non-standard interactions, depends on several main unknowns: i) the absolute neutrino mass, ii) the matter distribution (especially that of dark matter) in our galaxy, iii) the nature of neutrino masses (whether neutrinos are Dirac or Majorana particles), and iv) the true mass ordering. This last dependence is encoded in the |U ei | 2 factor in Equation (31) and it is only accessible if neutrinos are non-degenerate. A quantitative study on the PTOLEMY capabilities in determining the mass ordering has not been published yet, but a new Letter of Intent is in preparation [371] 26 .

VII. SUMMARY
Identifying the neutrino mass ordering is one of the major pending issues to complete our knowledge of masses and mixings in the lepton sector. The two possibilities, normal versus inverted, may result from very different underlying symmetries and therefore to single out the one realized in nature is a mandatory step to solve the flavor puzzle, i.e. to ensure a full theoretical understanding of the origin of particle masses and mixings. We have presented a comprehensive review on the current status and on future prospects of extracting the neutrino mass ordering via a number of different ongoing and upcoming observations. Furthermore, the most updated and complete result on the preference for a given neutrino mass ordering from a Bayesian global fit to all 2018 publicly available neutrino data has also been presented.
Currently, among the three available methods to extract the neutrino mass ordering (oscillations, neutrinoless double beta decay searches and cosmological observations), the leading probe comes from oscillations in matter, measured at long-baseline accelerator or atmospheric neutrino beams in combination with reactor experiments. The latest frequentists global data analysis results in a preference for normal mass ordering with ∆χ 2 = 11.7 (∼ 3.4σ), mostly arising from the combination of the long-baseline T2K and NOνA data with reactor experiments (Daya Bay, RENO and Double Chooz), plus the latest atmospheric neutrino results from Super-Kamiokande. Similar results for the preference in favor of the normal mass ordering arise from other global fit analyses [17].
Cosmological measurements are able to set indirect, albeit independent bounds on the neutrino mass ordering. Neutrinos affect Cosmic Microwave Background (CMB) primary anisotropies by changing the gravitational potential at the recombination period when they become non-relativistic. However, for sub-eV neutrino masses this effect is tiny and the most prominent effect on the CMB is via lensing, as neutrinos, having non-zero velocities, will reduce the lensing effect at small scales. Nevertheless, the largest impact of neutrinos in cosmology gets imprinted in the matter power spectrum. Once neutrinos become non-relativistic, their large velocity dispersions will prevent the clustering of matter inhomogeneities at all scales smaller than their free streaming length. At present, the cosmological constraints on the neutrino mass ordering come from the sensitivity to the total neutrino mass m ν and not via the effects induced in the CMB and matter power spectrum by each of the individual neutrino masses m i . Within the context of the minimal ΛCDM model with massive neutrinos, current cosmological probes cannot provide odds stronger than ∼ 3 : 1 in favor of normal ordering.
Neutrinoless double beta decay searches can also test the neutrino mass ordering if neutrinos are Majorana particles. However, present constraints on the so-called effective Majorana mass do not affect the overall Bayesian analyses.
All in all, the 2018 Bayesian global analysis, including all the neutrino oscillation data available before the Neutrino 2018 conference, results in a 3.2σ preference for the normal neutrino mass ordering which, in Bayesian words, implies a strong preference for such a scenario. One can then combine the oscillation data with 0νββ data from KamLAND-Zen, EXO-200 and Gerda and cosmological observations from Planck, SDSS BOSS, 6DF and SDSS DR7 MGS. Using this conservative cosmological data combination, the aforementioned preference becomes 3.4σ, which raises to 3.5σ if a prior on the Hubble parameter H 0 from local measurements is considered in addition. This clearly states the current power of oscillation results when dealing with neutrino mass ordering extractions.
While in the very near perspective an improved sensitivity (i.e. above the 3.5σ level) is expected mostly from more precise measurements of current long-baseline and atmospheric experiments, and, to a minor extent, from cosmological surveys (Planck, DES and eBOSS among others), there will be a number of planned experiments which will be crucial for extracting the neutrino mass ordering in the non-immediate future.
Of particular relevance are the upcoming neutrino oscillation facilities, as they will be able to measure the neutrino mass ordering with astonishing precision without relying on combinations of different data sets. Such is the case of the Deep Underground Neutrino Experiment (DUNE), that will be able to measure the neutrino mass ordering with a significance of 5σ with seven years of data. Atmospheric neutrino observatories as PINGU or ORCA will also mainly focus on the mass ordering measurement. Some of these future devices could also identify the neutrino mass ordering via the detection of matter effects in the neutrino fluxes emitted at the eventual explosion of a supernova in our galaxy or in its neighbourhood. On the other hand, medium baseline reactor neutrino detectors such as JUNO or RENO will also be able to extract the neutrino mass ordering despite matter effects are negligible for these two experiments. They will focus instead on an extremely accurate measurement of the survival electron antineutrino probability.
Improved masses and detection techniques in neutrinoless double beta decay future searches could go down the 10 meV region in the effective Majorana mass m ββ , and they could be able to discard at some significance level the inverted mass ordering scenario, in the absence of a positive signal. These limits, however, will apply only in case neutrinos have a Majorana nature. Moreover, the determination of the neutrino mass ordering may be complicated by the presence of a light sterile neutrino at the eV scale, as currently suggested by the NEOS and DANSS results.
Concerning future cosmological projects, the combination of different probes will still be required. Near-future CMB and large scale structure surveys will only be sensitive to the neutrino mass ordering via their achieved error on m ν . Furthermore, the accuracy in the extraction of the neutrino mass ordering will strongly depend on how far m ν lies from the minimum allowed value from oscillation probes. The future CMB mission CORE plus the DESI galaxy survey could provide odds of 9 : 1 for normal neutrino mass ordering assuming m ν = 0.056 eV. Even if very futuristic surveys, based on the observation of the 21 cm redshifted line in neutral hydrogen, may be able to extract the individual values of the neutrino masses, their precision on the m i values may not be enough to guarantee a direct determination of the neutrino mass ordering by these means, albeit they can achieve an accurate measurement of the ordering thanks to their unprecedented precision on m ν . Last, but not least, relic neutrino capture in tritium in a PTOLEMY-like experiment could also establish the neutrino mass ordering via an almost perfect energy reconstruction of the β-decay spectrum, ensured by the extremely large amount of tritium adopted. The detection is possible both from a kink in the β-decay spectrum which only appears if the ordering is inverted and from the peaks due to neutrino capture just above the endpoint.
All these future probes may either confirm or reject the current strong preference (∼ 3.5σ) in favor of the normal neutrino mass ordering. Such a preference has kept gaining significance in the recent years, thanks to the fact that current neutrino oscillation experiments have enormously improved our knowledge of neutrino flavor physics.