Updated fit to three neutrino mixing: exploring the accelerator-reactor complementarity

We perform a combined fit to global neutrino oscillation data available as of fall 2016 in the scenario of three-neutrino oscillations and present updated allowed ranges of the six oscillation parameters. We discuss the differences arising between the consistent combination of the data samples from accelerator and reactor experiments compared to partial combinations. We quantify the confidence in the determination of the less precisely known parameters $\theta_{23}$, $\delta_\text{CP}$, and the neutrino mass ordering by performing a Monte Carlo study of the long baseline accelerator and reactor data. We find that the sensitivity to the mass ordering and the $\theta_{23}$ octant is below $1\sigma$. Maximal $\theta_{23}$ mixing is allowed at slightly more than 90% CL. The best fit for the CP violating phase is around $270^\circ$, CP conservation is allowed at slightly above $1\sigma$, and values of $\delta_\text{CP} \simeq 90^\circ$ are disfavored at around 99% CL for normal ordering and higher CL for inverted ordering.


Introduction
Experiments measuring the flavor composition of solar neutrinos, atmospheric neutrinos, neutrinos produced in nuclear reactors and in accelerators have established that lepton flavor is not conserved in neutrino propagation, but it oscillates with a wavelength depending on distance and energy, because neutrinos are massive and the mass states are admixtures of the flavor states [1,2], see Ref. [3] for an overview.
With the exception of a set of unconfirmed "hints" of possible eV scale mass states (see Ref. [4] for a recent review), all the oscillation signatures can be explained with the three flavor neutrinos (ν e , ν µ , ν τ ), which can be expressed as quantum superpositions of three massive states ν i (i = 1, 2, 3) with masses m i . This implies the presence of a leptonic mixing matrix in the weak charged current interactions [5,6] which can be parametrized 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π]. Here P is a diagonal matrix which is the identity if neutrinos are Dirac fermions and it contains two additional phases if they are Majorana fermions, and plays no role in neutrino oscillations [7,8]. In this convention there are two non-equivalent orderings for the neutrino masses which can be chosen to be: normal ordering (NO) with m 1 < m 2 < m 3 , and inverted ordering (IO) with m 3 < m 1 < m 2 . Furthermore the data shows a relatively large 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 introduced in Ref. [9] and present our results in terms of the variable ∆m 2 3 , with = 1 for NO and = 2 for IO. Hence, ∆m 2 3 = ∆m 2 31 > 0 for NO and ∆m 2 3 = ∆m 2 32 < 0 for IO, i.e., it corresponds to the mass splitting with the largest absolute value.
In this article, we present an up-to-date (as of fall 2016) global analysis of neutrino data in the framework of three-neutrino oscillations. Alternative recent global fits have been presented in Refs. [10,11]. With current data from the accelerator long-baseline experiments MINOS, T2K, NOνA and modern reactor experiments like Daya-Bay, RENO, and Double-Chooz, their complementarity anticipated more than a decade ago [12][13][14] has become a reality, and the combined analysis starts to show some sensitivity to subtle effects like the θ 23 octant or the δ CP phase (though still at low statistical significance).
The outline of the paper is as follows: In Sec. 2.1 we describe the data samples included in our analysis (see also Appendix A for a schematic list). The presently allowed ranges of the six oscillation parameters are given in Sec. 2.2 assuming that ∆χ 2 follows a χ 2distribution, while Sec. 2.3 contains the corresponding measures of CP violation in terms of the leptonic Jarlskog invariant and the leptonic unitarity triangle. Deviations from the Gaussian approximation of the confidence intervals for θ 23 and δ CP and the confidence level for the mass ordering determination are quantified in Sec. 4. Several issues appearing in the present analysis are discussed in Sec. 3, in particular about the consistent combination of results from long baseline accelerator experiments with reactors results, now that both provide comparable precision in the determination of the relevant mass-squared difference. We also give the updated status on the ongoing tension in the ∆m 2 21 determination from solar experiments versus KamLAND, and comment on the stand-by in the analysis of the Super-Kamiokande atmospheric data. Sec. 5 contains the summary of our results.

