Neutrino masses and mixings: Status of known and unknown $3\nu$ parameters

Within the standard 3nu mass-mixing framework, we present an up-to-date global analysis of neutrino oscillation data (as of January 2016), including the latest available results from experiments with atmospheric neutrinos (Super-Kamiokande and IceCube DeepCore), at accelerators (first T2K anti-nu and NOvA nu runs in both appearance and disappearance mode), and at short-baseline reactors (Daya Bay and RENO far/near spectral ratios), as well as a reanalysis of older KamLAND data in the light of the"bump"feature recently observed in reactor spectra. We discuss improved constraints on the five known oscillation parameters (delta m^2, |Delta m^2|, sin^2theta_12, sin^2theta_13, sin^2theta_23), and the status of the three remaining unknown parameters: the mass hierarchy, the theta_23 octant, and the possible CP-violating phase delta. With respect to previous global fits, we find that the reanalysis of KamLAND data induces a slight decrease of both delta m^2 and sin^2theta_12, while the latest accelerator and atmospheric data induce a slight increase of |Delta m^2|. Concerning the unknown parameters, we confirm the previous intriguing preference for negative values of sin(delta) [with best-fit values around sin(delta) ~ -0.9], but we find no statistically significant indication about the theta_23 octant or the mass hierarchy (normal or inverted). Assuming an alternative (so-called LEM) analysis of NOvA data, some delta ranges can be excluded at>3 sigma, and the normal mass hierarchy appears to be slightly favored at 90% C.L. We also describe in detail the covariances of selected pairs of oscillation parameters. Finally, we briefly discuss the implications of the above results on the three non-oscillation observables sensitive to the (unknown) absolute nu mass scale: the sum of nu masses, the effective nu_e mass, and the effective Majorana mass.


Introduction
Yellow, blue and dark blue: this is the simple color palette used for painting and penning each of the two-sided Nobel Diplomas awarded to Takaaki Kajita [1] and Arthur B. McDonald [2]. On the left side, one can gaze at an artist's view-sketched with a few broad strokes-of the neutrino transformative trip from the bright yellow Sun, through the Earth's blue darkness, into a blue pool of water [3]. On the right side, one can read the-beautifully and precisely penned-Nobel laureate names and prize motivations, in ink colors that continuously change from deep blue to blue with yellow shades [4]. In a sense, the two sides of the Diplomas evoke the interplay between a broad-brush picture of ν masses and mixings (the pioneering era) and carefully designed measurements and theoretical descriptions (the precision era), in a continuous feedback between breakthrough and control, that may open the field to further fundamental discoveries [5].
In this paper, we aim at presenting both the broad-brush features and the fine structure of the current picture of neutrino oscillation phenomena, involving the mixing of the three neutrino states having definite flavor ν e,µ,τ with three states ν 1,2,3 having definite masses m i [6]. Information on known and unknown neutrino mass-mixing parameters is derived by a global analysis of neutrino oscillation data, that extends and updates our previous work [7] with recent experimental inputs, as discussed in Sec. 2 (see also [8,9] for previous global analyses by other groups). In Sec. 3, precise constraints (at few % level) are obtained on four well-known oscillation parameters, namely, the squared-mass differences δm 2 = m 2 2 − m 2 1 and ∆m 2 = m 2 3 − (m 2 1 + m 2 2 )/2, and the mixing angles θ 12 and θ 13 . Less precise constraints, including an octant ambiguity, are reported for the angle θ 23 . In this picture, we also discuss the current unknowns related to the neutrino mass hierarchy [sign(∆m 2 )] and to the possible leptonic CP-violating phase δ. The trend favoring negative values of sin δ appears to be confirmed, with best-fit values around δ 1.3-1.4 π (i.e., sin δ −0.9). More fragile indications, which depend on alternative analyses of specific data sets, concern the exclusion of some δ ranges at > 3σ, and a slight preference for normal hierarchy at 90% C.L. The covariances of selected parameter pairs, and the implications for non-oscillation searches, are presented in Sec. 4 and 5, respectively. Our conclusions are summarized in Sec. 6.

