Global analysis of three-flavour neutrino oscillations: synergies and tensions in the determination of theta_23, delta_CP, and the mass ordering

We present the results of a global analysis of the neutrino oscillation data available as of fall 2018 in the framework of three massive mixed neutrinos with the goal at determining the ranges of allowed values for the six relevant parameters. We describe the complementarity and quantify the tensions among the results of the different data samples contributing to the determination of each parameter. We also show how those vary when combining our global likelihood with the chi^2 map provided by Super-Kamiokande for their atmospheric neutrino data analysis in the same framework. The best fit of the analysis is for the normal mass ordering with inverted ordering being disfavoured with a Delta_chi^2 = 4.7 (9.3) without (with) SK-atm. We find a preference for the second octant of theta_23, disfavouring the first octant with Delta_chi^2 = 4.4 (6.0) without (with) SK-atm. The best fit for the complex phase is Delta_CP = 215_deg with CP conservation being allowed at Delta_chi^2 = 1.5 (1.8). As a byproduct we quantify the correlated ranges for the laboratory observables sensitive to the absolute neutrino mass scale in beta decay, m_nu_e, and neutrino-less double beta decay, m_ee, and the total mass of the neutrinos, Sigma, which is most relevant in Cosmology.


Introduction
Flavour transitions of neutrinos via the energy and distance dependent neutrino oscillation mechanism [1,2] is a well established phenomenon, which proves that at least two out of the three neutrinos in the Standard Model must have tiny but non-zero masses. In this work we revisit the status of three-flavour neutrino oscillations in view of latest global data.
To fix the convention, the three flavour neutrinos, ν e , ν µ , ν τ , are defined via the weak charged current. They are expressed as superposition of the three neutrino mass eigen-fields ν i (i = 1, 2, 3) with masses m i via a unitary leptonic mixing matrix [3,4] by U αi ν i (α = e, µ, τ ) . (1.1) The mixing matrix we parametrize as: where c ij ≡ cos θ ij and s ij ≡ sin θ ij . The angles θ ij can be taken without loss of generality to lie in the first quadrant, θ ij ∈ [0, π/2], and the phase δ CP ∈ [0, 2π]. Values of δ CP different from 0 and π imply CP violation in neutrino oscillations in vacuum [5][6][7]. P is a diagonal matrix which is the identity if neutrinos are Dirac fermions and it contains two additional phases, P = diag(e iα 1 , e iα 2 , 1), if they are Majorana fermions. The Majorana phases α 1 and α 2 play no role in neutrino oscillations [6,8].
In this convention there are two non-equivalent orderings for the neutrino masses, namely normal ordering (NO) with m 1 < m 2 < m 3 , and inverted ordering (IO) with m 3 < m 1 < m 2 . Furthermore the data show a hierarchy between the mass splittings, ∆m 2 21 |∆m 2 31 | |∆m 2 32 | with ∆m 2 ij ≡ m 2 i − m 2 j . In this work we follow the convention from Ref. [9] and present our results for both, NO and IO, using the smallest and largest mass splittings. The smallest one is always ∆m 2 21 , while the largest one we denote by ∆m 2 3 , with = 1 for NO and = 2 for IO. Hence, ∆m 2 3 = ∆m 2 31 > 0 for NO , ∆m 2 32 < 0 for IO .

(1.3)
Due to the wealth of experiments exploring neutrino oscillations, we are in the situation that a given parameter is determined by several measurements. Therefore, combined analyses such as the one presented below are an important tool to extract the full information on neutrino oscillation parameters. This is especially true for open questions, such as the octant of θ 23 , the type of the neutrino mass ordering, and the status of the complex phase δ CP , where some hints are emerging due to significant synergies between different experiments. However, also for parameters describing dominant oscillations, a significantly more accurate determination emerges by the combination of complementary data sets, such as for example for |∆m 2 3 |. We present below the global fit NuFIT-4.0, updating our previous analyses [9][10][11]. ∆χ 2 maps and future updates of this analysis will be made available at the NuFIT website [12]. For other recent global fits see [13,14].

Data samples analyzed
The analysis presented below uses data available up to fall 2018. A complete list of the used data including references can be found in appendix A. Here we give a brief overview of recent updates and mention changes with respect to our previous analysis [11].
We include latest data from the MINOS [15,16], T2K [17,18], and NOvA [19,20] long-baseline accelerator experiments from ν µ disappearance and ν µ → ν e appearance channels, both for neutrino and anti-neutrino beam modes. In particular, T2K and NOvA have presented updated results at the Neutrino18 conference, including also first data on anti-neutrinos from NOvA, whose impact will be discussed below.
Concerning reactor neutrino experiments, the fit of data with baselines in the km range (medium baseline, MBL) is completely dominated by modern experiments, most importantly by Daya Bay [21], with subleading contributions from RENO [22] and Double Chooz [23]. Moreover, those experiments are entirely based on relative spectra from detectors at different baselines, and are therefore largely independent of reactor neutrino flux predictions. In view of the unclear situation of reactor flux predictions and reactor data at very short baselines (see, e.g., Ref. [24] for a recent discussion), we decided to include only the modern MBL reactor experiments Daya Bay, RENO, and Double Chooz. For the analysis of KamLAND long-baseline reactor data we replaced predicted neutrino fluxes by the spectrum measured in Daya Bay near detectors [25], which makes also our KamLAND analysis largely independent of flux predictions.
Our solar neutrino data includes previous data from radio-chemical and the SNO experiments, as well as updated exposures from Super-Kamiokande and Borexino, see appendix A for the detailed list and references. 1 Atmospheric neutrino data generically are difficult to analyze outside the experimental collaborations. We present below two separate global analyses, depending on the used atmospheric neutrino data. Our default analysis makes use of IceCube/DeepCore 3-year data [27] which can be re-analyzed using the information provided by the collaboration [28]. Especially in the context of the mass ordering determination, atmospheric neutrino data from Super-Kamiokande 1-4 [29] seems to provide important information. Unfortunately there is not enough information available to reproduce these results outside the collaboration. However, Super-Kamiokande has published the results of their analysis in the form of a tabulated χ 2 map [30], which we can combine with our global analysis. We will show the results of this combination as an alternative global fit. A detailed discussion of atmospheric neutrino data, including also the potential impact of an alternative IceCube analysis [31] will be presented in section 3.3.

