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

CCTPU-PTC-21-02

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

Nishat Fiza, ∗ Mehedi Masud,

2, 3, † and Manimala Mitra

3, 4, ‡ Department of Physical Sciences, IISER Mohali,Knowledge City, SAS Nagar, Mohali - 140306, Punjab, India Center for Theoretical Physics of the Universe,Institute for Basic Science (IBS), Daejeon 34126, Korea Institute of Physics, Sachivalaya Marg, Bhubaneswar, Pin-751005, Odisha, India Homi Bhabha National Institute, BARC Training School Complex, Anushakti Nagar, Mumbai-400094, India (Dated: February 11, 2021)The various global analyses of available neutrino oscillation data indicate the presence of thestandard 3 + 0 neutrino oscillation picture. However, there are a few short baseline anomalies thatpoint to the possible existence of a fourth neutrino (with mass in the eV-scale), essentially sterilein nature. Should sterile neutrino exist in nature and its presence is not taken into considerationproperly in the analyses of neutrino data, the interference terms arising due to the additional CPphases 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) sterileneutrino and probe all the three CP phases ( δ , δ , δ ) in the context of the upcoming DeepUnderground Neutrino Experiment (DUNE) and also estimate how the results improve when datafrom NOvA, T2K and T2HK are added in the analysis. We illustrate the ∆ χ correlations of theCP phases among each other, and also with the three active-sterile mixing angles. Finally, we brieﬂyillustrate how the relevant parameter spaces in the context of neutrinoless double beta decay getmodiﬁed in light of the bounds in presence of a light sterile neutrino. ∗ [email protected] † [email protected] ‡ [email protected] a r X i v : . [ h e p - ph ] F e b I. INTRODUCTION

