Updated fit to three neutrino mixing: status of leptonic CP violation

We present a global analysis of solar, atmospheric, reactor and accelerator neutrino data in the framework of three-neutrino oscillations based on data available in summer 2014. We provide the allowed ranges of the six oscillation parameters and show that their determination is stable with respect to uncertainties related to reactor neutrino and solar neutrino flux predictions. We find that the maximal possible value of the Jarlskog invariant in the lepton sector is $0.0329 \pm 0.0009$ ($\pm 0.0027$) at the $1\sigma$ ($3\sigma$) level and we use leptonic unitarity triangles to illustrate the ability of global oscillation data to obtain information on CP violation. We discuss"tendencies and tensions"of the global fit related to the octant of $\theta_{23}$ as well as the CP violating phase $\delta_\mathrm{CP}$. The favored values of $\delta_\mathrm{CP}$ are around $3\pi/2$ while values around $\pi/2$ are disfavored at about $\Delta\chi^2 \simeq 6$. We comment on the non-trivial task to assign a confidence level to this $\Delta\chi^2$ value by performing a Monte Carlo study of T2K data.


Introduction
Thanks to remarkable discoveries by a number of neutrino oscillation experiments it is now an established fact that neutrinos have mass and leptonic flavors are not symmetries of Nature [1,2], see Ref. [3] for an overview. Ignoring controversial indications for the existence of neutrino mass states at the eV scale (see Ref. [4] and references therein) a consistent description of global data on neutrino oscillations is possible by assuming mixing among the three known 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 [7]: where c ij ≡ cos θ ij and s ij ≡ sin θ ij . In addition to the Dirac-type phase δ CP , analogous to that of the quark sector, there may also be two physical phases associated to a possible Majorana character of neutrinos, which however are not relevant for neutrino oscillations [8,9] and are therefore omitted in the present work. Given the observed hierarchy between the solar and atmospheric mass-squared splittings there are two possible non-equivalent orderings for the mass eigenvalues, which are conventionally chosen as with ∆m 2 ij ≡ m 2 i − m 2 j . As it is customary we refer to the first option, Eq. (1.2), as Normal Ordering (NO), and to the second one, Eq. (1.3), as Inverted Ordering (IO); in this form they correspond to the two possible choices of the sign of ∆m 2 31 . In this convention the angles θ ij can be taken without loss of generality to lie in the first quadrant, θ ij ∈ [0, π/2], and the CP phase δ CP ∈ [0, 2π]. In the following we adopt the (arbitrary) convention of reporting results for ∆m 2 31 for NO and ∆m 2 32 for IO, i.e., we always use the one which has the larger absolute value. Sometimes we will generically denote such quantity as ∆m 2 3 , with = 1 for NO and = 2 for IO.
In this article, we present an up-to-date (as of summer 2014) global analysis of solar, atmospheric, reactor and accelerator neutrino data in the framework of three-neutrino oscillations. Alternative recent global fits have been presented in Refs. [10,11]. In Sec. 2 we describe the data used in our analysis (listed also in Appendix A) and we present the results of the global analysis and the allowed ranges of the oscillation parameters. In Sec. 3 we focus on our knowledge on CP violation, discussing the present status of the leptonic Jarlskog invariant and displaying the results of our fit in terms of leptonic unitarity triangles. In Sec. 4 we comment on various "tensions and tendencies" in the global data, including the reactor anomaly, the tension in the ∆m 2 21 determination from solar experiments versus KamLAND, the determination of ∆m 2 31 , tendencies in fit results for θ 23 and δ CP , and statistical issues related to the determination of the CP violating phase δ CP . Finally in Sec. 5 we present our conclusions.
The numerical results of our analysis as well as figures are available at the website [12], where also one-and two-dimensional χ 2 tables are available for download. Furthermore, this website will be kept up-to-date when new data becomes available.
2 Oscillation parameters: results of the global analysis