Summary of global fit results
The results of our global fit are displayed in fig. 1 (one-dimensional ∆χ 2 curves) and fig. 2 (two-dimensional projections of confidence regions). In table 1 we give the best fit values as well as 1σ and 3σ confidence intervals for the oscillation parameters. We show two versions of the results. The default analysis is without Super-Kamiokande atmospheric neutrino data (SK-atm), and contains all the data for which a fit can be performed. For the alternative analysis, we add the pre-calculated ∆χ 2 table from SK-atm provided by the collaboration to our global fit, in order to illustrate the potential impact of these data. Let us summarize here the main features of the global fit result. More detailed discussions about how certain features emerge will be given in the following sections.  Figure 2. Global 3ν oscillation analysis. Each panel shows the two-dimensional projection of the allowed six-dimensional region after minimization with respect to the undisplayed parameters. The regions in the four lower panels are obtained from ∆χ 2 minimized with respect to the mass ordering. The different contours correspond to 1σ, 90%, 2σ, 99%, 3σ CL (2 dof). Coloured regions (black contour curves) are without (with) adding the tabulated SK-atm ∆χ 2 . Note that as atmospheric mass-squared splitting we use ∆m 2 31 for NO and ∆m 2 32 for IO.
Except for sin 2 θ 23 and δ CP the ∆χ 2 shapes are close to parabolic, indicating that the χ 2 approximation for the distribution should hold to good accuracy. The Monte Carlo studies performed in Refs. [11,32] indicate that also for sin 2 θ 23 , δ CP and the mass ordering the χ 2 approximation gives a reasonable estimate of the corresponding confidence level. Therefore, the ∆χ 2 values given below can be converted into an approximate number of standard deviations by the ∆χ 2 rule.
Defining the 3σ relative precision of the parameter by 2(x up − x low )/(x up + x low ), where x up (x low ) is the upper (lower) bound on a parameter x at the 3σ level, we obtain the following 3σ relative precisions (marginalizing over ordering): where the numbers between brackets show the impact of including SK-atm in the precision of that parameter determination. We notice that as ∆χ 2 shape for δ CP is clearly not gaussian this evaluation of its "precision" can only be taken as indicative.
Altogether the status of mass ordering discrimination, determination of sin 2 θ 23 , and the leptonic CP phase δ CP can be summarized as follows: • The best fit is for the normal mass ordering. Inverted ordering is disfavoured with a ∆χ 2 = 4.7 (9.3) without (with) SKatm.
• The best fit for the complex phase is at δ CP = 215 • . Compared to previous results (e.g., NuFIT 3.2 [12]), the allowed range is pushed towards the CP conserving value of 180 • , which now is only disfavoured with ∆χ 2 = 1.5 (1.8) without (with) SK-atm.
In table 1 we give the best fit values and confidence intervals for both mass orderings, relative to the local best fit points in each ordering. The global confidence intervals (marginalizing also over the ordering) are identical to the ones for normal ordering, which have also been used in eq. (2.1). The only exception to this statement is ∆m 2 3 in the analysis without SK-atm; in this case a disconnected interval would appear above 2σ corresponding to negative values of ∆m 2 3 (i.e., inverted ordering). Altogether we derive the following 3σ ranges on the magnitude of the elements of the leptonic mixing matrix:  Note that there are strong correlations between the elements due to the unitary constraint, see Ref. [33] for details on how we derive the ranges. The present status of leptonic CP violation is illustrated in figs. 2 and 3. In particular fig. 2 contains two projections of the confidence regions with δ CP on the vertical axis in which we observe the non-trivial correlations between δ CP and sin 2 θ 23 . In the left panel of fig. 3 we show the dependence of ∆χ 2 of the global analysis on the Jarlskog invariant which gives a convention-independent measure of CP violation [34], defined by: at 1σ (3σ) for both orderings. The preference of the present data for non-zero δ CP implies a best fit value J best CP = −0.019, which is favored over CP conservation with ∆χ 2 = 1.5 (1.8) without (with) SK-atm. These numbers can be compared with the size of the Jarlskog invariant in the quark sector, J quarks CP = (3.18 ± 0.15) × 10 −5 [35].