Data samples analyzed
In the analysis of solar neutrino data we consider the total rates from the radiochemical experiments Chlorine [15], Gallex/GNO [16] and SAGE [17], the results for the four phases of Super-Kamiokande [18][19][20][21][22], the data of the three phases of SNO included in the form of the parametrization presented in [23], and the results of both Phase-I and Phase-II of Borexino [24][25][26].
Results from long baseline (LBL) accelerator experiments include the final energy distribution of events from MINOS [27,28] in ν µ andν µ disappearance and ν e andν e appearance channels, as well as the latest energy spectrum for T2K in the same four channels [29,30] and for NOνA on the ν µ disappearance and ν e appearance neutrino modes [31].
Data samples onν e disappearance from reactor include the full results of the long baseline reactor data in KamLAND [32], as well as the results from medium baseline reactor experiments from CHOOZ [33] and Palo Verde [34]. Concerning running experiments we include the latest spectral data from Double-Chooz [35] and Daya-Bay [36], while for RENO we use the total rates obtained with their largest data sample corresponding to 800 days of data-taking [37].
For the analysis of atmospheric neutrinos we include the results from IceCube/DeepCore 3-year data [48].
The above data sets constitute the samples included in our NuFIT 3.0 analysis. For Super-Kamiokande atmospheric neutrino data from phases SK1-4 we will comment on our strategy in Sec. 3.3. A full list of experiments including the counting of data points in each sample can be found in Appendix A.

Results: oscillation parameters
The results of our standard analysis are presented in Figs. 1 and 2 where we show projections of the allowed six-dimensional parameter space. 1 In all cases when including reactor experiments we leave the normalization of reactor fluxes free and include data from shortbaseline (less than 100 m) reactor experiments. In our previous analysis [9,50] we studied the impact of this choice versus that of fixing the reactor fluxes to the prediction of the latest calculations [51][52][53]. As expected, the overall description is better when the flux normalization f flux is fitted against the data. We find χ 2 (f flux fix) − χ 2 (f flux fit) 6 which is just another way to quantify the well-known short baseline reactor anomaly to be ∼ 2.5σ. However, the difference in the resulting parameter determination (in particular for θ 13 ) between these two reactor flux normalization choices has become marginal, since data from the reactor experiments with near detectors such as Daya-Bay, RENO and Double-Chooz (for which the near-far comparison allows for flux-normalization independent analysis) is now dominant. Consequently, in what follows we show only the ∆χ 2 projections for our standard choice with fitted reactor flux normalization.
The best fit values and the derived ranges for the six parameters at the 1σ (3σ) level are given in Tab. 1. For each parameter x the ranges are obtained after marginalizing with respect to the other parameters 2 and under the assumption that ∆χ 2 marg (x) follows a χ 2 distribution. Hence the 1σ (3σ) ranges are given by the condition ∆χ 2 marg (x) = 1 (9). It is known that because of its periodic nature and the presence of parameter degeneracies 1 ∆χ 2 tables from the global analysis corresponding to all 1-dimensional and 2-dimensional projections are available for download at the NuFIT website [49].
2 In this paper we use the term "marginalization" over a given parameter as synonym for minimizing the the statistical distribution of the marginalized ∆χ 2 for δ CP and θ 23 (and consequently the corresponding CL intervals) may be modified [54,55]. In Sec. 4 we will discuss and quantify these effects. In Tab. 1 we list the results for three scenarios. In the first and second columns we assume that the ordering of the neutrino mass states is known a priori to be Normal or Inverted, respectively, so the ranges of all parameters are defined with respect to the minimum in the given scenario. In the third column we make no assumptions on the ordering, so in this case the ranges of the parameters are defined with respect to the global minimum (which corresponds to Normal Ordering) and are obtained marginalizing also over the ordering. For this third case we only give the 3σ ranges. In this case the range of ∆m 2 3 is composed of two disconnected intervals, one containing the absolute minimum (NO) and the other the secondary local minimum (IO).

Results: leptonic mixing matrix and CP violation
From the global χ 2 analysis described in the previous section and following the procedure outlined in Ref. [56] one can derive the 3σ ranges on the magnitude of the elements of the Note that there are strong correlations between the elements due to the unitary constraint. The present status of the determination of leptonic CP violation is illustrated in Fig. 3. In the left panel we show the dependence of ∆χ 2 of the global analysis on the Jarlskog invariant which gives a convention-independent measure of CP violation [57], defined as usual by: Im U αi U * αj U * βi U βj ≡ J max CP sin δ = cos θ 12 sin θ 12 cos θ 23 sin θ 23 cos 2 θ 13 sin θ 13 sin δ (2.2) where we have used the parametrization in Eq. (1.1). Thus the determination of the mixing angles yields at present a maximum allowed CP violation 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.033, which is favored over CP conservation with ∆χ 2 = 1.7. These numbers can be compared with the size of the Jarlskog invariant in the quark sector, which is determined to be J quarks CP = (3.04 +0. 21 −0.20 ) × 10 −5 [58]. In Fig. 4 we recast the allowed regions for the leptonic mixing matrix in terms of one leptonic unitarity triangle. Since in the analysis U is unitary by construction, any given pair of rows or columns can be used to define a triangle in the complex plane. In the figure we show the triangle corresponding to the unitarity conditions on the first and third columns which is the equivalent to the one usually shown for the quark sector.
NuFIT 3.0 (2016) NO IO Figure 4. Leptonic unitarity triangle for the first and third columns of the mixing matrix. After scaling and rotating the triangle so that two of its vertices always coincide with (0, 0) and (1, 0) we plot the 1σ, 90%, 2σ, 99%, 3σ CL (2 dof) allowed regions of the third vertex. Note that in the construction of the triangle the unitarity of the U matrix is always explicitly imposed. The regions for both orderings are defined with respect to the common global minimum which is in NO.
In this figure the absence of CP violation implies a flat triangle, i.e., Im(z) = 0. As can be seen, for NO the horizontal axis crosses the 1σ allowed region, which for 2 dof corresponds to ∆χ 2 ≤ 2.3. This is consistent with the present preference for CP violation, χ 2 (J CP = 0) − χ 2 (J CP free) = 1.7 mentioned above. We will comment on the statistical interpretation of this number in Sec. 4.

