The unfinished fabric of the three neutrino paradigm

In the current 3nu paradigm, flavor oscillations probe 3 mixing angles (theta_12, theta_23, theta_13), one CP phase delta, and two squared mass differences delta m^2>0 and Delta m^2, where sign(Delta m^2)=+ (-) for normal (inverted) ordering. Absolute nu masses can be probed by the effective m_beta in beta decay, by the total mass Sigma in cosmology and, if neutrinos are Majorana, by another effective m_{beta beta} in 0nu2beta decay. Within an updated global analysis of (non)oscillation data, we constrain these 3nu parameters, both separately and in selected pairs, and highlight the concordance or discordance among different constraints. Five oscillation parameters (delta m^2, Delta m^2, theta_12, theta_23, theta_13) are consistently measured, with an overall accuracy ranging from ~1% for Delta m^2 to ~6% for sin^2(theta_23) (due to its octant ambiguity). We find overall hints for normal ordering (at 2.5 sigma), as well as for theta_23<pi/4 and for sin(delta)<0 (both at 90% C.L.), and discuss some tensions among datasets. Concerning nonoscillation data, we include the recent KATRIN constraints on m_beta, and we combine the latest 76-Ge, 130-Te and 136-Xe bounds on m_{beta beta}, accounting for NME covariances. We also discuss some variants related to CMB anisotropy and lensing data, which may affect cosmological constraints on Sigma and hints on sign(Delta m^2). The default option, including all Planck results, irrespective of the lensing anomaly, sets upper bounds on Sigma at the level of ~10^-1 eV, and further favors normal ordering up to ~3 sigma. An alternative option, that includes recent ACT results + other independent results (from WMAP and selected Planck data) globally consistent with standard lensing, is insensitive to the ordering but prefers Sigma ~(few) x 10^-1 eV, with different implications for m_beta and m_{beta beta} searches. (Abridged)

Constraints on the oscillation parameters (δm 2 , ∆m 2 , θ ij , δ) and on the nonoscillation observables (Σ, m β , m ββ ) have been explored in several global neutrino data analyses, including our previous work [10] and the more recent papers [11][12][13], plus the preliminary contribution in [14]. In particular, the analyses in [12][13][14] are based on largely common datasets, based on updated information from the Conference Neutrino 2020 [15]. As a result, a solid fabric for the 3ν paradigm emerges from convergent measurements of five oscillation parameters (θ 12 , θ 23 , θ 13 , δm 2 , |∆m 2 |), with an overall accuracy ranging from ∼ 1% for |∆m 2 | to ∼ 6% for sin 2 θ 23 (dominated by the so-called θ 23 octant degeneracy). However, the fabric is still unfinished in ν oscillations, as far as the θ 23 octant, the mass ordering and the phase δ are concerned. In particular, a tension between recent long-baseline accelerator neutrino data (from T2K [17] and NOvA [18]) affects all these 3ν unknowns at the same time [12][13][14]16]. In addition, the Dirac-Majorana nature and the absolute neutrino mass scale remain undetermined in current nonoscillation searches, with Σ, m β and m ββ constrained at sub-eV scales but still consistent with null values [19].

arXiv:2107.00532v2 [hep-ph] 28 Sep 2021
In this work we discuss the status of the 3ν framework, including new data that have recently become available from both oscillation and nonoscillation searches, with particular attention to issues of concordance or discordance of various data sets, and to their implications on the unfinished fabric of the paradigm. By performing a global analysis of oscillation data, including the latest Super-Kamiokande atmospheric results made publicly available in 2021 [20], we find an indication for normal ordering at the level of 2.5σ, as well as 90% C.L. hints for θ 23 < π/4 and for sin δ < 0. We discuss the structure and interplay of such hints, especially in the light of the T2K and NOvA tension and of the complementarity among accelerator, atmospheric and reactor data. We surmise that further understanding of neutrino nuclear interactions may help to clarify some issues.
Concerning nonoscillation searches, we include the KATRIN 2021 data [21] that, for the first time, set sub-eV upper bounds on m β at 90% C.L. We analyze systematically all the latest 0νββ decay searches probing half lives T > 10 25 y (in 76 Ge,130 Te and 136 Xe), and translate them into m ββ bounds via correlated nuclear matrix elements. In the realm of cosmology and of its consensus ΛCDM model, increasing attention is also being paid to old and new data tensions, see e.g. [22][23][24]. In this context, we focus on the so-called A lens anomaly affecting Planck angular spectra (that show more lensing than expected in the ΛCDM model [25]), and we consider two possible options, leading to different implications for absolute mass observables. On the one hand, we revisit a previously considered "default" scenario [11], including all Planck results (irrespective of the A lens anomaly), that sets stringent upper bounds on Σ at the level of ∼ 10 −1 eV, and further favors normal ordering, raising its overall preference to ∼ 3σ. On the other hand, we discuss an "alternative" option that makes use of the recent ACT CMB polarization data release 4 (ACTPol-DR4) [26] that is consistent with standard lensing, also in combination with WMAP 9-year data (WMAP9) [27] and selected data from Planck [25,28,29]; such option is insensitive to the mass ordering and prefers Σ ∼ few × 10 −1 eV, with different implications for m β and m ββ searches.
Building upon our previous works [10,11] we elaborate upon these recent topics as follows: In Sections II and III we update and discuss the analysis of oscillation and nonoscillation data, respectively. We pay particular attention to relevant correlations among various observables, and to some emerging tensions among different data sets. We provide a brief synthesis of oscillation and nonoscillation results in Sec. IV.

II. OSCILLATION DATA, ANALYSIS METHODS AND RESULTS
In this section we introduce recent oscillation data that were not included in our previous work [11], together with the methodology used for their analysis, in the light of some emerging issues in precision oscillation physics. We then discuss the resulting constraints on the parameters (δm 2 , ∆m 2 , sin 2 θ ij , δ), both separately and in selected pairs, highlighting the concordance, discordance and complementarity of various datasets.