Status of comparison of results of solar experiments versus KamLAND
The analyses of the solar experiments and of KamLAND give the dominant contribution to the determination of ∆m 2 21 and θ 12 . We show in fig. 4 the present determination of these parameters from the global solar analysis in comparison with that of KamLAND data. The results of the solar neutrino analysis are shown for the two latest versions of the Standard Solar Model, namely the GS98 and the AGSS09 models [36] obtained with two different determinations of the solar abundances [37]. This clearly illustrates the independence of the results with respect to the solar modeling. There are two main differences compared to our previous published results in Ref. [11]. In what respects the KamLAND region it has shifted towards slightly smaller values of θ 12 . This effect arises mainly from the new reactor fluxes used in our analysis of the KamLAND data. As mentioned in section 2.1, in our calculation of the event rates in KamLAND we have replaced the predicted neutrino fluxes by the spectrum measured in Daya Bay near detectors [25] which is unfolded for detector and remaining oscillation effects. In Ref. [11] we used instead the unoscillated reactor determined by including in the fit the results from a compilation of short baseline reactor data. The net result is that the current unoscillated reactor fluxes are slightly lower and consequently a slightly higher survival probability is required to better fit the data. Since in the context of 3ν-oscillations P 3ν ee,KLAND = sin 4 θ 13 + cos 4 θ 13 1 − a larger survival probability implies smaller values of θ 12 . As a result the best-fit value of θ 12 determined by KamLAND, sin 2 θ 12,bf-Kam = 0.290, does not perfectly align with the corresponding best fit value from the solar neutrino analysis, sin 2 θ 12,bf-sol = 0.315. Statistically, however, this is a very small effect as the best fit value of sin 2 θ 12 = 0.315 lies at ∆χ 2 KamLAND

1.
In what respects the determination of ∆m 2 21 it has been a result of global analyses for several years already, that the value of ∆m 2 21 preferred by KamLAND is somewhat higher than the one from solar experiments. The tension arises from a combination of two effects: the well-known fact that none of the 8 B measurements performed by SNO, SK and Borexino shows any evidence of the low energy spectrum turn-up expected in the standard LMA-MSW [38,39] solution for the value of ∆m 2 21 favored by KamLAND; and the observation of a non-vanishing day-night asymmetry in SK, whose size is larger than the one predicted for the ∆m 2 21 value indicated of KamLAND. The new addition to this issue in the present analysis is the inclusion of the 2860-day energy spectrum of SK4 [40] (compared to the 2365-day energy spectrum used in [11]). For the day-night variation of the results we still use the SK4 2055-day day-night asymmetry [41] because SK has not presented any update concerning the day-night dependence of the observed rates. The inclusion of the new spectral data makes the lack of the turn-up effect slightly stronger (for example the best fit ∆m 2 21 of KamLAND was at ∆χ 2 solar = 4 in the analysis of [11] with the GS98 fluxes and it is now at ∆χ 2 solar = 4.7). For illustration of the relevance of the day-night variation results we plot in fig. 4 the corresponding results of the solar analysis without including the day-night asymmetry. . Left: Allowed parameter regions (at 1σ, 90%, 2σ, 99%, and 3σ CL for 2 dof) from the combined analysis of solar data for GS98 model (full regions with best fit marked by black star) and AGSS09 model (dashed void contours with best fit marked by a white dot), and for the analysis of KamLAND data (solid green contours with best fit marked by a green star) for fixed sin 2 θ 13 = 0.0224 (θ 13 = 8.6). We also show as orange contours the results of a global analysis for the GS98 model but without including the day-night information from SK. Right: ∆χ 2 dependence on ∆m 2 21 for the same four analyses after marginalizing over θ 12 .
3.2 θ 23 , δ CP and mass ordering from LBL accelerator and MBL reactor experiments The determination of the atmospheric parameters θ 23 and ∆m 2 3 is illustrated in fig. 5. We observe significant synergy from combining the various experiments, since the combined region is clearly smaller than any individual one. Moreover, the striking agreement of LBL accelerator and MBL reactor data in the determination of ∆m 2 3 within comparable accuracy is a non-trivial cross check of the 3-flavour oscillation paradigm. Let us now discuss in more detail how the indication of non-maximal mixing and preference for the second octant for θ 23 emerges.

