2020 Global reassessment of the neutrino oscillation picture

We present an updated global fit of neutrino oscillation data in the simplest three-neutrino framework. In the present study we include up-to-date analyses from a number of experiments. Namely, we have included all T2K measurements as of December 2019, the most recent NO$\nu$A antineutrino statistics, and data collected by the Daya Bay and RENO reactor experiments. Concerning the atmospheric and solar sectors, we have also updated our analyses of DeepCore and SNO data, respectively. All in all, these new analyses result in more accurate measurements of $\theta_{13}$, $\theta_{12}$, $\Delta m_{21}^2$ and $|\Delta m_{31}^2|$. The best fit value for the atmospheric angle $\theta_{23}$ lies in the second octant, but first octant solutions remain allowed at $\sim2\sigma$. Regarding CP violation measurements, the preferred value of $\delta$ we obtain is 1.20$\pi$ (1.54$\pi$) for normal (inverted) neutrino mass ordering. These new results should be regarded as extremely robust due to the excellent agreement found between our Bayesian and frequentist approaches. Taking into account only oscillation data, there is a preference for the normal neutrino mass ordering at the $2.7\sigma$ level. While adding neutrinoless double beta decay from the latest Gerda, CUORE and KamLAND-Zen results barely modifies this picture, cosmological measurements raise the significance to $3.1\sigma$ within a conservative approach. A more aggressive data set combination of cosmological observations leads to a stronger preference for normal with respect to inverted mass ordering, at the $3.3\sigma$ level. This cosmological data set provides $2\sigma$ upper limits on the total neutrino mass corresponding to $\sum m_\nu<0.13$ ($0.15$)~eV in the normal (inverted) neutrino mass ordering scenario. These bounds are among the most complete ones in the literature, as they include all currently available neutrino physics inputs.


I. INTRODUCTION
This paper updates the results from a long ongoing series of global fits to neutrino oscillation data [1][2][3][4][5][6][7][8] * . Neutrino flavor conversion was first observed in solar [11] and atmospheric neutrinos [12]. This discovery led to the Nobel prize in Physics in 2015 [13,14] and was confirmed by subsequent results from the KamLAND reactor experiment [15] as well as long baseline accelerator experiments. These were crucial to identify neutrino oscillations as the explanation of the solar  neutrino problem and the atmospheric neutrino anomaly † . In the simplest three-neutrino scenario, the probability for a neutrino to oscillate between flavors is described by six parameters, ∆m 2 21 , |∆m 2 31 |, θ 12 , θ 13 , θ 23 and δ. In addition, there are two possible mass orderings (MO) for neutrinos, according to the positive or negative sign of ∆m 2 31 . In the first case, we talk about normal ordering (NO), and in the latter, about inverted ordering (IO). In Tab. I we summarize the sensitivity of the various experiment types in probing each of the oscillation parameters. Since many of the parameters are measured by several classes of experiments, a combined or global fit of all data will give more precise results than a measurement of a single experiment on its own. Performing such global analysis is precisely the purpose of this study.
The paper is structured as follows: in Sec. II we present the analysis of each class of experiments, focusing on solar experiments and KamLAND, short-baseline reactor experiments, atmospheric experiments and, finally, long-baseline accelerator experiments. Next, we show the results from our global fit to neutrino oscillation data, following a frequentist approach in Sec. III, and a Bayesian approach in Sec. IV. In Sec. V we discuss the effects of the inclusion of non-oscillation data sets and present our final results on the neutrino mass ordering. Finally, we summarize all our results in Sec. VI. † Other mechanisms, such as magnetic moments [16][17][18] or non-standard interactions [19][20][21] could be present only at a sub-leading level [7], for recent analyses see, e.g. [22,23]. ‡ Here, we use the term short-baseline for baselines of the order of 1 km. We will not discuss the searches for light sterile neutrinos. We refer the interested reader to Refs. [24][25][26][27][28][29][30][31].

II. EXPERIMENTAL DATA
In this section we discuss the experimental results included in our global fit with more detail.
We dedicate one subsection to describe each class of experiments, discussing the main details of the data sets analyzed. The results of the oscillation analysis in each sector are presented as well.