Global analysis: methodology and updates
In this section we discuss methodological issues and input updates for the global analysis. Readers interested only in the fit results may jump to Sec. 3.
In general, no single oscillation experiment can currently probe, with high sensitivity, the full parameter space spanned by the mass-mixing variables (δm 2 , ±∆m 2 , θ 12 , θ 13 , θ 23 , δ). One can then group different data sets, according to their specific sensitivities or complementarities with respect to some oscillation parameters. We follow the methodology of Refs. [7,10] as summarized below . We first combine the data coming from solar and KamLAND reactor experiments ("Solar+KL") with those coming from long-baseline accelerator searches in both appearance and disappearance modes ("LBL Acc"). The former data set constrains the (δm 2 , θ 12 ) parameters (and, to some extent, also θ 13 [7,10]), which are a crucial input for the 3ν probabilities relevant to the latter data set. The combination "LBL Acc+Solar+KL data" provides both upper and lower bounds on the (δm 2 , |∆m 2 |, θ 12 , θ 13 , θ 23 ) parameters but, by itself, is not particularly sensitive to δ or to sign(±∆m 2 ) (+ for normal hierarchy, NH, and − for inverted hierarchy, IH).
The LBL Acc+Solar+KL data are then combined with short-baseline reactor data ("SBL Reactors"), that provide strong constraints on the θ 13 mixing angle via disappearance event rates, as well as on useful bounds on ∆m 2 via spectral data (when available). The synergy between LBL Acc+Solar+KL data and SBL Reactor data increases significantly the sensitivity to δ [7].
Finally, we add atmospheric neutrino data ("Atmos"), which probe both flavor appearance and disappearance channels for ν and ν, both in vacuum and in matter, with a very rich phenomenology spanning several decades in energy and path lengths. This data set is dominantly sensitive to the mass-mixing pair (∆m 2 , θ 23 ) and, subdominantly, to all the other oscillation parameters. Despite their complexity, atmospheric data may thus add useful pieces of information on subleading effects (and especially on the three unknown parameters), which may either support or dilute the indications coming from the previous data sets.
In all cases, the fit results are obtained by minimizing a χ 2 function, that depends on the arguments (δm 2 , ±∆m 2 , θ 12 , θ 13 , θ 23 , δ) and on a number of systematic nuisance parameters via the pull method [7,11]. Allowed parameter ranges at N σ standard deviations are defined via N 2 σ = χ 2 − χ 2 min [6]. The same definition is maintained in covariance plots involving parameter pairs, so that the previous N σ ranges are recovered by projecting the allowed regions onto each axis. Undisplayed parameters are marginalized away.
A final remark is in order. The definition N 2 σ = χ 2 −χ 2 min is based on Wilks' theorem [6], that is not strictly applicable to discrete choices (such as NH vs IH, see [12][13][14] and references therein) or to cyclic variables (such as δ, see [15,16]). Concerning hierarchy tests, it has been argued that the above N σ prescription can still be used to assess the statistical difference between NH and IH with good approximation [17]. Concerning CP violation tests, the prescription appears to lead (in general) to more conservative bounds on δ, as compared with the results obtained from numerical experiments [8,18,19]. In principle, one can construct the correct χ 2 distribution by generating extensive replicas of all the relevant data sets via Monte Carlo simulations, randomly spanning the space of the neutrino oscillation and systematic nuisance parameters. However, such a construction would be extremely time-consuming and is beyond the scope of this paper. For the sake of simplicity, we shall adopt the conventional N σ definition, supplemented by cautionary comments when needed.