Issues in present analysis
The 3ν fit results in the previous section provide a statistically satisfactory description of all the neutrino oscillation data considered. There are however some issues in the determination of some of the parameters which, although not of statistical significance at present, deserve some attention. 2 21 in 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 . 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. This tension arises from a combination of two effects which have not changed significantly over the last lustrum: a) 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 [59,60] solution for the value of ∆m 2 21 favored by KamLAND; b) the observation of a non-vanishing daynight asymmetry in SK, whose size is larger than the one predicted for the ∆m 2 21 value indicated of KamLAND (for which Earth matter effects are very small). In Ref. [9] we discussed the differences in the physics entering in the analyses of solar and KamLAND data which are relevant to this tension, and to which we refer the reader for details. Here . 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 θ 13 = 8.5 • . Right: ∆χ 2 dependence on ∆m 2 21 for the same three analyses after marginalizing over θ 12 .

Status of ∆m
for sake of completeness we show in Fig. 5 the quantification of this tension in our present global analysis. As seen in the figure, the best fit value of ∆m 2 21 of KamLAND lays at the boundary of the 2σ allowed range of the solar neutrino analysis.
Also for illustration of the independence of these results with respect to the solar modeling, the solar neutrino regions are shown for two latest versions of the Standard Solar Model, namely the GS98 and the AGSS09 models [61] obtained with two different determinations of the solar abundances [62]. Figure 6 illustrates the contribution to the present determination of ∆m 2 3 from the different data sets. In the left panels we focus on the determination from long baseline experiments, which is mainly from ν µ disappearance data. We plot the 1σ and 2σ allowed regions (2 dof) in the dominant parameters ∆m 2 3 and θ 23 . As seen in the figure, although the agreement between the different experiments is reasonable, some "tension" starts to appear in the determination of both parameters among the LBL accelerator experiments. In particular we see that the recent results from NOνA, unlike those from T2K, favor a non-maximal value of θ 23 . It is important to notice that in the context of 3ν mixing the relevant oscillation probabilities for the LBL accelerator experiments also depend on θ 13 (and on the θ 12 and ∆m 2 21 parameters which are independently well constrained by solar and KamLAND data). To construct the regions plotted in the left panels of Fig. 6 NuFIT 3.0 (2016) Figure 6. Determination of ∆m 2 3 at 1σ and 2σ (2 dof), where = 1 for NO (upper panels) and = 2 for IO (lower panels). The left panels show regions in the (θ 23 , ∆m 2 3 ) plane using both appearance and disappearance data from MINOS (green line), T2K (red lines), NOνA (light blue lines), as well as IceCube/DeepCore (orange lines) and the combination of them (colored regions). In these 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 lines), reactor data without Daya-Bay (violet lines), and their combination (colored regions). In all panels solar and KamLAND data are included to constrain ∆m 2 21 and θ 12 . Contours are defined with respect to the global minimum of the two orderings.