A. Solar neutrino experiments and KamLAND
Solar neutrinos are produced in thermonuclear reactions in the interior of the Sun when burning hydrogen into helium. The main nuclear chains producing neutrinos are the so-called protonproton (pp) chain and the CNO cycle. Neutrinos are produced in different reactions with energies ranging from 0.1 to 20 MeV. Our solar oscillation analysis includes data from all past and present solar neutrino oscillation experiments. We use the total rate measurements performed at the radiochemical experiments Homestake [32], GALLEX/GNO [33] and SAGE [34], the low-energy 7 Be neutrino data from Borexino [35,36], as well as the zenith-angle or day/night spectrum from phases I-IV in Super-Kamiokande [37][38][39][40]. Finally, we also include the last results from the Sudbury Neutrino Observatory (SNO), combining the solar neutrino data from the three phases of the experiment [41]. As in previous works, we have considered the low metallicity version of the standard solar model, labeled as AGSS09 [42]. The result of our combined analysis of solar neutrino oscillation data is shown in Fig. 1.
The solar neutrino oscillation parameters were also measured at the KamLAND experiment [43][44][45]. This long-baseline reactor neutrino experiment used a single detector to detect neutrinos from 56 nuclear reactors at an average distance of 180 km. This long distance made KamLAND sensitive to the values of the mass splitting ∆m 2 21 indicated by the solar data analysis. In our global fit, we include KamLAND data as presented in Ref. [44]. The result of our analysis is shown together with the result from the analysis of solar neutrino oscillation data in Fig. 1. As can be seen in the figure, the solar experiments provide a more precise measurement of the solar mixing angle, while KamLAND gives a better determination of the solar mass splitting. Note that, since KamLAND is mostly sensitive to sin 2 2θ 12 , using KamLAND data alone, we would obtain a second minimum in the upper octant of sin 2 θ 12 . This solution is excluded when combining with solar neutrino data, sensitive to sin 2 θ 12 through the observation of the adiabatic conversion in the solar medium. Note, however, that the upper-octant solution may emerge in the presence of non-standard interactions [21,46,47]. Regarding the determination of ∆m 2 21 , a mild tension appears in the combined analysis of solar and KamLAND results. The best fit value obtained by solar experiments, ∆m 2 21 = 4.8 × 10 −5 eV 2 , is excluded by KamLAND with a high confidence level. However, the regions overlap above 90% C.L.. Moreover, notice that the solar experiments and KamLAND show also a marginal sensitivity to θ 13 , which can be enhanced at the combined analysis [6,48]. In order to generate Fig. 1, we have marginalized over θ 13 , without taking any constraint from short baseline reactor data, which we discuss in the next subsection. (colored regions). The best fit values are indicated with dots for the independent analyses and with a star for the combined solar + KamLAND analysis. The reactor mixing angle θ 13 has been marginalized over without further constraint from short-baseline reactor experiments.

