Discovery Potential of T2K and NOvA in the Presence of a Light Sterile Neutrino

We study the impact of one light sterile neutrino on the prospective data expected to come from the two presently running long-baseline experiments T2K and NOvA when they will accumulate their full planned exposure. Introducing for the first time, the bi-probability representation in the 4-flavor framework, commonly used in the 3-flavor scenario, we present a detailed discussion of the behavior of the numu to nue and numubar to nuebar transition probabilities in the 3+1 scheme. We also perform a detailed sensitivity study of these two experiments (both in the stand-alone and combined modes) to assess their discovery reach in the presence of a light sterile neutrino. For realistic benchmark values of the mass-mixing parameters (as inferred from the existing global short-baseline fits), we find that the performance of both these experiments in claiming the discovery of the CP-violation induced by the standard CP-phase delta13 equivalent to delta, and the neutrino mass hierarchy get substantially deteriorated. The exact loss of sensitivity depends on the value of the unknown CP-phase delta14. Finally, we estimate the discovery potential of total CP-violation (i.e., induced simultaneously by the two CP-phases delta13 and delta14), and the capability of the two experiments of reconstructing the true values of such CP-phases. The typical (1 sigma level) uncertainties on the reconstructed phases are approximately 40 degree for delta13 and 50 degree for delta14.


Introduction and Motivation
More than fifteen years ago, pioneering observations of neutrinos originating from natural sources (the sun core and the earth atmosphere) led to the first evidence of neutrino oscillations establishing the massive nature of these fundamental particles. Such a discovery, recently awarded with the Nobel Prize [1], has been the first of a long series of milestones in our understanding of neutrinos, whose properties have been gradually clarified by the subsequent findings of several other experiments performed with man-made neutrino sources (reactors and accelerators).
The two apparently disjoint effective 2-flavor descriptions initially introduced to explain separately the solar and the atmospheric neutrino anomalies have been gradually recognized as two pieces of a single and more complex mosaic, which is currently accepted as the standard picture of neutrino oscillations. Such a 3-flavor framework involves two distinct mass-squared splittings (∆m 2 31 , ∆m 2 21 ), three mixing angles (θ 12 , θ 23 , θ 13 ), and one CPphase δ.
The latest fundamental step in the establishment of the 3-flavor scheme has been accomplished very recently (during 2012), with the determination of the long-sought third mixing angle θ 13 by means of dedicated reactor experiments [2][3][4]. This discovery has opened the way to the measurement of the last unknown 3-flavor parameters: the CPphase δ and the neutrino mass hierarchy (MH), i.e. the sign of ∆m 2 31 . Both properties are at the center of an intense world-wide research program, which will be carried out with new accelerator, reactor and atmospheric neutrino experiments (for a recent review see [5]).
In spite of its tremendous success and of its beautiful structure the standard 3-flavor framework may not constitute the ultimate description of neutrinos, which may reserve surprises. In fact, a few anomalies have been recorded at the short baseline experiments, which cannot be accommodated in the 3-flavor scheme (see [6][7][8] for a review of the topic). The two standard mass-squared splittings ∆m 2 21 ≡ m 2 2 − m 2 1 and ∆m 2 31 ≡ m 2 3 − m 2 1 are too small to produce observable effects in such setups and (at least) one new much larger masssquared difference O(eV 2 ) must be introduced. The hypothetical fourth mass eigenstate must be essentially sterile. A rich and diverse program of new more sensitive short-baseline experiments is underway in order to test such an intriguing hypothesis (see the review in [9]), which, if confirmed would represent a revolution in our understanding of neutrinos, as important as the discovery of neutrino oscillations.
At a phenomenological level, the existence of sterile neutrinos would make necessary to extend the standard 3-flavor framework. The enlarged scheme must be realized in such a way to preserve the very good description of all the other (non short-baseline) data. In the minimal extension, involving only one sterile species, the so-called 3 + 1 scheme, the new mass eigenstate ν 4 is assumed to be weakly mixed with the active neutrino flavors (ν e , ν µ , ν τ ) and it is separated from the standard mass eigenstates (ν 1 , ν 2 , ν 3 ) by a large O(eV 2 ) splitting, giving rise to the hierarchical pattern |∆m 2 21 | |∆m 2 31 | |∆m 2 41 |. The 3+1 scheme involves six mixing angles and three (Dirac) CP-violating phases. Hence, in case of discovery of a sterile neutrino, we would face the formidable challenge of identifying six more properties (3 mixing angles, 2 CP-phases and the sign of ∆m 2 41 ) in addition to those involved in the standard 3-flavor framework.
The 3 + 1 scheme naturally predicts sizable effects at the short baselines, where the oscillating factor ∆ 41 ≡ ∆m 2 41 L/4E (L being the baseline and E the neutrino energy) is of order one, and one expects an oscillating behavior with the characteristic L/E dependency. However, it must be emphasized that sterile neutrinos are observable in other (non-shortbaseline) types of experiments where they may manifest in a more subtle way. In the solar sector, for example, a non zero-value of the electron neutrino mixing with ν 4 (parametrized by the matrix element U e4 ) can be felt as a small deviation of the unitarity of the (ν 1 , ν 2 , ν 3 ) sub-system [10,11] (see also [12]). In the atmospheric sector, as first evidenced in [13], at very high [O(TeV)] energies one expects a novel MSW resonance, which may reveal as a distortion of the zenith angle distribution. To this regard we mention the dedicated analysis underway in the IceCube collaboration whose preliminary results have been shown very recently in [14]. Complementary information on sterile neutrino mixing using atmospheric neutrinos has been also extracted from Super-Kamiokande [15].
Sterile neutrino oscillations can influence also the long-baseline (LBL) accelerator ex-periments 1 . We recall that these setups, when working in the ν µ → ν e (andν µ →ν e ) appearance channel, can probe the 3-flavor CP-violating phenomena encoded by the CPphase δ. Their sensitivity is related to the fact that at long distances the ν µ → ν e transition probability develops a small interference term (which is completely negligible at SBL) between the oscillations induced by the atmospheric splitting and those driven by the (much smaller) solar splitting. As first evidenced in [19], in the presence of sterile neutrinos a new interference term appears in the transition probability, which depends on one additional CP-phase. Notably, the size of the new non-standard 4-flavor interference term is expected to be similar to that of the standard 3-flavor interference term. Therefore, the sensitivity to the new interference term can disclose the possibility to explore the new enlarged CPviolation (CPV) sector, which involves two additional CP-phases. It should be stressed that both in the 3-flavor and 4-flavor schemes the CP-phases cannot be observed in SBL experiments, since in such setups the two standard oscillating frequencies have negligible values. Therefore, the LBL and SBL experiments are complementary in the exploration of the 3+1 scheme (and of any scheme involving more than one sterile neutrino).
This basic observation provides the motivation for the study performed in the present paper, in which we explore the physics potential of the two currently running LBL experiments T2K and NOνA in the presence of a hypothetical light sterile neutrino. 2 The analyses performed in [19,37] with the existing data from T2K and NOνA have already shown that these two experiments 3 can probe one of the new CP-phases for realistic values of the mixing angles indicated by the global 3+1 fits [42,43]. In addition, in [37] it has been pointed out that the statistical significance of the current indications concerning the standard CP-phase δ and the neutrino mass hierarchy is modified (reduced) in the presence of sterile neutrinos. Therefore, it is timely and interesting to investigate if such features are expected to persist even when the full exposure will be reached in both experiments. With this aim, we perform a prospective study addressing in a quantitative way these questions.
The paper is organized as follows. In section 2, we present a detailed discussion of the behavior of the 4-flavor ν µ → ν e andν µ →ν e transition probabilities. For the first time we extend the bi-probability representation, commonly used in the 3-flavor framework, to the more general 3+1 scheme. Section 3 deals with the experimental details and also discusses the bi-events plots. In section 4, we describe the details of the statistical method that we use for the analysis. Section 5 is devoted to the presentation of the results of the sensitivity 1 In this paper we focus on the νµ → νe channel. However, information on sterile neutrinos can be obtained also from the LBL disappearance νµ → νµ searches and from the neutral current data. See the analyses performed by MINOS [16,17]. Also the appearance νµ → ντ channel can provide information, albeit currently the low statistics limits the sensitivity. See the analysis performed by OPERA [18]. 2 Old works on sterile neutrinos at LBL can be found in [20][21][22][23][24][25][26][27][28]. More recent studies focusing on the future LBNE/DUNE experiment [29] have been recently performed in [30][31][32]. In principle, the CERN-Pyhäsalmi baseline of 2290 km actively studied under the umbrella of the LBNO collaboration [33,34] can also be very sensitive to these issues. The same is also true for the future T2HK experiment [35,36], which is a bigger version of T2K. 3 In principle (see the 4-flavor analysis performed in [38]), the CP-phases can impact also the νµ → νe searches of ICARUS [39,40] and OPERA [41]. However, the very low statistics prevents to extract any information on the phases. study of T2K and NOνA. We draw our conclusions in section 6.
2 Conversion probability in the 3+1 scheme