Solar + KL data analysis and the reactor spectrum bump
With respect to [7], the solar neutrino analysis is unchanged. Concerning KamLAND (KL) reactor neutrinos, we continue to use the 2011 data release [20] as in [7]. We remark that the latest published KL data [21] are divided into three subsets, with correlated systematics that are difficult to implement outside the Collaboration. 1 In this work, we reanalyze the 2011 KL data for the following reason.
The KL analysis requires the (unoscillated) absolute reactor ν e spectra as input. In this context, a new twist has been recently provided by the observation of a ∼ 10% event excess in the range E ν ∼ 5-7 MeV (the so-called "bump" or "shoulder") [22], with respect to the expectations from reference Huber-Müller (HM) spectra [23,24], in each of the current high-statistics SBL reactor experiments RENO [25], Double Chooz [26] and Daya Bay [27].
This new spectral feature is presumably due to nuclear physics effects (see the recent review in [28]), whose origin is still subject to investigations and debate [29][30][31][32][33]. In principle, one would like to know in detail the separate spectral modifications for each reactor fuel component [34,35]. However, the only information available at present is the overall energy-dependent ratio f (E) between data and HM predictions, which we extract (and smooth out) from the latest Daya Bay results (see Fig. 3 in [27]) . We use the f (E) ratio as an effective fudge factor multiplying the unoscillated HM spectra for Kam-LAND, which are thus anchored to the absolute Daya Bay spectrum [27]. In our opinion, this overall correction can capture the main bump effects in the KL spectral analysis. More refined KL data fits will be possible when the bump feature(s) will be better understood and broken down into separate spectral components. Concerning the KL dominant oscillation parameters (δm 2 , θ 12 ), we find that the inclusion of the bump fudge factor induces a slight negative shift of their best-fit values, which persists in combination with solar data (see Sec. 3).
Finally, we recall that the 3ν analysis of solar+KL data is performed in terms of three free parameters (δm 2 , θ 12 , θ 13 ), providing a weak but interesting indication for nonzero θ 13 [7,10,36]. Tiny differences between transition probabilities in NH and IH [37,38] are negligible within the present accuracy. The hierarchy-independent function χ 2 (δm 2 , θ 12 , θ 13 ), derived from the solar+KL data fit, is then used in combination with the following LBL accelerator data.

LBL Accelerator data analysis
With respect to [7], we include the most recent results from the Tokai-to-Kamioka (T2K) experiment in Japan and from the NOνA experiment at Fermilab, in both appearance and disappearance modes. In particular, we include the latest T2K neutrino data [39] and the first T2K antineutrino data [40], as well as the first NOνA neutrino data as of January 2016 [41,42]. The statistical analysis of LBL experiments has been performed using a modified version of the software GLoBES [43,44] for the calculation of the expected number of events. 2 For each LBL data set, the χ 2 function takes into account Poisson statistics [6] and the main systematic error sources, typically related to energy-scale errors and to normalization uncertainties of signals and backgrounds, as taken from [39][40][41][42]. Concerning NOνA ν e appearance data, the collaboration used two different event selection methods for increasing the purity of the event sample: A primary method based on a likelihood identification (LID) selector, and a secondary one based on a library event matching (LEM) selector, leading to somewhat different results for the ν e signal and background [41]. We shall consider the LID data as a default choice for NOνA, but we shall also comments on the impact of the alternative LEM data.
We have reproduced with good approximation the allowed parameter regions shown by T2K [39,40] and by NOνA [42] (in both LID and LEM cases [41]), under the same hypotheses or restrictions adopted therein for the undisplayed parameters. We remark that, in our global analysis (see Sec. 4), all the oscillation parameters are left unconstrained. Note that we define the parameter ∆m 2 , driving the dominant LBL oscillations, as in both NH (∆m 2 > 0) and IH (∆m 2 < 0) [37,38]. A comparative discussion of this and alternative conventions in terms of ∆m 2 13 , ∆m 2 23 , ∆m 2 A , ∆m 2 µµ and ∆m 2 ee is reported in [45,46] and references therein. Although any such convention is immaterial (as far as the full 3ν oscillation probabilities are used), the adopted one must be explicitly declared, since the various definitions differ by terms of O(δm 2 ), comparable to the current ±1σ uncertainty of ∆m 2 .