∆m 2 3 determination in LBL accelerator experiments versus reactors
currently followed by the LBL accelerator experiments: we marginalize with respect to θ 13 , taking into account the information from reactor data by adding a Gaussian penalty term to the corresponding χ 2 LBL . This is not the same as making a combined analysis of LBL and reactor data as we will quantify in Sec. 3.2.1.
Concerning ν e disappearance data, the total rates observed in reactor experiments at different baselines can provide an independent determination of ∆m 2 3 [50,63]. On top of this, the observation of the energy-dependent oscillation effect due to θ 13 now allows to further strengthen such measurement. In the right panels of Fig. 6 we show therefore the allowed regions in the (θ 13 , ∆m 2 3 ) plane based on global data on ν e disappearance. The violet contours are obtained from all the medium-baselines reactor experiments with the exception of Daya-Bay; these regions emerge from the baseline effect mentioned above plus spectral information from Double-Chooz. 3 The black contours are based on the energy spectrum in Daya-Bay, whereas the colored regions show the combination.
By comparing the left and right panels of Fig. 6 we observe that the combined ν µ and ν e disappearance experiments provide a consistent determination of |∆m 2 3 | with similar precision. However when comparing the region for each LBL experiment with that of the reactor experiments we find some dispersion in the best fit values and allowed ranges. This is more clearly illustrated in the upper panels of Fig. 7, where we plot the one dimensional projection of the regions in Fig. 6 as a function of ∆m 2 3 after marginalization over θ 23 for each of the LBL experiments and for their combination, together with that from reactor data after marginalization over θ 13 . The projections are shown for NO(right) and IO(left). Let us stress that the curves corresponding to LBL experiments in the upper panels of Fig. 7 (as well as those in the upper panels of Figs. 8 and 9) have been obtained by a partial combination of the information on the shown parameter (∆m 2 3 or θ 23 or δ CP ) from LBL with that of θ 13 from reactors, because in these plots only the θ 13 constraint from reactors is imposed while the dependence on ∆m 2 3 is neglected. This corresponds to the 1-dim projections of the function: However, since reactor data also depends on ∆m 2 3 the full combination of reactor and LBL results implies that one must add consistently the χ 2 functions of the LBL experiment with that of reactors evaluated the same value of ∆m 2 3 , this is We discuss next the effect of combining consistently the information from LBL and reactor experiments in the present determination of θ 23 , δ CP and the ordering.
3.2.1 Impact on the determination of θ 23 , mass ordering, and δ CP We plot in the lower panels of Figs. 7-9 the one dimensional projections of ∆χ 2 LBL+REA for each of the parameters θ 23 , δ CP , ∆m 2 3 (marginalized with respect to the two undisplayed parameters) for the consistent LBL+REA combinations with both the information on θ 13 and ∆m 2 3 from reactors included, Eq. (3.2). As mentioned before, the curves in the upper panels for these figures show the corresponding 1-dimensional projections for the partial combination, in which only the θ 13 constraint from reactors is used, Eq. (3.1). For each experiment the curves in these figures are defined with respect to the global minimum of the two orderings, so the relative height of the minimum in one ordering vs the other gives a measure of the ordering favored by each of the experiments.    Comparing the upper and lower panels in Figs. 7, 8 and 9 one sees how the contribution to the determination of the mass ordering, the octant and non-maximality of θ 23 , and the presence of leptonic CP violation of each LBL experiment in the full LBL+REA combination (Eq. 3.2) can differ from those derived from the LBL results imposing only the θ 13 constraint from reactors (Eq. 3.1). This is due to the additional information on ∆m 2 3 from reactors, which is missing in this last case. In particular: • (NO) 0.4 (1.7) for LBL = NOνA (T2K). This is in agreement with the analyses shown by the collaborations for example in Refs. [29,31]. However, when consistently com- bining with the reactor data, we find that the preference for NO by T2K+REA is reduced, and NOνA+REA actually favors IO. This is due to the slightly lower value of |∆m 2 3 | favored by the reactor data, in particular in comparison with NOνA for both orderings, and also with T2K for NO. Altogether we find that for the full combination of LBL accelerator experiments with reactors the "hint" towards NO is below 1σ.
• Figure 8 illustrates how both NOνA and MINOS favor non-maximal θ 23 . From this figure we see that while the significance of non-maximality in NOνA seems more evident than in MINOS when only the information of θ 13 is included (upper panels), the opposite holds for the full combination with the reactor data (lower panels). In  • Regarding the octant of θ 23 , for IO all LBL accelerator experiments are better described with θ 23 > 45 • , adding up to a ∼ 1.8σ preference for that octant. Conversely, for NO θ 23 < 45 • is favored at ∼ 1σ.
• From Fig. 9 we see that the "hint" for a CP phase around 270 • is mostly driven by T2K data, with some extra contribution from NOνA in the case of IO. Within the present precision the favored ranges of δ CP in each ordering by the combination of LBL accelerator experiments are pretty independent on the inclusion of the ∆m 2 3 information from reactors.