Data included in our analysis
We include in our global analysis the results from Super-Kamiokande atmospheric neutrino data from phases SK1-4 [13] (with addition of the 1775 days of phase SK4 over their published results on phases SK1-3 [14]), divided into 70 energy and zenith angle bins. For what concerns disappearance results from long baseline accelerator experiments (LBL) we use the energy distribution of events from MINOS in both ν µ (ν µ ) disappearance with 10.71 (3.36) × 10 20 protons on target (pot) [15], which amounts to 39 (14) data points, and from T2K in ν µ disappearance [16] with 6.57 × 10 20 pot (16 data points). For LBL appearance results we include both the neutrino and antineutrino events from MINOS [17], with exposure 10.6 × 10 20 and 3.3 × 10 20 pot, respectively, and from T2K in ν e appearance [18] with 6.57 × 10 20 pot; each of these samples contributes 5 data points.
In the analysis of solar neutrino experiments we include the total rates from the radiochemical experiments Chlorine [19], Gallex/GNO [20] and SAGE [21]. For real-time experiments we include the results from on electron scattering (ES) from the four phases in Super-Kamiokande: the 44 data points of the phase I (SK1) energy-zenith spectrum [22], the 33 (42) data points of the full energy and day/night spectrum in phase II (III), SK2 [23] (SK3 [24]), and the 24 data points of the energy spectrum and day-night asymmetry of the 1669-day of phase IV, SK4 [25]. The results of the three phases of SNO are included in terms of the parametrization given in their combined analysis [26] which amount to 7 data points. We also include the main set of the 740.7 days of Borexino data [27] as well as their high-energy spectrum from 246 live days [28]. In the analysis of solar neutrino data we use the GS98 version of the solar standard model [29] (see Sec. 4.2).
For oscillation signals at reactor experiments we include data from the finalized experiments CHOOZ [30] (energy spectrum data, 14 data points) and Palo Verde [31] (total rate) together with the spectrum from Double Chooz with 227.9 days live time [32] (18 data points), and the 621-day spectrum from Daya Bay [33] (36 data points), as well as the near and far rates observed at RENO with 800 days of data-taking [34] (2 data points with free normalization). We also include the observed energy spectrum in KamLAND data sets DS-1 and DS-2 [35] with a total exposure of 3.49 × 10 32 target-proton-year (2135 days). Although reactor experiments with baselines 100 m do not contribute to oscillation physics, they play an important role in constraining the unoscillated reactor neutrino flux. For this purpose we consider also data from Bugey4 [36], ROVNO4 [37], Bugey3 [38], Krasnoyarsk [39,40], ILL [41], Gösgen [42], SRP [43], and ROVNO88 [44], to which we refer as reactor short-baseline experiments (RSBL). Details on the RSBL analysis can be found in [4].
For convenience a detailed list of all the data used in our global analysis can also be found in Appendix A.

Description of the results
The results of the global analysis are presented in Figs. 1 and 2 where we show different projections of the allowed six-dimensional parameter space. To account for the possible effect of the so-called reactor anomaly [45][46][47], we follow the approach of Refs. [48,49] and study the dependence of the determined value of the parameters on the assumptions about the reactor fluxes. To bracket the possible impact of the anomaly, the results in Figs. 1 and 2 are shown for two extreme choices. The first option is to leave the normalization of reactor fluxes free and include data from short-baseline (less than 100 m) reactor experiments. This corresponds to the colored regions in Fig. 1 and the solid curves in Fig. 2 (labeled "Free+RSBL"). The second option is not to include short-baseline reactor data but assume reactor fluxes as predicted in [45] (including their uncertainties). This corresponds to the black contours in Fig. 1 Figure 1. Global 3ν oscillation analysis. Each panel shows a two-dimensional projection of the allowed six-dimensional region after minimization with respect to the undisplayed parameters. The different contours correspond to 1σ, 90%, 2σ, 99% and 3σ CL (2 dof). Full regions correspond to the analysis with free normalization of reactor fluxes and data from short-baseline (less than 100 m) reactor experiments included. For void regions short-baseline reactor data are not included but reactor fluxes as predicted in [45] are assumed. Note that as atmospheric mass-squared splitting we use ∆m 2 31 for NO and ∆m 2 32 for IO. The regions in the lower 4 panels are based on a ∆χ 2 minimized with respect to NO and IO. For solid curves the normalization of reactor fluxes is left free and data from short-baseline (less than 100 m) reactor experiments are included. For dashed curves short-baseline data are not included but reactor fluxes as predicted in [45] are assumed. Note that as atmospheric mass-squared splitting we use ∆m 2 31 for NO and ∆m 2 32 for IO. In what follows we will consider our default analysis choice the one with "Free Fluxes + RSBL". It is for this choice of fluxes that 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 the ranges are obtained after marginalizing with respect to the other parameters. We show 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 Inverted Ordering) and are obtained marginalizing also over the ordering. For this third case we only give the 3σ ranges. Of course in this case the range of ∆m 2 3 is composed of two disconnected intervals, one one containing the absolute minimum (IO) and the other the secondary local minimum (NO).