Theoretical framework
In the presence of a sterile neutrino ν s , the mixing among the flavor and the mass eigenstates is described by a 4 × 4 matrix. A convenient parameterization of the mixing matrix is The parameterization in Eq. (2.1) has the following properties: i) When the mixing invoving the fourth state is zero (θ 14 = θ 24 = θ 34 = 0) it returns the 3-flavor matrix in its common parameterization. ii) With the leftmost positioning of the matrixR 34 the vacuum ν µ → ν e conversion probability is independent of θ 34 and of the related CP-phase δ 34 (see [19]).

Analytical Expressions in Vacuum and Matter
Let us now consider the transition probability relevant for T2K and NOνA. As shown in [19], the ν µ → ν e conversion probability can be written as the sum of three contributions (2.4) The first (positive-definite) term is driven by the atmospheric frequency and it gives the leading contribution to the probability. The second and third terms are related to the interference of two distinct frequencies and can assume both positive and negative values. The first of the two interference terms is connected to the standard solar-atmospheric interference, while the second one is driven by the atmospheric-sterile interference. The conversion probability depends on the three small mixing angles θ 13 , θ 14 , θ 24 , whose best estimates, derived from the global 3-flavor (for θ 13 ) analyses [44][45][46] and from the 3+1 fits [42,43]  where ∆ ≡ ∆m 2 31 L/4E is the atmospheric oscillating factor, which depends on the baseline L and the neutrino energy E. The two LBL experiments under consideration, T2K and NOνA, make use of an off-axis configuration, which leads to a narrow-band sharply-peaked energy spectrum of the emitted neutrinos. In theory, the off-axis angle should be tuned exactly to match (at the peak energy) the condition ∆ ∼ π/2, corresponding to the first oscillation maximum. In practice this condition holds only approximately. In T2K, the neutrino flux is peaked at E = 0.6 GeV and the condition ∆ = π/2 is exactly matched. In NOνA, the peak of flux is located at E = 2 GeV, while the oscillations maximum is at E = 1.5 GeV. In the following, when discussing the behavior of the conversion probability we will chose the peak value for both experiments, i.e. we will use E = 0.6 GeV for T2K and E = 2 GeV for NOνA. This will allow a better understanding of the subsequent discussion at the events level presented in section 3. In fact, in both experiments the total rate keeps the leading contribution from the energies close to the peak. Also, it should be stressed that the sensitivity to the spectral distortions is quite limited (due to systematic errors and the limited statistics) and the total rate suffices to understand the basic feature of the numerical results, albeit our analysis includes a full treatment of the spectrum (see section 3).
The presence of matter slightly modifies the transition probability through the MSW effect, which introduces a dependency on the ratio is the constant matter potential along the neutrino trajectory in the earth crust. Both in T2K and in NOνA the value of v is relatively small, being v ∼ 0.05 in T2K, and v ∼ 0.17 in NOνA, where we have taken as a benchmark value the peak energy (E = 0.6 GeV in T2K, E = 2 GeV in NOνA). Therefore, v can be treated as a small parameter of order . The ν µ → ν e conversion probability in matter can be obtained (see the appendix in [19] and the works [47][48][49]) by performing, in the leading term of the the vacuum probability, the following substitution P ATM m (1 + 2v)P ATM , (2.10) which incorporates the (third order) corrections due to matter effects. It can be shown that the two interference terms acquire corrections which are of the fourth order. In this work, we will limit the expansion at the third order in . Therefore, the interference terms will have the vacuum expression. 1.0 Not marginalized Table 1. Parameter values/ranges used in the numerical calculations. The second column reports the true values of the oscillation parameters used to simulate the "observed" data set. The third column depicts the range over which sin 2 θ 23 , δ 13 , and δ 14 are varied while minimizing the χ 2 to obtain the final results.
Before closing this section we recall that a swap in the neutrino mass hierarchy is parametrized by the replacements respect to the vacuum case in the NH (IH) case. NOνA is expected to be more sensitive than T2K to the MH because of the larger value of the ratio v. Finally, we recall that the transition probability for antineutrinos is obtained from that of neutrinos with a change in the sign of the MSW potential V and of all the CP-phases. This, for a given choice of the MH, corresponds to the substitutions In the NH case v > 0 for neutrinos and v < 0 for antineutrinos. According to Eq. (2.10), in the NH case the leading contribution to the transition probability will increase (decrease) for neutrinos (antineutrinos). In the IH case the opposite conclusion holds. Figures 1 and 2  igure 2. ν µ → ν e transition probability as a function of neutrino energy for NOνA after performing the averaging over the fast oscillations. and we have taken the mass-squared difference ∆m 2 41 = 1 eV 2 and fixed the mixing angle θ 13 and the two mixing angles θ 14 and θ 24 at the benchmark values indicated in the second column of table 1, where we also report all the other mass-mixing parameters involved in the calculations. For such high values of the mass-squared splitting, the oscillating factor ∆ 14 is very large and the sterile-induced oscillations are completely averaged out by the finite energy resolution of the detector. Hence, we report the transition probability obtained after that such an averaging process has been taken into account. In each plot of Figs. 1 and 2 the value of the standard CP-phase δ 13 is kept fixed at the value displayed in the legend, while the new phase δ 14 assumes four representative values. In each panel, the 3-flavor probability is represented by a thick black line, while the four 3+1 cases are displayed by thin colored lines. The magenta curve corresponds to δ 14 = −π/2, the blue one to δ 14 = π/2, the red one to δ 14 = 0 and the green one to δ 14 = π. For clarity, we will adopt such color convention in all the figures presented in the paper. From Figs. 1 and 2 it clearly emerges that the impact of the 4-flavor corrections induced by a nonzero value of the mixing angles θ 14 and θ 24 is sizable. Their amplitude and shape depend on the particular value of the new CP-phase δ 14 . The plots show that the most evident effect of the 4-flavor corrections is a change in the overall normalization of the transition probability with respect to the standard 3-flavor case. In addition, an appreciably different energy dependence is also present, which reflects the different dependency from the L/E ratio of the standard interference term [Eq. (2.6)] and the non-standard one [Eq. (2.7)]. The changes induced in the overall normalization are as big as the modifications induced by varying δ 13 (compare the excursion of the black curves between different panels with the excursion of the colored curve in a given panel). This confirms the analytical estimates made in the previous section.