Disappearance results and non-maximal θ 23
We focus first on LBL disappearance data. The ν µ survival probability is given to good accuracy by [42,43] where L is the baseline, E ν is the neutrino energy, and  3 ) plane using both appearance and disappearance data from MINOS (green), T2K (red), NOνA (brown), as well as IceCube/DeepCore (orange), and SK-atm (from the table provided by the experiment, marron line) and the combination of them (blue coloured region). In the left panels the constraint on θ 13 from the global fit (which is dominated by the reactor data) is imposed as a Gaussian bias. The right panels show regions in the (θ 13 , ∆m 2 3 ) plane using only Daya Bay (black), Reno (violet) and Double Chooz (magenta) reactor data, and their combination (blue coloured region). In all panels ∆m 2 21 , sin 2 θ 12 are fixed to the global best fit values. Contours are defined with respect to the global minimum of the two orderings.
Hence the survival probability is symmetric with respect to the octant of θ µµ , which implies symmetry around s 2 23 = 0.5/c 2 13 ≈ 0.51. This behaviour is visible in the left panels of fig. 6, which show the results of LBL accelerator disappearance data from MINOS, T2K, NOvA, separated into the neutrino and anti-neutrino data samples (for fixed value of θ 13 at the best fit and NO). While most of the shown data samples prefer maximal mixing (especially T2K and NOvA neutrino data), maximal mixing is disfavoured by MINOS neutrino data (∆χ 2 ≈ 2) and NOvA anti-neutrino data (∆χ 2 ≈ 6). This behaviour can be traced back to the number of events in the corresponding data samples observed at the dip of the survival probability: for maximal mixing the survival probability is zero at the dip and no events should be observed. Qualitatively similar behaviour are found for IO. . LBL accelerator ν µ disappearance data only, from MINOS, T2K, and NOvA, separated into neutrino and anti-neutrino data. Left panels correspond to LBL accelerator data with constraint on θ 13 from the global fit (which is dominated by the MBL reactor data) imposed as a Gaussian bias. In the right panels LBL data are consistently combined with MBL reactor data from Daya Bay, RENO, and Double Chooz. Upper panels show the ∆χ 2 as a function of sin 2 θ 23 , lower panels show confidence regions at 2σ (2 dof). All panels assume NO and ∆m 2 21 , sin 2 θ 12 are fixed to the global best fit values. Qualitatively similar behaviour is found in IO.
In the lower-left panel of fig. 6 we observe in addition a correlation between sin 2 θ 23 and ∆m 2 31 for the data which prefer non-maximal mixing: larger values of ∆m 2 31 imply more deviation from maximal mixing. As visible in fig. 5, also MBL reactor data provide an accurate determination of ∆m 2 3 , which, however, pushes slightly to larger values than LBL data. Because of the above mentioned correlation, this leads to an even stronger preference for non-maximal mixing, once LBL data are consistently combined with reactor data, as visible in the right panels of fig. 6: in combination with reactors, MINOS neutrino and NOvA anti-neutrino data disfavour maximal mixing with ∆χ 2 ≈ 7 and 9, respectively.

Appearance results, second θ 23 octant and δ CP
The preference for the second octant of θ 23 is driven by ν µ → ν e appearance channel in LBL experiments (available both for neutrinos and anti-neutrinos). Following Ref. [32], the appearance probability can be approximated by where V is the effective matter potential. In the above equations we have expanded in the small parameters s 13 , ∆m 2 21 L/E ν , and A, and used that for T2K and NOvA |∆m 2 3 |L/4E ν ≈ π/2. 2 Using the respective mean neutrino energies we find A ≈ 0.05 for T2K and an empirical value of A = 0.1 (for which this approximation works better) at NOvA. Correspondingly the number of observed appearance events in T2K and NOvA is approximately proportional to the oscillation probability: Taking all the well-determined parameters θ 13 , θ 12 , ∆m 2 21 , |∆m 2 3 | at their global best fit points, we obtain numerically C ≈ 0.28. The normalization constants N ν,ν calculated from our re-analysis of T2K and NOvA are given for the various appearance samples in table 2. Those values can be compared with the background subtracted observed number of events, which we also report in the table. Within this approximation, there are only the two parameters s 2 23 and sin δ CP , plus the discrete parameter o = ±1 encoding the mass ordering, to fit the appearance event numbers shown in table 2, with sin 2 θ 23 being constrained in addition from disappearance data. Note that C depends only on sin 2θ 23 , which varies by less than 2% for 0.42 < s 2 23 < 0.64, and can be taken as constant for our purposes. The general trends from eqs. (3.8) and (3.9) are the following: • Both neutrino and anti-neutrino events are enhanced by increasing s 2 23 .
• Values of sin δ CP +1 (−1) suppress (increase) neutrino events, and have the opposite effect for anti-neutrino events.
• For NO (IO) neutrino events are enhanced (suppressed) due to the matter effect, whereas anti-neutrino events are suppressed (enhanced).
• For NO (IO) the matter effect increases (decreases) the impact of δ CP for neutrinos, while the opposite happens for anti-neutrinos.  Table 2. Normalization coefficients N ν and Nν for eqs. (3.8) and (3.9) for approximations used to qualitatively describe the various appearance event samples used in our analysis for T2K and NOvA.
We also give the observed number of events, as well as the corresponding background subtracted event numbers, as reported in Refs. [44,45] The last two items are more important for NOvA than for T2K, due to larger matter effects in NOvA because of the longer baseline.
In fig. 7, the determination of s 2 23 from LBL data (including appearance) combined with reactor data is shown. In the upper panels only θ 13 is constrained by reactor data, whereas in the lower panels LBL and reactor data are combined consistently, including also ∆m 2 3 information. For the reasons explained above, lower panels show larger significance of non-maximality, but now the symmetry between the octants is broken by appearance data. Fig. 8 shows the ∆χ 2 dependence on δ CP for various data samples.
Let us consider first the T2K samples. We see from table 2 that in both neutrino samples (especially CC1π) the observed number of events after background subtraction is large compared to N ν , while the anti-neutrino number is low. Hence, we need to maximize the expression in eq. (3.8) and minimize eq. (3.9). Since neutrino data dominates over antineutrinos, a slight preference for s 2 23 > 0.5 appears (constrained by disappearance data), while at the same time sin δ CP ≈ −1 serves to maximize (minimize) neutrino (anti-neutrino) appearance, as visible in fig. 8.
For NOvA neutrino data, the coefficient N ν in eq. (3.8) is also somewhat low compared to the observed number of events minus background. For NO, the matter effect enhances neutrino events, and therefore, s 2 23 (around maximal mixing favoured in disappearance) and δ CP can be adjusted, such that the event numbers can always be fitted, so ∆χ 2 (δ CP ) from NOvA neutrino data alone is < 1 for NO, cf. fig. 8. For IO, however, the matter effect suppresses neutrino events, and therefore, preference for the second octant and sin δ CP ≈ −1 appears to maximize the term in the square-bracket in eq. (3.8). For NOvA antineutrino data, table 2 shows that the observed event number is of the order of Nν (only slightly higher). Consequently we observe for NO only a very mild preference for sin δ CP ≈ 1 just to enhance slightly the rate of anti-neutrinos. For IO, the matter effect enhances anti-neutrinos, and therefore, choosing the combinations (first θ 23 octant/sin δ CP ≈ 1) or (second θ 23 octant/sin δ CP ≈ −1) can fit the events, which leads to negligible ∆χ 2 (δ CP ) dependence for IO NOvA anti-neutrinos, cf. fig. 8. The combination of those effects for NO, leads to a disfavouring of sin δ CP ≈ −1 with ∆χ 2 ≈ 3.5 from NOvA, somewhat in contradiction of the T2K preferred region: with the non-maximality of θ 23 from antineutrinos plus the matter enhancement for neutrinos, sin δ CP ≈ −1 would predict too many neutrino events, and is therefore disfavoured.
The conclusion of those considerations lead to the preference of the second octant for θ 23 in the global analysis, as well as pushing the confidence interval for δ CP towards 180 • , which implies that CP conservation is allowed by the combined data with ∆χ 2 ≈ 1.5.