The successful discovery of the phenomena of neutrino oscillation [1, 2] has led to an impressive amount of researchwith an aim to establish the standard 3-ﬂavour (referred to as 3+0 hereafter) oscillation over a wide range of en-ergy ( E ) and neutrino propagation length ( L ). In the commonly used Pontecorvo-Maki-Nakagawa-Sakata (PMNS)parametrization [3] of the leptonic mixing matrix, this requires ﬁtting of the six standard oscillation parameters,namely three mixing angles ( θ , θ , θ ), two mass squared diﬀerences (∆ m = m − m , ∆ m = m − m ) and oneDirac CP phase ( δ ). A complete knowledge about the oscillation parameters will help to shed light on the followingstill unresolved issues in the neutrino sector: whether there exists CP violation (CPV) in the leptonic sector ( i.e. ,whether δ (cid:54) = 0 , π ), whether the neutrino mass eigenstates are arranged in normal ordering (NO, i.e. , ∆ m >

0) orinverted ordering (IO, i.e. , ∆ m < θ lies in the higher octant (HO, i.e. , θ > π/

4) or in the loweroctant (LO, i.e. , θ < π/ θ octant will help in understanding the origin of neutrinomass [6–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 Oﬀ-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 δ around π/ σ conﬁdence level(C.L.) irrespective of mass ordering. It also excludes CP conservation ( δ = 0 or π ) at 95% C.L. NO ν A, in its latestdataset [15] using both ν and ¯ ν running mode hint towards NO at 1 . σ C.L. and shows a weak preference for θ lying in HO at a C.L. of 1 . σ . Upcoming LBL experiments such as the Deep Underground Neutrino Experiment(DUNE) [16–18], Tokai to Hyper-Kamiokande (T2HK) [19], Tokai to Hyper-Kamiokande with a second detector inKorea (T2HKK) [20], European Spallation Source ν Super Beam (ESS ν SB) [21] as well as future data from T2K andNO ν A are expected to resolve the issues mentioned above with an unprecedented level of precision.Various global analyses [22–25], analyzing neutrino data from a diverse array of sources such as atmosphere, parti-cle accelerators, sun and nuclear reactors have consolidated the eﬀort to build up the 3+0 oscillation picture. Thereare however, a few short-baseline (SBL) anomalies [26–30] that hint towards the existence of oscillation governedby O (eV ) mass squared diﬀerence (∆ m = m − m ∼ eV ) that cannot be accommodated by the standard 3+0scenario. The eﬀort to explain the SBL anomalies (excess of electron-like events at low energy) has led to the modelswith possible presence of a sterile, fourth type of neutrino (this scenario is referred to as 3+1 hereafter), which canhave small mixing with the three active neutrinos. This fourth generation of neutrino, if exists, has to be sterile sincethe LEP experiment [31] restricts the number of active neutrino ﬂavours to three. In addition to the six standardoscillation parameters mentioned above, 3+1 scenario is parametrized by the mass squared diﬀerence ∆ m , threeactive-sterile mixing angles ( θ , θ , θ ), and two additional CP phases which we refer as δ and δ . Recently theobservation of a low energy electron-like event excess at a statistically signiﬁcant 4 . σ by the Mini Booster Neutrino Ex-periment (MiniBooNE) in its latest dataset [32, 33] has further boosted the search for sterile neutrinos. Various otherexisting and future facilities aim to search for sterile neutrinos with high precision using diﬀerent neutrino sources anddetection techniques. These facilities include IceCube [34], Karlsruhe Tritium Neutrino Experiment (KATRIN) [35],FermiLab’s Short Baseline Neutrino (SBN) programme [36], Neutrino Experiment for Oscillation at Short baseline(NEOS) [37], Short baseline neutrino Oscillations with a novel Lithium-6 composite scintillator Detector(SoLid) [38],Neutrino-4 [39] , Precision Reactor Oscillation and SPECTrum Experiment (PROSPECT) [46], Sterile Reactor Neu-trino Oscillations (STEREO) [47, 48], Detector of the reactor AntiNeutrino based on Solid Scintillator (DANSS) [49],J-PARC Sterile Neutrino Search at J-PARC Spallation Neutron Source(JSNS ) [50, 51]. Experiments looking forCoherent Elastic Neutrino-Nucleus Scattering (CE ν NS) [52] can also act as very useful probes of sterile neutrinos, asshown by the authors of [53].The eﬀect of sterile neutrino, should it exist with ∆ m ∼ , is most pronounced around L/E ∼ L/E ∼

500 km/GeV at the far detectors (FD), the high frequency induced by∆ m ∼ gets averaged out due to the ﬁnite energy resolution of the detector. Nevertheless, it has beenshown [54, 56–77] that even at the FDs of these LBL experiments, the interference eﬀects provided by the additionalCP phases play very signiﬁcant roles in spoiling the sensitivities to the crucial issues of CPV, MH and θ octant .For instance, as shown in [62], the CPV sensitivity becomes a wide band whose width depends on the unknownmagnitudes of the sterile phases δ and δ , - leading to serious confusion in interpreting the results as CP violationor CP conservation. The constraints on the active-sterile mixing angles θ i ( i = 1 , ,

3) do of course reduce such Neutrino-4 has recently claimed to observe active-sterile oscillation at 3 . σ around the vicinity of ∆ m ≈ and θ ≈ ◦ byanalyzing the its reactor antineutrino data accumulated since 2016 [40–42]. Though it has later been argued that considering a moreappropriate log-likelihood distribution [43, 44] or correct energy resolution [45], the statistical signiﬁcance of the active-sterile observationresults in reactor neutrino experiments actually come down to a much lower value. For a recent comprehensive status report of the impact of a light sterile neutrino in probing these issues at LBL see, for e.g. , [78–80]and the references therein. ambiguities in the interpretation of the results to some extent. But, a clear idea about how the LBL experiments areable to measure the sterile phases δ and δ given the constraints on the active-sterile mixing angles, will certainlyminimise the obfuscation when their data are analyzed with a view to tackle the unresolved physics issues. The issue ofmeasurement of one sterile phase ( δ or δ , depending on the parametrization used) has been addressed to some extentin literature. [62] analyzes how the sterile phases impact the standard CPV and mass ordering measurements and alsoprobe the joint parameter space δ − δ for DUNE. The authors of [60] use a slightly diﬀerent parameterization anddiscuss the CPV arising from the individual sterile CP phases and probe the joint parameter space δ − δ at DUNE.They further extend their analysis by combining simulated data from other LBL experiments such as T2HK [70] andESS ν SB [74]. Exploring the parameter space of δ − δ has also been addressed in [69] for DUNE, T2HK, T2HKKand their combinations. The authors of [69] further illustrate how the individual phase δ can be measured by theseexperiments at various C.L. by assuming four possible true values (0 , ± π/ , π ). [77] shows how the recent data fromT2K and NO ν A can help to probe the parameter space of δ − δ . Most recently, the authors of [81] estimate howthe diﬀerence ( δ − δ ) in the sterile phases can impact in constraining the standard δ − θ parameter space.It is noteworthy that the sterile CP phase δ and its correlation with the other phases has been little addressedin literature. In the present manuscript we tackle the very relevant issue of estimating the capability to reconstructall three CP phases ( δ , δ , δ ), taking into consideration their ∆ χ correlations with each other and also withthe active-sterile mixing angles ( θ , θ , θ ). We carry out this exercise in the context of DUNE and illustrate theimprovement when combined with T2K, NO ν A (both these currently running experiments are simulated upto theirpresent exposure) and T2HK. Apart from studying these CP phases in detail, another crucial aspect in which ouranalysis diﬀers from the existing studies mentioned above is that we have taken into consideration the current 3 σ hintof CP violation and the corresponding exclusion region of the CP phase δ by T2K data [14]. Moreover, in additionto illustrating how the 2-d parameter spaces for the CP phases can be probed, we also analyze how the individualCP phases can be reconstructed (after marginalizing all other relevant parameters) by the experiments, irrespectiveof the actual value they might have in nature. Additionally we have also considered the ν µ → ν τ channel (in additionto ν µ → ν e and ν µ → ν µ ) in our study and estimated the capability of the projected data to measure all three CPphases in detail. This enables us to probe the parameter spaces associated to θ and δ with better sensitivity.The present manuscript is organised as follows. In Sec. II we give a brief account of the constraints on the sterileneutrino parameters and how the CP phases aﬀect the relevant probabilities. In Sec. III we describe the methodologyof our statistical analyses. In Sec. IV we discuss the ∆ χ correlations among various CP phases, taking a pair of phasesat a time. We also discuss the potential to reconstruct the three CP phases for all possible true values by performingsimulations of DUNE, T2K, NOvA and T2HK. The role of individual oscillation channels in such reconstructions areanalyzed in Sec. V. We estimate the ∆ χ correlations among all the CP phases and the active-sterile mixing anglesin Sec. VI. Finally in Sec. VII we brieﬂy discuss how the relevant parameter spaces associated to Neutrinoless DoubleBeta Decay gets modiﬁed in light of the constraints on one eV scale sterile neutrino, followed by conclusion. II. BASICS

We ﬁrst discuss the oscillation probabilities for the three channels ( P ( ν µ → ν e ), P ( ν µ → ν µ ) and P ( ν µ → ν τ )) in3+1 scenario. Since the expressions become immensely complicated in matter, we show them in vacuum and thesewill act as useful templates for gaining the physics insights in explaining our subsequent sensitivity results. For themixing matrix we follow the parametrization scheme adopted in [54]: U = R ( θ , δ ) R ( θ , δ ) R ( θ ) R ( θ ) R ( θ , δ ) R ( θ ) , (1)where R ( θ ij , δ ij ) is a rotation in the ij − th plane with an associated phase δ ij such that, for e.g. , R ( θ , δ ) = θ e − iδ sin θ − e iδ sin θ cos θ . (2)Now we discuss brieﬂy on the allowed ranges of the sterile sector parameters, as estimated in great detail in the globalanalysis of various neutrino data [82]. | U e | is bounded by ν e and ¯ ν e disappearance searches and is equal to sin θ .The combined atmospheric neutrino data from IceCube, DeepCore and SK, at 99% C.L. (2 DOF) put the bound | U e | (cid:46) .

1. This implies θ (cid:46) . ◦ . It can be seen from all ν e and ¯ ν e searches that the best ﬁt value of | U e | isapproximately equal to 0 .

01, which gives θ ≈ . ◦ . The data from ν µ and ¯ ν µ disappearance searches put the following99% C.L. (2 D.O.F.) constraints: | U µ | (cid:46) .

01 and | U τ | (cid:46) .

17. Since, in our parametrization | U µ | = cos θ sin θ and | U τ | = cos θ cos θ sin θ , the corresponding bounds on θ and θ can easily be translated to θ (cid:46) . ◦ and θ (cid:46) . ◦ . The allowed values of ∆ m roughly lie in the range of 1 −

10 eV and we consider it to be 1 . in the present work as per the global analysis .Using the standard approach for deriving oscillation probability, we obtain for the ν α → ν β ( α, β = e, µ, τ, s and α (cid:54) = β ) transition probability, P αβ = 4 | U α U β | × . − U α U ∗ β U ∗ α U β ) sin ∆ + 2Im( U α U ∗ β U ∗ α U β ) sin 2∆ − U α U ∗ β U ∗ α U β ) sin ∆ + 2Im( U α U ∗ β U ∗ α U β ) sin 2∆ − U α U ∗ β U ∗ α U β ) sin ∆ + 2Im( U α U ∗ β U ∗ α U β ) sin 2∆ , (3)where ∆ ij = ∆ m ij L E . In deriving Eq. 3, we have used the unitarity of the 4 × U and applied theusual assumptions that the term containing mass square splitting between m and m i ( i = 1 , , i.e. , sin ∆ i andsin 2∆ i average out to 0.5 and 0 respectively at long baseline ( i = 1 , , i.e. , neglecting the oscillation eﬀects due to ∆ m ) and arrive at the following simpliﬁed expressionfor the dominant channel ν µ → ν e . P νµe ≈

12 sin θ νµe + ( a sin θ νµe −

14 sin θ sin θ νµe ) sin ∆ + cos( δ + δ ) a sin 2 θ νµe sin 2 θ νµe cos 2 θ sin ∆ + sin( δ ) ba sin 2 θ νµe sin θ sin 2∆ + 12 sin( δ + δ ) a sin 2 θ νµe sin 2 θ νµe sin 2∆ , (4)where we have followed the convention of [54] for the following quantities.sin 2 θ νµe = sin 2 θ sin θ , b = cos θ cos θ sin θ ,sin 2 θ νµe = sin 2 θ sin θ , a = cos θ cos θ .Eq. 4 tells us that in vacuum P ( ν µ → ν e ) is sensitive to both δ and δ , but not to δ . As explained in [54], smalldependence on δ creeps in when matter eﬀect is taken into account. The expressions for the less dominant channels P ( ν µ → ν µ ) and P ( ν µ → ν τ ) can similarly be derived from Eq. 3, but the expressions are quite lengthy. Interestedreaders can see [73, 83] for those probability expressions. Parameter Best-ﬁt-value 3 σ interval 1 σ uncertainty θ [Deg.] 34.3 31.4 - 37.4 2.9% θ (NH) [Deg.] 8.58 8.16 - 8.94 1.5% θ (IH) [Deg.] 8.63 8.21 - 8.99 1.5% θ (NH) [Deg.] 48.8 41.63 - 51.32 3.5% θ (IH) [Deg.] 48.8 41.88 - 51.30 3.5%∆ m [eV ] 7 . × − [6.94 - 8.14] × − m (NH) [eV ] +2 . × − [2.46 - 2.65] × − m (IH) [eV ] − . × − -[2.37 - 2.55] × − δ (NH) [Rad.] − . π [ − π, ∪ [0 . π, π ] − δ (IH) [Rad.] − . π [ − . π, − . π ] − TABLE I. Standard oscillation parameters and their uncertainties used in our study. The values were taken from the global ﬁtanalysis in [22]. 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 ) / x u + x l )% [17]. Using the widely used General Long Baseline Experiment Simulator (GLoBES) [84, 85] and the relevant plugin snu.c [86, 87] for implementing sterile neutrinos, we now illustrate how the probabilities for diﬀerent oscillation In our notation, m , , are the masses of the three active neutrinos, and m denotes the mass of the sterile neutrino. Also, ∆ m ij = m i − m j . δ ∈ [-π,π]δ ∈ [-π,π]δ ∈ [-π,π] P μ e P μ μ P μ τ FIG. 1.