Mixing matrix and leptonic CP violation
From the global χ 2 analysis described in the previous section and following the procedure outlined in Ref. [50] one can derive the 3σ ranges on the magnitude of the elements of the leptonic mixing matrix to be: By construction the derived limits in Eq. (3.1) are obtained under the assumption of the matrix U being unitary. In other words, the ranges in the different entries of the matrix are correlated due to the constraints imposed by unitarity, as well as the fact that, in general, the result of a given experiment restricts a combination of several entries of the matrix. As a consequence choosing a specific value for one element further restricts the range of the others. The present status of the determination of leptonic CP violation is illustrated in Fig. 3 where we show the dependence of the ∆χ 2 of the global analysis on the Jarlskog invariant which gives a convention-independent measure of CP violation [51], defined as usual by: Using the parametrization in Eq. (1.1) we get J max CP = cos θ 12 sin θ 12 cos θ 23 sin θ 23 cos 2 θ 13 sin θ 13 . at 1σ (3σ) for both orderings. The preference of the present data for non-zero δ CP implies a best fit J best CP = −0.033, which is favored over CP conservation at the ∼ 1.3σ level. These numbers can be compared with the size of the Jarlskog invariant in the quark sector, which is determined to be J quarks . In Fig. 4 we recast the allowed regions for the leptonic mixing matrix in terms of leptonic unitarity triangles, which are obtained as different combinations of the entries of the U matrix. 1 Since in our analysis U is unitary by construction, any given pair of rows or columns can be used to define a triangle in the complex plane. On the left (right) panels we show the triangles corresponding to the unitarity conditions In drawing these triangles we have rescaled and rotated their sides so that two of their vertices always coincide with (0, 0) and (1, 0) in the complex plane. To this aim we have defined a complex variable z as follows: and then we have plot the 1σ, 90%, 2σ, 99%, 3σ CL (2 dof) allowed regions of the third vertex of the triangle as the real and imaginary parts of z. For convenience in each panel NuFIT 2.0 (2014) Figure 4. Six leptonic unitarity triangles. After scaling and rotating each triangle so that two of its vertices always coincide with (0, 0) and (1, 0) (see text for details) we plot the 1σ, 90%, 2σ, 99%, 3σ CL (2 dof) allowed regions of the third vertex. Note that in the construction of the triangles the unitarity of the U matrix is always explicitly imposed.
we have chosen the normalization side (the one which lies on the horizontal (0, 0) → (0, 1) segment) as the best determined of the two longer sides of each triangle. In this way all the triangles have more or less the same size, and the uncertainty in the position of the third vertex is not too much affected by the uncertainty of the normalization side. Note that the most common unitarity triangle in the quark sector is the one based on the d-quark and b-quark columns [7], which corresponds to the 1st and 3rd column in the leptonic matrix, i.e., to the triangle in the middle-right panel in Fig. 4. In this kind of diagrams the absence of CP violation implies a flat triangle, i.e., Im(z) = 0. As can be seen, in all the panels the horizontal axis marginally 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.