Preference for normal ordering
An important result of the present global fit is the growing significance of the preference for the normal mass ordering. This indication emerges by a subtle interplay of various subsets of the global data. Sensitivity to the mass ordering is provided by the matter effect [38,39,46] in oscillations with ∆m 2 3 , observable in LBL accelerator and atmospheric neutrino experiments, as well as the comparison oscillations in the ν e and ν µ disappearance channels [43,47,48]. Let us first discuss the indication coming from LBL accelerator experiments. We find that T2K + the θ 13 constraint from reactors disfavours IO by ∆χ 2 ≈ 4, see upper panels of figs. 7, 8 and 9. This can be understood from the numbers in table 2 and eqs. (3.8) and (3.9), where the matter effect for NO helps to increase (decrease) events for neutrinos (anti-neutrinos). NOvA data + the θ 13 constraint also disfavours IO by about 2 units in χ 2 , driven by neutrino data, while anti-neutrinos are insensitive to the ordering, cf. fig. 8. Interestingly, by combining T2K, NOvA, and MINOS, decreases the ∆χ 2 of IO to about 2. An explanation for this effect is the slight tension between NOvA and T2K in the determination of δ CP for NO visible in fig. 8. This leads to a worse fit of NO compared to IO, where both experiments prefer the same region for δ CP .    An interesting additional effect sensitive to the mass ordering has been pointed out in Refs. [43,47]: the ν µ disappearance probability is symmetric with respect to the sign of ∆m 2 µµ given in eq. (3.4), while ν e disappearance is symmetric with respect to a slightly different effective mass-squared difference: ∆m 2 ee = cos 2 θ 12 ∆m 2 31 + sin 2 θ 12 ∆m 2 32 . (3.10) Hence, from a precise determination of the oscillation frequencies in ν µ and ν e disappear-ance experiments, information on the sign of ∆m 2 3 can be obtained. 3 Indeed, we observe in fig. 9 that this effect already contributes notably to the mass ordering discrimination in present data: the upper panels show the determination of ∆m 2 3 from the individual LBL accelerator experiments (ν µ disappearance) compared to the one from MBL reactors (ν e disappearance). We have verified that those curves are indeed symmetric with respect to the sign of ∆m 2 µµ and ∆m 2 ee , respectively, within excellent accuracy. When displaying them for common parameters (∆m 2 3 in fig. 9), we observe that the agreement is better for NO than for IO. The difference between the upper and lower panels in the ∆χ 2 for IO is largely due to this ∆m 2 3 effect. We see that the ∆χ 2 for the LBL combination is pushed from 2 to about 4.5, when combined consistently with reactor data taking into account the ∆m 2 3 dependence. In summary, we obtain from LBL+reactor data a preference for NO at about 2σ. As mentioned in section 2.2, this gets further enhanced by atmospheric neutrino data, with the main contribution from Super-Kamiokande, leading to the exclusion of IO at about 3σ, see fig. 1. In the following subsection we discuss in more detail various aspects of the atmospheric neutrino analyses from IceCube and Super-Kamiokande.

