Exploring the new physics phases in 3+1 scenario in neutrino oscillation experiments

The various global analyses of available neutrino oscillation data indicate the presence of the standard $3+0$ neutrino oscillation picture. However, there are a few short baseline anomalies that point to the possible existence of a fourth neutrino (with mass in the eV-scale), essentially sterile in nature. Should sterile neutrino exist in nature and its presence is not taken into consideration properly in the analyses of neutrino data, the interference terms arising due to the additional CP phases in presence of a sterile neutrino can severely impact the physics searches in long baseline (LBL) neutrino oscillation experiments. In the current work we consider one light (eV-scale) sterile neutrino and probe all the three CP phases ($\delta_{13}$, $\delta_{24}$, $\delta_{34}$) in the context of the upcoming Deep Underground Neutrino Experiment (DUNE) and also estimate how the results improve when data from NOvA, T2K and T2HK are added in the analysis. We illustrate the $\Delta \chi^2$ correlations of the CP phases among each other, and also with the three active-sterile mixing angles. Finally, we briefly illustrate how the relevant parameter spaces in the context of neutrinoless double beta decay get modified in light of the bounds in presence of a light sterile neutrino.


Introduction
The successful discovery of the phenomena of neutrino oscillation [1,2] has led to an impressive amount of research with an aim to establish the standard 3-flavour (referred to as 3+0 hereafter) oscillation over a wide range of energy (E) and neutrino propagation length (L). In the commonly used Pontecorvo-Maki-Nakagawa-Sakata (PMNS) parametrization [3] of the leptonic mixing matrix, this requires fitting of the six standard oscillation parameters, namely three mixing angles (θ 12 , θ 13 , θ 23 ), two mass squared differences (∆m 2 21 = m 2 2 − m 2 1 , ∆m 2 31 = m 2 3 − m 2 1 ) and one Dirac CP phase (δ 13 ). A complete knowledge about the oscillation parameters will help to shed light on the following still unresolved issues in the neutrino sector: whether there exists CP violation (CPV) in the leptonic sector (i.e., whether δ 13 = 0, π), whether the neutrino mass eigenstates are arranged in normal ordering (NO, i.e., ∆m 2 31 > 0) or inverted ordering (IO, i.e., ∆m 2 31 < 0), and, whether θ 23 lies in the higher octant (HO, i.e., θ 23 > π/4) or in the lower octant (LO, i.e., θ 23 < π/4). If leptonic CPV exists, it would provide a crucial missing ingredient [4] in resolving another very fundamental elusive puzzle that is baryon asymmetry in the observed universe via a mechanism called Leptogenesis [5]. Determination of mass ordering and θ 23 octant will help in understanding the origin of neutrino mass [6][7][8], its Dirac/Majorana nature via neutrinoless double beta decay [9] and in exploring a new symmetry called µ − τ symmetry [10,11]. Presently running long baseline (LBL) experiments such as Tokai to Kamioka (T2K) [12] and NuMI Off-axis ν e Appearance (NOνA) [13] have started uncovering a few of the open issues mentioned above. Latest T2K results [14] have been able to rule out a large range of values of δ 13 around π/2 at 3σ confidence level (C.L.) around L/E ∼ 1 km/GeV [56,57]. In LBL experiments where L/E ∼ 500 km/GeV at the far detectors (FD), the high frequency induced by ∆m 2 41 ∼ 1 eV 2 gets averaged out due to the finite energy resolution of the detector. Nevertheless, it has been shown [56, that even at the FDs of these LBL experiments, the interference effects provided by the additional CP phases play very significant roles in spoiling the sensitivities to the crucial issues of CPV, MH and θ 23 octant 1 . For instance, as shown in [64], the CPV sensitivity becomes a wide band whose width depends on the unknown magnitudes of the sterile phases δ 24 and δ 34 , -leading to serious confusion in interpreting the results as CP violation or CP conservation. The constraints on the active-sterile mixing angles θ i4 (i = 1, 2, 3) do of course reduce such ambiguities in the interpretation of the results to some extent. But, a clear idea about how the LBL experiments are able to measure the sterile phases δ 24 and δ 34 given the constraints on the active-sterile mixing angles, will certainly minimise the obfuscation when their data are analyzed with a view to tackle the unresolved physics issues. The issue of measurement of one sterile phase (δ 24 or δ 14 , depending on the parametrization used) has been addressed to some extent in literature. [64] analyzes how the sterile phases impact the standard CPV and mass ordering measurements and also probe the joint parameter space δ 13 −δ 24 for DUNE. The authors of [62] use a slightly different parameterization and discuss the CPV arising from the individual sterile CP phases and probe the joint parameter space δ 13 − δ 14 at DUNE. They further extend their analysis by combining simulated data from other LBL experiments such as T2HK [72] and ESSνSB [77]. Exploring the parameter space of δ 13 − δ 24 has also been addressed in [71] for DUNE, T2HK, T2HKK and their combinations. The authors of [71] further illustrate how the individual phase δ 24 can be measured by these experiments at various C.L. by assuming four possible true values (0, ±π/2, π). [80] shows how the recent data from T2K and NOνA can help to probe the parameter space of δ 13 − δ 14 . Most recently, the authors of [84] estimate how the difference (δ 14 − δ 24 ) in the sterile phases can impact in constraining the standard δ 13 − θ 23 parameter space.
It is noteworthy that the sterile CP phase δ 34 and its correlation with the other phases has been little addressed in literature. In the present manuscript we tackle the very relevant issue of estimating the capability to reconstruct all three CP phases (δ 13 , δ 24 , δ 34 ), taking into consideration their ∆χ 2 correlations with each other and also with the active-sterile mixing angles (θ 14 , θ 24 , θ 34 ). We carry out this exercise in the context of DUNE and illustrate the improvement when combined with T2K, NOνA (both these currently running experiments are simulated upto their present exposure) and T2HK. Apart from studying these CP phases in detail, another crucial aspect in which our analysis differs from the existing studies mentioned above is that we have taken into consideration the current 3σ hint of CP violation and the corresponding exclusion region of the CP phase δ 13 by T2K data [14]. Moreover, in addition to illustrating how the 2-d parameter spaces for the CP phases can be probed, we also analyze how the individual CP phases can be reconstructed (after marginalizing all other relevant parameters) by the experiments, irrespective of the actual value they might have in nature. Additionally we have also considered the ν µ → ν τ channel (in addition to ν µ → ν e and ν µ → ν µ ) in our study and estimated the capability of the projected data to measure all three CP phases in detail. This enables us to probe the parameter spaces associated to θ 34 and δ 34 with better sensitivity.
The present manuscript is organised as follows. In Sec. 2 we give a brief account of the constraints on the sterile neutrino parameters and how the CP phases affect the relevant probabilities. In Sec. 3 we describe the methodology of our statistical analyses. In Sec. 4 we discuss the ∆χ 2 correlations among various CP phases, taking a pair of phases at a time. We also discuss the potential to reconstruct the three CP phases for all possible true values by performing simulations of DUNE, T2K, NOvA and T2HK. The role of individual oscillation channels in such reconstructions are analyzed in Sec. 5. We estimate the ∆χ 2 correlations among all the CP phases and the active-sterile mixing angles in Sec. 6. Finally in Sec. 7 we briefly discuss how the relevant parameter spaces associated to Neutrinoless Double Beta Decay gets modified in light of the constraints on one eV scale sterile neutrino, followed by conclusion.

Basics
We first discuss the oscillation probabilities for the three channels (P (ν µ → ν e ), P (ν µ → ν µ ) and P (ν µ → ν τ )) in 3+1 scenario. Since the expressions become immensely complicated in matter, we show them in vacuum and these will act as useful templates for gaining the physics insights in explaining our subsequent sensitivity results. For the mixing matrix we follow the parametrization scheme adopted in [56]: where R(θ ij , δ ij ) is a rotation in the ij−th plane with an associated phase δ ij such that, for e.g., Now we discuss briefly on the allowed ranges of the sterile sector parameters, as estimated in great detail in the global analysis of various neutrino data [85]. |U e4 | 2 is bounded by ν e andν e disappearance searches and is equal to sin 2 θ 14 . The combined atmospheric neutrino data from IceCube, DeepCore and SK, at 99% C.L. (2 DOF) put the bound |U e4 | 2 0.1. This implies θ 14 18.4 • . It can be seen from all ν e andν e searches that the best fit value of |U e4 | 2 is approximately equal to 0.01, which gives θ 14 ≈ 5.7 • . The data from ν µ andν µ disappearance searches put the following 99% C.L. (2 D.O.F.) constraints: |U µ4 | 2 0.01 and |U τ 4 | 2 0.17. Since, in our parametrization |U µ4 | = cos θ 14 sin θ 24 and |U τ 4 | = cos θ 14 cos θ 24 sin θ 34 , the corresponding bounds on θ 24 and θ 34 can easily be translated to θ 24 6.05 • and θ 34 25.8 • . The allowed values of ∆m 2 41 roughly lie in the range of 1 − 10 eV 2 and we consider it to be 1.3 eV 2 in the present work as per the global analysis 2 . Using the standard approach for deriving oscillation probability, we obtain for the ν α → ν β (α, β = e, µ, τ, s and α = β) transition probability, In deriving Eq. 2.3, we have used the unitarity of the 4 × 4 mixing matrix U and applied the usual assumptions that the term containing mass square splitting between m 4 and m i (i = 1, 2, 3), i.e., sin 2 ∆ 4i and sin 2∆ 4i average out to 0.5 and 0 respectively at long baseline (i = 1, 2, 3). From Eq. 2.3 one can use the long baseline approximation (i.e., neglecting the oscillation effects due to ∆m 2 21 ) and arrive at the following simplified expression for the dominant channel ν µ → ν e . where we have followed the convention of [56] for the following quantities.
Eq. 2.4 tells us that in vacuum P (ν µ → ν e ) is sensitive to both δ 13 and δ 24 , but not to δ 34 . As explained in [56], small dependence on δ 34 creeps in when matter effect is taken into account. The expressions for the less dominant channels P (ν µ → ν µ ) and P (ν µ → ν τ ) can similarly be derived from Eq. 2.3, but the expressions are quite lengthy. Interested readers can see [75,86] for those probability expressions.
Using the widely used General Long Baseline Experiment Simulator (GLoBES) [87,88] and the relevant plugin snu.c [89,90] for implementing sterile neutrinos, we now illustrate how the probabilities for different oscillation channels depends on the three CP phases (δ 13 , δ 24 , δ 34 ) individually at the DUNE baseline of 1300 km. In Fig. 1, we plot the bands for P µe , P µµ , P µτ (in the three panels respectively) due to individual variation of δ 13 , δ 24 , δ 34 , each being varied in the range [−π, π]. The grey band shows the variation of the standard 2 In our notation, m1,2,3 are the masses of the three active neutrinos, and m4 denotes the mass of the sterile neutrino. Also, ∆m 2 ij = m 2 i − m 2 j .

Parameter
Best-fit-value 3σ interval 1σ uncertainty  Table 1. Standard oscillation parameters and their uncertainties used in our study. The values of 3+0 parameters were taken from the global fit analysis in [22] while the 3+1 parameter values were chosen from [85] (see also Sec. 2). If the 3σ upper and lower limit of a parameter is x u and x l respectively, the 1σ uncertainty is (x u − x l )/3(x u + x l )% [17]. For the active-sterile mixing angles, a conservative 5% uncertainty was used on sin 2 θ i4 (i = 1, 2, 3). and δ34 (red) in the whole range of [−π, π] at a baseline of 1300 km. The three panels correspond to the three channels P (νµ → νe), P (νµ → νµ) and P (νµ → ντ ). The insets in the second and third panels show magnified versions of the rectangular regions indicated. The three active-sterile mixing angles were taken as θ14 = 10 • , θ24 = 6 • , θ34 = 25 • . The values of the rest of the oscillation parameters were taken from Tab. 1. Normal hierarchy was assumed for generating this plot.
Dirac CP phase δ 13 . The variation due to δ 24 (δ 34 ) is shown with the blue (red) band. For the variation of each CP phase, the other two phases are kept fixed: δ 13 is kept fixed at the bf value of −0.8π while δ 24 , δ 34 were considered to be zero. The sterile phases are associated with the active-sterile mixing angles. We have a slightly high active-sterile mixing: θ 14 = 10 • , θ 24 = 6 • , θ 34 = 25 • . As expected 3 , the CP phases have larger impact on appearance channels rather than the disappearance channels. ν µ → ν e channel is most affected by the variation of δ 13 . Among the sterile CP phases, δ 24 has a visibly greater impact on P (ν µ → ν e ) than that of δ 34 . It is interesting to note this feature especially in light of the fact that the value of the active-sterile mixing angle θ 24 (taken as 6 • in Fig. 1) is almost 5 times smaller than θ 34 (25 • ). We also observe that with increase in energy the effect of δ 34 on P (ν µ → ν e ) further reduces. P (ν µ → ν τ ) is less prone to variation of the CP phases. Still it shows slight variation for all three CP phases with δ 24 and δ 34 having almost similar effects but less than that of δ 13 . P (ν µ → ν µ ) on the other hand, shows almost negligible variation with the CP phases. After a brief description of the simulation procedure we will estimate the capability of LBL experiments to reconstruct these CP phases.

Simulation details
We simulate the long baseline neutrino experiments DUNE, NOvA, T2K and T2HK using GLoBES [87,88]. DUNE is a 1300 km long baseline experiment employing a liquid argon far detector (FD) of 40 kt fiducial mass with a beam of power 1.07 MW and running 3.5 years each on ν andν mode (resulting in a total exposure of roughly 300 kt.MW.yr corresponding to 1.47 × 10 21 protons on target or POT). We have used the official configuration files [91] provided by the DUNE collaboration for its simulation. Following this, we have also taken into account the presence of a near detector (ND) at 459 m from the source. The ND helps in making a more precise measurement of the flux and cross-section, thereby reducing the relevant systematic uncertainties at the FD. We should mention here that we have done a ∆χ 2 analysis (discussed later) with the simulated data at FD alone, rather than a joint ∆χ 2 analysis using simulated data both at ND and FD 4 . Electron neutrino appearance signals (CC), muon neutrino disappearance signals (CC), as well as neutral current (NC) backgrounds and tau neutrino appearance backgrounds (along with the corresponding systematics/efficiencies etc.) are already included in the configuration files.
In the present analysis, we have additionally incorporated tau neutrino appearance as a separate signal following [75,76]. Charged current interaction of an incoming ν τ produces a τ lepton (requires a threshold energy of 3.4 GeV for the incoming ν τ ), which can decay hadronically (with a branching fraction ∼ 65%) or leptonically (with a branching fraction of ∼ 35%). The analysis of the hadronic decay channel involves the capability of the detector to study the resulting pions and kaons. More importantly, NC neutrino scattering constitutes the biggest background for the hadronic decay channel of τ . Following [75], we have used an efficiency to separate 30% hadronically decaying τ events (with about 1% NC events remaining). On the other hand, the leptonic decay channels of τ (τ − → e −ν e ν τ ; τ − → µ −ν µ ν τ ) are more difficult to analyse, due to the large background mainly consisting of ν e -CC and ν µ -CC respectively (along with backgrounds from NC and contaminations due to wrong sign leptons.). Following [76], we have taken the efficiency of the electron channel to be 15%. Due to the overwhelming background, we have taken a nominal efficiency of 5% in the muon channel. Naturally, the contribution of the leptonic decay channel of τ is very small. We should also mention here that the decay of τ at the detector will involve missing energy in the form of an outgoing ν τ , which in turn makes the energy reconstruction of the incoming ν τ difficult. From [75], we use a Gaussian energy reconstruction with a resolution of 20% which is a conservative estimate. We acknowledge that our implementation of the ν τ channel as a signal is conservative in nature. Nevertheless, this provides a small but non-negligible statistics in terms of events and ∆χ 2 sensitivity (see Sec. 5 and Appendix C). Using a much more sophisticated analysis of ν τ appearance channel at DUNE by implementing jet-clustering algorithms and machine learning techniques, as has been pioneered in [92], one certainly expects to exploit the rich physics capabilities hidden within this channel.
We have simulated NOvA with a baseline of 800 km employing an FD of fiducial mass of 14 kt and a beam of 742 kW. The simulation for NOvA was implemented according to [15] which generates 8.85 × 10 20 (12.33 × 10 20 ) POT in ν (ν) mode. T2K is a 295 km experiment with a 22.5 kt water cherencov FD. For T2K simulation we use the inputs from [14,93]. We have used a beam of 515 kW and simulating 1.97 × 10 21 (1.63 × 10 21 ) POT in ν (ν) mode. T2HK is an upgraded version of T2K with a higher beam of 1.3 MW and a much bigger fiducial mass of 187 kt of its water cherencov FD. For T2HK we simulate a total of 2.7 × 10 22 POT in 1:3 ratio of ν andν mode (with inputs taken from [19,94]). Note that, for the future experiments DUNE and T2HK we have used the full expected exposure, while for the currently running experiments T2K and NOvA we have simulated upto their current exposure.
To estimate the capability of LBL experiments to reconstruct the CP phases, we carry out a ∆χ 2 analysis using GLoBES and the relevant plugin snu.c [89,90] for implementing sterile neutrinos. In order to gain insight, let us examine the analytical form of the ∆χ 2 , is the set of true (test) oscillation parameters. The index i is summed over the energy bins of the experiment concerned (as discussed above). The indices j and k are summed over the oscillation channels (ν µ → ν e , ν µ → ν µ , ν µ → ν τ ) and the modes (ν andν) respectively. The term (N test − N true ) takes into account the algebraic difference while the log-term inside the curly braces considers the fractional difference between the test and true sets of events. The term summed over i, j, k and written inside the curly braces is the statistical part of ∆χ 2 . σ p l is the uncertainty in the prior measurement of the l-th oscillation parameter p l . The values of the true or best-fit oscillation parameters and their uncertainties as used in the present analysis are tabulated in Table 1. In the last term, σ ηm is the uncertainty on the systematic/nuisance parameter η m and the sum over m takes care of the systematic part of ∆χ 2 . This way of treating the systematics in the ∆χ 2 calculation is known as the method of pulls [95][96][97][98]. For DUNE, the ν e andν e signal modes have a normalization uncertainties of 2% each, whereas the ν µ andν µ signals have a normalization uncertainty of 5% each. The ν τ andν τ signals have a normalization uncertainties of 20% each. The background normalization uncertainties vary from 5% − 20% and include correlations among various sources of background (coming from beam ν e /ν e contamination, flavour misidentification, neutral current and ν τ ). The final estimate of ∆χ 2 obtained after a marginalization over the 3σ range of test parameters 5 p test and the set of systematics (η) is a function of the true values of the oscillation parameters. Technically this ∆χ 2 is the frequentist method of hypotheses testing [96,99].

Reconstruction of the CP phase in correlation with other phases
In Fig. 2, we illustrate how the CP phases can be reconstructed at 1σ C.L. in the plane of (test δ 13 -test δ 24 ), (test δ 13 -test δ 34 ) and (test δ 24 -test δ 34 ) in the three columns respectively. The red contour shows the reconstruction capability of DUNE. The green and blue contours illustrate the reconstruction when it is combined with (T2K + NOνA) and (T2K + NOνA + T2HK) respectively. The top (bottom) row of Fig. 2 shows the choice of the CP conserving (maximally CP violating) true values of δ 24 and δ 34 . Marginalisation has been carried over the test values of θ 13 , θ 23 , ∆m 2 31 (including both the mass hierarchies) with prior uncertainties mentioned in Tab. 1. the active-sterile mixing angles θ 14 , θ 24 , θ 34 (with a prior uncertainty of 5% on each sin 2 θ i4 ) and the third CP phase not shown along the axes of a particular panel has also been considered for marginalisation. Clearly the reconstruction of δ 13 is better (as evidenced by the narrowness of the contours along the test δ 13 axis) compared to the other two phases. This happens since the associated mixing angle θ 13 has been measured very precisely unlike the corresponding active-sterile mixing angles θ 24 and θ 34 (see Table 1 and Sec. 2). For the maximally CP violating choices of true δ 24 and true δ 34 , their reconstruction gets better at the cost of small degeneracies appearing for DUNE around test δ 24 ≈ 135 • , δ 34 ≈ 90 • . Adding data from other experiments lifts these degeneracies. We would also like to refer the reader to Appendices A and B, 5 The test parameters that are marginalised include those listed in Tab. 1 along with the specified ranges and priors. Note that, we have used a 5% prior for the sine of each of the active-sterile mixing angles while varying them over their entire range.  which show the effect of prior of θ 23 and different true choices of the standard oscillation parameters on the reconstruction of the CP phases.
In Fig. 3, we illustrate how efficiently the combination of LBL experiments can reconstruct the three CP phases δ 13 , δ 24 and δ 34 at 1σ C.L. in the three columns respectively given their true value lying anywhere in the whole parameter space of [−π, π]. In addition to the poorly measured 3+0 parameters (θ 23 , ∆m 2 31 ) and the active-sterile mixing angles (θ 14 , θ 24 , θ 34 ), in each panel we have also marginalised over the two other CP phases (∈ [−π, π]) not shown along the axes. The top (bottom) row depicts small (large) activesterile mixing with true θ 14 , θ 24 , θ 34 = 5.7 • , 5 • , 20 • (θ 14 , θ 24 , θ 34 = 10 • , 6 • , 25 • ). Note that, the latter choice of values represents the upper limits of the allowed active-sterile mixing (See Sec. 2). For each true value of the CP phases (∈ [−π, π]), the corresponding vertical width of the contours provide an estimate of the precision of reconstructing that true value. It is evident that in comparison to the reconstruction of the sterile phases, the standard Dirac phase δ 13 can be reconstructed much more efficiently by the LBL experiments. We note that the reconstruction of δ 13 does not noticeably depend upon the size of the activesterile mixing as assumed in the two rows. As the T2K data [14] suggests, if δ 13 indeed turns out to lie in the lower half plane of [−π, 0] with a best fit roughly around the maximal CPV (≈ −π/2), the future analyses with a combination of (DUNE + T2K + NOνA+ T2HK) with their projected runtime will be able to measure this phase in a narrow approx-  imate range of [−115 • , −75 • ] (at 1σ). The precision is marginally better if δ 13 turns out to be close to the CP conserving value (0). As far as the reconstruction of the sterile phases are concerned, δ 24 can be reconstructed with a better precision than δ 34 . For small activesterile mixing (top row of Fig. 3) and without considering the T2HK-projected data, δ 24 can be better reconstructed only in the lower half plane. But with T2HK-projected data, its reconstruction becomes much better throughout the entire parameter space. T2HK, due to its shorter baseline and much higher fiducial mass of its water cerenkov detector can offer very high statistics which helps to alleviate the degeneracy. On the other hand, the reconstruction of δ 34 is bad almost for the entire parameter space if T2HK is not considered. We have already seen in Fig. 1 that P (ν µ → ν e ) is most sensitive to δ 13 and least sensitive to δ 34 . This leads to a better precision in reconstructing δ 13 in comparison to the others. Fig. 4 illustrates the impact of different appearance and disappearance channels on the reconstruction of the CP phases (taken pairwise, like in Fig. 2) at 1σ C.L. at DUNE. The innermost red contours consisting of all the three channels for DUNE are the same as the red contours in the upper row of Fig. 2. We can clearly see the decrease in uncertainty in the measurement of the CP phases as we keep on adding ν µ → ν τ appearance and ν µ → ν µ disappearance channel to the ν µ → ν e appearance channel. The phase dependence of the ν µ → ν τ channel (see also Fig. 1) helps somewhat in this case. Note that, the improvement in result due to the addition of the ν τ appearance channel is not very significant and this is due to the very low statistics provided by this channel due to the difficulty in observing ν τ . On the other hand, though the phase dependence of ν µ → ν µ channel is small, it offers a very large number of events compared to the other two channels, thereby decreasing the uncertainty in reconstruction (see also, Appendix C). The improvement in reconstruction along the δ 13 direction is not significant with the addition of channels, underlining the pivotal role here played by the ν e appearance channel. In contrast, as observed from the first two panels of Fig. 4, the reconstruction of δ 24 and δ 34 is significantly better (more so in the former) particularly by the addition of ν µ disappearance channel. In the third panel, the role of ν µ disappearance is emphasized by the shrinking of the reconstruction contour in both δ 24 and δ 34 directions. Here we also note a slight improvement in reconstruction along the δ 34 axis (but not so much along δ 24 ) with the addition of ν τ appearance channel. In the next section we are going to estimate the correlations of the phases with the active-sterile mixing angles.

Reconstruction of the CP phase in correlation with active-sterile mixing angles
In Fig. 5 we show how efficiently the LBL experiments can reconstruct the CP phases in correlation with active-sterile mixing angles at 1σ C.L. (2 D.O.F.). The three phases (angles) are shown along the three rows (columns). In the second and third column the ranges of θ 24 and θ 34 are not shown beyond their allowed values. This results in the contours not being closed in these cases. As expected, the top row shows much less uncertainty along the test δ 13 direction for all the three angles involved. In all the panels, we note that combining NOνA and T2K data to DUNE gives mild improvement while a further addition of data projected by T2HK significantly improves the reconstruction. We have varied the test values of the angles upto their upper limit at 99% C.  The horizontal (vertical) span of the resulting contours indicate the uncertainties along the mixing angles (CP phases). We observe that the contours involving θ 24 and θ 34 do not close as long as we confine ourselves upto the allowed upper limit of θ 24 and θ 34 . It is interesting to note that even if θ 34 is allowed to vary in a larger range, the contours involving it show considerable uncertainty. This signifies the difficulty in constraining θ 34 in the experiments. It is mainly this reason which also hinders a good reconstruction of the associated phase δ 34 , as we have seen in Sec. 4. To have a quantitative idea about the potential to reconstruct the true values of the active-sterile mixing angle θ i4 (i = 1, 2, 3), we calculate how much the total horizontal span of each blue contour is. From this we estimate the maximum range of uncertainty (so as to obtain a conservative range) in reconstructing θ i4 and tabulate these below in Table 2

Impact of constraints on Neutrinoless Double Beta Decay
Before concluding, we include discussion of 3 + 1 framework on other observables, such as, neutrinoless double beta decay. The sterile neutrino, if a Majorana particle, gives non-zero contribution in the lepton number violating neutrinoless double beta decay (NDBD). In the presence of a non-zero θ 14 , the effective mass of NDBD process becomes [100], where α 2 , α 3 , α 4 are the relevant CP phases. Following the parametrisation given in Eq. 2.1, the relevant elements of the mixing matrix are as follows.
|U e1 | = c 12 c 13 c 14 , |U e2 | = s 12 c 13 c 14 , |U e3 | = s 13 c 14 , |U e4 | = sin θ 14 . (7. 2) The expression for half-life of 0νββ transition can be given as [101], where, In the above, m i is the mass of active neutrino and U ei is the PMNS mixing; whereas, θ ei is the mixing among the active and sterile and M i is the corresponding mass of heavy sterile. In the above, M ν and M N are the nuclear matrix elements (NME) for exchange of light and heavy neutrinos respectively. The values of NME and phase space factor G 0ν can The vertical light grey region is excluded by Planck data at 95% C.L. [102], while the horizontal dark region shows the 90% sensitivity from GERDA [103].
be find in Ref. [104]. The half-life of 0ν2β is can generally given as [105] 1 where j represents the number of light neutrino states and the additional heavy neutrino states. The parameters µ j and Θ ej represent the masses of the neutrino states and the mixing with SM neutrinos respectively. In the above, K 0ν = G 0ν (M N m p ) 2 and p 2 ≡ −m e m p M N Mν . Over the decade, many experiments on improving the lower limits on T 0ν 1/2 for 0νββ transition have been performed and best stringent bounds have been acquired from germanium-76, Xenon-136 and tellurium-130 [106] isotopes. The experiment GERDA-II, at 90% C.L. has obtained the corresponding lower limit as T 0ν 1/2 > 8.0 × 10 25 year for Ge (76) [107], whereas from KamLAND-Zen experiment Xe (136) , at 90% C.L. the lower limit on the half-life is obtained as T 0ν 1/2 > 1.07 × 10 26 year [108,109]. In presence of a sterile neutrino, the effective mass obtained from 3 + 0 framework will be significantly changed. Using Eq. 7.1, in Fig. 6, we show the effective mass for the standard 3 + 0 scenario, and for 3 + 1 scenario with the variation of the lightest mass. The left and right panels represent NH and IH, respectively. We consider a variation of ∆m 2 41 (for NH), and ∆m 2 43 (for IH) in between (1-10) eV 2 . The CP phases α 2 , α 3 and α 4 have been varied in between −π to π. Other oscillation parameters have been varied in their 3σ ranges as shown in Tab. 1 [22] as applicable for NH and IH. Overall we have considered 10 7 iterations at each value of the lightest mass in generating Fig. 6. The red and blue regions represent the variation of |m ef f | for the standard 3 + 0 scenario, and 3 + 1 scenario respectively. As can be seen from the figure the |m ef f | can be significantly large in the presence of a sterile neutrino with large mass value, and hence constrained from the experimental constraint. To make a connection with our oscillation analyses in the preceding sections, we shade the regions corresponding to ∆m 2 41 = 1.3 eV 2 with blue cross lines. Consequently we note a slight shrinking of the blue regions and it largely comes down to below the exclusion limit by Gerda. We note that for NH, there is a complete cancellation of m eff in the 3+0 case when the lightest mass approximately lies in the range of 10 −3 − 10 −2 eV. In 3+1 case, there is no complete cancellation in this region (due to the dominance of the m 4 -term in Eq. 7.1 in this range) and it only happens when the lightest mass becomes 10 −2 eV and beyond. For IH, the standard 3+0 case shows no total cancellation (dominant m 1 -term in Eq. 7.1), while there is a total cancellation in the 3+1 case in the range shown. Our results qualitatively agree with [100].

Summary and Conclusion
In this paper we have considered the presence of an eV-scale sterile neutrino (the so called 3+1 scenario which might turn out to be a possible resolution of the short baseline neutrino oscillation anomalies) and have analyzed how the present and future long baseline experiments T2K, NOvA, DUNE and T2HK can potentially probe the additional CP phases. We discuss how the three CP phases, namely δ 13 , δ 24 and δ 34 can individually affect the oscillation channels under consideration and appear in the probability expression. In light of the constraints on the active-sterile mixing from the global analysis, we estimate how the LBL experiments can probe the parameter spaces associated to the CP phases, by taking a pair of CP phases at a time. Though ν µ → ν e oscillation channel contributes the most in probing these parameter spaces, ν µ → ν µ and to a lesser extent ν µ → ν τ channel also help in exploring the δ 24 − δ 34 parameter space in particularly. By marginalizing over all other parameters we then show how the three individual CP phases can be reconstructed for all possible true values in the whole range of [−π, π]. We find that δ 24 and δ 34 cannot be reconstructed very efficiently by DUNE and also even after adding data from NOvA and T2K. But adding T2HK data removes much of the degeneracies and the uncertainties in reconstruction become much less. We found that if the active-sterile mixing angles turn out to be lying close to their current upper limits, the enhanced sensitivities to the associated phases make the reconstructions of δ 24 and δ 34 much better. In contrast the reconstruction of the standard CP phase δ 13 is much better even in presence of a light sterile neutrino and this conclusion is almost independent of the size of active-sterile mixing. We then analyze how efficiently the experiments can probe all the parameter spaces associated to one CP phase and one active-sterile mixing angle. It turns out that the parameter regions connected to the angle θ 14 can be probed relatively better that those related to the other two mixing angles. Finally, we briefly show how the relevant parameter spaces in 0νββ get modified in light of the active-sterile constraints used in this analysis.   In all the results so far, we have used a Gaussian prior of 3.5% on θ 23 . Since a true value lying in the higher octant of θ 23 has been used in the analyses, we also need to study the effect of using a non-Gaussian prior in order to properly consider both the octants when we marginalise over θ 23 . For this purpose, we take a prior uncertainty of 3.5% on sin 2 2θ 23 in the reconstruction of the CP phases and compare it with the case of considering a Gaussian prior in Fig. 7. We note that, taking the prior on sin 2 2θ 23 results in only a mild spread in the contours, especially along the axes of the sterile CP phases δ 24 or δ 34 . But this does not change our overall results qualitatively. For e.g., though the contours showing the degenerate solutions at DUNE (red) gets bigger with a prior on sin 2 2θ 23 , addition of data from other experiments still removes this degenerate fake solution as before.  The best fit values of the standard oscillation parameters used in our analyses (as mentioned in Tab. 1) were obtained under the 3+0 scenario. It has been shown in literature that some of these standard parameters are prone to significant changes in the presence of a light sterile neutrino. For e.g., the values of θ 23 , ∆m 2 31 , δ 13 can change so much that the issues of correct octant, mass hierarchy or CP Violation become very ambiguous (see for e.g., [63,64]). The reactor mixing angle θ 13 , on the other hand is quite robust even in the 3+1 scenario [110]. Though the marginalisation process during ∆χ 2 calculation partially takes care of the uncertainties of relevant parameters in the fit, a brief discussion regarding the effect of different true values of the standard parameters are in order. In Fig. 3, we had already illustrated our results for all possible true values of the CP violating phases.
Here we now show in Fig. 8, the impact of the true choice of lower octant (top row) and maximal mixing (bottom row) of θ 23 in probing the CP phases. This is similar to the bottom row of Fig. 7, but for two different possible true choices of θ 23 . We see that the contours remain qualitatively similar when different true octants of θ 23 are considered. For the case of maximal mixing, a slight elongation of the contours are observed, especially along the δ 24 and δ 34 directions. We have also checked that for inverted hierarchy (IH) of neutrino masses our results remain similar. This is due to the fact that even in presence of a sterile neutrino, the mass hierarchy sensitivity of DUNE, though deteriorates, remains high enough ( 5σ) [64].  Figure 9. Top row shows the bands for event spectra at DUNE due to variations of 3+1 CP phases for νµ → νe , νµ → νµ and νµ → ντ channels respectively. The band within blue (red) lines corresponds to the variation of δ24 (δ34) in the range [−π, π]. The bottom row shows the corresponding spread ∆ (square of difference between the two same coloured curves in the top row) over Ntrue (number of events when δ24, δ34 are fixed at zero). See text for details.
Here we make an attempt to understand the relative role of the three oscillation channels to probe the sterile CP phases δ 24 and δ 34 as observed in Fig. 4. For a simpler and more intuitive explanation we use the following Gaussian ∆χ 2 (the Poissonian definition in Eq. 3.1 reduces to the Gaussian version for sufficiently large number of events [3].): (C1) Here N true αβ (N test αβ ) is the true (test) set of events coming from the ν α → ν β oscillation channel. Note that, we have sketched here only the relevant statistical part of the ∆χ 2 and ignored the prior and systematics. In order to understand Fig. 4, we first estimate the band of N test αβ corresponding to individual variations of δ 24 ∈ [−π, π] and δ 34 ∈ [−π, π]. These spectra are illustrated for the three channels in Fig. 9 (top row). The ∆χ 2 for reconstructing δ 24 and δ 34 is governed by the spread/width of such event bands and we calculate the square of the vertical width at each energy bin and refer it by ∆. In accordance with Fig. 4 we then generate N true αβ corresponding to δ 24 = 0 and δ 34 = 0. Now, following Eq. C1, we plot ∆/N true αβ for the three channels in the bottom row of Fig. 9, both for δ 24 (blue) and δ 34 (red). This approximately gives us an idea of the relative contribution of the channels in the reconstruction of the sterile CP phases. As expected, the ν µ → ν e channel gives the dominant contribution and δ 24 -reconstruction is expected to be slightly better than δ 34 . ν µ → ν µ channel also plays a role and this mainly comes due to the large number of ν µ events. Finally ν µ → ν τ channel has a small but non-negligible contribution. Both ν µ → ν µ and ν µ → ν τ channels provide similar capability of reconstructing δ 24 and δ 34 . However, since the ν µ → ν µ channel events are also limited by systematic uncertainties, some remarks regarding the effect of systematics in Eq. C1 are in order. Extending our analysis in this appendix from the level of events to that of ∆χ 2 , we show in Fig. 10, a comparison between the ν µ → ν e and ν µ → ν µ channels in probing δ 24 or δ 34 . As in Fig. 9, here we also consider the true values of δ 24 or δ 34 to be 0 and then plot the ∆χ 2 as a function of the test values of δ 24 (blue) or δ 34 (red), after marginalising over all other relevant test parameters. The left (right) panel shows the contribution of the ν µ → ν e (ν µ → ν µ ) channel alone. As mentioned in Sec. 3, 2% (5%) systematics has been considered in the ν µ → ν e (ν µ → ν µ ) signal, following [91] (along with various background systematics as well). Due to this larger systematic uncertainty, we note from Fig. 10, that the relative role of ν µ → ν µ channel in probing the sterile phases with respect to ν µ → ν e channel further diminishes compared to that of Fig. 9. Nevertheless, even with such systematics the contribution of the ν µ → ν µ channel can still be roughly ∼ 15% − 20% of that of ν µ → ν e in the favourable region of the parameter space. Thus when all the oscillation channels are included to do a combined ∆χ 2 analysis, the interplay and complementarity among the channels improves the sensitivity to probe the CP phases. These findings are consistent with the full statistical analysis as discussed in Sec. 5.