A. Oscillation data update
Three-neutrino oscillations are currently constrained by experiments using long-baseline (LBL) accelerator, solar, long-baseline reactor (KamLAND), short-baseline (SBL) reactor and atmospheric neutrinos. With respect to [11], we update some of these datasets as follows. LBL accelerator data, in the form of neutrino and antineutrino energy spectra for flavor disappearance and appearance channels, are taken from the presentations of T2K [17] and NOvA [18] at Neutrino 2020 [15] (and subsequent conferences). Such spectra are endowed with statistical (Poisson) errors and systematic (normalization and energy scale) uncertainties, as well as with oscillation-independent backgrounds, in a modified version of GLOBES [30]. Solar neutrinos data from the Super-Kamiokande-IV 2970-days run (SK-IV energy spectrum and day-night asymmetry) are taken from the presentation at Neutrino 2020 [31], while the input solar model BP16-GS98 [32] is unchanged. Concerning SBL reactors, we update from Neutrino 2020 the RENO [33] and Double Chooz data [34], while the Daya Bay data [35] are unchanged. Note that Daya Bay and and RENO measure both θ 13 and ∆m 2 , while the latter parameter is not significantly constrained by Double Chooz. IceCube-DeepCore (IC-DC) atmospheric data are taken as in [11]; a new IC-DC data release is expected in the near future [37]. Finally, our oscillation dataset is completed by the recent SK-IV atmospheric results [31,36], included through the χ 2 map recently made available by the collaboration [20].

B. Analysis method and emerging issues
We adopt the methodology proposed in [38], see also [10,11]. In particular, we start with the combination of solar, KamLAND and LBL accelerator neutrino data, that represent the minimal dataset sensitive to all the oscillation parameters (δm 2 , ∆m 2 , θ ij , δ). We then add SBL reactor neutrino data, that sharpen the constraints on (|∆m 2 |, θ 13 ) and indirectly affect the parameters (θ 23 , δ) and sign(∆m 2 ) via correlations. We add atmospheric neutrino data at the end, for two reasons: (a) they provide rich but rather entangled information on the parameters (∆m 2 , θ 23 , θ 13 , δ); (b) their χ 2 (∆m 2 , θ 23 , δ) maps assume an input on (δm 2 , θ 12 , θ 13 ) from the combination of solar, KamLAND and SBL reactor data. Finally, a frequentist approach based on χ 2 functions is used for all datasets. Best fits are obtained by χ 2 minimization, while allowed regions around best fits are expanded in terms of "number of standard deviations" N σ = ∆χ 2 . In particular, two-dimensional contours are shown for N σ = 1, 2 and 3, which, for a χ 2 distribution with two degrees of freedom, correspond to C.L. of 39.35%, 86.47% and 98.89%, respectively. Their one-dimensional projections provide the N σ ranges for each parameter, corresponding to C.L. of 68.27%, 95.45% and 99.73%, respectively. The difference ∆χ 2 IO−NO between the minima in IO and NO may-or may not-be accounted for, when reporting fit results; these two options will be clearly distinguished in each context.
We briefly discuss some issues arising in global data analyses, in the era of increasingly precise measurements and of growing sensitivity to subleading effects. Data fits usually involve the comparison of experimental event rates R expt with their theoretical predictions R theo where, from left to right, the integrands represent the source flux of ν α , the probability of ν α → ν β oscillations, and the interaction cross section, detector resolution and efficiency for ν β events. Some factors may be differential functions that need multiple integrations or convolutions, as alluded by the cross product ⊗. Integrands are endowed with various uncertainties, that may be shared by (i.e., correlated among) various rates R in the same or different experiment(s). In solar neutrino searches, all these features can be accurately implemented to a large extent [39]. Also short-baseline reactor experiments (Daya Bay, RENO, Double Chooz) generally provide enough public information to allow reproducible analyses, although a more precise joint analysis by the different collaborations (accounting for minor correlated uncertainties) would be desiderable [40]. More relevant issues arise in the context of long-baseline accelerator searches, currently carried out by T2K [17] and NOvA [18]. Their event spectra are usually given in terms of a "reconstructed" (unobservable) neutrino energy E rec ν , that is processed from observable event energies at far detectors, through models (for some R theo integrands) constrained by near-detector data. In principle, both T2K and NOvA should share a common theoretical model for σ (and to some extent for Φ), leading to possible correlations among their uncertainties for E rec ν (affecting ∆m 2 ) and for the event rates R theo (E rec ν ) (affecting θ 23 and θ 13 ). In practice, the adopted models are different, and possible covariances are ignored. A joint analysis planned by the T2K and NOvA collaborations [41] might shed light on these issues. Finally, in current atmospheric neutrino searches at SK and IC-DC, the data processing and analysis are too complex to be reproducible by external users with an acceptable accuracy. Oscillation results are given in terms of public χ 2 maps that, when summed up, cannot account for known covariances, such as those related to the (common, in principle) input models for Φ and σ. Once again, joint analyses or in-depth comparisons of data by different atmospheric ν experiments would be desirable [42]. We surmise that global oscillation analyses could benefit from a better control of those systematics that are shared by different experiments (such as model uncertainties for Φ α and σ β ), but whose correlations are not yet properly accounted for. See also the remarks in Sec. II E.