Treatment of atmospheric results from Super-Kamiokande and Deep-Core
In what respects the atmospheric neutrino data, in our default analysis -full lines in figs. 1 and 3 (one-dimensional ∆χ 2 curves) and coloured regions in fig. 2 (two-dimensional projections of confidence regions) -we include the results of the Deep-Core 3-years data of Ref. [27,28] (which we refer here as DC16) for which the collaboration has provided enough information on their effective areas to allow for our own reanalysis. Its impact in the parameter determination obtained from the combination of solar, reactor and LBL data is very marginal, see fig. 10, which displays as example its contribution to the determination of ∆m 2 3 and the ordering where it adds about 0.3 units to χ 2 min of IO because of the slightly better matching between the ∆m 2 3 from reactor+LBL experiments with that of DC in NO. 4 In this respect it is interesting to notice that the ICECUBE collaboration has recently published the results of a dedicated analysis of another set of three-years data [31,49] leading to a better determination of ∆m 2 3 (which we refer to as DC17). Unfortunately we cannot reproduce this analysis because the corresponding effective areas have not been made public. The experiment has only made available the bi-dimensional χ 2 map (as a function of ∆m 2 3 and sin 2 θ 23 for a fixed value of sin 2 θ 13 = 0.0217 and δ CP = 0) corresponding to that analysis. Strictly this cannot be added in the global analysis without making some assumption about their possible θ 13 and δ CP dependence. Still, to illustrate the possible impact of using these results we show also in fig. 10 the corresponding contribution to the determination of ∆m 2 3 and the ordering obtained by naively adding their χ 2 map to our results of the global reactor+LBL experiments (neglecting any possible depen-   dence on the fixed parameters). As seen in the figure, using the DC17 results in the global combination disfavours IO by ∼ 1.2 additional units of χ 2 . One must notice, however, that the ICECUBE collaboration has recently performed a reanalysis of the same data sample which leads to similar precision but a somewhat shifted range for ∆m 2 3 [50]. In what respects to the results of Super-Kamiokande, in the last five years the collaboration has developed a more sophisticated analysis method for their atmospheric neutrino data with the aim of constructing ν e +ν e enriched samples which are then further classified into ν e -like andν e -like subsamples, thus increasing the sensitivity to subleading parameters such as the mass ordering and δ CP . The official results obtained with this method were published in Ref. [29] and show -once θ 13 is constrained to be within the range determined by reactor experiments -a preference for NO with a ∆χ 2 (IO) = 4.3, variation of χ 2 (δ CP ) with the CP phase at the level of ∼ 90% CL (with favouring δ CP ∼ 270 • ), and a slight favouring of the second octant of θ 23 (see fig.14 in Ref. [29]).
Unfortunately with the information at hand we are not able to reproduce the elements driving the main dependence on these subdominant oscillation effects in our own reanalysis of the data samples which can be simulated outside of the collaboration. However, Super-Kamiokande has also published the results of their analysis in the form of a tabulated χ 2 map [30] as a function of the four relevant parameters ∆m 2 3 , θ 23 , θ 13 , and δ CP which we can add to our global analysis χ 2 in the multidimensional parameter space in a fully consistent form and then perform the corresponding parameter marginalization to obtained the combined one-dimensional or two dimensional parameter ranges. The results of such combination are shown as dashed curves in in figs. 1 and 3 (one-dimensional ∆χ 2 curves) and void regions in fig. 2 (two-dimensional projections of confidence regions). As can be seen from fig. 1, adding the SK-atm χ 2 information results into: • Increase of the χ 2 min for IO by 4.6 units (from 4.7 to 9.3).
In other words, as the SK-atm tendencies for these subdominant effects are very well aligned with those from the combination of LBL experiments (currently dominated by T2K), their impact in the determination of δ CP and θ 23 in the global analysis is almost equivalent to just adding for each of those parameters their marginalized χ 2 (with fixed θ 13 at the reactor value) to that from the global analysis without SK-atm.

Projections on neutrino mass scale observables
Oscillation experiments provide information on the mass-squared splittings ∆m 2 ij and on the leptonic mixing angles U ij , but they are insensitive to the absolute mass scale for the neutrinos. Of course, the results of an oscillation experiment do provide a lower bound on the heavier mass in ∆m 2 ij , |m i | ≥ ∆m 2 ij for ∆m 2 ij > 0, but there is no upper bound on this mass. In particular, the corresponding neutrinos could be approximately degenerate at a mass scale that is much higher than ∆m 2 ij . Moreover, there is neither an upper nor a lower bound on the lighter mass m j .
Information on the neutrino masses, rather than mass differences, can be extracted from kinematic studies of reactions in which a neutrino or an anti-neutrino is involved. In the presence of mixing the most relevant constraint comes from the study of the end point (E ∼ E 0 ) of the electron spectrum in Tritium beta decay 3 H → 3 He + e − +ν e . This spectrum can be effectively described by a single parameter, m νe , if for all neutrino states E 0 − E m i : where the second equality holds if unitarity is assumed and m 0 = m 1 (m 3 ) in NO (IO) denotes the lightest neutrino mass. At present we only have an upper bound, m νe ≤ 2.2 eV at 95% CL [51], which is expected to be superseded soon by KATRIN [52] with about one order of magnitude improvement in sensitivity. Direct information on neutrino masses can also be obtained from neutrinoless double beta decay (A, Z) → (A, Z + 2) + e − + e − . This process violates lepton number by two units, hence in order to induce the 0νββ decay, neutrinos must be Majorana particles. In particular, for the case in which the only effective lepton number violation at low energies is induced by the Majorana mass term for the neutrinos, the rate of 0νββ decay is proportional to the effective Majorana mass of ν e :  [55].
Neutrino masses have also interesting cosmological effects. In general, cosmological data mostly give information on the sum of the neutrino masses, i m i , while they have very little to say on their mixing structure and on the ordering of the mass states.
Correlated information on these three probes of the neutrino mass scale can be obtained by mapping the results from the global analysis of oscillations presented previously. We show in fig. 11 the present status of this exercise. The relatively large width of the regions in the right panel are due to the unknown Majorana phases. Thus from a positive determination of two of these probes, in principle information can be obtained on the value of the Majorana phases and/or the mass ordering [56,57].