SBL Reactor data analysis
With respect to [7], we include herein the spectral data on the far-to-near detector ratio as a function of energy, as recently reported by the experiments Daya Bay (Fig. 3 of [47]) and RENO (Fig. 3 of [25]). Besides the statistical errors, we include a simplified set of pulls for energy-scale and flux-shape systematics, since the bin-to-bin correlations are not publicly reported in [25,47]. We neglect systematics related to the spectral bump feature, which affect absolute spectra (see Sec. 2.1) but largely cancel in the analysis of far/near ratios (see [25,47]) . We reproduce with good accuracy the joint allowed regions reported in [47] and [25] for the mixing amplitude sin 2 2θ 13 and their effective squared mass parameters ∆m 2 ee , for both NH and IH. 3 We then combine the Daya Bay and RENO analyses, in terms of our default parameters sin 2 θ 13 and ∆m 2 . The combined fit results are dominated, for both mass-mixing parameters, by the high-statistics Daya Bay data. While the reactor bounds on θ 13 are extremely strong, the current bounds on ∆m 2 are not yet competitive with those coming from LBL accelerator data in disappearance mode, although they help in reducing slightly its uncertainty (see Sec. 4).

Atmospheric ν data analysis
With respect to [7], we update our analysis of Super-Kamiokande (SK) atmospheric neutrino data by including the latest (phase I-IV) data as taken from [49,50]. We also include for the first time the recent atmospheric data released by the IceCube DeepCore (DC) collaboration [51][52][53]. We reproduce with good accuracy the joint bounds on the sin 2 θ 23 and ∆m 2 32 parameters shown by DC in [51], under the same assumptions used therein. In this work, the χ 2 functions for SK and DC have been simply added. In the future, it may be useful to isolate and properly combine possible systematics which may be common to SK and DC (related, e.g., to flux and cross section normalizations).

Global 3ν analysis: constraints on single oscillation parameters
In this Section we discuss the constraints on known and unknown oscillation parameters, coming from the global 3ν analysis of all the data discussed above. The impact of different data sets will be discussed in the next Section.  Figure 1 shows the bounds on single oscillation parameters, in terms of standard deviations N σ from the best fit. Linear and symmetric curves would correspond to gaussian uncertainties-a situation realized with excellent approximation for the (∆m 2 , sin 2 θ 13 ) mass-mixing pair and, to a lesser extent, for the (δm 2 , sin 2 θ 12 ) pair. The best fit of the sin 2 θ 23 parameter flips from the first to the second octant by changing the hierarchy from normal to inverted, but this indication is not statistically significant, since maximal mixing (sin 2 θ 23 = 1/2) is allowed at ∼ 1.6σ (∼ 90% C.L.) for NH and at ∼ 1σ for IH. In any case, all these parameters have both upper and lower bounds well above the 3σ level. If we define the average 1σ error as 1/6 of the ±3σ range, our global fit implies the following fractional uncertainties: δm 2 (2.4%), sin 2 θ 12 (5.8%), ∆m 2 (1.8%), sin 2 θ 13 (4.7%), and sin 2 θ 23 (9%).
The parameter δ is associated to a Dirac phase in the neutrino mixing matrix, which might induce leptonic CP violation effects for sin δ 0 [6]. Recent fits to global ν data [7][8][9] and partial (LBL accelerator) data [19,54] have consistently shown a preference for negative values of sin δ, as a result of the combination of LBL accelerator ν and ν data and of SBL reactor data. The reason is that the LBL appearance probability contains a CP-violating part proportional to − sin δ (+ sin δ) for neutrinos (antineutrinos) [6]. With respect to the CP-conserving case sin δ = 0, values of sin δ < 0 are then expected to produce a slight increase (decrease) of events in ν µ → ν e (ν µ → ν e ) oscillations for θ 13 fixed (by reactors), consistently with the appearance results of T2K (using both ν [39] and ν [39]) and in NOνA (using ν [41]), although within large statistical uncertainties. 4 This trend for δ is clearly confirmed by the results in Fig. 1, which show a best fit for sin δ −0.9 (δ 1.3-1.4π) in both NH and IH, while opposite values around sin δ +0.9 are disfavored at almost 3σ level. Although all values of δ are still allowed at 3σ, the emerging indications in favor of sin δ < 0 are intriguing and deserve further studies in T2K and NOνA, as well as in future LBL accelerator facilities. We remark that our bounds on δ are conservative, and that dedicated constructions of the χ 2 distributions via extensive numerical simulations might lead to stronger indications on δ, as discussed in Sec. 2. Table 1: Results of the global 3ν oscillation analysis, in terms of best-fit values and allowed 1, 2 and 3σ ranges for the 3ν mass-mixing parameters. See also Fig. 1 for a graphical representation of the results. We recall that ∆m 2 is defined as m 2 3 − (m 2 1 + m 2 2 )/2, with +∆m 2 for NH and −∆m 2 for IH. The CP violating phase is taken in the (cyclic) interval δ/π ∈ [0, 2]. The last row reports the (statistically insignificant) overall χ 2 difference between IH and NH.