B. Reactor neutrino experiments
Besides KamLAND, there are several other reactor neutrino oscillation experiments. Here, we use data coming from the reactor experiments RENO [49] and Daya Bay [50]. Unlike KamLAND, they lie quite close to the nuclear power plants. This makes them sensitive to θ 13 and ∆m 2 31 § .
Using current reactor neutrino data, it was shown that there is also some sensitivity to the solar parameters [52]. Note that these, however, are not competitive with the results coming from KamLAND and solar experiments and therefore we fix in our analyses the solar parameters to the ones measured by those experiments, as discussed in the previous section.
The Reactor Experiment for Neutrino Oscillation (RENO) is a neutrino oscillation experiment located at the Hanbit Nuclear Power Plant (South Korea), that has been taking data since August 2011. Two functionally identical 16 ton detectors placed at 294 m and 1383 m from the centerline of the antineutrino sources, detect electron antineutrinos produced by six pressurized water reactors (all equally distributed in space along a 3 km line), each with output thermal powers of 2.6 GW th § Actually, short-baseline reactor experiments are sensitive to the effective mass splitting ∆m 2 ee = cos 2 θ12∆m 2 31 + sin 2 θ12∆m 2 32 [51]. or 2.8 GW th . The average relative fission fractions for these reactor cores can be found in Ref. [53].
In the most recent publication [49], the RENO collaboration reported results that correspond to 2200 days of data taking. From the observation of electron antineutrino disappearance, RENO reported a value for the reactor mixing angle of sin 2 (2θ 13 ) = 0.0896 ± 0.0068, and a value of |∆m 2 ee | = (2.68 ± 0.14) × 10 −3 eV 2 for the observed neutrino mass squared difference. In our analysis, we consider antineutrino events (background subtracted) at the near and far detectors, as reported by RENO [49], distributed along 26 energy bins in prompt energy, ranging from 1.2 MeV to 8.0 MeV. A total of nine systematical uncertainties, accounting for reactor-flux uncertainties σ r = 0.9% (correlated between detectors), uncorrelated detection uncertainty σ du = 0.21% [53,54], and an overall normalization uncertainty σ o = 2%, have been included in the analysis. In the calculation of the signal events, a Gaussian energy smearing was assumed to account for the detector energy resolution with a width σ E /E ≈ 7%/ E[MeV] [53].
The Daya Bay Reactor Neutrino experiment analyzes the antineutrino flux produced by six reactor cores at the Daya Bay and Ling Ao nuclear power plants. The electron antineutrino oscillation probability is measured by eight identical antineutrino detectors (ADs). Two detectors are placed in each of the two near experimental halls of the experiment (EH1 and EH2), while the remaining four are located at the far experimental hall (EH3). Detailed studies on the antineutrino flux and spectra have been performed in order to determine the fission fractions (see Tab. 9 in Ref. [55]) as well as the thermal power (see Tab. I in Ref. [56]). Baseline distances range in ∼0.3 -1.3 km for the near experimental halls and ∼1.5 -1.9 km for the far hall. The Daya Bay collaboration analyzed data collected after 1958 days of running time [50] and reported the measurements sin 2 (2θ 13 ) = 0.0856 ± 0.0029 and |∆m 2 ee | = (2.522 +0.068 −0.070 ) × 10 −3 eV 2 . To obtain the oscillation parameters, our analysis uses the number of antineutrino events after background subtraction, considering the ratios of EH3 to EH1 and EH2 to EH1. Regarding the statistical methods, the Daya Bay collaboration has followed three different approaches (covariant approach, nuisance parameters and a hybrid approach). Consistent results can be obtained with the three methods and we have chosen to use nuisance parameters in our analysis. The uncertainties arising from the power and fission fractions at each of the 6 nuclear reactors are encoded in these nuisance parameters (σ r = 0.2% and σ f rac = 0.1%). In addition, characteristics of each detector, such as the differences in the running time or the efficiencies, have been accounted for in the simulation.
Other sources of uncertainties, such as shifts in the energy scale (σ scale = 0.6%), have also been included.
The results of our analyses of short-baseline reactor data are shown in Fig. 2. As can be seen in the figure, there is a total overlap between the parameter regions determined by RENO and Daya Bay, although the latter clearly dominates the measurement of the relevant oscillation parameters.
Note also that our results are almost identical for normal and inverted mass spectra, since these experiments are not sensitive to the mass ordering.

C. Atmospheric neutrino experiments
When cosmic rays collide with particles in the Earth's atmosphere, they start a particle shower which eventually creates the atmospheric neutrino flux. The energy of ν µ and ν e (and their antiparticles) produced in the atmosphere can range from a few MeV up to roughly 10 9 GeV, although only events up to ∼ 100 TeV are currently detectable. The energy of the atmospheric neutrinos relevant to oscillation studies, however, ranges from ∼ 0.1 GeV to ∼ 100 GeV. In our global fit we include data from Super-Kamiokande [57] and from IceCube DeepCore [58,59]. Since the largest part of the atmospheric neutrino flux is composed by ν µ and ν µ , and given that it is more difficult to identify electrons in the detector, the main channel used in current atmospheric neutrino experiments is ν µ → ν µ , which makes them mostly sensitive to the oscillation parameters θ 23 and ∆m 2 31 . Note, however, that the Super-Kamiokande experiment also detected a large sample of electron events from ν e appearance [57,60]. The results of this analysis, however, are not available outside of the collaboration. As a result, we do not analyze Super-Kamiokande atmospheric data ourselves, but only include the latest χ 2 -table made available by the collaboration [61].
For the current global fit, we update our analysis of DeepCore data. In addition to track-like events, the released experimental data now includes also shower-like events, increasing the number of events from roughly 6000 [62] to around 20000 [58,59]. The data analyzed correspond to 3 years of observations of the full sky, from April 2012 to May 2015. The details of the analysis are described in Ref. [59], and the full data set can be downloaded from Ref. [63]. Two data samples are provided: Sample A and Sample B, corresponding to the same data taking period but different cuts. For this analysis we have chosen Sample A. Several sources of systematic uncertainties are included in our analysis. They can be divided into detector-related and flux-related uncertainties. We account for neutrino scattering and absorption in the ice, and include several uncertainties related to the optical efficiencies. Concerning the atmospheric neutrino flux, we include systematic uncertainties on the ratio of neutrinos to antineutrinos, the ratio of electron to muon neutrinos, the spectral index, the ratio of vertically to horizontally incoming neutrinos and an overall normalization. The results of our analysis are depicted in Fig. 3, together with the ones from Super-Kamiokande. As in the reactor case, the regions allowed by the two experiments totally overlap. However, one can see that the mixing angle is slightly better measured by Super-Kamiokande, while DeepCore provides a more stringent result on the atmospheric mass splitting.