Impact of reactor flux uncertainties
Within the 3-flavor framework the so-called reactor anomaly leads to a "tension" of about 2.7σ between the predicted reactor neutrino fluxes [45,46] and the event rates observed in short-baseline reactor experiments. By adopting two extreme approaches in dealing with this tension we have shown in Sec. 2.2 that the impact on the determination of the oscillation parameters in the global fit is quite small, at the level of 0.5σ for sin 2 θ 13 (see Figs. 1 and 2). This is further illustrated in Fig. 5 where we show the allowed regions in the plane of θ 13 and the flux normalization f flux (relative to the one predicted in [45]) for several combinations of the reactor experiments. Short-baseline data (green contours) essentially determine the flux normalization. Adding also data from experiments at around 1 km without a dedicated near detector (red contours) provides already a signal for non-zero θ 13 , but such result is affected by significant correlation with the flux normalization. However, once the precise data on near-far comparison from Daya Bay and RENO are included (colored regions) no correlation is left between the determination of θ 13 and f flux . Thus in the 3ν analysis the unexplained reactor anomaly mostly translates in an overall increase of the χ 2 in the analysis with fluxes from Ref. [45] with χ 2 (f flux = 1) − χ 2 (f flux free) 7. Details of our analysis in this respect can be found in Ref. [4], where a discussion of a possible explanation in terms of sterile neutrinos is also given. 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 (see text for details). Right: ∆χ 2 dependence on ∆m 2 21 for the same four analysis after marginalizing over θ 12 .
4.2 Determination of ∆m 2 21 : solar and KamLAND We show in Fig. 6 the results of the analysis of the solar experiments and of KamLAND which give the dominant contribution to the determination of ∆m 2 21 and θ 12 . Here θ 13 is fixed to the present best fit value of the global analysis. For the sake of completeness the solar neutrino results are shown for two different versions of the Standard Solar Model, namely the GS98 and the AGSS09 models [29]. Let us remind that GS98 is based on the older solar abundances leading to high metallicity and which perfectly agreed with helioseismological data, whereas AGSS09 uses the new precise determination of the solar abundances which imply a lower metallicity and cannot reproduce the helioseismological data. This conflict constitutes the so-called "solar composition problem". Although it is a pretty serious problem in the context of solar physics, its impact in the determination of the relevant oscillation parameters is very small, as can be seen clearly from Fig. 6.
The left panel in Fig. 6 illustrates the complementarity of solar and KamLAND in the determination of the "12" parameters. Solar experiments provide the best precision of θ 12 while KamLAND gives a better determination of ∆m 2 21 . We remind the reader that the relevant survival probabilities for these experiments in the framework of three neutrino oscillations can be written as: P 3ν ee = sin 4 θ 13 + cos 4 θ 13 P 2ν ee (∆m 2 21 , θ 12 ) , where we have used the fact that L osc 31 = 4πE ν /∆m 2 31 is much shorter than the distance traveled by both solar and KamLAND neutrinos, so that the oscillations related to L osc 31 are averaged. In presence of matter effects P 2ν ee (∆m 2 21 , θ 12 ) should be calculated taking into account the evolution in an effective matter density n eff e = n e cos 2 θ 13 . For 10 −5 ∆m 2 /eV 2 10 −4 , P 2ν ee (∆m 2 21 , θ 12 ) presents the following asymptotic behaviors [55]: At present most of the precision of the solar analysis is provided by SNO and SK for which the relevant MSW survival probability [56,57] provides a direct measurement of sin 2 θ 12 , as seen in Eq. (4.3). In the MSW regime the determination of ∆m 2 21 in solar experiments comes dominantly from the ratio between the solar potential and the ∆m 2 21 term required to simultaneously describe the CC/NC data at SNO and the undistorted spectra of 8 B neutrinos as measured in both SK and SNO. Conversely KamLANDν e survival probability proceeds dominantly as vacuum oscillations and provides a most precise determination of ∆m 2 21 via the strong effect of the oscillating phase in the distortion of the reactor energy spectrum. On the contrary it yields a weaker constraint on θ 12 as the vacuum oscillation probability depends on the double-valued and "flatter" function sin 2 (2θ 12 ).
As seen in the left panel in Fig. 6 for either version of the solar model the best fit points of solar and KamLAND analysis lie at very similar values of θ 12 . As it was pointed out in Ref. [58] and widely discussed in the literature [59][60][61][62][63], the matching in the determination of θ 12 requires the presence of a non-zero value of θ 13 . With the present determination of θ 13 provided by the medium baseline reactor experiments, the agreement between the best fit point values of θ 12 is remarkable.
From the same figure, however, we see that the value of ∆m 2 21 preferred by KamLAND is higher than the one from solar experiments. At present this is about a 2σ effect, as can be seen in the right panel where we show the ∆χ 2 dependence as a function of ∆m 2 21 when marginalized over θ 12 . This tension has been present during the last two years and it arises from a combination of two effects: (a) the well-known fact that none of the 8 B measurement performed by SNO, SK and Borexino show any evidence of the spectrum low energy turnup expected in the standard LMA-MSW solution, and (b) the indication of a non-vanishing day-night asymmetry in SK, which disfavors the KamLAND ∆m 2 21 best fit value for which Earth matter effects are too small. The relevance of these effects is illustrated in Fig. 6 where we show the results of our analysis both with and without the inclusion of the SK day-night information. As can be seen, once the SK day-night information is removed the solar best-fit point shifts upwards and the solar allowed region extends to much larger values of ∆m 2 21 , as expected, so that the tension with KamLAND is reduced to about 1.4σ. Modified matter potential due to non-standard interactions [64,65] and super-light sterile neutrinos [66] have been proposed as extended scenarios which could relax this tension. 3 ) plane using both appearance and disappearance data from MINOS (green) and T2K (black), as well as SK atmospheric data (green) and a combination of them (colored regions). Here θ 13 is constrained to the 3σ range from the global fit. The right panels show regions in the (sin 2 θ 13 , ∆m 2 3 ) plane using data from Daya Bay (black), reactor data without Daya Bay (violet), 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 local minimum in each panel.