C. Results on single oscillation parameters
In this Section we present the constraints on the six oscillation parameters (δm 2 , ∆m 2 , sin 2 θ ij , δ) for increasingly rich data sets. We explicitly account for the χ 2 difference between NO and IO, in order to show its variations. Figure 1 shows the results for the combination of solar and KamLAND data (sensitive to δm 2 , sin 2 θ 12 , and sin 2 θ 13 ) with LBL accelerator data (mainly sensitive to ∆m 2 , sin 2 θ 23 , sin 2 θ 13 and δ), for both NO (blue) and IO (red). The latter mass ordering is slightly favored (by accelerator data) at the level of ∼ 1σ, as also discussed later. The parameters δm 2 and sin 2 θ 12 are rather precisely measured, with nearly linear and symmetrical (i.e., almost gaussian) uncertainties, and no significant difference between constraints in NO and IO. The parameters sin 2 θ 13 and sin 2 θ 23 are less accurately constrained. In particular, the two minima in sin 2 θ 23 reflect the θ 23 octant ambiguity in the ν µ → ν µ disappearance searches at LBL accelerators, inducing two correlated minima in sin 2 θ 13 via the leading amplitude of ν µ → ν e appearance (∝ sin 2 θ 23 sin 2 θ 13 ) [43]. The phase δ is poorly constrained, although it appears to be slightly favored around π in NO and around 3π/2 in IO, while it is disfavored around π/2 in both cases. Figure 2 shows the effect of adding SBL reactor data, which are sensitive to |∆m 2 | and sin 2 θ 13 . One can notice the strong reduction of the sin 2 θ 13 uncertainty, inducing also correlated changes on the relative likelihood of the lower and upper octant of θ 23 via ν µ → ν e appearance in LBL accelerators [43]. The synergy of SBL reactor and LBL accelerator data also helps to break mass-ordering degeneracies via independent measurements of ∆m 2 (see, e.g., [44]) and currently flips the fit preference from IO to NO (at the level of ∼ 1.3σ), together with an increase of the best-fit value of ∆m 2 with respect to Fig. 1. The preference for δ ∼ π (∼ 3π/2) in NO (IO) remains unaltered.   Figure 3 shows the effect of adding atmospheric ν data, which add further sensitivity to ∆m 2 (and to its sign), as well as to sin 2 θ 23 and δ. In particular, the inclusion of SK-IV data [20,31] corroborates the preference in favor of NO (at an overall level of ∼ 2.5σ), flips the θ 23 preference from the upper to the lower octant in NO (at ∼ 1.6σ) and also moves the best fit of δ slightly above the CP-conserving value π (disfavored at ∼ 1.6σ). The latter hints in favor of θ 23 < π/4 and on δ > π, currently emerging in NO at the statistical "threshold of interest" of 90% C.L., represent interesting updates with respect to previous global analyses not including SK-IV atmospheric data [12][13][14]. Table I reports a numerical summary of the same information shown in Fig. 3, for the separate cases of NO and IO (whose χ 2 difference is reminded in the last row). The two squared mass splittings ∆m 2 and δm 2 are measured with a formal 1σ accuracy of 1.1% and 2.3%, respectively. The mixing parameters sin 2 θ 13 , sin 2 θ 12 and sin 2 θ 23 are measured with an accuracy of ∼ 3%, 4.5%, and ∼ 6%, respectively. The latter uncertainty is largely affected by the θ 23 octant ambiguity; if one of the two quasi-degenerate θ 23 options could be removed, such uncertainty would be reduced by factor of ∼ 2 in both NO and IO.
Summarizing, five oscillation parameters are known with (few) percent accuracy, while only some hints emerge about the remaining three oscillation "unknowns". In particular, we find a preference for NO at ∼ 2.5σ and, in such ordering, we also find a preference at 90% C.L. for θ 23 in the lower octant (with respect to the secondary best fit in the upper octant) and for δ 1.24π (with respect to the CP-conserving value δ = π). Conversely, maximal θ 23 mixing is disfavored at ∼ 1.8σ and the range δ ∈ [0, 0.77π] is disfavored at > 3σ in NO. I: Global 3ν analysis of oscillation parameters: best-fit values and allowed ranges at Nσ = 1, 2 and 3, for either NO or IO, including all data. The latter column shows the formal "1σ fractional accuracy" for each parameter, defined as 1/6 of the 3σ range, divided by the best-fit value and expressed in percent. We recall that ∆m 2 = m 2 3 − (m 2 1 + m 2 2 )/2 and that δ ∈ [0, 2π] (cyclic). The last row reports the difference between the χ 2 minima in IO and NO.

D. Results on selected pairs of oscillation variables
By studying selected pairs of variables we can gain further insights about current unknowns (the mass ordering, the octant of θ 23 and the CP phase δ), and appreciate their interplay with known features of 3ν oscillations. We discuss the pairs (sin 2 θ 12 , δm 2 ), (sin 2 θ 23 , sin 2 θ 13 ), (sin 2 θ 23 , |∆m 2 |), (sin 2 θ 23 , δ), as well as pairs of total ν e and ν e events (bi-event plots) as observed in the appearance channel by T2K and NOvA. Figure 4 shows the regions separately allowed by solar and KamLAND neutrino data in the plane charted by (sin 2 θ 12 , δm 2 ), assuming fixed sin 2 θ 13 = 0.02 and NO. The two regions were somewhat displaced in the past, leading to a < 2σ tension between the best-fit δm 2 values [3] (see, e.g., the analogous Fig. 4 in [38]). The current regions in Fig. 4 appear to be in very good agreement, largely as a result of a slightly smaller day-night asymmetry in SK-IV 2970-day solar data, shifting the solar δm 2 best fit upwards and closer to the KamLAND one [31]. We find that this shift does not alter the combined solar and KamLAND constraints on θ 13 , namely, sin 2 θ 13 0.014 ± 0.015 (see, e.g., Fig. 5 in [38]). Results for IO (not shown) would be almost identical for all parameters (δm 2 , sin 2 θ 12 , sin 2 θ 13 ). In conclusion, solar and KamLAND data are not only in very good agreement about the (ν 1 , ν 2 ) oscillation parameters (δm 2 , sin 2 θ 12 ), but are also consistent with the measurement sin 2 θ 13 0.02 at SBL reactors. Figure 5 shows the covariance of the pair (sin 2 θ 23 , sin 2 θ 13 ) for increasingly rich data sets, in both NO (top) and IO (bottom), with the corresponding χ 2 functions separately minimized for each mass ordering. The θ 23 octant ambiguity leads to two quasi-degenerate solutions at 1σ, that generally merge at ∼ 2σ. The leading appearance amplitude in LBL accelerators, scaling as sin 2 θ 23 sin 2 θ 13 , induces an anticorrelation between the two angles in the left panels: the higher θ 23 , the smaller θ 13 . In the middle panels, the results from SBL reactors only (represented by ±2σ error bars) tend to prefer slightly the upper-octant solution (with lower values of θ 13 ) in both NO and IO, as confirmed by the combination with LBL accelerators (continuous curves). In the right panels, however, adding atmospheric data (that include SK-IV [20,31]) flips the octant preference in NO, while confirming it in IO. We conclude that current hints about the θ 23 octant are still rather fragile.