We show the probability bands due to individual variation of the CP phases δ (grey), δ (blue) and δ (red) in the wholerange 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 magniﬁed versions of the rectangular regions indicated. The three active-sterile mixingangles were taken as θ = 10 ◦ , θ = 6 ◦ , θ = 25 ◦ . The values of the rest of the oscillation parameters were taken from Tab. I. Normalhierarchy was assumed for generating this plot. channels depends on the three CP phases ( δ , δ , δ ) individually at the DUNE baseline of 1300 km. In Fig. 1, weplot the bands for P µe , P µµ , P µτ (in the three panels respectively) due to individual variation of δ , δ , δ , each beingvaried in the range [ − π, π ]. The grey band shows the variation of the standard Dirac CP phase δ . The variation dueto δ ( δ ) is shown with the blue (red) band. For the variation of each CP phase, the other two phases are kept ﬁxed: δ is kept ﬁxed at the bf value of − . π while δ , δ were considered to be zero. The sterile phases are associatedwith the active-sterile mixing angles. We have a slightly high active-sterile mixing: θ = 10 ◦ , θ = 6 ◦ , θ = 25 ◦ .As expected , the CP phases have larger impact on appearance channels rather than the disappearance channels. ν µ → ν e channel is most aﬀected by the variation of δ . Among the sterile CP phases, δ has a visibly greaterimpact on P ( ν µ → ν e ) than that of δ . It is interesting to note this feature especially in light of the fact that thevalue of the active-sterile mixing angle θ (taken as 6 ◦ in Fig. 1) is almost 5 times smaller than θ (25 ◦ ). We alsoobserve that with increase in energy the eﬀect of δ on P ( ν µ → ν e ) further reduces. P ( ν µ → ν τ ) is less prone tovariation of the CP phases. Still it shows slight variation for all three CP phases with δ and δ having almostsimilar eﬀects but less than that of δ . P ( ν µ → ν µ ) on the other hand, shows almost negligible variation with theCP phases. After a brief description of the simulation procedure we will estimate the capability of LBL experimentsto reconstruct these CP phases. III. SIMULATION DETAILS