D. Accelerator experiments
Long-baseline accelerator neutrino experiments measure neutrinos which are created in particle accelerators. They originate in meson decays. The mesons, typically pions and kaons, are created in the accelerator and then focused into a beam. Next, they decay into muon-neutrinos, while a beam dump absorbs the ones which do not decay. Using different polarities of the focusing horns one can separate mesons from antimesons, resulting in a mostly pure beam of neutrinos or antineutrinos. Note, however, that creating a really pure beam is not possible, and there will always be a background contamination of so-called "wrong-sign" neutrinos. The long-baseline This makes them sensitive to the oscillation parameters ∆m 2 31 , θ 23 , θ 13 , δ and, in principle, also to the neutrino mass ordering. In our global fit, we use data from several long-baseline experiments: NOνA [64], T2K [65,66], MINOS [67] and K2K [68].
The T2K collaboration has presented an updated analysis of neutrino and antineutrino data, These results improve their former ones [69], allowing them now to exclude CP-conserving values of δ at close to 3σ confidence level.
On the other hand, NOνA is now including also antineutrino data in their analysis. They can not exclude any value of δ in a significant way at the moment and, therefore, its preference for δ = 0 might be due to a statistical fluctuation. Note also that the sensitivity of both experiments to CP violation gets reduced when relaxing the reactor prior on θ 13 , used in the analyses of both experiments.
In order to perform our analysis, we extract the relevant data for each experiment from the corresponding reference. We simulate the signal and background rates using the GLoBES software [71,72]. For the energy reconstruction we assume Gaussian smearing. We include bin-to-bin efficiencies, which are adjusted to reproduce the best-fit spectra reported in the corresponding references. Finally, for our statistical analysis we include systematic uncertainties, related to the signal and background predictions, which we minimize over. The results of our analysis (without a prior on θ 13 ) are presented in Figs. 4 and 5. We find that T2K and NOνA measure the atmospheric parameters θ 23 and |∆m 2 31 | rather well and with similar sensitivity. Note, however, that T2K shows a better sensitivity to θ 13 and δ for normal neutrino mass ordering, as indicated by the 90% C.L. closed regions in the left panel Fig. 5. For inverted ordering, both experiments show similar sensitivity to δ, although T2K does somewhat better on the θ 13 measurement than NOνA.
In any case, these results are not competitive with short-baseline reactor experiments, discussed above. We also show the results from our analysis of MINOS data [73,74], which still contributes to the determination of |∆m 2 31 |, as seen in Fig. 4. Unfortunately, in this case there is no sensitivity to θ 13 and δ. The same applies to the pioneering K2K experiment [75], included in our global fit as well, but with a sensitivity to the oscillation parameters which has been overcome by the more recent long-baseline accelerator experiments.

III. RESULTS FROM THE GLOBAL FIT
In the previous section, we have presented the individual results of our neutrino data analysis, obtained sector-by-sector. In this section, we shall describe the results obtained by combining all previous data into our global neutrino oscillation fit. We will first briefly discuss the main contributions to the well-measured parameters, and then enter into more detail in the discussion of the remaining unknowns of the three-neutrino picture.

A. Well-measured oscillation parameters
So far the solar parameters θ 12 and ∆m 2 21 have only been measured by KamLAND and the solar neutrino experiments, and they have already been discussed in Sec. II A. After combining with data from other experiments the determination of the solar parameters improves further, due to a better determination of θ 13 , but the effect is not very visible. The future reactor experiment JUNO is expected to measure the solar parameters with great precision [76]. In contrast, the measurement of the remaining oscillation parameters emerge from the combinations of several data sets, as seen values, as preferred by reactor data. Note also that this combined analysis shifts the best fit value of sin 2 θ 23 to the second octant. A more distinctive feature appears when combining LBL with reactor data. As expected, this combination results in a much more restricted range for θ 13 , and therefore the partial breaking of the θ 23 -θ 13 degeneracy, arising from the LBL appearance data, see the red lines in Fig. 7. Finally, when combining all data, we obtain the green lines in Fig. 6