Determination of ∆m 2
3 : ν µ and ν e disappearance Fig. 7 illustrates the determination of ∆m 2 3 from different data sets. In the left panels we focus on long-baseline ν µ disappearance data. It is clear that in this case the final precision on |∆m 2 3 | emerges from the combination of T2K and MINOS data, while the determination of sin 2 θ 23 is dominated by T2K.
Concerning ν e disappearance data, Eq. (4.8) in Sec. 4.4 implies that the rates observed in reactor experiments at different baselines can provide an independent determination of ∆m 2 3 [49,67]. On top of this, the observation of the energy-dependent oscillation effect of θ 13 in Daya Bay [68] allows a rather precise determination of |∆m 2 3 |. In the right panels of Fig. 7 we show therefore the allowed regions in the (θ 13 , ∆m 2 3 ) plane based on global data on ν e disappearance. The blue contours are obtained from all the mediumbaselines reactor experiments with the exception of Daya Bay. Those regions emerge from the baseline effect mentioned above. The black contour are based on the energy spectrum in Daya Bay, whereas the colored regions show the combination.
By comparing the left and right panels we observe that ν µ and ν e disappearance experiments by now provide a consistent determination of |∆m 2 3 | with similar precision.

Mass ordering, θ 23 octant and CP phase: role of different data sets
As we have seen in Sec. 2, several 1-2σ "tendencies" appear in the global analysis in the determination of the mass ordering, the octant of θ 23 , and the CP violating phase. To illustrate the role of the different data sets on such tendencies, we show in Fig. 8 the ∆χ 2 as a function of ∆m 2 3 , θ 23 , and δ CP for different combinations of experiments. In each panel the results have been marginalized with respect to all undisplayed parameters except the mass ordering, which is fixed to Inverted (Normal) for the left (right) panels. Note, however, that for each combination of experiments the ∆χ 2 is defined with respect to the absolute minima between the two orderings. In this way the difference between the "height" of the minimum of the curve on the left and the corresponding one on the right gives the contribution of that set of observables to the determination of the mass ordering.
All the lines plotted in Fig. 8 include "by default" solar and reactor data, which take care of precisely determining the undisplayed parameters ∆m 2 21 , θ 12 and θ 13 . To this basic set we progressively add more and more data, to see how each new piece of information affects the results of the fit. Let us then start with the purple (and blue) curve, which shows the dependence of ∆χ 2 on the analysis of solar, reactor and MINOS (plus T2K) ν µ andν µ disappearance data. Being all disappearance experiments they provide very weak information on δ CP , as clearly visible in the bottom panels. Comparing the minima in the left and right panels we note a relative difference of χ 2 (NO) − χ 2 (IO) ∼ 0.2, which means that this combination of data "favors" Inverted Ordering by ∼ 0.5σ. More interestingly, from the central panels we see that MINOS disappearance data favors a non-maximal θ 23 with ∆χ 2 (θ 23 = 45 • ) = 2.8 (2.2) for IO (NO). Neglecting subleading ∆m 2 21 and matter effects, the relevant survival probability in MINOS is given by where L is the baseline and E ν is the neutrino energy. Hence, the probability is symmetric under θ dis → π/2 − θ dis . In the limit θ 13 = 0 the effective angle θ dis reduces to θ 23 , and a preference for non-maximal θ dis mixing leads to the appearance of two symmetric minima in the first and second octant of θ 23 . Such degeneracy persists also for θ 13 = 0, and is responsible for the presence of two quasi-degenerate minima at sin 2 θ 23 = 0.63 and 0.39. On the other hand, T2K disappearance data are better fitted with maximal θ dis , so once they are included in the analysis (blue line) the positions of the two minima move to values sin 2 θ 23 = 0.58 and 0.44 while the preference for non-maximal mixing reduces to ∆χ 2 (θ 23 = 45 • ) = 1. The comparison of the purple and blue curves also shows the impact of the inclusion of T2K disappearance data on the overall determination of θ 23 and ∆m 2 3 . The green line shows the effect of further adding to the analysis the T2K ν e appearance data. First, we see that the absolute minima now occurs for NO with ∆χ 2 (IO) = 0.6. In the central panels we see that the quasi-degeneracy of the octant of θ 23 is now broken and the second octant becomes favored with ∆χ 2 (θ 23 ≤ 45 • ) = 2.5 (1.5) for IO (NO). The lower panels show that after the inclusion of T2K ν e appearance data a minimum appears for δ CP = 270 • (300 • ) for IO (NO) with CP conservation disfavored at ∆χ 2 (sin δ CP = 0) = 2.5 (1.0). This can be understood from the relevant ν e appearance probability at T2K and MINOS, which, at the second order in the small parameters sin θ 13 and α ≡ ∆m 2 21 /∆m 2 31 and assuming a constant matter density, takes the form [69][70][71]: Here L is the baseline, E ν is the neutrino energy, and V is the effective matter potential [56] which for T2K yields |A| ∼ few %. The first term in Eq. (4.6) (which dominates for large θ 13 ) depends on sin 2 θ 23 and therefore is sensitive to the octant. Reactor experiments with L ∼ 1 km, on the other hand, provide a measurement of θ 13 independent of θ 23 At present the ν e appearance results from T2K points towards an excess with respect to what is expected for the best fit value of sin 2 θ 13 determined by the reactor experiments for maximal θ 23 (i.e., for 2 sin 2 θ 23 = 1), hence the tendency towards the θ 23 > 45 • minimum. The matter effects in Eq. (4.6) make this tendency different for NO and IO, while the last term introduces a δ CP modulation of the effect. For fixed θ 13 and θ 23 , P νµ→νe (δ CP ) − P νµ→νe (π) ≥ 0 (≤ 0) for δ CP ≥ π (≤ π). For the best fit values of θ 13 and θ 23 from the previous reactor and LBL ν µ disappearance results, the T2K ν e appearance signal is better fitted with δ CP values which enhance the corresponding appearance probability.
Conversely we see that adding the less significant MINOS ν e appearance data in the analysis (red curves) tends to slightly reduce the size of these effects for NO and it shifts the global minimum from NO to IO. Finally the orange curves show the impact of including the atmospheric data in the analysis. Comparing the orange and red curves we see that atmospheric data contributes positively to the significance of the tendency towards IO and δ CP > π. While for IO it does not affect the tendency towards second θ 23 octant, for NO it "shifts" this tendency to the first octant. The preference for θ 23 < 45 • for NO is related to an excess of sub-GeV e-like events, an effect which has already been discussed since many years (see, e.g., [72][73][74][75]). The fact that this preference is not visible for IO is probably related to multi-GeV data, which are affected by matter effects and therefore provides some sensitivity to the mass ordering. Identifying the relevant bins is difficult, given the large amount of data points entering the atmospheric fit. We stress that such effects happen at the level of 1-2 units in χ 2 and hence are not statistically significant. In order to highlight the pattern of correlations between δ CP and sin 2 θ 23 we show in Fig. 9 the allowed regions of the global analysis projected into the plane of these two parameters. Correlations between δ CP and other oscillation parameters are mostly trivial and are therefore omitted.