Parameter
Hierarchy  Table 1 shows the same results of Fig. 1 in numerical form, with three significant digits for each parameter. In the last row of the table we add a piece of information not contained in Fig. 1, namely, the ∆χ 2 difference between normal and inverted hierarchy. The NH is slightly favored over the IH at the (statistically insignificant) level of ∼ 1σ in the global fit. We remark that both Fig. 1 and Table 1 use the NOνA LID data set in appearance mode (see Sec. 2.2).
By adopting the alternative NOνA LEM data set, we find no variation for the δm 2 and sin 2 θ 12 parameters (dominated by Solar+KL data) and for the ∆m 2 parameter (dominated by LBL data in disappearance mode, in combination with atmospheric and reactor spectral data). We find slight variations for sin 2 θ 13 and sin 2 θ 23 , and a small but interesting increase of the bounds on δ above the 3σ level. Figure 2 shows the corresponding results for the sin 2 θ 23 and δ parameters, to be compared with the rightmost panels of Fig. 1.  In Table 2 we report the results of the global fit using NOνA LEM data, but only for those parameters bounds which differ from Table 1. Some intervals surrounding δ π/2 can be excluded at > 3σ. We also find an increased sensitivity to the hierarchy, with the NH slightly favored (at 90% C.L.) over the IH. These indications, although still statistically limited, deserve some attention, for reasons that will be discussed in more detail at the end of the next Section.