We simulate the long baseline neutrino experiments DUNE, NOvA, T2K and T2HK using GLoBES [84, 85]. DUNEis a 1300 km long baseline experiment employing a liquid argon far detector (FD) of 40 kt ﬁducial mass with a beamof power 1.07 MW and running 3.5 years each on ν and ¯ ν mode (resulting in a total exposure of roughly 300 kt.MW.yrcorresponding to 1 . × protons on target or POT). We have used the oﬃcial conﬁguration ﬁles [88] provided bythe DUNE collaboration for its simulation. Electron neutrino appearance signals (CC), muon neutrino disappearancesignals (CC), as well as neutral current (NC) backgrounds and tau neutrino appearance backgrounds (along with thecorresponding systematics/eﬃciencies etc. ) are already included in the conﬁguration ﬁles. In the present analysis, wehave additionally incorporated tau neutrino appearance signal (CC) including ﬂavour misidentiﬁcation background(due to the leptonic decay of tau lepton with a branching fraction ∼ ∼ . × (12 . × ) POT in ν (¯ ν ) mode. T2K is a 295 km experiment with a 22.5 kt water cherencov FD. For T2Ksimulation we use the inputs from [14, 90]. We have used a beam of 515 kW and simulating 1 . × (1 . × )POT in ν (¯ ν ) mode. T2HK is an upgraded version of T2K with a higher beam of 1.3 MW and a much bigger ﬁducial In Eq. 3 for α = β , the imaginary part in the RHS vanishes, diminishing the eﬀect of the CP phases. mass of 187 kt of its water cherencov FD. For T2HK we simulate a total of 2 . × POT in 1:3 ratio of ν and ¯ ν mode (with inputs taken from [19, 91]). Note that, for the future experiments DUNE and T2HK we have used thefull expected exposure, while for the currently running experiments T2K and NOvA we have simulated upto theircurrent exposure.To estimate the capability of LBL experiments to reconstruct the CP phases, we carry out a ∆ χ analysis usingGLoBES and the relevant plugin snu.c [86, 87] for implementing sterile neutrinos. In order to gain insight, let usexamine the analytical form of the ∆ χ ,∆ χ ( p true ) = Min p test ,η (cid:34) mode (cid:88) k channel (cid:88) j bin (cid:88) i (cid:40) N test ijk ( p test ; η ) − N true ijk ( p true ) + N true ijk ( p true ) ln N true ijk ( p true ) N test ijk ( p test ; η ) (cid:41) + (cid:88) l ( p true l − p test l ) σ p l + (cid:88) m η m σ η m (cid:35) , (5)where N true ( N test ) is the true ( test ) set of events, while p true ( p test ) 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. Theterm ( N test − N true ) takes into account the algebraic diﬀerence while the log-term inside the curly braces considersthe fractional diﬀerence between the test and true sets of events. The term summed over i, j, k and written insidethe curly braces is the statistical part of ∆ χ . σ p l is the uncertainty in the prior measurement of the l -th oscillationparameter p l . The values of the true or best-ﬁt oscillation parameters and their uncertainties as used in the presentanalysis are tabulated in Table I. 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 ∆ χ . This way of treating the systematics in the ∆ χ calculation is known as the method of pulls [92–95]. For DUNE, the ν e and ¯ ν e signal modes have a normalizationuncertainties of 2% each, whereas the ν µ and ¯ ν µ signals have a normalization uncertainty of 5% each. The ν τ and T e s t δ [ D e g . ] −180−135−90−4504590135180 T e s t δ [ D e g . ] −180−135−90−4504590135180 T r u e ( δ , δ , δ ) = ( - . π , , ) DUNE DUNE + T2K + NOvA DUNE + T2K + NOvA + HK T e s t δ [ D e g . ] −180−135−90−4504590135180 T e s t δ [ D e g . ] −180−135−90−4504590135180 Test δ [Deg.]−180 −90 −45 0 45 90 135 180 T e s t δ [ D e g . ] −180−135−90−4504590135180 Test δ [Deg.]−180 −90 −45 0 45 90 135 180 T r u e ( δ , δ , δ ) = ( - . π , - . π , - . π ) T e s t δ [ D e g . ] −180−135−90−4504590135180 Test δ [Deg.]−180 −90 −45 0 45 90 135 180 FIG. 2.

Reconstruction of the CP phases, taken pairwise at a time, at a C.L. of 1 σ (2 D.O.F.) at DUNE (red), DUNE + T2K + NO ν A(green), and DUNE + T2K + NO ν A + T2HK (blue). The top (bottom) row corresponds to the choice δ = δ = 0 ( − π/ δ is ﬁxed at − . π . The black dot indicates the true values assumed. Other oscillation parameterswere taken from Tab. I. This is the

Poissonian deﬁnition of ∆ χ , which in the limit of large sample size, reduces to the Gaussian form. ¯ ν τ signals have a normalization uncertainties of 20% each. The background normalization uncertainties vary from5% −

20% and include correlations among various sources of background (coming from beam ν e / ¯ ν e contamination,ﬂavour misidentiﬁcation, neutral current and ν τ ). The ﬁnal estimate of ∆ χ obtained after a marginalization overthe 3 σ range of test parameters ( p test ) and the set of systematics ( η ) is a function of the true values of the oscillationparameters. Technically this ∆ χ is the frequentist method of hypotheses testing [93, 96]. IV. 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 δ -test δ ), (test δ -test δ ) and (test δ -test δ ) in the three columns respectively. The red contour shows the reconstruction capabilityof 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 δ and δ . Clearly the reconstruction of δ is better (as evidenced bythe narrowness of the contours along the test δ axis) compared to the other two phases. This happens since theassociated mixing angle θ has been measured very precisely unlike the corresponding active-sterile mixing angles θ and θ (see Table I and Sec. II). For the maximally CP violating choices of true δ and true δ , their reconstructiongets better at the cost of small degeneracies appearing for DUNE around test δ ≈ ◦ , δ ≈ ◦ . Adding datafrom other experiments lifts these degeneracies. T e s t δ [ D e g . ] −180−135−90−4504590135180 T e s t δ [ D e g . ] −180−135−90−4504590135180 T r u e ( θ , θ , θ ) = ( . o , o , o ) T e s t δ [ D e g . ] DUNE DUNE + T2K + NOvA DUNE + T2K + NOvA + HK−180−135−90−4504590135180 T e s t δ [ D e g . ] −180−135−90−4504590135 True δ [Deg.]−180 −90 −45 0 45 90 135 180 T e s t δ [ D e g . ] −180−135−90−4504590135 True δ [Deg.]−180 −90 −45 0 45 90 135 180 T r u e ( θ , θ , θ ) = ( o , o , o ) T e s t δ [ D e g . ] −180−135−90−4504590135 True δ [Deg.]−180 −90 −45 0 45 90 135 180 FIG. 3.

Reconstruction of the CP phases δ , δ and δ , for all the choices of their true values in [ − π, π ] at a C.L. of 1 σ (1 D.O.F.) atDUNE (red), DUNE + T2K + NO ν A (green), and DUNE + T2K + NO ν A + T2HK (blue). The top (bottom) row corresponds to thechoice of the active-sterile mixing angles as true θ , θ , θ = 5 . ◦ , ◦ , ◦ (10 ◦ , ◦ , ◦ ). The true values of the phases not shown in apanel are ﬁxed at: δ , δ and δ = − . π, σ allowed valuesmeasured by T2K [14]. In Fig. 3, we illustrate how eﬃciently the combination of LBL experiments can reconstruct the three CP phases δ , δ and δ at 1 σ C.L. in the three columns respectively given their true value lying anywhere in the wholeparameter space of [ − π, π ]. In addition to the poorly measured 3+0 parameters ( θ , ∆ m ) and the active-sterilemixing angles ( θ , θ , θ ), in each panel we have also marginalised over the two other CP phases ( ∈ [ − π, π ]) notshown along the axes. The top (bottom) row depicts small (large) active-sterile mixing with true θ , θ , θ =5 . ◦ , ◦ , ◦ ( θ , θ , θ = 10 ◦ , ◦ , ◦ ). Note that, the latter choice of values represents the upper limits of theallowed active-sterile mixing (See Sec. II). For each true value of the CP phases ( ∈ [ − π, π ]), the correspondingvertical width of the contours provide an estimate of the precision of reconstructing that true value. It is evidentthat in comparison to the reconstruction of the sterile phases, the standard Dirac phase δ can be reconstructedmuch more eﬃciently by the LBL experiments. We note that the reconstruction of δ does not noticeably dependupon the size of the active-sterile mixing as assumed in the two rows. As the T2K data [14] suggests, if δ indeedturns out to lie in the lower half plane of [ − π,

0] with a best ﬁt roughly around the maximal CPV ( ≈ − π/ ν A+ T2HK) with their projected runtime will be able tomeasure this phase in a narrow approximate range of [ − ◦ , − ◦ ] (at 1 σ ). The precision is marginally better if δ turns out to be close to the CP conserving value (0). As far as the reconstruction of the sterile phases are concerned, δ can be reconstructed with a better precision than δ . For small active-sterile mixing (top row of Fig. 3) andwithout considering the T2HK-projected data, δ can be better reconstructed only in the lower half plane. But withT2HK-projected data, its reconstruction becomes much better throughout the entire parameter space. T2HK, due toits shorter baseline and much higher ﬁducial mass of its water cerenkov detector can oﬀer very high statistics whichhelps to alleviate the degeneracy. On the other hand, the reconstruction of δ is bad almost for the entire parameterspace if T2HK is not considered. For e.g. , if the true value of δ in nature turns out to be around [ − ◦ , − ◦ ],the range of values reconstructed at 1 σ by the combination of DUNE + T2K + NO ν A can be anywhere between − ◦ to 180 ◦ . Large active-sterile mixing (bottom row of Fig. 3) signiﬁcantly improves the sensitivities of the LBLexperiments to the sterile CP phases, resulting in signiﬁcant improvement of the reconstruction of δ and δ . Wehave already seen in Fig. 1 that P ( ν µ → ν e ) is most sensitive to δ and least sensitive to δ . This leads to a betterprecision in reconstructing δ in comparison to the others. V. ROLE OF INDIVIDUAL CHANNELS IN RECONSTRUCTION

Fig. 4 illustrates the impact of diﬀerent appearance and disappearance channels on the reconstruction of the CPphases (taken pairwise, like in Fig. 2) at 1 σ C.L. at DUNE. The innermost red contours consisting of all the threechannels for DUNE are the same as the red contours in the upper row of Fig. 2. We can clearly see the decreasein 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 alsoFig. 1) helps somewhat in this case. Note that, the improvement in result due to the addition of the ν τ appearancechannel is not very signiﬁcant and this is due to the very low statistics provided by this channel due to the diﬃcultyin observing ν τ . On the other hand, though the phase dependence of ν µ → ν µ channel is small, it oﬀers a very largenumber of events compared to the other two channels, thereby decreasing the uncertainty in reconstruction. Theimprovement in reconstruction along the δ direction is not signiﬁcant with the addition of channels, underliningthe pivotal role here played by the ν e appearance channel. In contrast, as observed from the ﬁrst two panels of Fig.4, the reconstruction of δ and δ is signiﬁcantly 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 δ [ D e g . ] −150−100−50050100150 δ [Deg.]−150 −100 −50 0 50 100 150 δ [ D e g . ] −150−100−50050100150 δ [Deg.]−150 −100 −50 0 50 100 150ν e app. channel ν e + ν τ app. channel All channels δ [ D e g . ] −150−100−50050100150 δ [Deg.]−150 −100 −50 0 50 100 150 FIG. 4.