Remarks on confidence levels for δ CP
In order to study the information from data on the CP phase we consider the quantity where the first term on the right hand side is minimized with respect to all oscillation parameters except δ CP (x = θ 12 , θ 13 , θ 23 , ∆m 2 21 , ∆m 2 31 ) and the last term is the χ 2 minimum with respect to all oscillation parameters. We have shown ∆χ 2 (δ CP ) for various data sets in the lower panels of Fig. 8, as well as in the corresponding panel in Fig. 2 for the global data. The standard way to derive confidence intervals for δ CP is to assume that ∆χ 2 (δ CP ) follows a χ 2 -distribution with 1 dof, and then apply cuts corresponding to, e.g., ∆χ 2 = 0.99, 2.71, 3.84, 6.63 for 68%, 90%, 95%, 99% CL, respectively. This procedure relies on Wilks theorem to hold [76]. However, in the case of δ CP some of the hypothesis of this theorem may be violated [77,78]. One reason for this is the complicated non-linear dependence of the event rates on δ CP . Present sensitivity is so poor, that those non-linearities (as well as the periodic character of δ CP ) become relevant already at very low CL. Furthermore, Hence, although our results represent the most up-to-date analysis of the atmospheric neutrino data which can be performed outside the collaboration, such an analysis has unavoidable limitations. For details on our simulation of the data samples and the statistical analysis see the Appendix of Ref. [3]. parameter degeneracies (especially with θ 23 , see discussion in the previous sub-section) affect the distribution of the test statistic ∆χ 2 from Eq. (4.9). In order to address such concerns we have performed a Monte Carlo study of T2K data (appearance and disappearance). We consider a test statistic similar to the ∆χ 2 given in Eq. (4.9); however, in order to keep calculation time manageable we fix all oscillation parameters except δ CP and θ 23 to their best fit values from the global fit assuming normal mass ordering. Hence, in the notation of Eq. (4.9) we have only x = θ 23 . In particular, since we keep also θ 13 fixed, the main feature of the complementarity of long-baseline appearance and medium-baseline reactor data is maintained. We have checked that allowing θ 13 to vary imposing the constraint from Daya Bay data has a negligible impact on ∆χ 2 (δ CP ) compared to fixing it to the best fit value. The resulting ∆χ 2 (δ CP ) is shown as blue curve in Fig. 10 (identical in all three panels). It differs somewhat from the global result displayed in Figs. 2 or 8, which include more data, but it captures the essential features and suffices for the purpose of studying the statistical properties of the test statistic.
In order to estimate the probability distribution for ∆χ 2 (δ CP ) we proceed as follows. We scan the parameter space of δ CP and sin 2 θ 23 and for a given point of assumed true values we generate a large number of pseudo-data samples for T2K. For each data set we calculate the value of the test statistic ∆χ 2 (δ CP ) and in this way we obtain a distribution for it. We scan 41 points in δ CP and 3 points for sin 2 θ 23 , and for each of those points we generate 5000 pseudo-data samples. The black curves in Fig. 10 show the values of ∆χ 2 (δ CP ) which are larger than 68%, 90%, 95% and 99% of all generated data samples. We observe quite large deviations from the corresponding values based on the χ 2 -distribution for 1 dof, shown by the dashed lines in the figure. Interestingly we find also a rather strong dependence on the assumed true value of θ 23 .
The behavior of the curves can be understood qualitatively. Due to the non-linearity of δ CP (its cyclic nature) and the poor sensitivity mentioned above it actually counts as less than 1 full degree of freedom, which implies distributions more concentrated at lower values than the χ 2 -distribution for 1 dof, as observed in Fig. 10. The rather strong variations for non-maximal values of θ 23 , including a flipped behavior for δ CP smaller or larger π between sin 2 θ 23 = 0.4 and 0.6 can be understood in terms of a degeneracy. For θ 23 < π/4 and δ CP ∼ 3π/2 as well as for for θ 23 > π/4 and δ CP ∼ π/2 there is a degeneracy between the two octants of θ 23 which effectively enhances the number of degrees of freedom in the fit. 3 Now we can compare ∆χ 2 (δ CP ) obtained from the observed data to the expected distribution. If the observed ∆χ 2 (δ CP ) is larger than the values obtained for x% of the pseudo-data samples for that true value of δ CP we exclude this value of δ CP at the x% CL. In Fig. 10 we show as an example the resulting 90% confidence interval for δ CP as brown shaded area. This corresponds to the confidence interval according to the Feldman-Cousins (FC) prescription [79]. It has to be compared to the corresponding interval based on the χ 2 -approximation, indicated by the gray area in the plot.
We can draw the following conclusions from the exercise shown in Fig. 10: 1. the confidence intervals based on the Monte Carlo simulation are smaller than the ones based on the χ 2 -approximation. Hence, the latter is conservative; 2. for confidence levels 90% the confidence intervals are similar, whereas for higher confidence levels differences become significant. In particular, at 99% CL all values of δ CP are allowed using the χ 2 -approximation, whereas a region around δ CP ∼ π/2 remains excluded by the 99% CL FC interval; 3. the CL with which δ CP ∼ π/2 can be disfavored depends strongly on the unknown true value of θ 23 . For sin 2 θ 23 = 0.6, δ CP π/2 is excluded at about 99% CL, whereas for sin 2 θ 23 = 0.4 it is excluded at very high CL. In all cases, the CL based on the Monte Carlo is higher than in the χ 2 -approximation which again can be considered conservative.
Let us conclude this section by commenting that ideally such a simulation should be performed also for the global analysis. Unfortunately this is currently out of question, in particular due to atmospheric neutrino data, which is very computational intensive and does play a non-negligible role in the global fit for ∆χ 2 (δ CP ), see Fig. 8. However, we believe that the above results based on T2K are approximately representative also for the global fit. One may expect that, with more statistics, distributions become more close to the expected χ 2 -distribution. However, preliminary estimates indicate that parameter degeneracies may lead to deviations also in a high-statistics scenario.