Conclusions
We have presented the results of the updated (as of fall 2018) analysis of relevant neutrino data in the framework of mixing among three massive neutrinos. We have shown our results for two analyses. The first contains our own statistical combination of all the experimental data for which we are able to reproduce the results of the partial analysis performed by the different experiments, and therefore does not include the information of the Super-Kamiokande atmospheric neutrino data. In the second analysis we combine the likelihood of the first one with the four-dimensional χ 2 map provided by Super-Kamiokande for the analysis of their atmospheric data. Quantitatively the present determination of the two mass differences, three mixing angles and the relevant CP violating phase for the two analysis is listed in table 1, and the corresponding leptonic mixing matrix is given in Eq. (2.2). In both analysis the maximum allowed CP violation in the leptonic sector parametrized by the Jarlskog determinant is J max CP = 0.0333 ± 0.0006 (±0.0019) at 1σ (3σ). We have performed a detail study of the role of the different data samples and their correct combination in the determination of the less known parameters, θ 23 , δ CP and the ordering in section 3. We can summarize the main conclusions in this section as follows: • The long standing tension between the best ∆m 2 21 determined in the solar neutrino analysis and that from KamLAND persists. The inclusion of latest spectral data from SK4 and the use of Daya Bay near detector data for reactor flux normalization in KamLAND has made this tension slightly stronger, but it is still a ∼ 2σ effect.
• We obtain preference for the second octant of θ 23 in the global analysis with a best fit at sin 2 θ 23 = 0.58. There are two effects contributing to this results: -While most data samples in ν µ andν µ disappearance at LBL prefer close to maximal mixing (especially T2K and NOvA neutrino data), maximal mixing is disfavoured by MINOS neutrino data (∆χ 2 ≈ 2) and NOvA anti-neutrino data (∆χ 2 ≈ 6). This disfavouring increases when fully combining with the reactor neutrino determination of ∆m 2 3 to ∆χ 2 ≈ 7 and 9, respectively. -The appearance results both in T2K and NOvA (SK-atm adds in the same direction) favour the second octant. The final value of the best fit θ 23 results of this effect in combination with the substantial non-maximality favoured by NOvA anti-neutrino and MINOS neutrino disappearance data.
• The determination of δ CP is mostly driven by T2K neutrino and anti-neutrino appearance results which favour δ CP ∼ 3π/2 and disfavours δ CP ∼ π/2 for both NO and IO. NOvA neutrino appearance data align with this behaviour, and are more statistically significant in IO. On the contrary NOvA anti-neutrino appearance data are better (worse) described with δ CP ∼ π/2 (δ CP ∼ 3π/2) in NO. This slight tension results into a shift of the best fit to δ CP = 215 • in NO. So the allowed range is pushed towards the CP conserving value of 180 • , which now is only disfavoured with ∆χ 2 2.
• Regarding the mass ordering: -Both T2K and NOvA prefer NO individually: T2K (NOvA) + the θ 13 constraint from reactors disfavours IO by ∆χ 2 ≈ 4 (2), but combining T2K, NOvA (and MINOS) decreases the ∆χ 2 of IO to about 2. This is a consequence of the slight tension between NOvA and T2K in the determination of δ CP for NO.
-Additional sensitivity to the ordering is found by the precise determination of the oscillation frequency in ν µ and ν e disappearance data at LBL and reactors, respectively. This effect increases ∆χ 2 of IO from 2 to about 4.5.
-Inclusion of the atmospheric neutrino results (mainly from SK) further increases ∆χ 2 of IO to the 3σ level.
Future updates of this analysis will be provided at the NuFIT website quoted in Ref. [12].

B Technical details and validation cross checks
This this appendix we provide some details on our analysis of the most recent data from accelerator and reactor experiments, and show that we can reproduce the results of the experimental collaborations with good accuracy, when using the same assumptions in the analysis.

B.1 T2K
The predicted number of events in the T2K far detector in a given energy bin i and for a given channel α can be calculated as • N bkg,i is the number of background events in that bin, which we have extracted from Ref. [71] and consistently re-scaled to the latest exposure. If there is a neutrino component, its oscillation has to be consistently included.
• [E i , E i+1 ] are the bin limits.
• E rec is the reconstructed neutrino energy.
• E ν is the true neutrino energy.
• R(E rec , E ν ) is the energy reconstruction function, that we take to be Gaussian.
• ε is the detection efficiency, which is adjusted to reproduce the observed spectra in Ref. [72].
For the antineutrino channel, one has to switch ν byν. If we assume a Poissonian χ 2 with the data points in Ref. [72], add an overall normalisation systematic uncertainty 5 and combine all the data, we get the contours in fig. 12. The other oscillation parameters are fixed as specified in Ref. [72] and the reactor uncertainty on θ 13 is included as a Gaussian bias and marginalised over. Finally, following the ad-hoc procedure described in section 8.4.2 in Ref. [72], the disappearance ∆m 2 32 contours are manually Gaussian-smeared with a standard deviation σ = 4.1 · 10 −5 eV 2 .