Reconstruction of the CP phases, taken pairwise at a time for three diﬀerent channels at DUNE at a C.L. of 1 σ (2 D.O.F.). for ν e appearance channel (cyan), ν e + ν τ appearance channel (grey) and all channels (red). The three panels depict the test values of thethree CP phases. The true values for the CP phases were assumed as δ , δ , δ = − . π, , θ , θ , θ = 5 . ◦ , ◦ , ◦ . The black dot indicates the true values assumed. reconstruction contour in both δ and δ directions. Here we also note a slight improvement in reconstruction alongthe δ axis (but not so much along δ ) with the addition of ν τ appearance channel. In the next section we are goingto estimate the correlations of the phases with the active-sterile mixing angles. VI. RECONSTRUCTION OF THE CP PHASE IN CORRELATION WITH ACTIVE-STERILE MIXING ANGLES

In Fig. 5 we show how eﬃciently 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). Inthe second and third column the ranges of θ and θ are not shown beyond their allowed values. This results in thecontours not being closed in these cases. As expected, the top row shows much less uncertainty along the test δ direction for all the three angles involved. In all the panels, we note that combining NO ν A and T2K data to DUNEgives mild improvement while a further addition of data projected by T2HK signiﬁcantly improves the reconstruction.We have varied the test values of the angles upto their upper limit at 99% C.L.: θ (cid:46) . ◦ , θ (cid:46) . ◦ , θ (cid:46) . ◦ ,and the true values to be reconstructed were θ = 5 . ◦ , θ = 5 ◦ , θ = 20 ◦ (as discussed in Sec. II). T e s t δ [ D e g . ] −180−135−90−4504590135180 DUNEDUNE + T2K + NOvADUNE + T2K + NOvA + HK T e s t δ [ D e g . ] −180−135−90−4504590135 T e s t δ [ D e g . ] −180−135−90−4504590135 Test θ [Deg.]0 2 4 6 8 10 12 14 Test θ [Deg.]0 1 2 3 4 5 Test θ [Deg.]0 5 10 15 20 25 FIG. 5.