Global 3ν analysis: Parameter covariances and impact of different data sets
In this section we show and interpret the joint N σ contours (covariances) for selected pairs of oscillation parameters. We also discuss the impact of different data sets on such bounds.
We start with the analysis of the (δm 2 , sin 2 θ 12 , sin 2 θ 13 ) parameters, which govern the oscillations phenomenology of solar and KamLAND neutrinos. Figure 3 shows the corresponding bounds derived by a fit to Solar+KL data only (solid lines). By themselves, these data provide a ∼ 1σ hint of θ 13 > 0 [36], with a best fit (sin 2 θ 13 0.16) close to current SBL reactor values.  The (δm 2 , sin 2 θ 12 ) parameters in Fig. 3 appear to be slightly anticorrelated, with a best-fit point at (7.37 × 10 −5 eV 2 , 0.297). These values are slightly lower than those reported in our previous work (7.54 × 10 −5 eV 2 , 0.308) [7], as a result of altering the absolute KL spectra to account for the bump feature (see Sec. 2.1). Statistically, these deviations amount to about −1σ for δm 2 and −0.6σ for sin 2 θ 12 , and thus are not entirely negligible. A better understanding of the absolute reactor spectra (in both normalization and shape) is thus instrumental to analyze the KamLAND data with adequate precision.
Finally, Fig. 3 shows the joint bounds on the (δm 2 , sin 2 θ 12 , sin 2 θ 13 ) parameters from the global fit including all data (dashed lines). The bounds on the pair (δm 2 , sin 2 θ 12 ) are basically unaltered, while those on sin 2 θ 13 are shrunk by more than an order of magnitude, mainly as a result of SBL reactor data.
Let us consider now the interplay between sin 2 θ 13 and the mass-mixing parameters ∆m 2 and sin 2 θ 23 , which dominate the oscillations of LBL accelerator neutrinos. Figure 4 shows the covariance plot for the (sin 2 θ 13 , ∆m 2 ) parameters. Starting from the leftmost panels, one can see that the LBL Acc.+Solar+KL data, by themselves, provide both upper and lower bounds on sin 2 θ 13 at 3σ level. The best-fit values of sin 2 θ 13 lie around ∼ 0.02 in either NH or IH, independently of SBL reactor data. The best-fit values of ∆m 2 are slightly higher than in our previous work [7], mainly as a result of the recent NOνA data. The joint (sin 2 θ 13 , ∆m 2 ) contours appear to be somewhat bumpy, as a result of the octant ambiguity discussed below. 5 In the middle panels, the inclusion of SBL reactor data improves dramatically the bounds on sin 2 θ 13 and, to a small but nonnegligible extent, also those on ∆m 2 . Finally, in the rightmost panels, atmospheric data induce a small increase of its central value (mainly as a result of DeepCore data), and a further reduction of the ∆m 2 uncertainty. In comparison with [7], the ∆m 2 value is shifted by ∼ 1σ upwards in the global fit. Figure 5 shows the covariance plot for the (sin 2 θ 23 , sin 2 θ 13 ) parameters. The leftmost panels show a slight negative correlation and degeneracy between these two variables, which is induced by the dominant dependence of the LBL appearance channel on the product sin 2 θ 13 sin 2 θ 23 , as also discussed in [7,10]. The overall LBL acc+Solar+KL preference for relatively low values of sin 2 θ 13 (∼ 0.02) breaks such a degeneracy and leads to a weak preference for the second octant.  refer to the analysis of LBL Acc+Solar+KL data (left panels), plus SBL reactor data (middle panels), plus Atmospheric data (right panels), with best fits marked by dots. The three upper (lower) panels refer to NH (IH). SBL reactor data (middle panels of Fig. 5) shrink the sin 2 θ 13 range for both NH and IH. For IH, however, they do not significantly change the central value of sin 2 θ 13 , nor the correlated best-fit value of sin 2 θ 23 , which stays in the second octant. Conversely, for NH, the SBL reactor data do shift the central value of sin 2 θ 13 upwards (with respect to the left panel), and best-fit value of sin 2 θ 23 is correspondingly shifted into first octant. Finally, the inclusion of atmospheric data (rightmost panels) alters the N σ contours, but does not change the qualitative preference for the first (second) octant of θ 23 in NH (IH). Figure 6 shows the octant ambiguity in terms of bounds on the mass-mixing parameters (∆m 2 , sin 2 θ 23 ). The fragility of current octant indications stems from the data themselves rather than on analysis details: nearly maximal mixing is preferred by T2K (accelerator) and DeepCore (atmospheric) data, while nonmaximal mixing is preferred by MINOS and NOνA (accelerator) and by SK (atmospheric) data. The combined results on θ 23 appear thus still fragile, as far as the long-standing octant degeneracy [57] is concerned.  Fig. 4, but for the (sin 2 θ 13 , δ) parameters.
A very recent example of the (non)maximal θ 23 issue is provided by the NOνA data in disappearance mode, which entailed a preference for maximal mixing with preliminary data [56] and for nonmaximal mixing with definitive data [42]. We trace this change to the migration of a few events among reconstructed energy bins in the final NOνA data (not shown).
Let us complete the covariance analysis by discussing the interplay of the CP-violating phase δ with the mixing parameters sin 2 θ 13 and sin 2 θ 23 . Figure 7 shows the N σ bounds in the (sin 2 θ 13 , δ) plane, which is at the focus of LBL accelerator searches in appearance mode [40,41,55]. The leftmost panels show the wavy bands allowed by LBL Acc.+Solar+KL data, with a bumpy structure due to the octant ambiguity (which was even more evident in older data fits [10]). In the middle panels, SBL Reactor data select a narrow vertical strip, which does not alter significantly the preference for δ ∼ 3π/2 stemming from LBL Acc.+Solar+KL data alone. In this context, it is sometimes asserted that the current preference for δ ∼ 3π/2 emerges from a "tension" between LBL accelerator and SBL reactor data on θ 13 ; however, Fig. 7 clearly show that these data are currently highly consistent with each other about θ 13 , and that their interplay should be described in terms of synergy rather than tension. Finally, the inclusion of atmospheric data (rightmost panels) corroborates the previous indications for δ, with a global best fit around 1.3-1.4π and a slight reduction of the allowed ranges at 1 and 2σ (at least for NH).
Note that, with NOνA LEM data, the wavy bands in the leftmost panels of Fig. 7 would be slightly shifted to the right (not shown), leading to slightly stronger bounds on δ in combination with SBL reactor and atmospheric data, as reported in Fig. 2 of the previous Section; see also the official NOνA LID and LEM results in [41]. In this case, one might invoke a slight "tension" between LBL accelerator and and SBL reactor data, but only at the level of ∼ 1σ differences on θ 13 in the worst case (IH). Figure 8 shows the N σ bounds in the (sin 2 θ 23 , δ) plane, which is gaining increasing attention from several viewpoints, including studies of degeneracies among these parameters and θ 13 [58], of the interplay between LBL appearance and disappearance channels [59], and of statistical issues in the interpretation of N σ bounds [19]. The bounds in Fig. 8 appear to be rather asymmetric in the two half-ranges of both θ 23 and δ, and also quite different in NH and IH. This is not entirely surprising, since this is the only covariance plot (among Figs. 3-8) between two unknowns: the θ 23 octant (in abscissa) and the CP-violating phase δ (in ordinate). Therefore, the contours of Fig. 8 may evolve significantly as more data are accumulated, especially by oscillation searches with atmospheric and LBL accelerator experiments.
We conclude this Section by commenting the ∆χ 2 I−N values reported in Sec. 3, which differ only by the inclusion of NOνA LID data (Table 1) vs LEM data ( Table 2). In the first case, the ∆χ 2 I−N takes the value −1.2 for a fit to LBL Acc.+Solar+KL data, becomes −0.88 by including SBL reactor data, and changes (also in sign) to +0.98 by including atmospheric data. 6 Since these ∆χ 2 I−N values are both small (at the level of ±1σ) and with unstable sign, we conclude that there is no significant indication about the mass hierarchy, at least within the global fit including default (LID) NOνA data. By replacing NOνA LID with LEM data (Table 2), the same exercise leads to the following progression of ∆χ 2 I−N values: +0.61 (LBL Acc.+Solar+KL), +2.2 (plus SBL Reactor), +2.8 (plus Atmospheric). In this case, a weak hint for NH (at ∼ 1.6σ, i.e., ∼ 90% C.L.) seems to emerge from consistent (same-sign) indications coming from different data sets, which is the kind of "coherent" signals that one would hope to observe, at least in principle. Time will tell if these fragile indications about the hierarchy will be strengthened or weakened by future data with higher statistics.