Analysis of Super-Kamiokande atmospheric data
In all the results discussed so far we have not included information from Super-Kamiokande atmospheric data. The reason is that our oscillation analysis cannot reproduce that of the collaboration presented in their talks in the last two years (see for example Ref. [66] for their latest unpublished results). Already since SK2 the Super-Kamiokande collaboration has been presenting its experimental results in terms of a growing number of data samples. The rates for some of those samples cannot be predicted (and therefore included in a statistical analysis) without a detailed simulation of the detector, which can only be made by the experimental collaboration itself. Our analysis of Super-Kamiokande data has been always based on the "classical" set of samples for which our simulations were reliable enough: sub-GeV and multi-GeV e-like and µ-like fully contained events, as well as partially contained, stopping and through-going muon data, each divided into 10 angular bins for a total of 70 energy and zenith angle bins (details on our simulation of the data samples and the statistical analysis  Table 2. Three-flavor oscillation parameters from our fit to global data, including also our reanalysis of SK1-4 (4581 days) atmospheric data. The numbers in the 1st (2nd) column are obtained assuming NO (IO), i.e., relative to the respective local minimum, whereas in the 3rd column we minimize also with respect to the ordering. The omitted parameters are identical to Tab. 1.
are given in the Appendix of Ref. [3]). Despite the limitations, until recently our results represented the most up-to-date analysis of the atmospheric neutrino data which could be performed outside the collaboration, and we were able to reproduce with reasonable precision the oscillation results of the full analysis presented by SK -both for what concerns the determination of the dominant parameters ∆m 2 3 and θ 23 , as well as their rather marginal sensitivity to the subdominant ν e appearance effects driven by θ 13 (and consequently to δ CP and the ordering). Thus we confidently included our own implementation of the Super-Kamiokande χ 2 in our global fit.
However, in the last two years Super-Kamiokande has developed a new analysis method in which a set of neural network based selections are introduced, some of them with the aim of constructing ν e +ν e enriched samples which are then further classified into ν e -like andν elike subsamples, thus increasing the sensitivity to subleading parameters such as the mass ordering and δ CP [65,67]. The selection criteria are constructed to exploit the expected differences in the number of charged pions and transverse momentum in the interaction of ν e versusν e . With this new analysis method Super-Kamiokande has been reporting in talks an increasing sensitivity to the ordering and to δ CP : for example, the preliminary results of the analysis of SK1-4 (including 2520 days of SK4) [66] in combination with the reactor constraint of θ 13 show a preference for NO with a ∆χ 2 (IO) = 4.3 and variation of χ 2 (δ CP ) with the CP phase at the level of ∼ 1.7σ.
Unfortunately, with publicly available information this analysis is not reproducible outside the collaboration. Conversely our "traditional" analysis based on their reproducible data samples continues to show only marginal dependence on these effects. This is illustrated in Fig. 10 and Tab. 2 where we show the impact of inclusion of our last re-analysis of SK atmospheric data using the above mentioned 70 bins in energy and zenith angle. 4 We only show the impact on the determination of sin 2 θ 23 , δ CP , and the mass ordering as the effect on all other parameters is negligible. We observe that ∆χ 2 for maximal mixing and the second θ 23 octant receive an additional contribution of about 1 unit in the case of NO, whereas the θ 23 result for IO is practically unchanged. Values of δ CP 90 • are slightly more disfavoured, whereas there is basically no effect on the mass ordering discrimination.
In summary, with the information at hand we are not able to reproduce the elements driving the main dependence on the subdominant effects of the official (though preliminary and unpublished) Super-Kamiokande results, while the dominant parameters are currently well determined by LBL experiments. For these reasons we have decided not to include our re-analysis of Super-Kamiokande data in our preferred global fit presented in the previous section. Needless to say that when enough quantitative information becomes available to allow a reliable simulation of the subdominant ν e -driven effects, we will proceed to include it in our global analysis.
4 Monte Carlo evaluation of confidence levels for θ 23 , δ CP and ordering At present the three least known neutrino oscillation parameters are the Dirac CP violating phase δ CP , the octant of θ 23 and the mass ordering (which in what follows we will denote by "O"). In order to study the information from data on these parameters one can use two ∆χ 2 test statistics [55,68]: where the minimization in the first equation is performed with respect to all oscillation parameters except δ CP and the ordering (x 1 = {θ 12 , θ 13 , θ 23 , ∆m 2 21 , |∆m 2 3 |}), while in the second equation the minimization is over all oscillation parameters except θ 23 and the ordering (x 2 = {θ 12 , θ 13 , δ CP , ∆m 2 21 , |∆m 2 3 |}). Here χ 2 min indicates the χ 2 minimum with respect to all oscillation parameters including the mass ordering.
We have plotted the values of these test statistics in the lower right and central left panels in Fig. 2. We can use them not only for the determination of δ CP and θ 23 , respectively, but also of the mass ordering. For instance, using Eq. (4.1) we can determine a confidence interval for δ CP at a given CL for both orderings. However, below a certain CL no interval will appear for the less favored ordering. In this sense we can exclude that ordering at the CL at which the corresponding interval for δ CP disappears. Note that a similar prescription to test the mass ordering can be built for any other parameter as well, e.g., for θ 23 using Eq. (4.2). 5 In Sec. 2 we have presented confidence intervals assuming that the test statistics follow a χ 2 -distribution with 1 dof, relying on Wilks theorem to hold [70] (this is what we call the Gaussian limit). However, the test statistics in Eqs. (4.1) and (4.2) are expected not to follow Wilks' theorem because of several reasons [68]: • Sensitivity of current data to δ CP is still limited, as can be seen in Fig. 2: all values of δ CP have ∆χ 2 < 14, and for NO not even ∆χ 2 = 6 is attained.
• Regarding θ 23 , its precision is dominated by ν µ disappearance experiments. Since the relevant survival probability depends dominantly on sin 2 2θ 23 , there is both a NuFIT 3.0 (2016) Figure 11. Allowed regions from the global data at 1σ, 90%, 2σ, 99% and 3σ CL (2 dof). We show projections onto different planes with δ CP on the vertical axis after minimizing with respect to all undisplayed parameters. The lower (upper) panels correspond to IO (NO). Contour regions are derived with respect to the global minimum which occurs for NO and is indicated by a star. The local minimum for IO is shown by a black dot.
physical boundary of their parameter space at θ 23 = 45 • (because sin 2θ 23 < 1), as well as a degeneracy related to the octant.
• The mass ordering is a discrete parameter.
• The dependence of the theoretical predictions on δ CP is significantly non-linear, even more considering the periodic nature of this parameter. Furthermore, there are complicated correlations and degeneracies between δ CP , θ 23 , and the mass ordering (see Fig. 11 for illustration).
Therefore, one may expect deviations from the Gaussian limit of the ∆χ 2 distributions, and confidence levels for these parameters should be cross checked through a Monte Carlo simulation of the relevant experiments. We consider in the following the combination of the T2K, NOνA, MINOS and Daya-Bay experiments, which are most relevant for the parameters we are interested in this section. For a given point of assumed true values for the parameters we generate a large number (10 4 ) of pseudo-data samples for each of the experiments. For each pseudo-data sample we compute the two statistics given in Eqs. (4.1) and (4.2) to determine their distributions numerically. In Ref. [68] it has been shown that the distribution of test statistics for 2-dimensional parameter region (such as for instance the middle panels of Fig. 11) are more close to Gaussianity than 1-dimensional ones such as Eqs. (4.1) and (4.2). Therefore we focus here on the 1-dimensional cases. First, let us note that in order to keep calculation time manageable one can fix all parameters which are known to be uncorrelated with the three we are interested in (i.e., θ 23 , δ CP , O). This is certainly the case for ∆m 2 21 and θ 12 which are determined independently by solar and KamLAND data. As for θ 13 , presently the most precise information arises from reactor data whose results are insensitive to δ CP and θ 23 . Consequently, marginalizing over θ 13 within reactor uncertainties or fixing it to the best fit value gives a negligible difference in the simulations. Concerning |∆m 2 3 | we observe that there are no strong correlations or degeneracies with δ CP (see Fig. 11), and we assume that the distributions of the test statistics do not significantly depend on the assumed true value. Therefore we consider only the global best fit values for each ordering as true values for |∆m 2 3 | to generate pseudo-data. However, since the relevant observables do depend non-trivially on its value, it is important to keep |∆m 2 3 | as a free parameter in the fit and to minimize the χ 2 for each pseudo-data sample with respect to it. Hence, we approximate the test statistics in Eqs. (4.1) and (4.2) by using with the other oscillation parameters kept fixed at their best fit points: ∆m 2 21 = 7.5 × 10 −5 eV 2 , sin 2 θ 12 = 0.31, and sin 2 θ 13 = 0.022.