Bi-Probability Plots
In the 3-flavor framework, the behavior of the transition probability is often represented with the CP-phase trajectory diagrams in bi-probability space, first introduced in [50]. Such plots, commonly dubbed as bi-probability plots, represent the parametric curves of the two transition probabilities (ν µ → ν e andν µ →ν e ) where the varying parameter is the CP-phase δ 13 . Since the two transition probabilities are cyclic functions of the phase δ 13 the resulting contours form a closed curve. This representation is particularly advantageous as it gives a bird-eye view of the salient features of a given experimental setups, in particular its sensitivity to MH and CPV. Here we attempt to generalize the biprobability representation to the more general 4-flavor scheme. This makes sense because also in the 4-flavor scheme the probability remains a cyclic function of the (more numerous) CP-phases. Also in this case, as we will show, this representation is very useful for the interpretation of the numerical results. In the following we first recall the basic features of the bi-probability plots in the standard 3-flavor framework, then we generalize our study to the 4-flavor scheme.

The 3-flavor case
In the 3-flavor case the neutrino and antineutrino transition probabilities can be written as (2.14) In general, due to the presence of matter effects one has P 0 =P 0 and A =Ā. As discussed in the previous section the matter effects shift 4 P 0 by an amount proportional to the small the interference term. The amplitude of the interference term is also modified with respect to the vacuum case, its relative change being proportional to v. Since the amplitude A is of order O( 3 ), the corrections are of order O( 4 ). Therefore, truncating the expansion of the probability at the third order corresponds to consider A =Ā, with A 8s 13 s 12 c 12 s 23 c 23 (α∆) sin ∆ . (2.15) The relations (2.13)-(2.14) represent the parametric equations of an ellipse of center (P 0 ,P 0 ). Under the assumption A =Ā, as already shown in [50], the major (minor) axis of the ellipse is proportional to sin ∆ (cos ∆) and has an inclination of −π/4 (π/4). To see this one can perform a counter-clockwise rotation R of the parametric curve around its center (P 0 ,P 0 ) by the angle ω = π/4 obtaining for the rotated probabilities From these relations one arrives at the equation of an ellipse in the canonical form with the two semi-axes having lengths The combination of signs in the parametric equations (2.17)-(2.18) implies that for the NH case the "chirality" of the ellipse is positive, i.e. the trajectory winds in the counterclockwise sense as the phase δ 13 increases. The chirality is opposite (negative) in the IH case since the coefficient A changes its sign under a swap of the mass hierarchy [see Eq. (2.15)]. The energy spectrum of the neutrino beams employed in typical LBL experiments is peaked around the first oscillation maximum, where ∆ ∼ π/2 and therefore one expects a b, i.e. the major axis much bigger than the minor one. This feature is particularly pronounced in T2K since, as already remarked in the previous section, the peak energy (E = 0.6 GeV) almost exactly matches the condition ∆ = π/2. In this case the ellipse becomes almost degenerate with a line. This behavior can be observed in Fig. 3, where in all four panels the two black curves correspond to the 3-flavor limit for the two cases of NH (solid line) and IH (dashed line). The colored curves correspond to four representative 4-flavor cases that will be discussed later. In the NOνA experiment, at the peak energy (E = 2 GeV) we have ∆ = 0.4π and the ratio of the major over the minor axis is given by a/b = tan ∆ ∼ 3, as one can appreciate from the plots in Fig. 4, where again like for T2K, we display the two cases of NH (black solid line) and IH (black dashed line).
The bi-probability representation is particularly useful because it neatly shows that the presence of matter effects tend to split the two ellipses corresponding to the two mass hierarchies, thus giving a qualitative bird-eye view of the sensitivity of a given experiment to the MH. In addition, the ellipse curves show pictorially the effect of the genuine (or intrinsic) CPV due to sin δ 13 , disentangling it from the fake (or extrinsic) CPV induced by the matter effects. In particular, for δ 13 = (−π/2, π/2) the representative point in the bi-probability space (respectively a circle and a square) lies on the intercepts of the ellipse with the major axis and one has the maximal (intrinsic) CPV. Conversely, the effect from the CP conserving cos δ 13 term is proportional to the length of the minor axis. For δ 13 = (0, π) the representative point on the ellipse (respectively a triangle and an asterisk) basically coincide.
From the comparison of Fig. 3 and Fig. 4, it emerges that the splitting between the NH and IH curves is less pronounced in T2K than in NOνA. This is a consequence of the fact that, as already discussed in the previous section, the matter effects are larger in the second experiment. It is useful to recall [see Eq. (2.10)] that the matter effects induce modifications proportional to the dimensionless quantity v = 2V E/∆m 2 31 , which is v 0.05 at the T2K peak energy and v 0.17 at the NOνA peak energy.