Reconstruction of the CP phases in correlation with the mixing angles θ , θ , θ at a C.L. of 1 σ (2 D.O.F.) at DUNE (red),DUNE + T2K + NO ν A (green), and DUNE + T2K + NO ν A + T2HK (blue). The three columns (rows) depict the test values of thethree active-sterile mixing angles (three CP phases). The true values for the CP phases were assumed as δ , δ , δ = − . π, , θ , θ , θ = 5 . ◦ , ◦ , ◦ . The black dot indicates the true values assumed. The horizontal (vertical) span of the resulting contours indicate the uncertainties along the mixing angles (CPphases). We observe that the contours involving θ and θ do not close as long as we conﬁne ourselves upto theallowed upper limit of θ and θ . It is interesting to note that even if θ is allowed to vary in a larger range, thecontours involving it show considerable uncertainty. This signiﬁes the diﬃculty in constraining θ in the experiments.It is mainly this reason which also hinders a good reconstruction of the associated phase δ , as we have seen in Sec.IV. To have a quantitative idea about the potential to reconstruct the true values of the active-sterile mixing angle θ i ( i = 1 , , θ i and tabulate these belowin Table II.0 Angle Value to be reconstructed Reconstructed range[Degree] [Degree] θ . (cid:46) θ (cid:46) . θ . (cid:46) θ (cid:46) θ

20 9 . (cid:46) θ (cid:46) θ , θ , θ as estimated from Fig. 5 for DUNE + T2K + NOvA + HK. VII. IMPACT OF CONSTRAINTS ON NEUTRINOLESS DOUBLE BETA DECAY

Before concluding, we include discussion of 3 + 1 framework on other observables, such as, neutrinoless doublebeta decay. The sterile neutrino, if a Majorana particle, gives non-zero contribution in the lepton number violatingneutrinoless double beta decay (NDBD). In the presence of a non-zero θ , the eﬀective mass of NDBD processbecomes, m eﬀ = m U e + m U e e iα + m U e e iα + m U e , (6)where α and α are the majorana phases. Following the parametrisation given in Eq. 1, the relevant elements of themixing matrix are as follows. U e = c c c , U e = s c c , U e = s c e − iδ , U e = sin θ . (7)The expression for half-life of 0 νββ transition can be given as [97],1 T ν / = G ν | M ν η ν + M N η N | , (8)where, η ν = U e i m i m e , η N = V e i m p M i . (9)In the above, m i is the mass of active neutrino and U ei is the PMNS mixing; whereas, θ ei is the mixing among theactive and sterile and M i is the corresponding mass of heavy sterile. In the above, M ν and M N are the nuclear matrixelements (NME) for exchange of light and heavy neutrinos respectively. The values of NME and phase space factor G ν can be ﬁnd in Ref. [100]. The half-life of 0 ν β is can generally given as [101]1 T / = K ν (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Θ ej µ j (cid:104) p (cid:105) − µ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (10)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 ν = G ν ( M N m p ) and (cid:104) p (cid:105) ≡ − m e m p M N M ν . Over the decade, many experiments on improving the lower limitson T ν / for 0 νββ transition have been performed and best stringent bounds have been acquired from germanium-76,Xenon-136 and tellurium-130 [102] isotopes. The experiment GERDA-II, at 90% C.L. has obtained the correspondinglower limit as T ν / > . × year for Ge (76) [103], whereas from KamLAND-Zen experiment Xe (136) , at 90% C.L.the lower limit on the half-life is obtained as T ν / > . × year [104, 105].In presence of a sterile neutrino, the eﬀective mass obtained from 3 + 0 framework will be signiﬁcantly changed. InFig. 6, we show the eﬀective mass for the standard 3 + 0 scenario, and for 3 + 1 scenario with the variation of thelightest mass. The left and right panels represent NH and IH, respectively. We consider a variation of ∆ m (for NH),and ∆ m (for IH) in between (1-10) eV , and the mixing angle θ in between 0 ◦ -18 ◦ . The Majorana phases α and α have been varied in between − π to π . The Dirac CP phase δ and other oscillation parameters θ , θ , θ and∆ m have been varied in their 3 σ ranges [22] as applicable for NH and IH The red and blue regions represent thevariation of | m eff | for the standard 3 + 0 scenario, and 3 + 1 scenario respectively. As can be seen from the ﬁgurethe | m eff | can be signiﬁcantly large in the presence of a sterile neutrino with large mass value, and hence constrainedfrom the experimental constraint. To make a connection with our oscillation analyses in the preceding sections, weshade the regions corresponding to ∆ m = 1 . with blue cross lines. Consequently we note a slight shrinking ofthe blue regions and it largely comes down to below the exclusion limit by Gerda.1 NH P l a n c k % Gerda 90% | m e ff | [ e V ] −5 −4 −3 [eV]10 −5 −4 −3 P l a n c k % SMSterile m [eV]10 −5 −4 −3 FIG. 6.

Eﬀective mass of NDBD versus the smallest neutrino mass. The left (right) panel shows the case of NH (IH). The red (blue)region correspond to 3+0 (3+1) scenario, labelled as SM (Sterile). The region hatched with blue cross lines represents ∆ m = 1 . .The vertical light grey region is excluded by Planck data at 95% C.L. [98], while the horizontal dark region shows the 90% sensitivity fromGERDA [99]. VIII. SUMMARY AND CONCLUSION

In this paper we have considered the presence of an eV-scale sterile neutrino (the so called 3+1 scenario which mightturn out to be a possible resolution of the short baseline neutrino oscillation anomalies) and have analyzed how thepresent and future long baseline experiments T2K, NOvA, DUNE and T2HK can potentially probe the additional CPphases. We discuss how the three CP phases, namely δ , δ and δ can individually aﬀect the oscillation channelsunder consideration and appear in the probability expression. In light of the constraints on the active-sterile mixingfrom the global analysis, we estimate how the LBL experiments can probe the parameter spaces associated to the CPphases, by taking a pair of CP phases at a time. Though ν µ → ν e oscillation channel contributes the most in probingthese parameter spaces, ν µ → ν µ and to a lesser extent ν µ → ν τ channel also help in exploring the δ − δ parameterspace in particularly. By marginalizing over all other parameters we then show how the three individual CP phasescan be reconstructed for all possible true values in the whole range of [ − π, π ]. We ﬁnd that δ and δ cannot bereconstructed very eﬃciently by DUNE and also even after adding data from NOvA and T2K. But adding T2HK dataremoves much of the degeneracies and the uncertainties in reconstruction become much less. We found that if theactive-sterile mixing angles turn out to be lying close to their current upper limits, the enhanced sensitivities to theassociated phases make the reconstructions of δ and δ much better. In contrast the reconstruction of the standardCP phase δ is much better even in presence of a light sterile neutrino and this conclusion is almost independent ofthe size of active-sterile mixing. We then analyze how eﬃciently the experiments can probe all the parameter spacesassociated to one CP phase and one active-sterile mixing angle. It turns out that the parameter regions connectedto the angle θ can be probed relatively better that those related to the other two mixing angles. Finally, we brieﬂyshow how the relevant parameter spaces in 0 νββ get modiﬁed in light of the active-sterile constraints used in thisanalysis. ACKNOWLEDGMENTS

N.F. is grateful for a visit to IOP, Bhubaneswar where this project was initialized. N.F. acknowledges the sup-port of Dr. Satyajit Jena of IISER Mohali. M. Masud acknowledges Dr. S.K. Agarwalla of IOP, Bhubaneswarfor providing the ﬁnancial support from the Indian National Science Academy (INSA) Young Scientist Project[INSA/SP/YS/2019/269]. M. Masud is supported by IBS under the project code IBS-R018-D1. M. Mitra acknowl-edges the ﬁnancial support from DST INSPIRE Faculty research grant (IFA-14-PH-99), and thanks Indo-French2Centre for the Promotion of Advanced Research for the support (grant no: 6304-2). [1] Y. Fukuda et al. (Super-Kamiokande Collaboration), Phys.Rev.Lett. , 1562 (1998), hep-ex/9807003.[2] Q. Ahmad et al. (SNO), Phys. Rev. Lett. , 011301 (2002), nucl-ex/0204008.[3] P. Zyla et al. (Particle Data Group), PTEP , 083C01 (2020).[4] A. Sakharov, Sov. Phys. Usp. , 392 (1991).[5] M. Fukugita and T. Yanagida, Phys. Lett. B174 , 45 (1986).[6] R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. , 912 (1980).[7] J. Schechter and J. Valle, Phys. Rev. D , 2227 (1980).[8] S. Petcov, Phys. Lett. B , 245 (1982).[9] W. Haxton and G. Stephenson, Prog. Part. Nucl. Phys. , 409 (1984).[10] C. Lam, Phys. Lett. B , 214 (2001), hep-ph/0104116.[11] P. Harrison and W. Scott, Phys. Lett. B , 219 (2002), hep-ph/0210197.[12] K. Abe et al. (T2K), Phys. Rev. Lett. , 061802 (2014), 1311.4750.[13] D. S. Ayres et al. (NOvA) (2004), hep-ex/0503053.[14] K. Abe et al. (T2K), Nature , 339 (2020), [Erratum: Nature 583, E16 (2020)], 1910.03887.[15] M. A. Acero et al. (NOvA), Phys. Rev. Lett. , 151803 (2019), 1906.04907.[16] R. Acciarri et al. (DUNE) (2015), 1512.06148.[17] B. Abi et al. (DUNE) (2020), 2002.03005.[18] B. Abi et al. (DUNE) (2020), 2006.16043.[19] K. Abe et al. (Hyper-Kamiokande Proto-Collaboration), PTEP , 053C02 (2015), 1502.05199.[20] K. Abe et al. (Hyper-Kamiokande), PTEP , 063C01 (2018), 1611.06118.[21] E. Baussan et al. (ESSnuSB), Nucl. Phys. B , 127 (2014), 1309.7022.[22] P. de Salas, D. Forero, S. Gariazzo, P. Mart´ınez-Mirav´e, O. Mena, C. Ternes, M. T´ortola, and J. Valle (2020), 2006.11237.[23] Valencia-Globalﬁt, http://globalfit.astroparticles.es/ (2020).[24] F. Capozzi, E. Lisi, A. Marrone, and A. Palazzo, Prog. Part. Nucl. Phys. , 48 (2018), 1804.09678.[25] I. Esteban, M. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou (2020), 2007.14792.[26] A. Aguilar-Arevalo et al. (LSND), Phys. Rev. D , 112007 (2001), hep-ex/0104049.[27] A. Aguilar-Arevalo et al. (MiniBooNE), Phys. Rev. Lett. , 101802 (2009), 0812.2243.[28] T. Mueller et al., Phys. Rev. C , 054615 (2011), 1101.2663.[29] G. Mention, M. Fechner, T. Lasserre, T. Mueller, D. Lhuillier, M. Cribier, and A. Letourneau, Phys. Rev. D , 073006(2011), 1101.2755.[30] P. Huber, Phys. Rev. C , 024617 (2011), [Erratum: Phys.Rev.C 85, 029901 (2012)], 1106.0687.[31] D. Decamp et al. (ALEPH), Phys. Lett. B , 519 (1989).[32] A. Aguilar-Arevalo et al. (MiniBooNE), Phys. Rev. Lett. , 221801 (2018), 1805.12028.[33] A. Aguilar-Arevalo et al. (MiniBooNE) (2020), 2006.16883.[34] M. Aartsen et al. (IceCube) (2020), 2005.12943.[35] A. Osipowicz et al. (KATRIN) (2001), hep-ex/0109033.[36] M. Antonello et al. (MicroBooNE, LAr1-ND, ICARUS-WA104) (2015), 1503.01520.[37] Y. Ko et al. (NEOS), Phys. Rev. Lett. , 121802 (2017), 1610.05134.[38] Y. Abreu et al. (SoLid), JINST , P04024 (2017), 1703.01683.[39] A. Serebrov et al. (NEUTRINO-4), Pisma Zh. Eksp. Teor. Fiz. , 209 (2019), 1809.10561.[40] A. Serebrov and R. Samoilov, Pisma Zh. Eksp. Teor. Fiz. , 211 (2020), 2003.03199.[41] A. Serebrov et al. (2020), 2005.05301.[42] A. Serebrov (Neutrino-4), Observation of sterile antineutrino oscillation in Neutrino-4 experiment at SM-3 reactor , talkat Neutrino2020, June, 2020. ” https://indico.fnal.gov/event/43209/contributions/187878/attachments/129237/158638/Serebrov_Neutrino-4_25july.pdf ” (2020).[43] C. Giunti, Phys. Rev. D , 095025 (2020), 2004.07577.[44] P. Coloma, P. Huber, and T. Schwetz (2020), 2008.06083.[45] C. Giunti, Y. F. Li, C. A. Ternes, and Y. Y. Zhang (2021), 2101.06785.[46] J. Ashenfelter et al. (PROSPECT), Phys. Rev. Lett. , 251802 (2018), 1806.02784.[47] N. Allemandou et al. (STEREO), JINST , P07009 (2018), 1804.09052.[48] H. Almaz´an et al. (STEREO), Phys. Rev. D , 052002 (2020), 1912.06582.[49] I. Alekseev et al. (DANSS), Phys. Lett. B , 56 (2018), 1804.04046.[50] C. Rott (JSNS2), J. Phys. Conf. Ser. , 012176 (2020).[51] T. Maruyama (JSNS2), JSNS Experiment , talk at Neutrino2020, June, 2020. ” https://indico.fnal.gov/event/43209/contributions/187854/attachments/129166/159526/Neutrino2020_JSNS2_v3.pdf ” (2020).[52] D. Akimov et al. (COHERENT), Science , 1123 (2017), 1708.01294.[53] O. Miranda, D. Papoulias, O. Sanders, M. T´ortola, and J. Valle (2020), 2008.02759.[54] R. Gandhi, B. Kayser, M. Masud, and S. Prakash, JHEP , 039 (2015), 1508.06275. [55] S. B¨oser, C. Buck, C. Giunti, J. Lesgourgues, L. Ludhova, S. Mertens, A. Schukraft, and M. Wurm, Prog. Part. Nucl.Phys. , 103736 (2020), 1906.01739.[56] N. Klop and A. Palazzo, Phys. Rev. D , 073017 (2015), 1412.7524.[57] J. M. Berryman, A. de Gouvˆea, K. J. Kelly, and A. Kobach, Phys. Rev. D , 073012 (2015), 1507.03986.[58] A. Palazzo, Phys. Lett. B , 142 (2016), 1509.03148.[59] S. K. Agarwalla, S. S. Chatterjee, A. Dasgupta, and A. Palazzo, JHEP , 111 (2016), 1601.05995.[60] S. K. Agarwalla, S. S. Chatterjee, and A. Palazzo, JHEP , 016 (2016), 1603.03759.[61] S. K. Agarwalla, S. S. Chatterjee, and A. Palazzo, Phys. Rev. Lett. , 031804 (2017), 1605.04299.[62] D. Dutta, R. Gandhi, B. Kayser, M. Masud, and S. Prakash, JHEP , 122 (2016), 1607.02152.[63] J. Rout, M. Masud, and P. Mehta, Phys. Rev. D , 075035 (2017), 1702.02163.[64] K. J. Kelly, Phys. Rev. D , 115009 (2017), 1703.00448.[65] M. Ghosh, S. Gupta, Z. M. Matthews, P. Sharma, and A. G. Williams, Phys. Rev. D , 075018 (2017), 1704.04771.[66] S. Choubey, D. Dutta, and D. Pramanik, Phys. Rev. D , 056026 (2017), 1704.07269.[67] P. Coloma, D. V. Forero, and S. J. Parke, JHEP , 079 (2018), 1707.05348.[68] J. Tang, Y. Zhang, and Y.-F. Li, Phys. Lett. B , 217 (2017), 1708.04909.[69] S. Choubey, D. Dutta, and D. Pramanik, Eur. Phys. J. C , 339 (2018), 1711.07464.[70] S. K. Agarwalla, S. S. Chatterjee, and A. Palazzo, JHEP , 091 (2018), 1801.04855.[71] S. Gupta, Z. M. Matthews, P. Sharma, and A. G. Williams, Phys. Rev. D , 035042 (2018), 1804.03361.[72] S. Choubey, D. Dutta, and D. Pramanik, Eur. Phys. J. C , 968 (2019), 1811.08684.[73] A. De Gouvˆea, K. J. Kelly, G. Stenico, and P. Pasquini, Phys. Rev. D , 016004 (2019), 1904.07265.[74] S. Kumar Agarwalla, S. S. Chatterjee, and A. Palazzo, JHEP , 174 (2019), 1909.13746.[75] R. Majhi, C. Soumya, and R. Mohanta, J. Phys. G , 095002 (2020), 1911.10952.[76] M. Ghosh, T. Ohlsson, and S. Rosauro-Alcaraz, JHEP , 026 (2020), 1912.10010.[77] S. S. Chatterjee and A. Palazzo (2020), 2005.10338.[78] C. Giunti and T. Lasserre, Ann. Rev. Nucl. Part. Sci. , 163 (2019), 1901.08330.[79] A. Diaz, C. Arg¨uelles, G. Collin, J. Conrad, and M. Shaevitz, Phys. Rept. , 1 (2020), 1906.00045.[80] A. Palazzo, Universe , 41 (2020).[81] A. Chatla and B. A. Bambah, in (2020), 2010.06321.[82] M. Dentler, A. Hern´andez-Cabezudo, J. Kopp, P. A. Machado, M. Maltoni, I. Martinez-Soler, and T. Schwetz, JHEP ,010 (2018), 1803.10661.[83] B. Yue, W. Li, J. Ling, and F. Xu, Chin. Phys. C , 103001 (2020), 1906.03781.[84] P. Huber, M. Lindner, and W. Winter, Comput. Phys. Commun. , 195 (2005), hep-ph/0407333.[85] P. Huber, J. Kopp, M. Lindner, M. Rolinec, and W. Winter, Comput. Phys. Commun. , 432 (2007), hep-ph/0701187.[86] J. Kopp, Int. J. Mod. Phys. C19 , 523 (2008), physics/0610206.[87] J. Kopp, M. Lindner, T. Ota, and J. Sato, Phys. Rev.

D77 , 013007 (2008), 0708.0152.[88] T. Alion et al. (DUNE) (2016), 1606.09550.[89] J. Rout, S. Roy, M. Masud, M. Bishai, and P. Mehta, Phys. Rev. D , 116018 (2020), 2009.05061.[90] K. Abe et al. (T2K), Nucl. Instrum. Meth. A , 106 (2011), 1106.1238.[91] K. Abe et al. (Hyper-Kamiokande) (2018), 1805.04163.[92] P. Huber, M. Lindner, and W. Winter, Nucl. Phys.

B645 , 3 (2002), hep-ph/0204352.[93] G. L. Fogli, E. Lisi, A. Marrone, D. Montanino, and A. Palazzo, Phys. Rev.

D66 , 053010 (2002), hep-ph/0206162.[94] M. Gonzalez-Garcia and M. Maltoni, Phys.Rev.

D70 , 033010 (2004), hep-ph/0404085.[95] R. Gandhi, P. Ghoshal, S. Goswami, P. Mehta, S. U. Sankar, and S. Shalgar, Phys. Rev.

D76 , 073012 (2007), 0707.1723.[96] X. Qian, A. Tan, W. Wang, J. J. Ling, R. D. McKeown, and C. Zhang, Phys. Rev.

D86 , 113011 (2012), 1210.3651.[97] M. Mitra, G. Senjanovic, and F. Vissani, Nucl. Phys.