FIG. 5:
Regions allowed in the plane (sin 2 θ 23 , sin 2 θ 13 ) for increasingly rich data sets: Solar + KamLAND + LBL accelerator data (left panels), plus SBL reactor data (middle panels), plus Atmospheric data (right panels). Top and bottom panels refer, respectively, to NO and IO as taken separately (i.e., without any relative ∆χ 2 offset). The error bars in the middle panels show the ±2σ range for θ 13 arising from SBL reactor data only. FIG. 6: As in Fig. 5, but in the plane (sin 2 θ 23 , |∆m 2 |). The error bars in the middle panels show the ±2σ range for |∆m 2 | arising from SBL reactor data only. Figure 6 shows the covariance of the pair (sin 2 θ 23 , |∆m 2 |). In this case, there is a marked preference of SBL reactor data for relatively "high" values of |∆m 2 | (±2σ error bars), as compared with LBL accelerator data. A compromise is more easily reached for NO, featuring a smaller difference between the ∆m 2 values derived from SBL reactor and LBL accelerator data. This explains why the preferred mass ordering flips from inverted to normal, when passing from Fig. 1 to Fig. 2; see also [12,13]. For the same reason, maximal values of θ 23 (corresponding to the lowest values of |∆m 2 | allowed by LBL accelerator data) are slightly more disfavored by adding SBL reactor data. The overall preferences for NO and for nonmaximal θ 23 are confirmed by atmospheric data that, however, move the best fit in NO from the upper to the lower octant. We emphasize that SBL reactor data, despite having no direct sensitivity to sign(∆m 2 ) and θ 23 , contribute to constrain (via covariances) these two variables, in combination with other datasets. Figure 7 shows the covariance of the pair (sin 2 θ 23 , δ). The octant ambiguity leads to two quasi-degenerate best fits, surrounded by allowed regions that merge at 2σ or 3σ. In IO there is rather stable preference for the CP-violating case δ 3π/2 in all data combinations, with no significant correlation with θ 23 . In NO the allowed δ range is always larger, and includes the CP-conserving case δ π at 2σ; moreover, a slight negative correlation between δ and sin 2 θ 23 emerges when adding SBL reactor data. It is difficult to trace the origin of these null or small covariances, since the interplay between δ and sin 2 θ 23 (and with sin 2 θ 13 ) is rather subtle, see e.g. [45,46]. In any case, the negative correlation emerging in NO slightly amplifies the effect of adding atmospheric neutrino data, that prefer both δ ∼ 3π/2 and the lower octant of θ 23 , thus disfavoring δ π in a synergic way. Quantitatively, we find that the CP-conserving value δ = π is disfavored at 90% C.L. (or ∼ 1.6σ, see Fig. 3), while recent analyses not including SK-IV atmospheric data allowed this value at < 1σ [12,13]. Although these covariance effects are admittedly small in current data, they are expected to grow with increasing statistics and accuracy in LBL accelerator experiments, whose results we comment in more detail through the so-called bi-event plots, derived from bi-probability plots.  The slanted ellipses represent the theoretical expectations for NO (blue) and IO (red), and for two representative values of sin 2 θ 23 : 0.45 (lower octant, thin ellipses) and 0.57 (upper octant, thick ellipses). The CP-conserving value δ = π and the CP-violating value δ = 3π/2 are marked as a circle and a star, respectively. Each gray band represents one datum with its ±1σ statistical error (from [17,18]); the combination of any two data provides a (black dashed) 1σ error ellipse, whose center is marked by a cross. See the text for details.
Bi-probability plots, charted by the ν µ → ν e and ν µ → ν e appearance probabilities in LBL accelerator experiments at fixed neutrino energy, display the cyclic dependence on δ through ellipses [47] and help to understand parameter degeneracies [48]. See [16] for a related discussion, in the context of recent T2K and NOvA data. After integration over energy (weighted by ν fluxes and cross sections), the probabilities can be converted into total number of appearance events and thus into bi-event plots, preserving elliptic shapes [49]. Such theoretical ellipses can be directly compared with the measured number of events; see, e.g., the presentations at Neutrino 2020 [15] by T2K [17] and NOvA [18]. Although we use the full energy spectra (and not their integrals) in our LBL accelerator data analysis, we think that bi-event plots can help to highlight some issues emerging in the comparison of current T2K and NOvA data. Figure 8 show the ν and ν appearance events for T2K and NOvA, in four possible combinations. The grey bands and the black ellipses represent the data with their 1σ statistical errors, while the colored ellipses represent the theoretical expectations for NO (blue) and IO (red), for two representative values of θ 23 in the lower octant (thin) or upper octant (thick). Two representative values of δ (π and 3π/2) are also marked on each ellipse.
We first consider the two experiments separately, as shown in the upper left panel for T2K (ν vs ν) and in the lower right panel for NOvA (ν vs ν). In T2K, the best agreement of theory and data is reached for NO; in this ordering, there is a clear preference for δ = 3π/2, and a slight preference for the upper octant of θ 23 . In NOvA, all the four theoretical ellipses are close to the experimental one, but the overlap is larger in NO; in this ordering, there is a preference for δ = π/2 with respect to 3π/2, with no significant distinction of the θ 23 octants.
We then rearrange exactly the same information (both data and predictions) by combining T2K and NOvA separately in the ν and ν channels, as shown in the upper right panel for ν, and in the lower left panel for ν. In both plots, the best agreement of data and theory is now reached for IO. In such ordering, there is a clear preference for δ = 3π/2, as well as for the upper (lower) octant in the ν (ν) channel.
In conclusion, Fig. 8 shows that, as far as the three oscillation unknowns are concerned (mass ordering, θ 23 octant, CP symmetry), separate and combined T2K and NOvA data provide us with different indications, signalling a "tension" between such results; see also the discussion in [12,13,16]. Ultimately, the tension reflects the fact that NOvA (T2K) observes relatively symmetric (asymmetric) rates of ν and ν in their current appearance data.

E. Remarks
The current hints about the three oscillation unknowns are less converging and more fragile than in the recent past [11], due to the T2K-NOvA tension. As for its origin, Fig. 8 shows that possible statistical fluctuations of the data (at the level of one or two standard deviations) might play a role. However, it makes sense to speculate if there is more than just statistics behind the tension. One possibility is to invoke nonstandard neutrino interactions, that would induce different effects along the T2K and NOvA baselines [50,51]. Barring new physics beyond the 3ν paradigm, we surmise that standard interactions of neutrinos in nuclei might also play a systematic role.
It is widely recognized that both the total and the differential neutrino cross sections in nuclei are not known accurately enough for the purposes of LBL accelerator experiments [52,53]. Roughly speaking, normalization and energy reconstruction uncertainties in ν µ cross sections affect the measurement of θ 23 and ∆m 2 , respectively, while systematics on the relative ν µ /ν e and ν/ν cross sections affect the LBL experimental sensitivity to θ 13 , δ and the mass ordering. Note, e.g., that a 1% systematic error on the reconstructed neutrino energy E rec ν is transferred to ∆m 2 via the leading ∆m 2 /E dependence of the oscillation phase. The formal ∼ 1% error on ∆m 2 emerging from the global fit (Table I) implicitly posits that energy reconstruction errors are known at sub-percent level and are independent in different experiments, which may be an optimistic representation of the current uncertainties.
With increasingly higher statistics and accuracy, cross-section systematics may start to affect both known and unknown oscillation parameters extracted from detailed energy spectra. Possible parameter biases have been shown to arise by swapping different cross section models in simulations of prospective LBL data [54][55][56]. We remark that a percent-level bias on ∆m 2 as measured at LBL accelerators, in comparison with the ∆m 2 measurement at reactors, might alter the current combined preference for NO (see Fig. 6 and related discussion). Although all these effects can be reduced by tuning interactions models to cross-section data from near detectors [17,18], as a matter of fact T2K and NOvA use two different such models, while no single model or neutrino generator can be currently tuned to agree with world cross-section data [53,57].
Summarizing, the global analysis of current data shows a subtle interplay between the known oscillation parameters (|∆m 2 |, sin 2 θ 23 ) and the three unknowns (δ, sign(θ 23 − π/4), sign(∆m 2 )), as discussed through Figs. 5-7. Although there are overall hints in favor of NO (at 2.5σ), CP violation (at 1.6σ) and lower θ 23 octant (at 1.6σ), the T2K-NOvA tension (Fig. 8) warrants some caution. Neutrino interaction systematics might affect all these parameters, in a way that escapes control in global fits by external users, since the complexity of the near-to-far analysis chain can be handled only by the experimental collaborations. It is thus encouraging that T2K and NOvA are planning a joint analysis [41]. In this context, we practically suggest that these two experiments try to swap interaction models or neutrino generators in their separate simulations, so as to gauge the relative size of cross-section systematics and tuning effects, before attempting a combined fit. In general, we suggest that experiments sharing potential relevant systematics (e.g., neutrino fluxes Φ α and interactions in water σ β for atmospheric neutrinos) collaborate on a detailed comparison of such uncertainties and possibly towards joint data analyses. Of course, experimental developments on Φ α and σ β should be accompanied by corresponding advances in nuclear theory.