Implications for absolute neutrino masses
Let us discuss the implications of the previous oscillation results on the three observables sensitive to the (unknown) absolute ν mass scale: the sum of ν masses Σ (probed by precision cosmology), the effective ν e mass m β (probed by β decay), and the effective Majorana mass m ββ (probed by 0νββ decay if neutrinos are Majorana fermions). Definitions and previous constraints for these observables can be found in [10,38,60]; here we just remark that the following discussion is not affected by the current uncertainties on θ 23 or δ. Figure 9 shows the constraints induced by our global 3ν analysis at 2σ level, for either NH (blue curves) or IH (red curves), in the planes charted by any two among the parameters m β , m ββ and Σ. The allowed bands for NH and IH, which practically coincide in the (so-called degenerate) mass region well above O(10 −1 ) eV, start to differ significantly at relatively low mass scales of O( √ ∆m 2 ) and below. At present, βand 0νββ-decay data probe only the degenerate region of m β and m ββ , respectively [6], while cosmological data are deeply probing the sub-eV scale, with upper bounds on Σ as low as ∼ 0.1-0.2 eV, see e.g. [61][62][63][64][65][66][67] and references therein. Taken at face value, the cosmological bounds would somewhat disfavor the IH case, which entails Σ values necessarily larger than ∼ 2 √ ∆m 2 ∼ 0.1 eV (see Fig. 9). Interestingly, these indications are consistent with a possible slight preference for NH from the global 3ν analysis (with NOνA LEM data), as discussed at the end of the previous Section. We do not attempt to combine cosmological and oscillation data, but we remark that the evolution of such hints will be a major issue in neutrino physics for a long time, with challenging implications for β-decay and 0νββ-decay searches [68].