C. The CP phase δ
We now discuss the measurement of the CP-violating phase, δ. This phase induces opposite shifts in the ν µ → ν e and ν µ → ν e oscillation probabilities and, therefore, information on this parameter can be obtained by analyzing neutrino and antineutrino oscillation data in the appearance channels. Note, however, that the separate analysis of neutrino and antineutrino channels can not provide, at present, a sensitive measurement of δ [77]. The CP phase can therefore be measured by the long-baseline accelerator experiments T2K and NOνA, and also by Super-Kamiokande atmospheric neutrino data, see the black and blue lines in the lower right panel of Fig. 6. In addition to Fig. 6, in Fig. 8 we show the ∆χ 2 profiles for the CP-violating phase δ as obtained from the analysis of data from T2K (blue) and NOνA (red), the combination of all long-baseline data (black) and the result from the global fit (green). For normal neutrino mass ordering (left panel), a new tension arises between the determinations of δ obtained from T2K and NOνA data. Indeed, the analysis of NOνA results shows a preference for δ ≈ 0, disfavoring the region around δ ≈ 1.5π, where we encounter the best fit value for T2K. This does not happen for inverted ordering (right panel), for which NOνA shows better sensitivity to δ and also an excellent agreement with T2K.
Note that this behavior is due to the antineutrino data sample recently released by NOνA, and it is the reason why our sensitivity to δ in the current global fit is worse than it was in Ref. [1]. The inclusion of reactor data can help to improve the determination of δ, due to the existing correlation between the CP phase and θ 13 . This is illustrated in Figs. 6 and 9. From the global combination, we obtain the best fit value for the CP phase at δ = 1.20π (1. denote the best fit value obtained from the global fit for normal and inverted ordering, respectively.

D. The neutrino mass ordering
Finally, in this subsection, we present the results of our present analysis of the neutrino mass ordering issue. Combining all neutrino oscillation data, we obtain a preference for normal mass ordering with respect to inverted with ∆χ 2 = 9.1. This corresponds to a 3σ preference in favor of NO. This preference comes from several contributions, which we will discuss in the following. Our independent analyses of NOνA and T2K data both slightly prefer NO with ∆χ 2 ≈ 1.5. Such a small value is expected, due to the rather small matter effects present in the neutrino propagation over the corresponding baselines. However, if we combine all long-baseline accelerator data we obtain ∆χ 2 = 0.5 in favor of NO. This reduction appears as a consequence of the tension in the measurement of δ by T2K and NOνA, as discussed in the previous subsection. Since the tension appears only in normal ordering, the minimum χ 2 from the combined long-baseline analysis for this ordering is worse than the sum of the individual T2K and NOνA fits. This results in a lower value for ∆χ 2 = χ 2 min (IO) − χ 2 min (NO). If we now perform a combined analysis of accelerator and reactor data, we obtain a preference for NO with ∆χ 2 = 4.5. This increase comes from the difference in the measurements of sin 2 θ 13 in accelerator and reactor experiments. As shown in Figs. 7 and 9, the values of θ 13 preferred by accelerator experiments are slightly different for normal and inverted ordering: sin 2 θ 13 = 0.0275 and sin 2 θ 13 = 0.029, respectively. The slightly larger values preferred in IO are then further penalized upon combining with reactor data (that prefer sin 2 θ 13 0.022) than the NO counterparts, what results in a larger value for ∆χ 2 .
On the other hand, the atmospheric neutrino results from the Super-Kamiokande and Deep-Core experiments show some sensitivity to the neutrino mass ordering on their own. From Super-Kamiokande data alone (neither imposing a prior on θ 13 nor combining with data from reactor experiments), there is already a preference for normal mass ordering with ∆χ 2 ≈ 3.5, while Deep-Core gives ∆χ 2 ≈ 1.0. Combining these atmospheric results with long-baseline accelerator data, this preference grows up to ∆χ 2 = 6.7. From Fig. 7 we also see that, after this combination, the measurement of sin 2 θ 13 agrees better with the reactor measurement than in the case of longbaseline data alone. However, while the best fit values for normal ordering nearly coincide, there is still a small tension in inverted ordering. Therefore, after the global combination with data from reactor experiments, we obtain the final preference of ∆χ 2 = 9.1, corresponding to a significance of 3σ. As for the CP phase δ, the current preference for normal mass ordering is lower than reported in Ref. [1]. The explanation is the same as before, namely the tension in the combined analysis of T2K and NOνA for normal mass ordering, due to the different preferred values for δ. Therefore, any development on this tension will affect the sensitivity of neutrino oscillation data to the mass ordering.

IV. BAYESIAN ANALYSIS OF NEUTRINO OSCILLATION DATA
In this section we turn to the discussion of our Bayesian analysis of neutrino oscillation data.

A. The Bayesian method
In order to perform a Bayesian analysis of neutrino oscillation data, we convert the χ 2 functions described in the previous sections into a likelihood, using the expression The analysis is performed using MontePython  and frequentist (dashed) determinations obtained for normal (blue) and inverted (magenta) ordering. Note that inverted ordering results are normalized with respect to the minimum χ 2 of inverted ordering.

B. Oscillation parameter results
In order to obtain a Bayesian comparison of the NO and IO spectra, we have to perform numerical analyses which we also use to produce Bayesian neutrino oscillation parameter determinations.
While the likelihood is the same, converted from the χ 2 discussed in the previous sections according to Eq. (1), minor differences appear between the frequentist and Bayesian analyses, which are shown in Fig. 10. In the figure, we show the frequentist one-dimensional ∆χ 2 profiles (dashed lines) and the marginalized Bayesian posterior probabilities P (x) (solid lines), converted into an effective χ 2 using where x represents any one of the six oscillation parameters. In the figure, we show NO (blue) and IO (magenta), normalizing in both cases with respect to the best fit for the same ordering of the spectrum. Apart for the normalization in the IO case, the dashed lines are the same we show in the global fit summary in Fig. 15. As we can see, most of the posterior distributions are exactly the same as the frequentist profiles. Minor differences only appear in sin 2 θ 23 and δ, but none of the conclusions of the paper are changed.
In the first line of Tab. II and in Fig. 12, we report the significance of the Bayesian comparison of NO and IO. As we can see, the significance decreased slightly with respect to the previous results obtained in Ref. [84], due to the mismatch in the determination of δ by T2K and NOνA, as already explained. Neutrino oscillation data alone give now ln B NO,IO = 5.05 ± 0.04, corresponding to a 2.7σ probability for a Gaussian variable.

V. ABSOLUTE SCALE OF NEUTRINO MASSES
Since neutrino oscillations depend only on the mass splittings between the neutrino mass eigenstates, in order to probe the absolute scale of the neutrino mass, other experiments are required.
In this section, we discuss the status of current probes of the absolute neutrino mass: kinematic measurements through the observation of the energy spectrum of tritium β decay, neutrinoless double β decay, plus cosmological constraints.
At the moment, the strongest limits on the effective electron antineutrino mass m β are set by the KATRIN experiment [86], which obtained the upper limit m β < 1.

B. Neutrinoless double β decay
If neutrinos are Majorana particles, one expects that a neutrinoless variety of double beta decay in which no neutrinos are emitted as real particles should take place. This is called neutrinoless double β decay (0νββ) and, if it is ever detected, it implies the Majorana nature of neutrinos [90].
The non-observation of 0νββ can then be used to set complementary limits on the neutrino mass scale. The decay amplitude is given as where m j are the Majorana masses of the three light neutrinos. Notice the absence of complex conjugation of the lepton mixing matrix elements. These contain the new CP phases [91,92] characteristic of the Majorana neutrinos (note that the Dirac phase that appears in neutrino oscillations does not appear in the 0νββ amplitude). This is manifest within the original symmetrical parametrization of the lepton mixing matrix [91] in which the Majorana phases are treated symmetrically, each one associated to the corresponding mixing angles * * .
One finds that 0νββ probes constrain the half-life T 0ν 1/2 (N ) of the isotope involved in the decay (see e.g. [94] for a recent review). Assuming that the dominant mechanism responsible for the 0νββ events is light neutrino exchange, one finds constraints on T 0ν 1/2 (N ) which can be translated into bounds on the effective Majorana mass m ββ .

The conversion between half-life and effective Majorana mass is
where m e is the electron mass, G N 0ν is the phase space factor and M N 0ν is the nuclear matrix element for the decay. The latter two terms depend on the isotope under consideration.
The latter expression was proposed in [98], while the former two have been obtained using the information from [95] and [96], respectively, and using the general fitting formula proposed in [98].
When performing analyses including constraints from neutrinoless double β decay probes, we marginalize over the two Majorana phases in their allowed range. Moreover, in order to take into account the theoretical uncertainties in the calculation of the nuclear matrix elements, we vary them within the ranges M 130 Te 0ν M 136 Xe 0ν which correspond to the proposed 1σ range from [99] (see Tab. 6 in that reference).
This way, we find that the bounds on the half-life T 0ν 1/2 (N ) from the three experiments under consideration imply the following upper limits on the effective mass: m ββ < 104 − 228 meV by Gerda [95], m ββ < 75 − 350 meV by CUORE [96] and m ββ < 61 − 165 meV by KamLAND-Zen [97], respectively, where the lower (upper) values correspond to the most aggressive (conservative) choices for the nuclear matrix elements.

C. Cosmological probes
Stronger, though model-dependent (see e.g. [100]), are the limits on the sum of the neutrino masses provided by cosmological observations. They arise mainly from the combination of Cosmic Microwave Background (CMB) and Baryon Acoustic Oscillation (BAO) measurements (see e.g. [101]). Due to the anticorrelation between the sum of the neutrino masses and the Hubble parameter, it is also interesting to consider constraints on the latter quantity.
In our analyses, we consider the most recent observations of the CMB spectrum by Planck [102,103], which measures the temperature and polarization spectra in a wide range of multipoles through the respective two-point correlation functions [104] and the lensing [105] potential through the four-point correlation function. We include BAO observations from the 6dF [106], SDSS DR7 Main Galaxy Sample (MGS) [107] and BOSS DR12 [108] galaxy redshift surveys. Bounds on the expansion of the universe quantified by the Hubble parameter H(z) also come from measurements at z = 0.45 [109]. Constraints from observations of Type Ia Supernovae are also taken into account, by means of the Pantheon sample [110]. Here we will indicate by "Cosmo" the combination that includes Planck CMB temperature, polarization and lensing spectra, BAO measurements, H(z) observations and Supernovae luminosity distance data. Finally, the most recent determination of the Hubble parameter today, H 0 = 74.03 ± 1.42 km/s/Mpc from [111], is also included in some cases in order to illustrate the impact of the Σm ν − H 0 degeneracy.
The calculation of predicted cosmological observables is performed using the Boltzmann solver code CLASS [112][113][114]. Our fiducial cosmological model is a minimal extension of the ΛCDM model, which is described by the usual six free parameters. Namely, the baryon and cold dark matter physical densities Ω b h 2 and Ω c h 2 , the angular size of the sound horizon at last-scattering θ s , the optical depth to reionization τ and the amplitude and tilt of the primordial scalar power spectrum A s and n s . We fix the number of ultra-relativistic species to zero and we add three massive neutrinos, each with its own mass. Such masses are derived from the lightest neutrino mass m lightest and the two mass splittings before performing the cosmological calculations. The total neutrino mass in the two orderings reads as Σm NO ν = m lightest + m 2 1 + ∆m 2 21 + m 2 1 + ∆m 2 31 , D. Global results on the neutrino mass scale and mass ordering Once we introduce in our analyses the constraints from β decay, neutrinoless double-β decay and cosmology, we are able to obtain upper bounds on the absolute scale of neutrino masses. Here we present the results in terms of m lightest , which is the quantity we can compare in an easier way when discussing the various probes. In Fig. 11 we report the constraints on m lightest in a priorindependent way, using the method of Refs. [115][116][117] and recently revived in [118]. The plotted function R(m lightest , 0) ≡ p(m lightest )/π(m lightest ) p(m lightest = 0)/π(m lightest = 0) (14) shows the ratio between the posterior p(m lightest ) and the prior π(m lightest ) normalized with respect to the same ratio computed at m lightest = 0, for different data sets, and comparing normal (blue) to inverted (magenta) mass ordering. Such quantity, which has the property of being independent of the shape and normalization of the prior π(m lightest ), is statistically equivalent to a Bayes factor between a model where m lightest has been fixed to some value and one where m lightest is equal to zero. Since the considered measurements are insensitive to the value of m lightest when it is very small, the function R is expected to be equal to one for small m lightest † † , while it decreases when large values of m lightest become disfavored. In the same way as a Bayes factor, we can compare the constraining power of different data sets by means of the Jeffreys' scale. The horizontal lines in Fig. 11 show the values ln R = −1, −3, −6, which separate regions where the significance is none, weak, moderate and strong, according to the Jeffreys' scale we adopt.
Analyzing the figure, we can see that the results obtained from the β-decay data (solid) are completely insensitive to the mass ordering, and provide the weakest constraints on m lightest nowadays.
We must remember, however, that β decay measurements provide the most robust constraints on the absolute scale of neutrino masses, as they are completely model independent. Neutrinoless double-β decay bounds (dashed), which only apply to Majorana neutrinos, provide stronger Refs. [119][120][121] after considering very similar cosmological observations. Namely, in Ref. [103] it is quoted Σm ν < 0.12 eV from CMB temperature, polarization, lensing and BAO observations. This is due to the different lower prior assumed in our global fit: while the Planck collaboration just assumes a physical prior on the sum of the neutrino masses, i.e. Σm ν > 0, here the lower value of the prior is determined by neutrino oscillation experiments assuming m lightest = 0, and it is therefore different for NO and IO. See Ref. [120] and Ref. [100] for an assessment of the changes in the upper bounds on the total neutrino mass after taking into account neutrino oscillation information and the uncertainty on the underlying cosmological model, respectively.
Considering absolute neutrino mass measurements also affects the preference in favor of NO  the ∼ 0.06 eV that apply for NO, which is therefore preferred.
Finally, we report in Fig. 13 the allowed regions at 1, 2σ for the mass parameters m lightest , m β , m ββ and Σm ν , obtained considering the OSC+Cosmo data set. We do not show the regions allowed by β-decay and 0νββ experiments as they are outside the scale of the respective effective parameter: in other words, the constraints that we obtain considering cosmological data are much tighter than those obtained from terrestrial experiments (see Secs. V A and V B).

VI. SUMMARY OF THE GLOBAL FIT
In this study, we have first analyzed global data coming only from neutrino oscillation experiments. These hold in complete generality. In a second step, we have combined the oscillation results with direct neutrino mass probes such as β decay, neutrinoless double β decay and cosmological observations. The latter has its caveats, for example, 0νββ restrictions apply only to Majorana neutrinos, while cosmological considerations suffer from a higher degree of model-dependence.   the new best fit value slightly smaller. While the determination of the solar mass splitting ∆m 2 21 remains unchanged, the best fit value is also smaller than that obtained before. Due to new shortbaseline reactor data from Daya Bay and RENO, we obtain an improved measurement of sin 2 θ 13 and a larger best fit value of sin 2 θ 13 = 0.0225. Also the best fit value of the atmospheric mass splitting ∆m 2 31 is now larger, although the size of the allowed interval remains roughly the same for both mass orderings. Regarding the atmospheric angle, we obtain the best fit value sin 2 θ 23 = 0.566, in the second octant for both mass orderings. Indeed, the preference for the second octant obtained here is stronger than in Ref. [1], and lower octant solutions are now disfavored with ∆χ 2 ≥ 4.3 (5.1) for normal (inverted) ordering.
However, for the case of the CP-violating phase δ, we obtain a weaker result in comparison with our previous global fit [1], due to the mismatch in the value of δ extracted by T2K and NOνA. The best fit is obtained for δ = 1.2π (1.54π) for normal (inverted) ordering. Concerning the CP-conserving values, δ = 0 is disfavored with ∆χ 2 = 9.1 (15.0) for NO (IO), while δ = π, remains allowed with ∆χ 2 = 2.1 for NO and it is excluded with ∆χ 2 = 17.4 in IO. This is due to the fact that the aforementioned mismatch in the extracted value of δ by the T2K and NOνA experiments only occurs for normal ordering. Indeed, for inverted mass ordering both experiments prefer values close to maximal CP violation, with δ ≈ 1.5π. Finally, the very same mismatch reduces the statistical significance of the preference for NO from the 3.4σ obtained in Ref. [1] to the 3.0σ derived in the frequentist analysis presented in this work.
Our results from the Bayesian analysis are summarized in Fig. 12   ordering from neutrino oscillation data alone. This would correspond to a Gaussian preference of 2.7σ. The determination of neutrino oscillation parameters shows also an excellent agreement among the frequentist and Bayesian analyses.
While the inclusion of β-decay data does not change the former Bayes factor, data from 0νββ experiments mildly increase that figure to ln B = 5.46 ± 0.22, still indicating moderate preference for normal neutrino mass ordering. Note, however, that this improvement would only apply if neutrinos are Majorana particles. The combination with data from cosmological observations is independent of the neutrino nature and leads to a Bayes factor of ln B = 6.09 ± 0.30, strongly preferring normal neutrino mass ordering and corresponding to a statistical significance of 3.1σ.
Finally, when we also include in the cosmological observations a prior on the Hubble constant we obtain ln B = 7.05 ± 0.36, corresponding to a preference for NO of 3.3σ.
Concerning the cosmological limits on the sum of the neutrino masses, the tightest 2σ bound we obtain here is m ν < 0.13 (0.15) eV for NO (IO) taking into account CMB temperature, polarization and lensing measurements from the Planck satellite, BAO observations, H(z) information and Supernovae Ia data. These limits are slightly weaker than those existing in the literature due to our prior, that takes into account neutrino oscillation results as an input.
Overall, we have seen that the determination of some of the neutrino parameters has improved thanks to new oscillation data, while the determination of δ and the neutrino mass ordering has worsened, due to a new tension in current long-baseline accelerator measurements. We have also seen that the inclusion of non-oscillation data, especially from cosmological observations, enhances the preference for normal neutrino mass ordering from "moderate" to "strong".
In summary, neutrino oscillation parameters are currently measured with very good precision.
In the upcoming years these accurate measurements will further improve, allowing for better sensitivities to New Physics effects, which may show up as sub-leading effects in the neutrino oscillation probabilities.