Summary
We have presented the results of an updated (as of summer 2014) global analysis of solar, atmospheric, reactor and accelerator neutrino data in the framework of three-neutrino oscillations. Quantitatively the present determination of the oscillation parameters is listed in Table 1, and the corresponding leptonic mixing matrix is given in Eq. (3.1). From the present analysis we have derived the maximum allowed CP violation in the leptonic sector as parametrized by the Jarlskog determinant, J CP = 0.033 ± 0.010 (± 0.027) at 1σ (3σ). All these results have also been shown in terms of unitarity triangles in Fig. 4 which further illustrate the ability of global oscillation data to obtain information on leptonic CP violation.
The global analysis presents a series of tensions between data sets as well as some 1-2σ effects in the determination of less known parameters (θ 23 , mass ordering, and δ CP ) which we denote as "tendencies" and we discuss in Sec. 4. We can summarize these results as follows: • due to the very precise determination of the flux-independent near-far ratio from Daya Bay and RENO, the so-called reactor neutrino anomaly (i.e., the tension between the predicted reactor fluxes in Refs. [45,46] and the event rates observed in short-baseline reactor experiments) results only in a 0.5σ uncertainty on the determination of θ 13 ; • the long-standing ∼ 2σ tension between the best fit values of ∆m 2 21 as determined from the analysis of KamLAND and solar data is still unresolved. This tension is driven by both the indication of a non-zero day-night effect at SK, and by the lack of evidence of a low energy turn-up in the 8 B energy spectrum as measured by SNO, SK4 and Borexino. In both cases the ∆m 2 21 value favored by KamLAND is in disagreement with the expectations from the standard LMA-MSW solution; • the uncertainty on the determination of ∆m 2 21 and θ 12 due to the choice of Standard Solar Model associated with the "solar composition problem" is negligible; • at present the precision on the determination of |∆m 2 3 | from ν µ disappearance experiments (mainly T2K and MINOS) is comparable to that from ν e disappearance experiments (i.e. reactor experiments including, in particular, the spectral information from Daya Bay); • for Inverted Ordering, the "tendency" towards non-maximal mixing and second octant of θ 23 is driven mainly by two effects: (a) the non-maximality favored by MINOS ν µ disappearance, and (b) the "mismatch" between the best fit θ 13 obtained fromν e disappearance at reactors and from ν µ → ν e at T2K. Atmospheric results do not alter this; • for Normal Ordering, such preference for non-maximal θ 23 mixing is considerably weaker than for IO; also, in this case the global best-fit occur in the first θ 23 octant, mostly driven by atmospheric data; • the "mismatch" between reactor and T2K results is the driving effect in the present dependence of the global ∆χ 2 on the CP violating phase with a best fit value close to δ CP = 3 2 π. Inclusion of the atmospheric results adds positively to this effect for both orderings; • the tendency towards IO or NO in the present analysis does not seem to result from any consistent effect and it shifts in sign depending on the data sets considered.
Finally in Sec. 4.5 we have addressed the issue of the "gaussianity" of the confidence levels attributed to ∆χ 2 (δ CP ) by performing a Monte Carlo study of T2K data, and we have compared the resulting probability distribution to that of a χ 2 -distribution as usually assumed. Deviations are expected due to the cyclic nature of δ CP and to the presence of parameter degeneracies. The conclusion is that, within the present data, the use of the χ 2 -distribution approximation is slightly conservative in the determination of the excluded range of δ CP at confidence levels 90%. The differences however are not very significant as illustrated in Fig. 10.
Future updates of this analysis will be provided at the website quoted in Ref. [12].