B.2 NOvA
For NOvA the predicted number of events is given, as for T2K, by Eq. (B.1). Fluxes and detector response are extracted from Refs. [74,75]. We take the data points and backgrounds from Ref. [75] and with those construct a Poissonian χ 2 including an normalisation systematic uncertainty 6 and we get the contours in fig. 13 when the undisplayed oscillation parameters are fixed as specified in Ref. [75]. In particular the reactor uncertainty on θ 13 is included as a Gaussian bias and marginalised over. . Different projections to our fit to latest T2K data (dashed lines) compared to the corresponding results of the experimental collaboration (solid curves), as presented in Ref. [72], when adopting the same assumptions. All unshown parameters are marginalised over.

B.3 Daya Bay
The Daya Bay 3ν oscillation analysis is based on [21], the data is taken from the supplementary material provided on arXiv.
In each experimental hall in Daya Bay there are more than one detector, so the number of events in hall H, N H i is a sum of the contributions in all the detectors in the hall. The predicted numbers of events in a detector d in an energy bin i is computed as (B. 2) The indices i, r, d, iso refer to the energy bin, reactor, detector, and fissible isotope, respectively. d are the detector efficiencies, including ε µ × ε m multiplied by the life time days (taken both from table I in [21]) as well as the relative difference of target protons ∆N p in each detector, obtained from table VI in [76]. L rd are the baselines between the reactor  Figure 13. Different projections to our fit to latest NOvA data (dashed lines) compared to the corresponding results of the experimental collaboration (solid curves), as presented in Ref. [45], when adopting the same assumptions. The upper plot only includes disappearance data, whereas the bottom left plot corresponds to ν e appearance and the bottom right plot toν e appearance. r and detector d, obtained from table I in [76]. E ν and E rec are the true and reconstructed neutrino energy, which are related by the detector response function R(E rec , E prompt ), provided in the complementary material to [21]. The relation between the prompt energy E prompt and the neutrino energy E ν is E ν = m n − m p − m e + E prompt . σ(E ν ) is the Inverse Beta Decay cross section computed performing the integral over cos θ of the differential cross section in [77]. φ iso (E ν ) are the Huber-Mueller flux predictions [78,79] and f iso are the fission fractions. For each isotope, f iso is computed as the average of the fission fractions in table 9 of Ref. [25]. Following section 2.6 of [25], we apply non-equilibrium corrections by adding the relative correction from table VII of [79] to the fluxes [78,79]. Since Daya Bay has run for a long period we take the row corresponding to 450 days. P rd νe→νe (E ν ) is the oscillation probability. The global constant N will cancel when taking ratios of event numbers.
Our DayaBay χ 2 is based on the ratios of the observed spectra in experimental halls 3 and 1 as well as 2 and 1: 3), η is the vector of pull parameters and V η is the pull correlation matrix which accounts for the systematic uncertainties and their correlations. In order to reproduce the Daya Bay results [21], systematics in detection efficiency, relative energy scale and cosmogenic Li-He background have to be taken into account. The detection efficiency and relative energy scale uncertainties are given by 0.13% and 0.2%, respectively [21,76], and the Li-He background uncertainty is given by 30% [21]. We include also the uncertainties in accidental (1%) and fast neutron (13% (17%) in EH1 and EH2 (EH3)) backgrounds [76], which however play only a subleading role. We take into account the 5 systematics in each of the 3 experimental hall as uncorrelated, so V η is a diagonal 15 × 15 matrix. In order to use the correct uncertainties in each experimental hall we divide the detector uncertainties by √ 2 in EH1 and EH2 and by √ 4 in EH3, since there are 2 and 4 detectors, respectively In the left panel of fig. 14, our re-analysis is compared to the one published in [21].

B.4 RENO
The RENO 3ν oscillation analysis is based on [22]. The number of events in the near and far detectors are computed as in eq. (B.2), using the Daya Bay response function. The average fission fractions are taken from [80], the baselines from [81], and life time days can be found in [22]. In order to compute the total relative efficiency between the near and far detectors, a normalization to the total number of predicted events without oscillations in the far detector is performed. The RENO χ 2 is based on the far/near spectral ratio is implemented as follows:  Δχ 2 Figure 14. Our fit to DayaBay and RENO (black-dashed lines) compared to the results of the experimental collaborations (solid lines), as published in Refs. [21] and [22], respectively, when adopting the same assumptions. Following the collaborations, on the vertical axes, the parameter ∆m 2 ee defined in eq. (3.10) is used.
Here, O . η is the vector of the pull parameters and V η is the pull correlation matrix which accounts for the systematic uncertainties associated to each pull parameter and their correlations.
Ref. [22] quotes systematic uncertainties in the relative detection efficiency and relative energy scale of 0.21% and 0.15%, respectively. The cosmogenic Li-He background plays no significant role, but is included as well with a relative uncertainty of 5% (8%) for the near (far) detector, cf. table I of [22]. In order to match precisely the results shown in fig. 3 of [22], an extra factor of 0.984 to the Far/Near ratio has to be included and the relative detection efficiency uncertainty has to be increased by a factor of 1.4. Our re-analysis of RENO data is compared to the official one in the right panel of fig. 14.