δ CP and the mass ordering
The value of the test statistics (4.1) is shown in Fig. 12 for the combination of T2K, NOνA, MINOS and Daya-Bay as a function of δ CP for both mass orderings. In the generation of the pseudo-data we have assumed three representative values of θ 23,true as shown in the plots. The broken curves show, for each set of true values, the values of ∆χ 2 (δ CP , O) which are larger than 68%, 95%, and 99% of all generated data samples.
From the figure we read that if the ∆χ 2 from real data (solid curve, identical in the three panels) for a given ordering is above the x% CL lines for that ordering for a given value of δ CP , that value of δ CP and the mass ordering can be rejected with x% confidence. So if the minimum of the ∆χ 2 curve for one of the orderings (in this case IO is the one with non-zero minimum) is above the x% CL line one infers that that ordering is rejected at that CL.
For the sake of comparison we also show in Fig. 12 the corresponding 68%, 95% and 99% Gaussian confidence levels as horizontal lines. There are some qualitative deviations from Gaussianity that have already been reported [68]: • For θ 23 < 45 • , δ CP = 90 • , and IO as well as for θ 23 > 45 • , δ CP = 270 • and NO, the confidence levels decrease. This effect arises because at those points in parameter space the ν µ → ν e oscillation probability has a minimum or a maximum, respectively. Therefore, statistical fluctuations leading to less (or more) events than predicted cannot be accommodated by adjusting the parameters. ∆χ 2 is small more often and the confidence levels decrease. This is an effect always present at boundaries in parameter space, usually referred to as an effective decrease in the number of degrees of freedom in the model.
• Conversely for δ CP ∼ 90 • for θ 23 > 45 • , and δ CP ∼ 270 • for θ 23 < 45 • , the confidence levels increase. This is associated with the prominent presence of the octant degeneracy. Degeneracies imply that statistical fluctuations can drive you away from the true value, ∆χ 2 increases, and the confidence levels increase. This is usually referred to as an effective increase in the number of degrees of freedom in the model due to degeneracies.
• Overall we find that with present data confidence levels are clearly closer to Gaussianity than found in Refs. [9,68], where similar simulations have been performed with less data available. For those data sets confidence levels were consistently below their Gaussian limit. This was mainly a consequence of the limited statistics and the cyclic nature of δ CP which lead to an effective decrease in the number of degrees of  freedom. We now find that when the full combination of data currently available is included this effect is reduced, as expected if experiments become more sensitive.
• For all true values considered, IO is not rejected even at 1σ. In particular we find IO disfavored at 30% − 40% for sin 2 θ 23 = 0.44 − 0.60.
Quantitatively we show in Tab. 3 the CL at which CP conservation (δ CP = 0, 180 • ) is disfavored as well as the 90% and 95% confidence intervals for δ CP . We find that the CL of rejection of CP conservation as well as the allowed ranges do not depend very significantly on θ 23,true . This can be understood from Fig. 12: the dependence on θ 23,true occur mostly for δ CP ∼ 90 • and IO, a region discarded with a large CL, and for δ CP ∼ 270 • and NO, a region around the best fit.
Note that in the table the intervals for δ CP are defined for both orderings with respect to the global minimum (which happens for NO). Hence the intervals for IO include the effect that IO is slightly disfavored with respect to NO. They cannot be directly compared to the intervals given in Tab. 1, where we defined intervals relative to the local best fit point for each ordering.
A similar comment applies also to the CL quoted in the table to reject CP conservation. For IO this is defined relative to the best fit point in NO. We find that for NO, CP conservation is allowed at 70% CL, i.e., slightly above 1σ (with some deviations from the Gaussian result of 80% CL), while for IO the CL for CP conservation is above 2σ. Note that values of δ CP 90 • are disfavored at around 99% CL for NO, while for IO the rejection is at even higher CL: the ∆χ 2 with respect to the global minimum is around 14, which would correspond to 3.7σ in the Gaussian limit. Our Monte Carlo sample of 10 4 pseudo-data sets is not large enough to confirm such a high confidence level.  Figure 13. 68%, 95% and 99% confidence levels (broken curves) for the test statistics (4.2) along with its value (solid curves) for the combination of T2K, NOνA, MINOS and reactor data. The value of δ CP above each plot corresponds to the assumed true value chosen to generate the pseudoexperiments and for all panels we take ∆m 2 3 ,true = −2.53 × 10 −3 eV 2 for IO and +2.54 × 10 −3 eV 2 for NO. The solid horizontal lines represent the 68%, 95% and 99% CL predictions from Wilks' theorem.