III. NONOSCILLATION DATA, ANALYSIS METHODS AND RESULTS
In this section we introduce recent nonoscillation data that were not included in our previous work [11], together with the methodology used for their analysis in terms of the three absolute mass observables (m β , m ββ , Σ). We include the latest m β upper bounds from the KATRIN experiment [21] and introduce a method to combine the latest 0νββ constraints in terms of upper bounds on m ββ , including correlated uncertainties on their nuclear matrix elements. We also enlarge the ensemble of cosmological datasets presented in [11] in the light of the current lively discussion on tensions in cosmology [22][23][24], that might suggest possible inconsistencies among different data (or alterations of the standard ΛCDM model). In particular, we consider an "alternative" dataset, that is exempt from the Planck lensing anomaly, at the price of being restricted to recent cosmic microwave background (CMB) anisotropy observations from ACTPol-DR4 [26] plus WMAP9 [27] and selected Planck data [25,28,29]. The alternative option provides a nonzero best fit and more conservative upper bounds on Σ, as compared with the "default" dataset described in [11]. We shall highlight the different sensitivities to Σ and to the mass ordering in the default and alternative options, as examples of admissible cosmological variants with rather different impact on global neutrino data analyses.

A. Single beta decay and constraints on m β
The KATRIN β-decay experiment has recently released the results of the second campaign of measurements [21]. In combination with the results of the first campaign [58], they constrain at 1σ the effective squared mass m 2 β as: with an approximately gaussian distribution around the best fit, currently in the physical region m 2 β > 0 [21] ( while it was negative in the first campaign [58]). The upper bound at 90% C.L. (∼ 1.6σ) corresponds to m 2 β < 0.6 eV 2 or m β < 0.8 eV, representing the first constraint on the effective β-decay mass in the sub-eV range [21]. Note that variants of the statistical analysis may lead to small differences in the second significant digit of m 2 β or m β [21], not considered herein. We implement the datum in Eq. (2) via a contribution ((x − 0.1)/0.3) 2 ∈ χ 2 , where x = (m β /eV) 2 .