The 4-flavor case
In the 3+1 scheme the transition probabilities have the general form where, neglecting O( 4 ) corrections, we have In the expression of the coefficients A and B in Eqs. (2.24)-(2.25) we have left evident the dependency from the oscillation factor ∆ and from the sign of the ratio α. This will be useful when discussing the role of the neutrino mass hierarchy. We note that, under a swap of the MH (implying ∆ → −∆ and α → −α), both coefficients A and B change sign and therefore their product AB remains unaltered and positive definite. The equations (2.22)-(2.23) can be re-expressed in the form where the new coefficients C, D,C,D depend on ∆ and on the CP-phase δ 14 as follows . We can understand the basic behavior of a 4-flavor ellipse using the following relation for its inclination, 6 valid under the assumption that the perturbations induced by the matter effects on the interference terms are negligible (i.e. A =Ā and B =B), (2.37) The angle ω ∈] − π/4, π/4[ represents the inclination (with respect to the axis of the abscissas) of the major (minor) axis depending on the negative (positive) sign of the denominator in Eq. (2.37) (i.e. the sign of sin δ 14 since the product AB is positive definite). In the limit sin δ 14 → 0 the inclination of the major axis is |ω| = π/4. In this case, the sign of ω can be determined by looking at the sign of the numerator in Eq. (2.37). If the numerator is positive, one has ω = −π/4, if it is negative one has ω = π/4. 5 If one makes explicit the dependency on δ14 (instead of δ13) and treats δ14 as the varying parameter one still obtains a (different) ellipse. In this case, the center of the ellipse depends on the value of δ13. 6 For the derivation one has to write the equation of the ellipse by eliminating the parameter δ13 and than use the general formulae available in textbooks. In T2K we have ∆ π/2 and Eq. (2.37) takes the simpler form This result is independent of the mass hierarchy and therefore the inclination of the ellipses will be identical in the two cases of NH and IH. This is confirmed by Fig. 3. cases of NH and IH. We also observe that the centers of the 4-flavor ellipses almost coincide with those of the 3-flavor ones since, as discussed in section 2 [see Eq. (2.10)], the matter effects enter in a similar way in the two schemes. The very small differences in the location of the centers of the 3-flavor and 4-flavor ellipses is imputable to corrections of order O( 4 ), which are neglected in our treatment. For the two values δ 14 = (0, π), the inclination of the major axis is ω = −π/4 since the numerator in Eq. (2.38) is positive in both cases. This is confirmed by the first (red curves) and second (green curves) panel of Fig. 3. For δ 14 = ±π/2, one has tan 2ω = ±0. 22, approximately corresponding to ω ±0.11 (or ±6 0 ). In the case δ 14 = −π/2, the sign of the denominator in Eq. (2.38) is negative and the inclination of −6 0 is that of the major axis. In the case δ 14 = π/2, the sign of the denominator in Eq. (2.38) is positive and the inclination of +6 0 is that of the minor axis. This behavior is corroborated by the third panel (magenta curves) and fourth panel (blue curves) of Fig. 3.
In NOνA we have ∆ 0.4π and Eq. (2.37) takes the form tan 2ω = k 1 1 ∓ k 2 cos δ 14 sin δ 14 . (2.39) where the two constants k 1 , k 2 are given by The minus (plus) sign in Eq. (2.39) refers to the case of NH (IH). So at the NOνA peak energy, which corresponds to a value of ∆ different from π/2, differently from T2K, we expect a dependency of the orientation of the ellipse from the mass hierarchy. For the two values δ 14 = (0, π) the denominator goes to zero so the absolute inclination of the ellipses is |ω| = π/4. The sign of ω is determined by the sign of the numerator, which in the normal hierarchy case is negative for δ 14 = 0 and positive for δ 14 = π. Therefore, in the NH case the inclination of the major axis is π/4 for δ 14 = 0 and −π/4 for δ 14 = π. In the IH case the situation is reversed, since the sign in the numerator in Eq. (2.39) is opposite. This behavior is basically confirmed by the first two panels of Fig. 4. Coming now to the two cases δ 14 = ±π/2, one has tan 2ω = k 1 / sin δ 14 , which is a relation independent of the neutrino mass hierarchy. Due to the small value of the coefficient k 1 , the value of ω is approximately zero. The inclination refers to the major (minor) axis for δ 14 = −π/2 (δ 14 = π/2). This behavior is confirmed by the numerical results displayed in the third and fourth panel in Fig. 4. Hence, one can see that for both experiments the relatively simple formulae illustrated above allow us to explain analytically all the properties of the ellipses displayed in Figs. 3 and 4, which are obtained by a full numerical calculation. In the case of T2K the formula for the inclination of the ellipse is accurate at the level of less than one degree. In NOνA the accuracy, in some cases, is at the level of a few degrees, due to the larger impact of the fourth order corrections related to matter effects. For clarity, in table 2 we report the approximated properties of the ellipses for the 3-flavor and the 4-flavor cases.
The bi-probability plots shown in Figs. 3 and 4 are obtained for fixed values of the CP-phase δ 14 . Since the value of δ 14 is unknown, it is interesting to ask what happens if one superimposes all the (theoretically infinite) ellipses corresponding to all the possible choices of δ 14 . The result of this exercise is shown in figure 5, which has been produced by drawing the convolution of all the ellipses 7 obtained with a dense grid for the parameter δ 14 in its range of variability [−π, π]. Alternatively, Fig. 5 may be seen as a dense scatter plot obtained by varying simultaneously both CP-phases δ 13 and δ 14 . This plot provides a bird-eye view of the degree of separation of the two neutrino mass hierarchies in the 3+1 δ 14 (true) MH Chirality Inclination (T2K) Inclination (NOνA) Table 2. Geometrical properties of the ellipses for the 3-flavor and 4-flavor schemes. The first column reports the value of the CP-phase δ 14 (not defined in the 3-flavor case). The second column reports the neutrino mass hierarchy. The third column reports the chirality of the ellipse (which is the same for T2K and NOνA). The plus (minus) sign means that the trajectory winds in the counter-clockwise (clockwise) sense as the phase δ 13 increases. The fourth and fifth columns report the inclination of the major axis of the ellipse for T2K and NOνA, respectively. The values of the inclinations are those found with the third order expansion of the transition probabilities.
scheme. We see that a separation persists also in such an enlarged scheme. This means that there will exist some combinations of the two CP-phases δ 13 and δ 14 (corresponding to those points which do not lie in the superposition area of the blue and orange regions) for which it will be possible to distinguish between the two hierarchies at some non-zero confidence level. The numerical analysis of the section 5 will allow us to determine such specific combinations of the two CP-phases and the exact confidence level of the separation of the two hierarchies.
3 Experimental features and discussion at the events level

The off-axis experiments: T2K and NOνA
In this section, we briefly mention the key experimental features of the currently running T2K [52,53] and NOνA [54][55][56][57] experiments that go into carrying out the simulation. The T2K experiment in Japan is collecting the data since 2010. Neutrinos are being produced at the J-PARC accelerator facility in Tokai, and are being observed in the 22.5 kton (fiducial) Super-Kamiokande waterČerenkov detector at Kamioka, at a distance of 295 km from the source at an off-axis angle of 2.5 • [52]. Due to the off-axis nature of the beam [58], it peaks sharply at the first oscillation maximum of 0.6 GeV. Another major benefit of using the offaxis technique is that it helps to reduce the intrinsic ν e contamination in the beam and also the background coming from neutral current events, improving the signal-to-background ratio by great extent. As a result, the T2K experiment has already been able to provide an important breakthrough to establish the three-flavor paradigm by observing the θ 13driven appearance signal in ν µ → ν e oscillation channel [59]. In May 2014, T2K started its operation in the antineutrino mode, and after collecting 10% of their expected antineutrino data set, recently they have announced the first appearance results in the antineutrino channel [60,61], clearly taking a first step towards probing the CP symmetry in a direct fashion. As mentioned earlier, in this paper, we consider the full projected exposure of 7.8 × 10 21 protons on target (p.o.t.) which the T2K experiment plans to achieve during their entire run with a proton beam power of 750 kW and with a proton energy of 30 GeV. We also assume that the T2K experiment would use half of its full exposure in the neutrino mode which is 3.9 × 10 21 p.o.t. and the remaining half would be used during antineutrino run. We follow the recent publication by the T2K collaboration [62] in great detail to simulate the signal and background event spectra and their total rates to obtain our final results. Following the same reference [62], we assume an uncorrelated 5% normalization error on signal and 10% normalization error on background for both the appearance and disappearance channels to analyze the prospective data from the T2K experiment. We use the same set of systematics for both the neutrino and antineutrino channels which are also uncorrelated.
The US-based long-baseline experiment NOνA is currently taking data. It uses a 14 kton liquid scintillator far detector at Ash River, Minnesota to detect the oscillated NuMI 8 muon neutrino beam produced at Fermilab [56,57,63]. The NOνA far detector is placed 810 km away from the source at an off-axis angle of 14 mrad (0.8 • ) with respect to the beam line, and sees a narrow-band beam which peaks around 2 GeV. Based on the exposure of 2.74 × 10 20 p.o.t., recently, the NOνA experiment has released their first ν e appearance data providing a solid evidence of ν µ → ν e oscillation over a baseline of 810 km which is the longest baseline in operation now [64][65][66][67]. In this work, we take the full projected exposure of 3.6×10 21 p.o.t. which the NOνA experiment aims to use during their full running time with a NuMI beam power of 700 kW and 120 GeV proton energy [56]. In our simulation, we assume that NOνA would also use 50% of its full exposure in the neutrino mode which is 1.8 × 10 21 p.o.t. and the remaining 50% would be utilized to collect the data in the anti-neutrino mode. Following references [57,68,69], we estimate the signal and background event spectra and their total rates in our calculations. We use a simplified systematic treatment for NOνA: an uncorrelated 5% normalization uncertainty on signal and 10% normalization uncertainty on background for both the appearance and disappearance channels. This is true for both the neutrino and antineutrino modes which are also assumed to be uncorrelated.

Event spectra
We devote this section to discuss the expected event spectra in 3ν and 3+1 schemes for both the T2K and NOνA setups using their full projected exposures as mentioned in the previous section. The number of expected appearance electron events 9 in the i-th energy bin in the detector is estimated using the following well known expression where φ(E) is the neutrino flux, T is the total running time, n n is the number of target nucleons in the detector, is the detector efficiency, σ νe is the neutrino interaction crosssection, and R(E, E A ) is the Gaußian energy resolution function of the detector. The quantities E and E A are the true and reconstructed (anti-)neutrino energies respectively, and L is the baseline. In Fig. 6, we show the expected signal event spectra for the ν e appearance channel as a function of reconstructed neutrino energy for both the experiments under consideration. As expected due to their off-axis nature, we see a narrow peak in the projected event spectrum around 0.6 GeV for the T2K experiment (see the left panel), and for the NOνA experiment (see the right panel), the events mainly occur around 2 GeV where the flux is maximum. In both panels, the thick black lines correspond to the 3ν case assuming δ 13 = 0 • . The other colored histograms (red, green, magenta, and blue) are drawn in the 3+1 scheme assuming different values of δ 14 which are mentioned in the figure legends. Next, we discuss the bi-events plots to get more physics insight.

Bi-events plots
The bi-probability plots presented in section 2.3.2 give a very clear idea of the behavior of the transition probability at the specific value of the energy corresponding to the peak of the spectrum and allow us to approximately predict the behavior of a given off-axis experiment, since the dominant contribution to the total rate comes from the energies close to the peak. In this section we present, for completeness, also the bi-events plots, where on the two axes it is represented the theoretical value of number of events (ν e on the x-axis,ν e on the y-axis) expected in a given experiment. Such plots provide a more precise information on the behavior of a given experiment, because the event rates take into account the complete energy spectrum and provide the information on the statistics involved in the experiment. Figure 7 shows the bi-events plots for T2K, where we have used the same contour style convention of the bi-probability plots. We still have elliptical curves, since we are just replacing the coefficients in the parametric equations (2.32)-(2.35) with appropriate weighted averages. A quick comparison of the T2K bi-events plot in Fig. 7 with the corresponding bi-probability one (Fig. 3) shows that the geometrical properties of the ellipses are slightly different from those obtained for the probabilities. Apart from an obvious deformation factor due to the different scale used for the events, we can appreciate other differences, which are introduced by the contribution of several energies in the integration. In particular, appreciable differences are now visible between the two cases of NH and IH. Most importantly, figure 7 gives a clear feeling on the number of events expected in the experiment T2K. Figure 8 shows the bi-events plots for NOνA. The comparison with the corresponding bi-probability plot in Fig. 4 shows that the geometrical properties of the ellipses are quite similar to those obtained for the probabilities. This is due to the fact the the energy spectrum of NOνA is more sharp than the T2K one (see Fig. 6). As a consequence the peak energy is more important in determining the global behavior of the total rate. Finally, in Fig. 9 we show the convolution plot in the bi-event space, which gives a visual information on the degree of separation of the two neutrino mass hierarchies. igure 7. Bi-events plots for T2K for four fixed values of the CP-phase δ 14 . In each panel, we also show the 3-flavor ellipses for the sake of comparison. In both the 3-flavor and 4-flavor ellipses, the running parameter is the CP-phase δ 13 varying in the range [−π, π]. The solid (dashed) curves refer to NH (IH). We have assumed that half of the full T2K exposure will be used in the neutrino mode and the other half in the antineutrino mode.

Details of the Statistical Method
This section deals with the numerical technique and analysis procedure which we follow to compute our main results. We use the GLoBES software [70,71] along with its new physics tools to obtain our results. We include the 4-flavor effects both in the ν µ → ν e appearance channel and in the ν µ → ν µ disappearance channel. We have found that, for the ν µ → ν µ disappearance channel, the survival probability is very close to the 3-flavor case, in agreement with the analytical considerations made in [19]. We consider the true value of sin 2 2θ 13 to be 0.085 (see table 1) to generate the data and keep it fixed in the fit expecting that the Daya Bay experiment would be able to measure θ 13 with a very In each panel, we also show the 3-flavor ellipses for the sake of comparison. In both the 3-flavor and 4-flavor ellipses, the running parameter is the CP-phase δ 13 varying in the range [−π, π]. The solid (dashed) curves refer to NH (IH). We have assumed that half of the full NOνA exposure will be used in the neutrino mode and the other half in the antineutrino mode.
high-precision (∼ 3% relative precision at 1σ C.L.) by the end of 2017 [72]. As far as the atmospheric mass-squared splitting is concerned, we take the true value of ∆m 2 32 = 2.4 × 10 −3 eV 2 (−2.4 × 10 −3 eV 2 ) for NH (IH). Accordingly, we take ∆m 2 31 = 2.475 × 10 −3 eV 2 (−2.4 × 10 −3 eV 2 ) for NH (IH). We also do not marginalize over this parameter in the fit since the present precision on this parameter is already quite good [46], and the future data from the running T2K and NOνA experiments would certainly improve this further [62,73]. This should remain true also in the presence of sterile neutrino oscillations because the value of ∆m 2 31 is extracted from the ν µ → ν µ searches which, for the small values of the mixing angles θ 14 and θ 24 considered in the present analysis (see table 1), are almost unaffected by the 4-flavor effects. For θ 23 , we consider the maximal mixing For both T2K and NOνA we have assumed that half of the full exposure will be used in the neutrino mode and the other half in the antineutrino mode.
(π/4) as the true choice, and in the fit, we marginalize over the range given in table 1. We marginalize over both the choices of hierarchy in the fit for all the analyses, except for the mass hierarchy discovery studies where our aim is to exclude the wrong hierarchy in the fit. We vary the true value of δ 13 in its allowed range of [−π, π], and it has been marginalized over its full range in the fit if the performance indicator demands so. We take the line-averaged constant Earth matter density 10 of 2.8 g/cm 3 for both the baselines.
We take the mass-squared splitting ∆m 2 41 = 1 eV 2 , which is the value currently suggested by the SBL anomalies. However, we stress that our results would remain unaltered for different choices of such parameter, provided that ∆m 2 41 0.1 eV 2 . For such values, the fast oscillations induced by the new large frequency get completely averaged because of the finite resolution of the detector. For the same reason, the LBL setups are insensitive to the sign of ∆m 2 41 and we can safely assume positive sign for it. Concerning the active-sterile mixing angles, we take the true value of 0.025 for both the sin 2 θ 14 and sin 2 θ 24 and keep them fixed in the fit. These values are close to the best fit obtained by the global 3+ fits [42]. We vary the true value of δ 14 in its allowed range of [−π, π], and it has been marginalized over its full range in the fit as needed. We assume sin 2 θ 34 = 0 and δ 34 = 0 in our simulation. 11 In our analysis, we do not explicitly consider the near detectors of 10 The line-averaged constant Earth matter density has been computed using the Preliminary Reference Earth Model (PREM) [74]. 11 We recall that the vacuum νµ → νe transition probability is independent of θ34 (and δ34). In matter, a tiny dependence appears which is more appreciable in NOνA than in T2K (see the appendix of [19] for T2K and NOνA which may shed some light on θ 14 and θ 24 , but certainly, the near detector data are not sensitive to the CP-phases which is the main thrust of this work. In our simulation, we have performed a full spectral analysis using the binned events spectra for both experiments. In the statistical analysis, the Poissonian ∆χ 2 is marginalized over the uncorrelated systematic uncertainties (as mentioned in section 3) using the method of pulls as discussed in Refs. [75,76]. When showing the results, we display the 1, 2, 3σ confidence levels for 1 d.o.f. using the relation nσ ≡ ∆χ 2 . In [77], it was shown that the above relation is valid in the frequentist method of hypothesis testing.

CP-violation Searches in the Presence of Sterile Neutrinos
In this section we explore the impact of sterile neutrinos in the CPV searches of T2K and NOνA. As a first step we consider the discovery potential of the CPV induced by the standard 3-flavor CP-phase δ 13 , which is proportional to sin δ 13 . The discovery potential is defined as the confidence level at which one can reject the test hypothesis of no CPviolation, i.e. the cases δ 13 = 0 and δ 13 = π. We have taken the best fit values of all the parameters at the values specified in the second column of table 1. In the 3-flavor scheme, we marginalize over θ 23 and over the hierarchy. In the 3+1 scheme, in addition, we marginalize over the unknown value of δ 14 .
In Figure 10 we display the results of the numerical analysis. The upper panels refer to T2K, the middle ones to NOνA, and the lower ones to their combination. In the left (right) panels, we consider NH (IH) as the true hierarchy choice. In each panel, we present the results obtained for the 3-flavor case (black solid curve) and for the 3+1 scheme, in which case we select four different values of the true value of δ 14 (while its test value is left free to vary and is marginalized away). The values of the phase δ 14 and the colors of the corresponding curves are the same of the previous plots. The 3-flavor sensitivities (black curves) are in agreement with those shown in the official analyses [62]. We see that for all values of the new CP-phase δ 14 the discovery potential of the two experiments decreases with respect to that of the 3-flavor case. The loss of sensitivity is imputable to the degeneracy between the two CP-phases δ 13 and δ 14 . Similar to the 3-flavor case, the discovery potential has a maximum for δ 13 = −90 0 (δ 13 = 90 0 ) for NH (IH). Abrupt changes in the sensitivity are evident in the range [45 0 , 135 0 ] for the NH case and in the range [−135 0 , −45 0 ] for the IH case. This behavior can be traced to the degeneracy among the two CP-phases and the mass hierarchy. In fact, in these ranges the best fit is obtained for the false hierarchy. In the bi-events plots these ranges correspond to points where the ellipses of the two hierarchies tend to overlap.
Until now we have considered only four selected values of the CP-phase δ 14 . It is interesting to see what happens for a generic choice of such a parameter. To this purpose we have generalized the analysis by treating δ 14 as a free parameter. In Fig. 11 we show the iso-contour lines of the discovery potential of the CP-violation induced by δ 13 as a function a detailed discussion).
of the true values of the two phases δ 13 and δ 14 . Inside the red regions the discovery potential is larger than 2σ. Inside the beige regions it is larger than 1σ. The plots refer to the combination T2K + NOνA for the two cases of NH (left) and IH (right). One can easily check that horizontal cuts of the contour plots made in correspondence of the four particular values of the phase δ 14 considered in Fig. 10 return the 1σ and 2σ intervals derivable from the last two panels of Fig. 10.
In the 3+1 scheme also the new CP-phase δ 14 can be a source of CP-violation. Hence it is interesting to determine the discovery potential of the CPV induced by sin δ 14 alone and the total CP-violation induced simultaneously by sin δ 13 and sin δ 14 . Our numerical analysis shows that the discovery potential of non-zero sin δ 14 is always below the 2σ level so we do not show the corresponding plot. Instead in Fig. 12 we show the results for the total CPV discovery since it is appreciably different (larger) from that induced by sin δ 13 alone. The plots refer to the combination T2K + NOνA for the two cases of NH (left) and IH (right). Inside the small green regions the discovery potential is larger than 3σ. Inside the red regions it is larger than 2σ. Inside the beige regions it is larger than 1σ. As expected the regions in Fig. 12 contain as sub-regions those of Fig. 11 where the sole CPV induced by sin δ 13 is considered. In Fig. 12 we can recognize nine white regions where the discovery potential is below 1σ. We stress that these nine regions correspond to four physical regions because both CP-phases are cyclic variables. Geometrically one can view the square represented in Fig. 12 as an unwrapped torus, provided one takes into account that the upper edge is connected with the lower edge, and that left edge with the right one 12 . These four regions contain the points where the total CPV is zero i.e. sin δ 13 = sin δ 14 = 0. This condition is verified for the four (inequivalent) CP-conserving cases [δ 13 , δ 14 ] = [0, 0], [π, 0], [0, π], [π, π] and all the other five combinations obtainable by a change of sign of one of (or both) the phases equal to π.  Figure 10. Discovery potential of CP-violation induced by sin δ 13 . Upper panels refer to T2K. Middle panels to NOνA. Lower panels to T2K and NOνA combined. In the left (right) panels, we consider NH (IH) as the true hierarchy choice. In each panel, the black curve corresponds to the 3-flavor case. The colored curves are obtained in the 3+1 scheme for four different true values of δ 14 . We marginalize over θ 23 and δ 14 over their allowed ranges in the fit, and also over the hierarchy.   in the 3+1 scheme for the T2K + NOνA combined setup. Inside the green regions the discovery potential is ≥ 3σ. Inside the red regions it is ≥ 2σ. Inside the beige regions it is ≥ 1σ.

Reconstruction of the CP phases
The CP-violation discovery potential tells us how much one experiment will be able to rule out the case of CP conservation given a positive observation of CPV corresponding to a true value of the phases δ 13 and δ 14 . While this is certainly a very important feature, a complementary information is provided by the capability of reconstructing the values of the two CP-phases, independent of the amount of CP-violation (if any). Figure 13 gives a quantitative answer to such a different kind of question. The four plots represent the regions reconstructed around four representative points in the plane [δ 13 , δ 14 ]. In all cases we have taken the NH as the true hierarchy in the data and then marginalized over NH and IH in theory. Similar results (not shown) were obtained for the IH case. The two upper panels refer to the CP-conserving cases [0, 0] and [π, π] respectively. The third and fourth panels refer to the two (maximally) CP-violating cases [−π/2, −π/2] and [π/2, π/2]. The two confidence levels refer to 1σ and 2σ (1 d.o.f.). We see that in all cases we obtain a unique reconstructed region at the 1σ level. Note that this is true also in the second panel, because the four corners of the square form a connected region due to the cyclic properties of the two CP-phases. At the 2σ level we obtain a unique region only in the case [−π/2, −π/2] (bottom left). Small spurious islands start to appear in the other cases.
We have checked that these islands disappear if one assumes the prior knowledge of the correct mass hierarchy. In all cases the typical 1σ uncertainty is about 40 0 (50 0 ) for δ 13 (δ 14 ). As recently shown in the 4-flavor analysis performed in [37], the present data seem to indicate a slight preference for the combination [δ 13 , δ 14 ] = [−π/2, −π/2]. If this trend gets confirmed in a few years (assuming the existence of a sterile neutrino), the picture should resemble that of the left bottom panel of Fig. 13.
(IH) as the true hierarchy choice. In each panel, we give the results for the 3-flavor case (thick black curve) and for the 3+1 scheme (colored curves) for four different values of the true δ 14 (that is −90 0 , 90 0 , 0 0 and 180 0 ). The color convention is the same adopted in the rest of the paper.
We observe that in T2K (upper panels) the discovery potential is quite limited both in the 3-flavor framework and in the 3+1 scheme. This is due to the fact that the matter effects are small in T2K. Comparing the results of the 3+1 scheme (colored curves) with those of the 3-flavor one (black curve) we observe that, apart for the case δ 14 = −90 0 for NH (δ 14 = 90 0 for IH), in all the other cases the discovery potential is smaller than the 3-flavor one. The overall behavior of the 4-flavor curves is similar to that of the 3-flavor one. In particular, the sensitivity presents a maximum at δ 13 = −90 0 for NH. This similar behavior can be understood by observing that in the bi-events plots the point δ 13 = −90 0 (the squares in the four panels of Fig. 7) always provides the maximal separation from the cloud generated by the convolution of all the possible IH ellipses (see Fig. 9). A similar observation can be done for the specular case of δ 13 = +90 0 and IH.
Concerning NOνA (middle panels) we can make the following observations. Similarly to T2K the maximal discovery potential is obtained for δ 13 = −90 0 for NH and δ 13 = 90 0 for IH. For such two values the representative points on the ellipse (the squares for NH and the circles for IH in Fig. 8) always provide the maximal separation from the convolution of all the ellipses of the opposite MH (see again Fig. 9). However, there are also important differences with respect to T2K. First of all, we observe that the maximal discovery potential is much larger than that of T2K. This is imputable to the fact that the matter effects are much bigger in NOνA (see the discussion in section 2). Second, we can see that in the NH case (left middle panel) there is a good sensitivity not only for δ 14 = −90 0 (magenta curve) but also for δ 14 = 180 0 (green curve). In the IH case (right middle panel) there is a good sensitivity for δ 14 = 90 0 (blue curve) and δ 14 = 0 0 (red curve). This different behavior with respect to T2K can be traced to the fact that in NOνA the peak energy is not centered exactly at the first oscillation maximum but at ∆ = 0.4π. Finally, we notice that the combination of the two experiments (lower panels) is dominated by NOνA.
The study of the discovery potential of the neutrino mass hierarchy can be generalized to the case in which the CP-phase δ 14 can assume any value in its variability range. Figure 15 shows the results of such more general analysis, where we have treated δ 14 as a free parameter. We display the iso-contour lines of the discovery potential as a function of the true values of the two phases δ 13 and δ 14 . Inside the red regions the discovery potential is larger than 3σ. Inside the blue regions it is larger than 2σ. The plots refer to the combination T2K + NOνA for the two cases of NH (left) and IH (right). One can easily check that horizontal cuts of the contour plots made in correspondence of the four particular values of the phase δ 14 considered in Fig. 14 return the 1σ and 2σ intervals derivable from the last two panels of Fig. 14.   Figure 15. Discovery potential for excluding the wrong hierarchy for the combination of T2K and NOνA as a a function of the two CP-phases δ 13 and δ 14 . Inside the red regions the discovery potential is ≥ 3σ. Inside the blue regions it is ≥ 2σ.

Conclusions and Outlook
We have considered the impact of light sterile neutrinos on the prospective data expected to come from the two long-baseline experiments T2K and NOνA when the planned full exposure will be reached. We have presented a detailed discussion of the behavior of the 4-flavor ν µ → ν e andν µ →ν e transition probabilities, extending for the first time the biprobability representation, commonly used in the 3-flavor framework, to the 3+1 scheme.
We have also performed a comprehensive sensitivity study of the two experiments (taken alone and in combination) in order to assess their discovery potential in the presence of a sterile neutrino species. We have considered realistic benchmark values of the 3+1 massmixing parameters as inferred from the existing global short-baseline fits. We found that the performance of both the experiments in claiming the discovery of the CP-violation induced by the standard CP-phase δ 13 ≡ δ, and the neutrino mass hierarchy get substantially deteriorated. The degree of loss of sensitivity depends on the value of the unknown CPphase δ 14 . We have also assessed the discovery potential of total CP-violation (i.e., induced simultaneously by the two CP-phases δ 13 and δ 14 ) and the capability of the two experiments of reconstructing the true values of such CP-phases. The typical (1σ level) uncertainty on the reconstructed phases is approximately 40 0 for δ 13 and 50 0 for δ 14 .
In the eventuality of a discovery of a sterile neutrino at the new short-baseline experiments, we will face two challenges. First, we will have to reassess the status of the 3-flavor parameters whose best fit values will change in the 3+1 scheme. Among the 3-flavor properties, the most sensitive to the perturbations induced by the sterile neutrino oscillations are the CP-phase δ, and the neutrino mass hierarchy. In both cases, their determination is based on the observation of very tiny effects which, in the LBL setups, can be appre-ciably perturbed by new interference phenomena induced by the sterile neutrinos. Our study gives the first quantitative assessment of the discovery potential of the 3-flavor CPV and of the MH for the two LBL experiments T2K and NOνA. The second, perhaps more stimulating challenge, will be that of determining all the new parameters that govern the enlarged 3+1 scheme. According to our study, T2K and NOνA may be able to give the first indications on one of the new CP-phases involved in the 3+1 scheme. The future LBL experiments (DUNE, LBNO and T2HK) will be needed to extract more robust information on the enlarged CP-violation sector. We hope that the comprehensive analysis presented in this paper may play an important role in exploring light sterile neutrinos at the long baseline facilities.