θ 23 and the mass ordering
Moving now to the discussion of θ 23 , we show the value of the test statistics (4.2) in Fig. 13 for the combination of T2K, NOνA, MINOS and Daya-Bay experiments as a function of θ 23 , for both mass orderings. For the generation of the pseudo-data we have assumed three example values δ CP,true = 0, 180 • , 270 • . We do not show results for δ CP,true = 90 • , since this value is already quite disfavored by data, especially for IO. 6 The broken curves show for each set of true values, the values of ∆χ 2 (θ 23 , O) which are larger than 68%, 95%, and 99% of all generated data samples. From the figure we see that the deviations from Gaussianity are not very prominent and can be understood as follows: • The confidence levels decrease around maximal mixing because of the boundary on the parameter space present at maximal mixing for disappearance data.
• There is some increase and decrease in the confidence levels for δ CP = 270 • , in the same parameter region as the corresponding ones in Fig. 12.
In Tab  for sin 2 θ 23 for both orderings with respect to the global best fit. We observe from the table that the Gaussian approximation is quite good for both, the CL of maximal mixing as well as for the confidence intervals. We conclude that present data excludes maximal mixing at slightly more than 90% CL. Again we note that the intervals for sin 2 θ 23 for IO cannot be directly compared with the ones from Tab. 1, where they are defined with respect to the local minimum in each ordering. In Tab. 5 we show the CL at which a certain combination of mass ordering and θ 23 octant can be excluded with respect to the global minimum in the NO and 1st θ 23 octant. We observe that the CL of the second octant for NO shows relatively large deviations from Gaussianity and dependence on the true value of δ CP . In any case, the sensitivity is very low and the 2nd octant can be reject at most at 70% CL (1σ) for all values of δ CP . The first octant for IO can be excluded at between 83% and 91% CL, depending on δ CP . As discussed above, the exclusion of the IO/2nd octant case corresponds also to the exclusion of the IO, since at that point the confidence interval in IO would vanish. Also in this case we observe deviations from the Gaussian approximation and the CL of at best 32% is clearly less than 1σ (consistent with the results discussed in the previous subsection), showing that the considered data set has essentially no sensitivity to the mass ordering.