Conclusions
We have presented the results of a state-of-the-art global analysis of neutrino oscillation data, performed within the standard 3ν framework. Relevant new inputs (as of January 2016) include the latest data from the Super-Kamiokande and IceCube DeepCore atmospheric experiments, the long-baseline accelerator data from T2K (antineutrino run) and NOνA (neutrino run) in both appearance and disappearance mode, the far/near spectral ratios from the Daya Bay and RENO short-baseline reactor experiments, and a reanalysis of KamLAND data in the light of the "bump" feature recently observed in reactor antineutrino spectra.
The status of the three unknown oscillation parameters is as follows. The θ 23 octant ambiguity remains essentially unresolved: The central value of sin 2 θ 23 is somewhat fragile, and it can flip from the first to the second octant by changing the data set or the hierarchy. Concerning the CP-violating phase δ, we confirm the previous trend favoring sin δ < 0 (with a best fit at sin δ −0.9), although all δ values are allowed at 3σ. Finally, we find no statistically significant indication in favor of one mass hierarchy (either NH or IH). Some differences arise by changing the NOνA appearance data set, from the default (LID) sample to the alternative (LEM) sample. A few known parameters are slightly altered, as described in Fig. 2 and Table 2. There is no significant improvement on the octant ambiguity, while the indications on δ are strengthened, and some ranges with sin δ > 0 can be excluded at 3σ level. Concerning the mass hierarchy, the NH case appears to be slightly favored (at ∼ 90% C.L.).
We have discussed in detail the parameter covariances and the impact of different data sets through Figs. 3-8, that allow to appreciate the interplay among the various (known and unknown) parameters, as well as the synergy between oscillation searches in different kinds of experiments. Finally, we have analyzed the implications of the previous results on the non-oscillation observables (m β , m ββ , Σ) that can probe absolute neutrino masses (Fig. 9). In this context, tight upper bounds on Σ from precision cosmology appear to favor the NH case. Further and more accurate data are needed to probe the hierarchy and absolute mass scale of neutrinos, their Dirac or Majorana nature and CP-violating properties, and the θ 23 octant ambiguity, which remain as missing pieces of the 3ν puzzle.