B. Neutrinoless double beta decay and constraints on half-lives and m ββ
Neutrinoless double beta decay [9] can be considered as the process of creation of two matter particles (electrons) [59], occurring if neutrinos are of Majorana type [1]. Within the 3ν paradigm, the decay half-life T is given by where the index i labels the 0νββ nuclide, characterized by a phase space G i and a nuclear matrix element (NME) M i , while m ββ is the effective Majorana mass. The inverse half life S i = 1/T i represents, up to a constant factor, the observable decay rate or signal strength. Current experiments are consistent with null signal (S = 0), placing lower limits on T and upper limits on m ββ [19] via Eq. (3). In deriving separate and combined limits on m ββ , two issues arise: (1) experimental results are often given in terms of 90% C.L. bounds on T (say, T > T 90 ), with little or no information on the probability distribution of T (or of S); (2) theoretical uncertainties on the NME are rather large (and correlated among nuclides, e.g., via the axial coupling) [60]. See e.g. [61] for a recent discussion. We describe below our approach to these issues, in order to build first a probability distribution for half-lives T i and then, including NME's, for the effective mass m ββ .
We limit our analysis to experiments placing limits T 90 > 10 25 y, namely: GERDA [62] and MAJORANA [63] for 76 Ge; CUORE [64] for 130 Te; KamLAND-Zen 400 [65] and 800 (preliminary) [66] and EXO-200 [67] for 136 Xe. For each experiment we need not only a limit (T 90 ) but the probability distribution of T i or, equivalently, a function ∆χ 2 (S i ). Unfortunately, such detailed information is not provided by current experiments in an explicit or user-friendly way; see [61,68] for recent attempts to parametric reconstructions.
We adopt a ∆χ 2 (S i ) parametrization inspired by [68] and based on the following considerations. In 0νββ searches with zero background and nearly null results, the likelihood L of a signal S > 0 should be a poissonian with a scaling coefficient µ (L ∼ exp(−µS)), leading to a linear dependence on S (∆χ 2 ∼ ln L ∝ S) [61]. In 0νββ searches with nonnegligible background subtraction, the dependence is expected to be nearly gaussian (∆χ 2 ∝ (S − S 0 ) 2 ), where the best-fit signal S 0 may fall either in the physical region (S 0 ≥ 0) or in the unphysical one (S 0 < 0). All these limiting cases can be covered by a quadratic form as previously advocated in [68] on an empirical basis. Note that the offset c i is set by the condition that ∆χ 2 ≥ 0 in the physical region S i ≥ 0. In particular, for a i > 0, a ∆χ 2 minimum in the physical region implies b i ≤ 0 and c i = b 2 i /4a i , while in the unphysical one it implies b i > 0 and c i = 0. For a i = 0 one recovers the linear limit, that implies b i > 0 and c i = 0. The case a i < 0 is never realized.
In order to assess the coefficients (a i , b i , c i ), we have carefully sifted the information contained in the experimental publications [62][63][64][65][66][67] and in available PhD theses conducted within EXO-200 [69,70], KamLAND-Zen 400 [71] and KamLAND-Zen 800 (preliminary) [72]. We find that the linear approximation advocated in [61] for GERDA, can be roughly applied also to MAJORANA (up to subleading corrections at small S i , neglected herein). The other experimental bounds on T i can be reasonably approximated by parabolic ∆χ 2 curves, setting the various coefficients (a i , b i , c i ). We also require that our 90% C.L. limit ∆χ 2 (S 90 ) = 2.706 reproduces the T 90 limit reported by each experiment, up to their quoted significant digits. We have further checked that, by shifting the minimum of each ∆χ 2 function from S = S 0 to S = 0, the corresponding 90% C.L. limits are in reasonable agreement with the reported sensitivities for the null hypothesis. We are thus confident that current experimental results are fairly well represented by our ∆χ 2 's. Finally, we combine different experiments probing the same nuclide, by adding up their ∆χ 2 functions. II: Neutrinoless double beta decay: Details of the adopted parametrization ∆χ 2 (Si) = ai S 2 i + bi Si + ci for the signal strength Si = 1/Ti, expressed in units of 10 −26 y −1 . The first two columns report the nuclide and the name of the experiment(s). The next three columns report our evaluation of the coefficients (ai, bi, ci), for the various experiments, taken either separately (upper six rows) or in combination for the same nuclide (lower three rows). The sixth column reports our 90% C.L. (∆χ 2 = 2.706) half-life limits T90 in units of 10 26 y, to be compared with the experimentally quoted ones in the seventh column (in the same units). Pertinent references are listed in the last column.  [We formally keep up to four significant digits, to avoid accumulation of round-off errors in the analysis.] By combining GERDA and MAJORANA, we evaluate a 90% C.L. limit as high as T > 2.07 × 10 26 y for the 76 Ge half life, about 15% higher than from GERDA alone. A similar improvement is obtained for 136 Xe, reaching T > 1.27 × 10 26 y in the combination of KamLAND-Zen and EXO data. For 130 Te, CUORE alone sets the bound T > 0.22 × 10 26 y. Figure 9 shows our parametrized ∆χ 2 functions in terms of 1/T = S (bottom abscissa) and of T (top abscissa). The left and right panels refer to separate experiments and to combinations for the same nuclide, respectively. The dotted horizontal lines intersect all curves at the 90% C.L. limit T 90 . Note that the hierarchy of bounds at 90% C.L. is not necessarily preserved at different statistical levels, since some curves cross each other. This is another reason to suggest that the experimental collaborations explicitly provide their probability profiles for T , rather than focusing on a single C.L. limit (T 90 ), that provides a poor summary of the data and a limited comparison of different results.
In order to translate constraints from T i to m ββ , one needs to know the phase space factors G i and the nuclear matrix elements |M i | in Eq. (3). The G i can be accurately computed, see [73,74] for recent calculations. The |M i | embed complex nuclear physics, currently treated in a variety of approaches that, unfortunately, still carry significant uncertainties despite the theoretical progress, see e.g. the recent reviews in [60,[75][76][77][78][79][80][81]. It is common practice to select a set of published values {M i } in order to obtain a set of upper bounds {m ββ }, whose spread may be taken as indicative of the theoretical uncertainties; see [61] for a very recent application. However, this procedure overlooks significant correlations among the NME uncertainties of different nuclides [82][83][84][85] that, as shown below, are as important as the uncertainties themselves in obtaining conservative bounds. In a given nuclear model, estimating NME covariances requires massive numerical experiments to generate many M i variants. To our knowledge, this task has been performed only in [82] within the quasiparticle random phase approximation (QRPA), by varying the axial coupling g A , the short-range correlations, the model basis, and the renormalization procedure, while requiring consistency with available 2νββ data. Apart from the subsequent papers [83][84][85], and despite the potential relevance of NME covariance issues [86], we are aware of just one independent (but preliminary) correlation matrix estimate in a different model [87] and of a single statistical analysis [88] based on the correlations in [82]. In the absence of novel estimates of NME covariances, we shall use the only available results of [82] at face value. We have checked that, in any case, the NME uncertainties estimated therein are conservative enough to cover within ±2σ most (and within ±3σ all) of the M i values compiled in [60,[75][76][77][78][79][80][81] for the three nuclides.
We summarize and use the results in [82] as follows. In order to deal with large |M i | variations, possibly hitting the unphysical range |M i | < 0, the (adimensional) |M i | values are parametrized in terms of the logarithms η i , where the overlined symbols represent central values, |M i | = e η i . The index i = 1, 2, 3 runs over 76 Ge,130 Te, 136 Xe. The NME variations ∆η i are endowed with variances σ 2 i and a covariance matrix σ ij = ρ ij σ i σ j , whose inverse defines the weight matrix w ij = (σ ij ) −1 . For any choice of m ββ and of ∆η i , the signal strength in the i-th nuclide is given by where q i = G i |M i | 2 . The strengths S i carry two ∆χ 2 contributions: an experimental one coming from ∆χ 2 (S i ), and a theoretical one coming from NME covariances. These contributions are coupled by-and must be minimized over-the three variations {∆η i } where the first line holds in general, while the second line refers to the specific parametrization in Eq. (4). Minimization yields three coupled equations in the three ∆η i unknowns, to be solved numerically. Neglecting NME correlations (or errors) amounts to setting ρ ij = δ ij (or ∆η i = 0). Bounds for a single nuclide are obtained by selecting a specific index i. Subtraction of a χ 2 min offset (if any) yields the desired ∆χ 2 (m ββ ). Table III reports the values of σ ij and q i used herein. Figure 10 shows the resulting bounds on m ββ for the three nuclides taken separately (left panel) and in combination (right panel). The solid and dotted curves refer to our estimates with and and without NME uncertainties, respectively. Currently, the most constraining results are obtained by combining 76 Ge data, followed by weaker constraints from 136 Xe and 130 Te data. In the right panel, the case with uncorrelated NME uncertainties is also shown (dashed line). The effect of correlations is noticeable and leads to more conservative bounds; in fact, when the NME of different nuclides are positively correlated, they are more likely to become all smaller at the same time (with respect to the uncorrelated case), allowing larger values of m ββ . Including correlated errors, we obtain the overall bound m ββ < 0.11 eV at 2σ; the same bound was previously estimated (although in a less refined way) as m ββ < 0.14 eV in [11] and as m ββ < 0.18 eV in [10], reflecting the steady experimental progress in the last few years. III: Neutrinoless double beta decay: NME covariance matrix σ ij and auxiliary coefficients q i = G i |M i | 2 , as derived from the results in [82]; see the text for details. The q i are given in units of 10 −26 y −1 eV −2 , for m ββ expressed in eV. i Nuclide σij qi 1 76   √ ∆χ 2 . The left and right panels refer, respectively, to separate and combined bounds from the three nuclides, with (solid) or without (dotted) NME uncertainties. In the right panel, the case with uncorrelated uncertainties is also shown (dashed).
We emphasize that the above methodology can be applied to 0νββ data consistent with either a null signal or a positive detection. It may include generic likelihoods for the half-lives T i (hopefully provided by the experimental collaborations) and alternative evaluations of the NME covariances (possibly computed in different nuclear models).

C. Cosmology and constraints on Σ
In this section we discuss various choices for cosmological data combinations, enlarging the set of cases previously considered in [11], in order to deal with known tensions about lensing data. We then focus on two specific cases, dubbed as default and alternative options, leading to different implications for Σ and the sensitivity to mass ordering.
We remind that in [11] we analyzed the following data in various combinations: the complete Planck 2018 data (Planck) on the angular CMB temperature power spectrum (TT) plus polarization power spectra (TE, EE) [25,28] and lensing reconstruction power spectrum (lensing) [29]; a compilation of Baryon Acoustic Oscillation (BAO) measurements, given by data from the 6dFGS [89], SDSS MGS [90], and BOSS DR12 [91] surveys; and the Hubble constant datum [H 0 (R19)] from HST observations of Cepheids in the Large Magellanic Cloud measurements [92]. We assume the standard 6-parameter ΛCDM model augmented with nonzero neutrino masses (ΛCDM+Σ) and, in some cases, we add an extra empirical parameter A lens to marginalize over the excess amount of gravitational lensing emerging in the fits of the Planck data [25]; see also [93] for a similar approach.
Here we also consider an alternative way to deal with the Planck lensing problem, in the light of the wider debate about data tensions in the standard ΛCDM model [22,24]. In particular, rather than introducing an additional parameter to account for unknown systematics, one may consider restricting the analysis to alternative CMB datasets not affected by internal or mutual tensions. A possibility is offered by the combination of the recent CMB anisotropy data from ACTPol-DR4 [26] with WMAP9 data [27], together with a Planck-derived conservative gaussian prior on the optical depth to reionization τ [τ prior = 0.065 ± 0.015]. The ACT Collaboration explored this combination in [26], finding no evidence for a lensing anomaly, and obtaining a relaxed upper bound on the total neutrino mass, Σ < 1.2 eV at 2σ. Using the same data combination, we get an upper bound in excellent agreement, Σ < 1.21 eV; interestingly, we also find that the best fit is at Σ 0.70 eV, with no significant difference between NO and IO. If we replace the gaussian prior on τ with the actual Planck polarization data at large angular scale (dubbed Planck LowE) [28], that directly constrain the value of τ , we obtain a slightly stronger upper bound Σ < 1.12 eV (with a best fit at Σ 0.58 eV). Finally, in order to achieve a more complete combination of CMB data not in tension with each other, we further add to ACT and WMAP the independent Planck lensing reconstruction results [29], that do not suffer of the lensing anomaly. In this case the upper bound becomes Σ < 0.96 eV, with a best fit at Σ 0.54 eV, as obtained with CMB-only data. Of course, by including additional data, such as BAO measurements, the upper bound on Σ would be pushed back to ∼ 10 −1 eV (or even below as shown in [94]), at the price of a significant mutual tension among different datasets; these fits, that would be more comprehensive but also more discordant, are not further discussed herein. Results of the global 3ν analysis of cosmological data within the standard ΛCDM + Σ model (possibly augmented with the A lens parameter). The inputs numbered from 0 to 9 are the same as in [11], and refer to various combinations of the Planck 2018 angular CMB temperature power spectrum (TT) plus polarization power spectra (TE, EE), lensing potential power spectrum (lensing), Barion Acoustic Oscillations (BAO), and the Hubble constant from HST observations of Cepheids in the Large Magellanic Cloud, H0(R19). The inputs numbered from 10 to 12 are new and refer to ACTPol-DR4 and WMAP9 data, in combination with a prior on optical depth to reionization (τ prior ), Planck polarization data at large angular scale (lowE), and lensing data. For each case we report the 2σ upper bound on the sum of ν masses Σ (marginalized over NO and IO), together with the ∆χ 2 difference between IO and NO, using cosmology only. In the last two columns, we report the same information as in the previous two columns, but using cosmological data plus m β and m ββ constraints. The specific cases numbered 3 and 12 are dubbed as default and alternative, see the text for details.  Table IV reports, for convenience, the bounds on Σ for both the cosmological inputs considered in [11] (cases 0-9) and the new ones discussed above (cases [10][11][12]. These inputs can be roughly divided into two categories: a first group (cases 0-6) where one includes at face value Planck CMB results (plus other data) despite the Planck lensing anomaly, obtaining relatively strong upper bounds on Σ and a noticeable sensitivity to mass ordering; and a second group (cases 7-12), where one "solves" the lensing anomaly (by either adding an extra model parameter or by considering alternative CMB data), with weaker upper bounds on Σ and no significant sensitivity to NO vs IO. Hereafter we focus on just two representative cases for these two different categories, namely, case #3 (dubbed "default" as in [11]) and case #12 (dubbed "alternative"). Figure 11 shows the ∆χ 2 curves for the default and alternative cases. In the default case, cosmological data push Σ towards its lowest physical values in both IO and NO, and favors the latter at the level of ∼ 1.5σ. In the alternative case, there is a preference (< 2σ) for higher Σ values, with a best fit at Σ 0.54 eV and an upper limit Σ < 0.96 eV at 2σ, while there are only minor differences between NO and IO. Roughly speaking, the alternative case corresponds to a putative cosmological "signal" for neutrino masses, equivalent to Σ 0.54 ± 0.22 eV (with symmetrized 1σ errors).  Table IV as #3 (left, solid) and #12 (right, dashed), taken as representative of default and alternative options, respectively. The cases #10 and #11 (not shown) would be qualitatively similar to case #12.

Cosmological inputs for nonoscillation data analysis
We remark that, at fixed Σ, the small differences between NO and IO fits are due to: a slight sensitivity of cosmological data to the different ordering of ν masses at small Σ; the conversion from fit probability densities P (Σ) to χ 2 (Σ) functions, as the P normalization covers different physical Σ ranges in NO and IO; and, to a lesser extent (∆χ 2 < 0.1), to numerical fit inaccuracies for P → 0 (i.e., for high χ 2 ). See also the comments in Sec. II C of [10].
Summarizing, the two cases in Fig. 11 correspond to two qualitatively different outcomes, that might persist in future cosmological data analyses, as opposite examples of the tradeoff between completeness and consistency of inputs. On the one hand, by combining various datasets at face value (regardless of possible tensions), one may obtain strong upper limits, Σ < O(10 −1 ) eV, with some sensitivity to mass ordering (typically in favor of NO, where Σ attains its lowest possible values). On the other hand, by combining selected (and mutually consistent) data sets, one may end up with more relaxed upper bounds on Σ, possibly shifting the best fit towards Σ ∼ (few×)10 −1 eV, at the price of a reduced sensitivity to mass ordering. We think that, at present, both options deserve to be explored, especially because they imply rather different outcomes in combination with other (non)oscillation neutrino data.

D. Results on pairs of nonoscillation observables
The nonoscillation observables (m β , m ββ , Σ) are strongly and positively correlated, via their common dependence on the absolute neutrino mass scale. Contrary to the case of oscillation parameters, it is useful to show first the results on pairs of observables, and then the projections on single ones. With respect to our previous works [10,11], we shall use linear (rather than logarithmic) coordinates as, e.g., advocated for the pair (Σ, m ββ ) in [80,95]. Since cosmological data play a major role in constraining the nonoscillation parameter space, we shall discuss the two different options (default and conservative) defined in the previous Section. Figure 12 shows the correlation bands at 2σ for the pairs (Σ, m β ) and (Σ, m ββ ) in linear scales, including only the constraints from oscillation data, for NO and IO taken separately (i.e., without the offset ∆χ 2 IO−NO ). In the top panel, the bands have a tiny width, reflecting the small fractional errors on the oscillation parameters (δm 2 , ∆m 2 , θ 12 , θ 13 ) relevant for the pair (m β , Σ). In the bottom panel, the widening of the bands is almost entirely due to the unknown Majorana phases in m ββ . See also [11] and  Figure 13 shows how the data from β decay, 0νββ decay and cosmology further constrain the 2σ bands in Fig. 12, for the two cosmological input options considered. In the default case, cosmological bounds on Σ dominate-via correlations-the constraints on m β and m ββ , which are squeezed to the relatively small 2σ regions around the best fits, located close to the lowest possible values for Σ in both NO and IO. In the alternative case, there is an interplay between cosmological and 0νββ data: the first would prefer Σ 0.58 eV, implying relatively high values for the Majorana mass (m ββ > 0.06 eV, see Fig. 12); however, such values are disfavored by 0νββ data at > 1σ in Fig. 10. A best-fit compromise is reached for intermediate values, Σ ∼ 0.4 eV and m ββ 0.05 eV, surrounded by large 2σ allowed regions. In the right panel, note that both cosmology and 0νββ data constrain the correlations bands from above, leading to a joint 2σ bound Σ < 0.85 eV, stronger than the bound from cosmology only (Σ < 0.96 eV, see also Table IV). In all cases, current β-decay data play a minor role in the overall fit.
The implications of Fig. 13 can be summarized as follows. In the default case, it appears that the current KATRIN experiment (probing m β > 0.2 eV) is not expected to find any signal, while planned 0νββ experiments are expected to probe at least the region covered by both NO and IO (m ββ > 0.02 eV). The region covered only by NO (m ββ < 0.02 eV) is more difficult to probe, and becomes eventually prohibitive as m ββ vanishes, see e.g. [77,96,97]. In the alternative case, a much larger phase space is amenable to β decay and 0νββ decay searches. Cosmological searches may find a signal for Σ in a wide sub-eV range. Neutrinoless double beta decay data might find a signal for m ββ anywhere below the current bounds. The KATRIN experiment might find a signal in its sensitivity region (m β > 0.2) eV, or at least strengthen significantly the upper bounds on m β . E. Results on single nonoscillation observables Figures 14 and 15 show the projected bounds on single nonoscillation parameters for the default and alternative cases, respectively. In both cases we account for the NO-IO offset coming from the combination of all nonoscillation data (i.e., the rightmost numbers in rows #3 and #12 of Table IV (that were omitted in the previous Figs. 12 and 13). A vertical rise of N σ occurs when the lower physical limits are reached. These two figures quantify previous considerations about the default and alternative options: the first exemplifies the case of strong upper bounds on Σ from cosmology, accompanied by some sensitivity to the mass ordering, and by hard-to-probe phase spaces for m β and m ββ ; the second represents the case of weaker upper bounds (and a possible signal) for Σ, with scarce sensitivity to the mass ordering but more optimistic expectations for m β and m ββ signals. Together with Fig. 3, the above Figs. 14 and 15 provide a neat summary of what we (do not) know in the standard 3ν paradigm.  Table IV, as indicated by the horizontal lines (the thick one corresponding to the default case).

IV. SYNTHESIS
We conclude our work by merging the information coming from the analysis of oscillation and nonoscillation data, that have one important observable in common: the mass ordering. Figure 16 shows a histogram with separate and combined contributions to the ∆χ 2 IO−NO . The first bin adds up the contributions from oscillation data, starting from the negative one in the combination of LBL accelerator, solar and KamLAND data, that becomes positive by adding SBL reactor data, and further increases with atmospheric data. The second bin shows the range spanned by all the cases considered in Table IV, for the fit to cosmological data only. The third bin shows the slight change induced by adding current constraints on m β and m ββ , that provide an extra upward shift (see Tab. IV). The fourth bin adds up the contents of the first and third bins, providing an overall indication in favor of NO in the range ∼ 2.5-3.2σ. Although none of the single oscillation or nonoscillation data sets provides compelling evidence for normal ordering yet, their current combination clearly favors this option at a global ∼ 3σ level.
In conclusion, the main results of our global analysis can be essentially summarized in terms of bounds N σ = N σ (p), as shown in the following figures: Fig. 3 for the neutrino oscillation parameters p = δm 2 , |∆m 2 |, θ 12 , θ 23 , θ 13 , δ; Figs. 14 and 15 for the nonoscillation observables p = m β , m ββ , Σ, in two representative cases for cosmological inputs; and Fig. 16 for the discrete mass ordering parameter, p = sign(∆m 2 ) = NO(+)/IO(−). Finishing the fabric of the 3ν paradigm amounts to having convergent, narrow and linear N σ bounds, for one surviving mass ordering, in terms of any continuous 3ν oscillation parameter and nonoscillation observable p (with the possible exception of m ββ , if neutrinos have a Dirac nature). At present, this goal has been reached for p = δm 2 , |∆m 2 |, θ 12 , θ 13 and, to some extent, for p = θ 23 (up to an octant ambiguity). The current results for p = δ, m β , m ββ , Σ and NO/IO may instead be considered as initial, shaky steps of a long march towards the characterization of the neutrino-antineutrino differences and of the absolute neutrino mass spectrum. On the way, we shall learn a lot about neutrino properties in many different contexts, clarify the origin of old and new data tensions, and possibly find obstacles that, tearing away the fabric of the 3ν paradigm, may reveal hidden new physics.