Conclusions
We have presented the results of the updated (as of fall 2016) analysis of relevant neutrino data in the framework of mixing among three massive neutrinos. Quantitatively the present determination of the two mass differences, three mixing angles and the relevant CP violating phase obtained under the assumption that their log-likelihood follows a χ 2 distribution is listed in Table 1, and the corresponding leptonic mixing matrix is given in Eq. (2.1). We have found that the maximum allowed CP violation in the leptonic sector parametrized by the Jarlskog determinant is J max CP = 0.0329 ± 0.0007 ( +0.0021 −0.0024 )) at 1σ (3σ). We have studied in detail how the sensitivity to the least-determined parameters θ 23 , δ CP and the mass ordering depends on the proper combination of the different data samples (Sec. 3.2). Furthermore we have quantified deviations from the Gaussian approximation in the evaluation of the confidence intervals for θ 23 and δ CP by performing a Monte Carlo study of the long baseline accelerator and reactor results (Sec. 4). We can summarize the main conclusions in these sections as follows: • At present the precision on the determination of |∆m 2 3 | from ν µ disappearance in LBL accelerator experiments NOνA, T2K and MINOS is comparable to that from ν e disappearance in reactor experiments, in particular with the spectral information from Daya-Bay. When comparing the region for each LBL experiment with that of the reactor experiments we find some dispersion in the best fit values and allowed ranges.
• The interpretation of the data from accelerator LBL experiments in the framework of 3ν mixing requires using information from the reactor experiments, in particular about the mixing angle θ 13 . But since, as mentioned above, reactor data also constrain |∆m 2 3 |, the resulting CL of presently low confidence effects (in particular the non-maximality of θ 23 and the mass ordering) is affected by the inclusion of this information in the combination.
• We find that the mass ordering favored by NOνA changes from NO to IO when the information on ∆m 2 3 from reactor experiments is correctly included in the LBL+REA combination, and the ∆χ 2 of NO in T2K is reduced from around 2 to 0.5 (see Fig. 7). Our MC study of the combination of LBL and reactor data shows that for all cases generated, NO is favored but with a CL of less than 1σ.
• About the non-maximality of θ 23 , we find that when the information on ∆m 2 3 from reactor experiments is correctly included in the LBL+REA combination, it is not NOνA but actually MINOS which contributes most to the preference for non-maximal θ 23 (see Fig. 8). Quantitatively our MC study of the combination of LBL and reactor data shows that for all the cases generated the CL for rejection of maximal θ 23 is about 92% for NO. As seen in Fig. 13 and Tab. 4, the CL of maximal mixing as well as confidence intervals for sin 2 θ 23 derived with MC simulations are not very different from the corresponding Gaussian approximation.
The CL for rejection of the disfavored octant depends on the true value of δ CP assumed in the MC study and it is generically lower than the one obtained in the Gaussian limit (see Tab. 5). For example, for NO the second octant is disfavored at a confidence level between 0.9σ and 1.3σ depending on the assumed true value of δ CP .
• The present sensitivity to δ CP is driven by T2K with a minor contribution from NOνA for IO (see Fig. 9). The dependence of the combined CL of the "hint" towards leptonic CP violation and in particular for δ CP 270 • on the true value of θ 23 is shown in Fig. 12, from which we read that for all cases generated CP conservation is disfavored only at 70% (1.05σ) for NO. Values of δ CP 90 • are disfavored at around 99% CL for NO, while for IO the rejection is at higher CL (∆χ 2 14 with respect to the global minimum).
Finally we comment that the increased statistics in SK4 and Borexino has had no major impact in the long-standing tension between the best fit values of ∆m 2 21 as determined from the analysis of KamLAND and solar data, which remains an unresolved ∼ 2σ effect.
Future updates of this analysis will be provided at the NuFIT website quoted in Ref. [49].