# A New Approach to Probe Non-Standard Interactions in Atmospheric Neutrino Experiments

Anil Kumar, Amina Khatun, Sanjib Kumar Agarwalla, Amol Dighe

PPrepared for submission to JHEP

IP/BBSR/2020-8, TIFR/TH/20-49

A New Approach to Probe Non-Standard Interactionsin Atmospheric Neutrino Experiments

Anil Kumar, a,b,c

Amina Khatun, d Sanjib Kumar Agarwalla, a,c,e

Amol Dighe f a Institute of Physics, Sachivalaya Marg, Sainik School Post, Bhubaneswar 751005, India b Applied Nuclear Physics Division, Saha Institute of Nuclear Physics, Block AF, Sector 1, Bid-hannagar, Kolkata 700064, India c Homi Bhabha National Institute, Anushakti Nagar, Mumbai 400094, India d Comenius University, Mlynsk´a dolina F1, SK842 48 Bratislava, Slovakia e International Centre for Theoretical Physics, Strada Costiera 11, 34151 Trieste, Italy f Tata Institute of Fundamental Research, Homi Bhabha Road, Colaba, Mumbai 400005, India

E-mail: [email protected] (ORCID: 0000-0002-8367-8401) , [email protected] (ORCID: 0000-0003-3493-607X) , [email protected] (ORCID: 0000-0002-9714-8866) , [email protected] (ORCID: 0000-0001-6639-0951) Abstract:

We propose a new approach to explore the neutral-current non-standard neu-trino interactions (NSI) in atmospheric neutrino experiments using oscillation dips andvalleys in reconstructed muon observables, at a detector like ICAL that can identify themuon charge. We focus on the ﬂavor-changing NSI parameter ε µτ , which has the maximumimpact on the muon survival probability in these experiments. We show that non-zero ε µτ shifts the oscillation dip locations in L/E distributions of the up/down event ratios ofreconstructed µ − and µ + in opposite directions. We introduce a new variable ∆ d repre-senting the diﬀerence of dip locations in µ − and µ + , which is sensitive to the magnitude aswell as the sign of ε µτ , and is independent of the value of ∆ m . We further note that theoscillation valley in the ( E , cos θ ) plane of the reconstructed muon observables bends in thepresence of NSI, its curvature having opposite signs for µ − and µ + . We demonstrate theidentiﬁcation of NSI with this curvature, which is feasible for detectors like ICAL havingexcellent muon energy and direction resolutions. We illustrate how the measurement ofcontrast in the curvatures of valleys in µ − and µ + can be used to estimate ε µτ . Using theseproposed oscillation dip and valley measurements, the achievable precision on | ε µτ | at 90%C.L. is about 2% with 500 kt · yr exposure. The eﬀects of statistical ﬂuctuations, systematicerrors, and uncertainties in oscillation parameters have been incorporated using multiplesets of simulated data. Our method would provide a direct and robust measurement of ε µτ in the multi-GeV energy range. Keywords:

Atmospheric Neutrinos, NSI,

L/E

Analysis, Oscillation Dip, Oscillation Val-ley, ICAL, INO

ArXiv ePrint: a r X i v : . [ h e p - ph ] J a n ontents ε µτ on the oscillation dip in the L ν /E ν ε µτ on the oscillation valley in the ( E ν , cos θ ν ) plane 7 L rec µ /E rec µ distributions 103.2 Distributions in the ( E rec µ , cos θ rec µ ) plane 11 d for determining ε µτ ε µτ from the measurement of ∆ d ∆ m in the Absence of NSI 186 Identifying NSI through Oscillation Valley 20 ε µτ from oscillation valley 22 The phenomenon of neutrino oscillations has now been well-established, from measurementsat the solar, atmospheric, reactor as well as accelerator experiments with short and longbaselines [1]. Neutrino oscillations are the consequences of mixing among diﬀerent neutrinoﬂavors and non-degenerate values of neutrino masses, with at least two neutrino massesnonzero. However, neutrinos are massless in the Standard Model (SM) of particle physics,and therefore, physics beyond the SM (BSM) is needed to accommodate nonzero neutrinomasses and mixing. Many models of BSM physics suggest new non-standard interactions(NSI) of neutrinos [2], which may aﬀect neutrino production, propagation, and detection ina given experiment. The possible impact of these NSI at neutrino oscillation experimentshave been studied extensively, for example see Refs. [3–13]. In this paper, we propose a newmethod for identifying NSI at atmospheric neutrino experiments which can reconstruct theenergy, direction, as well as charge of the muons produced in the detector due to charged-current interactions of ν µ and ¯ ν µ . – 1 –e shall focus on the neutral-current NSI, which may be described at low energies viaeﬀective four-fermion dimension-six operators as [2] L NC − NSI = − √ G F ε fαβ,C (¯ ν α γ ρ P L ν β ) ( ¯ f γ ρ P C f ) , (1.1)where G F is the Fermi constant. The dimensionless parameters ε fαβ,C describe the strengthof NSI, where the superscript f ∈ { e, u, d } denotes the matter fermions ( e : electron, u :up-quark, d : down-quark), and the indices α, β ∈ { e, µ, τ } refer to the neutrino ﬂavors.The subscript C ∈ { L, R } represents the chiral projection operator P L = (1 − γ ) / P R = (1 + γ ) /

2. The hermiticity of the interactions demands ε fβα,C = ( ε fαβ,C ) ∗ .The eﬀective NSI parameter relevant for the neutrino propagation through matter is ε αβ ≡ (cid:88) f = e,u,d (cid:16) ε fLαβ + ε fRαβ (cid:17) N f N e ≡ (cid:88) f = e,u,d ε fαβ N f N e , (1.2)where N f is the number density of fermion f . In the approximation of a neutral andisoscalar Earth, the number densities of electrons, protons, and neutrons are identical,which implies N u ≈ N d ≈ N e . Thus, ε αβ ≈ ε eαβ + 3 ε uαβ + 3 ε dαβ . (1.3)In the presence of NSI, the modiﬁed eﬀective Hamiltonian for neutrino propagation throughmatter is H eﬀ = 12 E U m

00 0 ∆ m U † + V CC ε ee ε eµ ε eτ ε ∗ eµ ε µµ ε µτ ε ∗ eτ ε ∗ µτ ε ττ , (1.4)where U is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix that parametrizes neu-trino mixing, and ∆ m ij ≡ m i − m j are the mass-squared diﬀerences. The quantity V CC ≡ √ G F N e is the eﬀective matter potential due to the coherent elastic forward scat-tering of neutrinos with electrons inside the medium via the SM gauge boson W . Forantineutrinos, V CC → − V CC , U → U ∗ , and ε αβ → ε ∗ αβ .In the present study, we suggest a novel approach to unravel the presence of ﬂavor-changing neutral current NSI parameter ε µτ , via its eﬀect on the propagation of multi-GeVatmospheric neutrinos and antineutrinos through Earth matter. We choose ε µτ as the NSIparameter to focus on, since it can signiﬁcantly aﬀect the evolution of the mixing angle θ and the mass-squared diﬀerence ∆ m in matter [14], which in turn would alter the survivalprobabilities of atmospheric muon neutrinos and antineutrinos substantially [15, 16]. Ingeneral, ε µτ can be complex, i.e. ε µτ ≡ | ε µτ | e iφ µτ . However, in the disappearance channels ν µ → ν µ and ¯ ν µ → ¯ ν µ that dominate in our analysis , ε µτ appears only as | ε µτ | cos φ µτ at the leading order [14]. Thus, a complex phase only changes the eﬀective value of ε µτ to a real number between −| ε µτ | and + | ε µτ | at the leading order. We take advantage of The appearance channels ν e → ν µ and ¯ ν e → ¯ ν µ also contribute to the muon events in our analysis,however these channels are not aﬀected by ε µτ to leading order [14]. – 2 –his observation, and restrict ourselves to real values of ε µτ in the range − . ≤ ε µτ ≤ . ε µτ with | ε µτ | ≤ . | ε µτ | turns out tobe | ε µτ | < σ C.L.. The existing bounds on ε µτ from various neutrino oscillationexperiments are listed in Table 1. An important point to note is the energies of neutrinosinvolved in these measurements: the IceCube results are obtained using high energy events( >

300 GeV) [18], the energy threshold of DeepCore is around 10 GeV [19], while theSuper-K experiment is more eﬃcient in the sub-GeV energy range [20].Experiment 90% C.L. BoundIceCube [18] − . < ε µτ < . − . < ε µτ < . | ε µτ | < . Table 1 : Existing bounds on ε µτ at 90% conﬁdence level.The proposed 50 kt magnetized Iron Calorimeter (ICAL) detector at the India-basedNeutrino Observatory (INO) [21, 22] would be sensitive to multi-GeV neutrinos, since it caneﬃciently detect muons in the energy range 1–25 GeV. Note that the MSW resonance [2, 23,24] due to Earth matter takes place for neutrino energies around 4–10 GeV, so ICAL wouldalso be in a unique position to detect any interplay between the matter eﬀects and NSI.Another important feature is that ICAL can explore physics in neutrinos and antineutrinosseparately, unlike Super-K and IceCube/DeepCore. The studies of physics potential ofICAL for detecting NSI have shown that, using the reconstructed muon momentum, itwould be possible to obtain a bound of | ε µτ | < .

015 at 90% C.L. [25] with 500 kt · yrexposure. When information on the reconstructed hadron energy in each event is alsoincluded, the expected 90% C.L. bound improves to | ε µτ | < .

010 [26]. The results in [25,26] are obtained using a χ analysis with the pull method [15, 27, 28]The wide range of neutrino energies and baselines available in atmospheric neutrinoexperiments oﬀer an opportunity to study the features of “oscillation dip” and “oscillationvalley” in the reconstructed µ − and µ + observables, as demonstrated in [29]. These featurescan be clearly identiﬁed in the ratios of upward-going and downward-going muon eventsat ICAL. If the muon neutrino disappearance is solely due to non-degenerate masses andnon-zero mixing of neutrinos, then the valley in both µ − and µ + is approximately a straightline. The location of the dip, and the alignment of the valley, can be used to determine∆ m [29]. These features may undergo major changes in the presence of NSI, and can actas smoking gun signals for NSI.The novel approach, which we propose in this paper, is to probe the NSI parameter ε µτ based on the elegant features associated with the oscillation dip and valley. For the oscil-lation dip feature, we note that non-zero ε µτ shifts the oscillation dip location in opposite– 3 –irections for µ − and µ + . We demonstrate that this opposite shift in dip location due toNSI can be clearly seen in the µ − and µ + data by reconstructing L rec µ /E rec µ distributions,thanks to the excellent energy and direction resolutions for muons at ICAL. We develop awhole new analysis methodology to extract the information on ε µτ using the dip locations.For this, we deﬁne a new variable exploiting the contrast between the shifts in reconstructeddip locations, which eliminates the dependence of our results on the actual value of ∆ m .For the oscillation valley feature, we notice that the valley becomes curved in the presenceof non-zero ε µτ , and the direction of this bending is opposite for neutrino and antineutrino.We then demonstrate that this opposite bending can indeed be observed in expected µ − and µ + events. We propose a methodology to extract the information on the bending ofthe valley in terms of reconstructed muon variables, and use it for identifying NSI.In Sec. 2, we discuss the oscillation probabilities of neutrino and antineutrino in thepresence of non-zero ε µτ , and discuss the shifts in the dip locations as well as the bendingof the oscillation valleys in the survival probabilities of ν µ and ¯ ν µ . In Sec. 3, we investigatethe survival of these two striking features in the reconstructed L rec µ /E rec µ distributions andin the ( E rec µ , cos θ rec µ ) distributions of µ − and µ + events separately at ICAL. In Sec. 4,we propose a novel variable for identifying the NSI, which is based on the contrast in theshifts of dip locations in µ − and µ + . This variable leads to the calibration of ε µτ , and isused to ﬁnd the expected bound on ε µτ from a 500 kt · yr exposure of ICAL. In Sec. 5, wecome up with a new procedure for determining the alignment of the oscillation valley andestimating the value of ∆ m in the absence of NSI, which we extend to the NSI analysis inSec. 6. Here, we measure the contrast in the curvatures of the oscillation valleys in µ − and µ + in the presence of NSI, and use it for determining the expected bound on ε µτ from thevalley analysis at ICAL. Finally, in Sec. 7, we summarize our ﬁndings and oﬀer concludingremarks. In the limit of θ →

0, and the approximation of one mass scale dominance scenario[∆ m L/ (4 E ) (cid:28) ∆ m L/ (4 E )] and constant matter density, the survival probability of ν µ when traveling a distance L ν is given by [15] P ν µ → ν µ = 1 − sin θ eﬀ sin (cid:20) ξ ∆ m L ν E ν (cid:21) , (2.1)where sin θ eﬀ = | sin 2 θ + 2 β η µτ | ξ , (2.2) ξ = (cid:113) | sin 2 θ + 2 β η µτ | + cos θ , (2.3)and η µτ = 2 E ν V CC ε µτ ∆ m , (2.4)where β ≡ sgn(∆ m ). That is, β = +1 for normal mass ordering (NO), and β = − θ = 45 ◦ ), Eq. 2.1 reduces to the following simpleexpression [16]: P ν µ → ν µ = cos (cid:20) L ν (cid:18) ∆ m E ν + ε µτ V CC (cid:19)(cid:21) . (2.5)We further ﬁnd that the correction due to the deviation of θ from maximality, χ ≡ θ − π/

4, is second order in the small parameter χ , and hence can be neglected. Here, onenotes that the NSI parameter ε µτ primarily modiﬁes the wavelength of neutrino oscillations.It also comes multiplied with V CC , which increases at higher baselines inside the Earth’smatter. Thus, the modiﬁcation in ν µ survival probabilities due to NSI varies with thebaseline L ν (or cos θ ν ). ε µτ on the oscillation dip in the L ν /E ν In Fig. 1, we present the survival probabilities of ν µ and ¯ ν µ as functions of L ν /E ν , for twoﬁxed values of cos θ ν ( i.e. cos θ ν = − . , − . ε µτ ( i.e. ε µτ = +0 . , . , − . θ sin θ sin θ | ∆ m | (eV ) ∆ m (eV ) δ CP Mass Ordering0.855 0.5 0.0875 2 . × − . × − Table 2 : The benchmark values of oscillation parameters that we use in the analysis.Normal mass ordering corresponds to m < m < m .Figure 1 also indicates the sensitivity ranges for the atmospheric neutrino experimentsSuper-K, ICAL, and IceCube. Note that these ranges are diﬀerent for cos θ ν = − . − .

8, which correspond to L ν around 5100 km and 10200 km, respectively Whilecalculating these L ν /E ν ranges, the energy ranges chosen are those for which the detectorsperform very well: we use the E ν range of 100 MeV – 5 GeV for Super-K [30], 1–25 GeVfor ICAL [21], and 100 GeV – 10 PeV for IceCube [31]. Note that these ranges are onlyindicative. The following observations may be made from the ﬁgure. • For log [ L ν /E ν ] in the range of [0 – 1.5], the survival probabilities of both ν µ and¯ ν µ are observed to be suppressed in the presence of non-zero ε µτ . This is because,although the oscillations due to neutrino mass-squared diﬀerence do not develop forsuch small value of L ν /E ν , i.e. ∆ m L ν /E ν (cid:28)

1, the disappearance of ν µ is possibledue to the L ν ε µτ V CC term, which can be high for large baselines (for core passingneutrinos, V CC is large). For example, log [ L ν /E ν ] = 1 may correspond to L ν = 5000km and E ν = 500 GeV, so that ∆ m L ν /E ν ≈ . ρ ≈ . ε µτ = 0 .

1, we have L ν ε µτ V CC ≈ .

93, which takes the oscillation probability away from unity. The eﬀectof NSI at such small L ν /E ν is energy-independent, and can be seen at detectors likeIceCube [18] due to its better performance at high energy. • As we go to higher value of log [ L ν /E ν ] ( > Lε µτ V CC term in Eq. 2.5, and the competition– 5 – (GeV)] n (km)/E n [L Log ) mn ﬁ mn P ( ICALIceCube Super-K = -0.4 n q NO, cos = 0.1 tm e = 0.0 tm e = -0.1 tm e (GeV)] n (km)/E n [L Log ) mn ﬁ mn P ( ICALIceCube Super-K = -0.4 n q NO, cos = 0.1 tm e = 0.0 tm e = -0.1 tm e (GeV)] n (km)/E n [L Log ) mn ﬁ mn P ( ICALIceCube Super-K = -0.8 n q NO, cos = 0.1 tm e = 0.0 tm e = -0.1 tm e (GeV)] n (km)/E n [L Log ) mn ﬁ mn P ( ICALIceCube Super-K = -0.8 n q NO, cos = 0.1 tm e = 0.0 tm e = -0.1 tm e Figure 1 : The survival probabilities of ν µ and ¯ ν µ as functions of log ( L ν /E ν ) in left andright panels, respectively. The black lines correspond to ε µτ = 0 (SI), whereas red and bluelines are with ε µτ = 0 .

1, and − .

1, respectively. The neutrino direction taken in upper(lower) panels is cos θ ν = − . − . L ν /E ν ranges that these detectors are well-suited for. We consider normal mass ordering,and the benchmark oscillation parameters given in Table 2.between these two terms leads to an energy dependence. When the oscillations dueto the mass splitting are the dominating contribution, the change in oscillation wave-length due to non-zero ε µτ results in a shift of the oscillation minima towards left orright side (lower or higher values of L ν /E ν , respectively), depending on the amplitudeof ε µτ and its sign. The direction of shift in the dip location depends on whetherit is neutrino or antineutrino (since the sign of V CC for them are opposite), and onthe neutrino mass ordering. This eﬀect on the shift in the dip location is discussedin next section in detail. This region is relevant for Super-K, INO, and most of thelong-baseline experiments. – 6 –e now discuss the modiﬁcation of the oscillation dip due to non-zero ε µτ . First, usingthe approximate expression of ν µ survival probability from Eq. 2.5, we obtain the value of L ν /E ν at the ﬁrst oscillation minimum (or the dip): L ν E ν (cid:12)(cid:12)(cid:12)(cid:12) dip = 2 π | ∆ m + 4 ε µτ V CC E ν | . (2.6)The above expression may be written by expressing V CC in terms of the line-averagedmatter density ρ and taking care of units, L ν [km] E ν [GeV] (cid:12)(cid:12)(cid:12)(cid:12) dip = π (cid:12)(cid:12) . · ∆ m [eV ] ± . × − · ρ [g / cm ] · Y e · ε µτ · E ν [GeV] (cid:12)(cid:12) . (2.8)Here a positive sign in the denominator corresponds to neutrinos, whereas a negative sign isfor antineutrinos. This approximation is useful for understanding the shift of dip positionwith non-zero ε µτ . Let us take the case of neutrino with normal ordering. With positive ε µτ , the denominator of Eq. 2.6 increases, thus the oscillation minimum appears at a lowervalue of L ν /E ν than the SI. On the other hand, with negative ε µτ , the oscillation dip wouldoccur at a higher value of L ν /E ν . These two features can be seen clearly in the left panelsof Fig. 1. For antineutrinos and the same mass ordering, since the matter potential has theopposite sign, the shift of oscillation dip is in the opposite direction as compared that forneutrinos, given the same ε µτ .Expanding Eq. 2.8 to ﬁrst order in ε µτ , one can write L ν [km] E ν [GeV] (cid:12)(cid:12)(cid:12)(cid:12) dip = π (cid:12)(cid:12) . × ∆ m [eV ] (cid:12)(cid:12) ∓ π × . × − · ρ [g / cm ] · Y e · E ν [GeV] β (∆ m ) [eV ] · ε µτ , (2.9)where β ≡ sgn(∆ m ), as deﬁned earlier. Here, the negative sign corresponds to neutrinos,whereas the positive sign is for antineutrino. This indicates that, for small values of ε µτ ,the shift in the dip location will be linear in ε µτ , and will have opposite sign for neutrinosand antineutrinos, as well as for the two mass orderings. We have checked that for E ν < ε µτ < . ε µτ on the oscillation valley in the ( E ν , cos θ ν ) plane In Fig. 2, we show the oscillograms of survival probabilities for ν µ and ¯ ν µ , in the planeof neutrino energy and cosine of neutrino zenith angle (cos θ ν ), with NO as the neutrinomass ordering. We show only upward-going neutrinos (cos θ ν < θ ν > L ν is very small and neutrino oscillations do We can write V CC approximately as a function of matter density ρV CC ≈ ± . × Y e × ρ g / cm eV . (2.7)Here, Y e = N e N p + N n is the relative electron number density. In an electrically neutral and isoscalar medium, Y e = 0 .

5. The positive (negative) sign is for neutrino (antineutrino). The zenith angle θ ν is related to the baseline L ν by L ν = (cid:113) ( R + h ) − ( R − d ) sin θ ν − ( R − d ) cos θ ν , (2.10) – 7 – igure 2 : The oscillogram for the survival probabilities in ( E ν , cos θ ν ) plane for neutrinoand antineutrino in left and right panels, respectively, with normal mass ordering and thebenchmark oscillation parameters from Table 2. The value of ε µτ is taken as 0 (SI), +0.1,and -0.1 in top, middle and bottom panels, respectively. where R , h , and d are the radius of Earth, the average height from the Earth surface at which neutrinosare created, and the depth of the detector underground, respectively. In this study, we use R = 6371 km, h = 15 km, and d = 0 km. A small change in h and d does not aﬀect the ν µ oscillation probabilities because R (cid:29) h (cid:29) d . – 8 –ot develop, making P µµ ≈ .

0. The diﬀerences observed among the upper, middle, andlower panels are due to diﬀerent values of ε µτ . A major impact of NSI is observed onthe feature corresponding to the ﬁrst oscillation minima, seen in the ﬁgure as a broadblue/black diagonal band. We refer to this broad band as the oscillation valley [29]. Inthe SI case, the oscillation valley appears like a triangle with straight edges, while in thepresence of NSI, its edges seem to acquire a curvature. The sign of this curvature is oppositefor neutrinos and antineutrinos, as can be seen by comparing the left and right plots ofupper and lower panels of Fig. 2.The modiﬁcation in the oscillation valley due to NSI may be explained from the fol-lowing relation between E ν and cos θ ν at the ﬁrst oscillation minima. We rewrite Eq. 2.6with L ν ≈ R | cos θ ν | : E ν | valley ≈ | ∆ m | ( π/ | R cos θ ν | ) − β ε µτ V CC . (2.11)Taking care of units and expressing V CC in terms of line-averaged constant matter density ρ , the above expression may be written as E ν [GeV] | valley = | ∆m [eV ] | ( π/ | . · R[km] · cos θ ν | ) ∓ . × − · β · ε µτ · Y e · ρ [g / cm ] . (2.12)Here, the negative sign in the denominator corresponds to neutrinos, whereas the positivesign is for antineutrinos. Putting ε µτ = 0 in Eq. 2.12 gives the condition of the ﬁrstoscillation minima to be E ν = | (5 . /π ) · ∆ m · R [km] · cos θ ν | , and this relation clearlyshows that in the SI case, the minimum of the oscillation valley is a straight line in the( E ν , cos θ ν ) plane.Now in the normal mass ordering scenario, if ε µτ > E ν increases for a given cos θ ν as compared to the SI case. As a result, the oscillation valley(the broad blue/black band) bends towards higher values of energies in top left panel.On the other hand, for antineutrino, the oscillation valley tilts towards lower energies for ε µτ >

0, as can be seen in the top right panel. For negative values of ε µτ , the oscillationvalley bends in the opposite direction, for both neutrinos and antineutrinos, as shown inthe bottom panels of Fig. 2. In the inverse mass ordering scenario, the bending of theoscillation valley will be in the direction opposite to that in the normal mass orderingscenario.From Eqs. 2.9 and 2.12, the shift in the oscillation minima and the bending in theoscillation valley are in opposite directions for normal and inverted mass ordering. Theeﬀect of ε µτ therefore depends crucially on the mass ordering. In this paper, we presentour analysis in the scenario where the mass ordering is known to be NO. The analysis forIO may be performed in an exactly analogous manner. Our results are presented withthe exposure corresponding to 10 years of ICAL data-taking. It is expected that the massordering will be determined from the neutrino experiments (including ICAL) by this time.– 9 – Impact of NSI on the Event Distribution at ICAL

While the eﬀect of NSI on the neutrino and antineutrino survival probabilities was discussedin the last section, it is important to conﬁrm whether the features present at the probabilitylevel can survive in the observables at a detector, and if they can be reconstructed. Hereis where the response of the detector plays a crucial role. Neutrinos cannot be directlydetected in an experiment, however the charged leptons produced from their charged-current interactions in the detector contains information about their energy, direction, andﬂavor, which can be recovered depending on the nature of the detector.The upcoming ICAL detector at the India-based Neutrino Observatory (INO) will becomposed of 151 layers of magnetized iron plates; the 50 kt iron acting as the target, andaround 30,000 glass resistive plate chambers as detection elements. The magnetic ﬁeld ofstrength 1.5 Tesla and the time resolution of less than 1 ns [32–34] will enable ICAL to dis-tinguish µ − from µ + events in the multi-GeV energy range. In this study, we generate thecharged-current interactions of ν µ and ¯ ν µ , similar to Ref. [29], using neutrino and antineu-trino ﬂuxes calculated at the Theni site by Honda et.al. [35, 36] and the neutrino generatorNUANCE [37]. The reconstructed energy ( E rec µ ) and zenith angle ( θ rec µ ) of muons producedin neutrino and antineutrino interactions will be used in the analysis. The quantity L rec µ associated with the reconstructed muon direction cos θ rec µ is deﬁned (similar to Eq. 2.10) as L rec µ ≡ (cid:113) ( R + h ) − ( R − d ) sin θ rec µ − ( R − d ) cos θ rec µ . (3.1)Note that L rec µ is simply a proxy observable, and there is no need to associate it directlywith the distance traveled by the incoming neutrino. L rec µ /E rec µ distributions Observable Range Bin width Number of binslog (cid:104) L rec µ (km) E rec µ (GeV) (cid:105) [0, 1] 0.2 5 Table 3 : The binning scheme of log [ L rec µ /E rec µ ] for µ − and µ + events.Fig. 3 presents the log [ L rec µ /E rec µ ] distribution of the upward-going µ − and µ + events(U, cos θ rec µ <

0) with SI ( ε µτ = 0) and with SI + NSI ( ε µτ = ± .

1) expected with 10 yrexposure at ICAL. We include the statistical ﬂuctuation in the number of events by simulat-ing 100 independent data sets, and calculating the mean and root-mean-square deviation The values we mention for log [ L rec µ /E rec µ ] throughout the paper are calculated with L rec µ in the unitsof km and E rec µ in GeV. – 10 – .8 2 2.2 2.4 2.6 2.8 3 3.2 (GeV)] rec - m (km)/E rec - m [L Log U = 0.1 tm e = 0.0 tm e = -0.1 tm e - m yr (cid:215)

500 kt (GeV)] rec + m (km)/E rec + m [L Log U = 0.1 tm e = 0.0 tm e = -0.1 tm e + m yr (cid:215)

500 kt

Figure 3 : The log [ L rec µ /E rec µ ] distributions of µ − and µ + events in left and right panels,respectively, expected in 10 years at ICAL. The black lines correspond to SI ( (cid:15) µτ = 0),whereas red and blue lines correspond to ε µτ = − . [ L rec µ /E rec µ ], we use the same binning scheme as used in Ref. [29],which is shown in Table 3. We have a total of 34 non-uniform bins of log [ L rec µ /E rec µ ] inthe range 0–4. The downward-going events (D, cos θ rec µ >

0) are not aﬀected signiﬁcantlyby oscillations or by additional NSI interactions. The log [ L rec µ /E rec µ ]-distribution of µ − and µ + events is identical to the one shown in the upper panel of Fig. [4] in Ref. [29], andwe do not repeat it here.Fig. 3 clearly shows that the log [ L rec µ /E rec µ ] range of 2.4–3.0 would be important forNSI studies, since non-zero ε µτ has the largest eﬀect in this range. With positive ε µτ ,the number of µ − events is lower as compared to that of SI case, whereas the numberof µ + events is higher. If ε µτ is negative, then the modiﬁcations of the number of µ − and µ + events are the other way around. As a consequence, if a detector is not able todistinguish between the µ − and µ + events, the diﬀerence between SI and NSI would getdiluted substantially. The charge identiﬁcation capability of the magnetized ICAL detectorwould be crucial in exploiting this observation. E rec µ , cos θ rec µ ) plane In order to study the eﬀect of NSI on the distribution of events in the ( E rec µ , cos θ rec µ ) plane,we bin the data in reconstructed observables, E rec µ and cos θ rec µ . We have a total of 16non-uniform E rec µ bins in the range 1 – 25 GeV, and the whole range of − ≤ cos θ rec µ ≤ E rec µ and cos θ rec µ are binned with thesame binning scheme as used in [29], and is shown in Table 4. For demonstrating eventdistributions, we scale the 1000-year MC sample to an exposure of 10 years, and the– 11 –

1 -1 -1 -1 1 1 1 1 -1 1 1 1 -1 -1 -1 1 1 -1 1 1 -1 1 -1 -1 -1 -1 1 3 2 -2 -2 -1 1 3 1 1 -1 -2 -2 2 3 3 1 -1 -2 -2 -2 4 3 1 -1 -3 -4 -3 -2 2 1 -1 -2 -1 -2 -1 -1 - - - - - rec - m q cos ( G e V ) r e c -m E - - - - = 0.0) tm e U( - = 0.1) tm e , U( - m -1 -1 -1 -1 1 -1 -1 -1 1 -1 1 1 -1 2 1 -1 -1 1 1 -1 -1 -1 1 1 1 -1 -1 -2 2 1 1 -1 -1 1 1 1 1 1 1 1 1 2 1 1 1 2 2 2 2 2 2 2 1 1 1 2 2 1 1 1 1 - - - - - rec + m q cos ( G e V ) r e c + m E - - - - = 0.0) tm e U( - = 0.1) tm e , U( + m

1 -1 1 1 1 1 1 1 1 1 1 1 -1 1 -1 -1 -1 1 1 1 1 -1 -1 -1 1 1 1 1 -1 -1 -1 1 1 -1 -2 -2 1 2 3 3 1 1 1 3 3 2 1 1 3 2 3 4 4 4 3 2 1 4 5 5 5 4 3 2 2 1 2 2 2 2 2 1 1 1 - - - - - rec - m q cos ( G e V ) r e c -m E - - - - = 0.0) tm e U( - = -0.1) tm e , U( - m -1 -1 -1 -1 -1 -1 1 -1 1 -1 1 1 1 -1 1 1 -1 -1 1 -1 1 1 -1 1 2 1 -1 -1 1 2 1 -1 -1 2 2 1 -1 -2 -1 2 2 1 -2 -1 -1 -1 1 1 -1 -1 -1 -1 - - - - - rec + m q cos ( G e V ) r e c + m E - - - - = 0.0) tm e U( - = -0.1) tm e , U( + m Figure 4 : The ( E rec µ , cos θ rec µ ) distributions of diﬀerence of events between including non-zero NSI and only SI ( ε µτ = 0) expected with 500 kt · yr of ICAL. Left and right panels arefor µ − and µ + events, respectively, whereas upper and lower panels correspond to (cid:15) µτ = 0 . − .

1, respectively. We consider normal mass ordering, and the benchmark oscillationparameters given in Table 2.diﬀerence in the number of events with SI + NSI ( | ε µτ | = 0 .

1) and the number of eventswith SI, for µ − and µ + events are shown in Fig. 4.Observable Range Bin width Number of bins E rec µ [1, 5] 0.5 8 θ rec µ [-1.0, 1.0] 0.1 20 Table 4 : The binning scheme adopted for E rec µ and cos θ rec µ of µ − and µ + events.– 12 –rom this ﬁgure, we can make the following observations. • The events in the bins with with E rec µ > θ rec µ < − . µ − and µ + events. • There is almost no diﬀerence in the number of events between SI+NSI and SI at E rec µ < L ν ε µτ V CC (cid:46) ∆ m / (2 E ν ) at small energies. • As we go to higher energies, the mass-squared diﬀerence-induced neutrino oscillationsdie down, while the NSI-induced matter eﬀect term L ν ε µτ V CC increases. Therefore,large eﬀects of NSI are observed at high energies. • The eﬀect of NSI is also larger at higher baselines, since in this case, neutrinos passthrough inner and outer cores that have very high matter densities (around 11 g/cm and 13 g/cm , respectively), and hence higher values of L ν ε µτ V CC . • In many of the ( E rec µ , cos θ rec µ ) bins, the diﬀerence in the number of events in µ − and µ + events has the opposite sign. Therefore, muon charge information is a crucialingredient in the search for NSI. We have also seen this feature in the log [ L rec µ /E rec µ ]distributions of µ − and µ + events in Sec. 3.1. To study the eﬀect of NSI in oscillation dip, we focus on the ratio of upward-going (U)and downward-going (D) muon events as the observable to be studied. The up/down eventratio U/D is deﬁned for cos θ rec µ <

0, as [29]U/D( E rec µ , cos θ rec µ ) ≡ N ( E rec µ , −| cos θ rec µ | ) N ( E rec µ , | cos θ rec µ | ) , (4.1)where the numerator corresponds to the number of upward-going events and the denomina-tor is the number of downward-going events, with energy E rec µ and direction | cos θ rec µ | . Weassociate the U/D ratio with cos θ rec µ < L rec µ of upward-goingevents. This ratio is expected to serve as the proxy for the survival probabilities of ν µ and¯ ν µ at ICAL, since around 98% muon events at ICAL arise from the ν µ → ν µ oscillationchannel [29].One advantage of the use of the U/D ratio would be to nullify the eﬀects of systematicuncertainties like ﬂux normalization, cross sections, overall detection eﬃciency, and energydependent tilt error, as can be seen later. Note that the up-down symmetry of the detectorgeometry and detector response play an important role in this.Fig. 5 presents the L rec µ /E rec µ distributions of U/D with ε µτ = 0 (SI) and ± .

1, for anexposure of 10 years at ICAL. The mass ordering is taken as NO. The statistical ﬂuctuationsshown in the ﬁgure are the root-mean-square deviation obtained from the simulations of 100independent sets of 10 years data. Around log [ L rec µ /E rec µ ] ∼ . − .

6, we see a signiﬁcant– 13 – (GeV)] rec - m (km)/E rec - m [L Log U / D = 0.1 tm e = 0.0 tm e = -0.1 tm e - m yr (cid:215)

500 kt - d 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4 (GeV)] rec + m (km)/E rec + m [L Log U / D = 0.1 tm e = 0.0 tm e = -0.1 tm e + m yr (cid:215)

500 kt + d Figure 5 : The log ( L rec µ /E rec µ ) distributions of ratio of upward-going and downward-going µ − (left panel) and µ + (right panel) events using 10 years data of ICAL. The blacklines correspond to SI ( ε µτ = 0), whereas red and blue lines are for ε µτ = 0 . ε µτ ( ± . i.e. to lower orhigher values of L rec µ /E rec µ , respectively) from that of the SI case with non-zero ε µτ . Thisshift for µ − and µ + events is observed to be in opposite directions, as is expected fromthe discussions in Sec. 2.1. For example, for ε µτ >

0, the location of oscillation dip in µ − events gets shifted towards smaller L rec µ /E rec µ values. On other hand, in µ + events, theoscillation dip shifts to higher values of L rec µ /E rec µ . For ε µτ <

0, this shift is in the oppositedirections. Moreover, as the oscillation dip shifts towards the higher value of L rec µ /E rec µ ,the dip becomes shallower due to the eﬀect of rapid oscillation at high L ν /E ν . This eﬀectis visible in both, µ − and µ + . It can be easily argued that the modiﬁcation in oscillationdips of µ − and µ + due to non-zero ε µτ is the reﬂection of how the survival probabilitiesof neutrino and antineutrino changes in the presence of non-zero ε µτ , as expected fromEq. 2.6 in Sec. 2.1.Note that for a neutrino detector that is blind to the charge of muons, the dip itselfwould get diluted since diﬀerent average inelasticities of neutrino and antineutrino eventsat these energies would lead to slightly diﬀerent dip locations in µ − and µ + distributions.Moreover, the shift of oscillation dip towards opposite directions in µ − and µ + events dueto non-zero ε µτ would further contribute to the dilution of NSI eﬀects on the dip location.On the other hand, a detector like ICAL that has the charge identiﬁcation capability, notonly provides independent undiluted measurements of dip locations in µ − and µ + events,but also provides a unique novel observable that can cleanly calibrate against the value of ε µτ , as we shall see in the next section. – 14 – .65 2.7 2.75 2.8 2.85 2.9 2.95 3 - d - - tme eV -3 · = 2.37 m D eV -3 · = 2.46 m D eV -3 · = 2.57 m D + d - - tme eV -3 · = 2.37 m D eV -3 · = 2.46 m D eV -3 · = 2.57 m D - - - d D - - tme eV -3 · = 2.37 m D eV -3 · = 2.46 m D eV -3 · = 2.57 m D Figure 6 : Upper panels: The calibration of ε µτ with the dip location in U/D ratio of µ − (left panel) and µ + (right panel) events separately. Lower panel: The calibration of ε µτ with diﬀerence in dip location of µ − and µ + . Note that for − . < ε µτ < − .

07, the diplocation is simply the minimum U/D ratio since the ﬁtting is not stable for these values.We consider normal mass ordering, and the benchmark oscillation parameters given inTable 2. ∆ d for determining ε µτ Since the oscillation dip locations in the U/D distributions of both, µ − and µ + events,shift to lower / higher values of L rec µ /E rec µ depending on the value of ε µτ , the dip locationin either of these distributions may be used to estimate ε µτ , independently. We havealready developed a dip-identiﬁcation algorithm [29], where the dip location is not simplythe lowest value of U/D ratio, but takes into account information from the surroundingbins that have a U/D ratio below a threshold value. The algorithm essentially identiﬁes thecluster of contiguous bins that have a U/D ratio lower than any surrounding bins, and ﬁtsthese bins with a quadratic function whose minimum corresponds to the dip. The value oflog [ L rec µ /E rec µ ] corresponding to this minimum is termed as the dip location. We denotethe dip location for µ − and µ + events as d − and d + , respectively.– 15 –e obtain the calibration of ε µτ with the dip locations for µ − and µ + events separately,using 1000-year MC data. In the top panels of Fig. 6, we present the calibration curvesof ε µτ with d − (left panel) and d + (right panel), for three diﬀerent values of ∆ m . Thestraight line nature of the calibration curves can be understood qualitatively using Eq. 2.6.For a majority of events ( E r ecµ <

25 GeV), the ε µτ V CC E ν term is much smaller than ∆ m term, and therefore the linear expansion as shown in Eq. 2.9 is valid.In Eq. 2.9, the dominating ∆ m dependence of the dip location is through the ﬁrstterm. We can get rid of this dependence by introducing a new variable∆ d = d − − d + , (4.2)which is expected to be proportional to ε µτ , and can be used to calibrate for ε µτ . As may beseen in the bottom panel of Fig. 6, this indeed removes the dominant ∆ m -dependence inthe calibration of ε µτ . Note that some ∆ m -dependence still survives, which may be seenin the slightly diﬀerent slopes of the calibration curves for diﬀerent ∆ m values. However,this dependence is clearly negligible, as it appears in Eq. 2.9 as a multiplicative correctionto ε µτ . We have thus obtained a calibration for ε µτ that is almost independent of ∆ m .Note that ε µτ = 0 does not imply ∆ d = 0, since diﬀerent inelasticities and mattereﬀects in the neutrino and antineutrino channels mean that the dip locations in these twolocations are not identical even in the absence of NSI. ε µτ from the measurement of ∆ d We saw in the last section that one can infer the value of ε µτ in the observed data fromthe calibration curve between ∆ d and ε µτ , independent of the actual value of ∆ m . Weobtain the calibration curve using the 1000-year MC-simulated sample, and show it withthe solid blue line in the panels of Fig. 7.For estimating the expected constraints on ε µτ , we simulate 100 independent sets ofdata, each for 10 years exposure of ICAL with the benchmark oscillation parameters asgiven in Table 2, and ε µτ = 0. In Fig. 7, the range indicated along the x-axis by solid linesis the expected 90% C.L. range of the measured value of ∆ d , while the range indicatedalong the y-axis by the dark gray band is the expected 90% C. L. limit on ε µτ . From theﬁgure, we ﬁnd that one can expect to constrain ε µτ to − . < ε µτ < .

020 at 90% C.L.with 10 years of data.Since the calibration was found to be almost independent of the actual value of ∆ m ,we expect that even if the actual value of ∆ m were not exactly known, the resultswill not change. We have conﬁrmed this by generating 100 independent sets of 10-yrexposure, where the ∆ m values used as inputs follow a Gaussian distribution ∆ m =(2 . ± . × − eV , in accordance to the range allowed by the global ﬁt [38–40] ofavailable neutrino data.We also checked the impact of uncertainties in the values of the other oscillation pa-rameters. We ﬁrst simulated 100 statistically independent unoscillated data sets. Then foreach of these data sets, we take 20 random choices of oscillation parameters, according tothe Gaussian distributions∆ m = (7 . ± . × − eV , ∆ m = (2 . ± . × − eV , (4.3)– 16 – .3 - - - d D - - tme Fixed param., no syst.Varied param., with syst.

NO, 90% C.L.

Figure 7 : The blue stars and the blue solid line indicate the calibration curve of ε µτ against∆ d , obtained using 1000-year MC data with normal mass ordering. The gray bands showthe expected 90% C.L. limits on ∆ d (vertical bands), and hence on ε µτ (horizontal bands),with 500 kt · yr of ICAL exposure, when the actual value of ε µτ = 0. The light (dark)gray bands show the limits when the errors on oscillation parameters, and the impactof systematic uncertainties are included (excluded), as discussed in Sec. 4.2. For ﬁxed-parameter case, we use the benchmark oscillation parameters given in Table 2.sin θ = 0 . ± . , sin θ = 0 . ± . , sin θ = 0 . ± . , (4.4)guided by the present global ﬁt [41]. We keep δ CP = 0, since its eﬀect on ν µ survivalprobability is known to be highly suppressed in the multi-GeV energy range [14]. Thisprocedure eﬀectively enables us to consider the variation of our results over 2000 diﬀerentcombinations of oscillation parameters, to take into account the eﬀect of their uncertainties.We do not observe any signiﬁcant dilution of results due to the variation over oscillationparameters. This is expected, since (i) solar oscillation parameters ∆ m and θ do notcontribute signiﬁcantly to the ν µ survival probability in the multi-GeV range, (ii) thereactor mixing angle θ is already measured to a high precision, (iii) the mixing angle θ does not aﬀect the location of the oscillation minimum at the probability level, and (iv)our procedure of using ∆ d (see Eq. 4.2) to determine ε µτ minimizes the impact of ∆ m uncertainty, as shown in Sec. 4.1.In addition, we take into account the ﬁve major systematic uncertainties in the neutrinoﬂuxes and cross sections that are used in the standard ICAL analyses [21]. These ﬁveuncertainties are (i) 20% in overall ﬂux normalization, (ii) 10% in cross sections, (iii) 5%– 17 –n the energy dependence, (iv) 5% in the zenith angle dependence, and (iii) 5% in overallsystematics. For each of the 2000 simulated data sets, we modify the number of events ineach ( E rec µ , cos θ rec µ ) bin as N = N (0) (1 + δ )(1 + δ )( E rec µ /E ) δ (1 + δ cos θ rec µ )(1 + δ ) , (4.5)where N (0) is the theoretically predicted number of events, and E = 2 GeV. Here( δ , δ , δ , δ , δ ) is an ordered set of random numbers, generated separately for each sim-ulated data set, with the Gaussian distributions centered around zero and the 1 σ widthsgiven by (20% , , , , ε µτ to − . < ε µτ < .

024 at 90% C.L.. The results denotingthe eﬀects of these uncertainties have been shown in Fig. 7. ∆ m in the Absence of NSI The reconstruction of the oscillation valley using the distributions of the upward-goingand downward-going muons in ( E rec µ , cos θ rec µ ) plane is discussed in detail in Ref. [29]. Ithas been shown that the alignment of the valley can provide a measurement of ∆ m . Inthis paper, we use an alternative procedure for getting the alignment of the oscillationvalley, which is more robust and more well-motivated as compared to what was used inRef. [29], and determine the expected precision that the ICAL detector can achieve on theatmospheric oscillation parameter ∆ m , using this method.Fig. 8 shows the distributions of the ratio of upward-going and downward going eventsat ICAL in the ( E rec µ , cos θ rec µ ) plane for 10 years of ICAL data. The left and right panelsshow the U/D ratios for µ − and µ + events, respectively. For the sake of demonstration(in order to reduce the statistical ﬂuctuations in the ﬁgure and emphasize the physics),we show the binwise average of 100 independent data sets, each with 10 years exposure ofICAL. The blue bands in Fig. 8 correspond to the reconstructed oscillation valley. Notethat a good reconstruction of the oscillation valley would be an evidence for the ﬁdelity ofthe ICAL detector.To get the alignment of the oscillation valley, we ﬁt the U/D distribution for µ − or µ + independently with a functional form F ( E rec µ , cos θ rec µ ) = Z + N cos (cid:18) m cos θ rec µ E rec µ (cid:19) , (5.1)where Z , N , and m are the independent parameters to be determined from the ﬁttingof U/D distributions. The parameter Z quantiﬁes the minimum depth of the ﬁtted U/Dratio, N is the normalization constant, whereas m is the slope of the oscillation valley.– 18 – , SI - m U/D for - - - - - rec - m q cos ( G e V ) r e c - m E , SI - m U/D for , SI + m U/D for - - - - - rec + m q cos ( G e V ) r e c + m E , SI + m U/D for

Figure 8 : The distributions of U/D ratios in ( E rec µ , cos θ rec µ ) plane, for µ − and µ + events,in left and right panels, respectively, This is the average of 100 independent simulated datasets, each for 10 years exposure at ICAL, with ε µτ = 0. The white solid and dashed linescorrespond to the contours with ﬁtted U/D ratio 0.4 and 0.5, respectively, obtained afterthe ﬁtting of oscillation valley with the function F as given in Eq. 5.1. The black solidand dashed lines correspond to log [ L rec µ /E rec µ ] = 2.2 and 3.0. We consider normal massordering, and the benchmark oscillation parameters given in Table 2.Since more than 95% of the events at ICAL are contributed from the ν µ → ν µ and ¯ ν µ → ¯ ν µ survival probabilities, we expect that the function F , that resembles Eq. 2.5, would besuitable for ﬁtting the oscillation valley in the ( E rec µ , cos θ rec µ ) plane. Here, the slope m canbe used directly to calibrate ∆ m . The parameters N and Z also contain informationabout the atmospheric mixing parameters, however, they cannot be connected easily tothe determinations of these parameters.Figure 8 also shows the contour lines for the ﬁtted U/D ratio of 0.4 and 0.5. The contourlines for µ + are seen to reproduce the data better than those for µ − . The reason behindthis is the matter eﬀects, which appear in µ − data since the mass ordering is assumed to beNO here. Note that the ﬁtting function F is designed for only the vacuum oscillation, andwe focus only on the region where vacuum oscillation is expected to dominate. In orderto be safe from matter eﬀects, we remove the events with log [ L rec µ /E rec µ ] > .

0, wherethe matter eﬀects are signiﬁcant. We also discard the events with log [ L rec µ /E rec µ ] < . < . µ − and µ + events.We calibrate for ∆ m by ﬁtting the U/D ratio of 1000-year MC data with input | ∆ m | values in the range (1 . − . × − eV , and obtaining the corresponding value of m . The statistical ﬂuctuations are estimated by simulating 100 independent data sets of10 years for a given input value of | ∆ m | = 2 . × − eV , and ﬁtting for the U/D ratio– 19 – m - · m D Fixed param., no syst.Varied param., with syst.

NO90% C.L. - m

24 26 28 30 32 34 36 38 m - · m D Fixed param., no syst.Varied param., with syst.

NO90% C.L. + m Figure 9 : The blue stars and the blue lines indicate the calibration curves of ∆ m against m , obtained using 1000-year MC data with normal mass ordering. The graybands represent the 90% C.L. allowed ranges of m (vertical bands), and hence of ∆ m (horizontal bands), with a given input value of ∆ m = 2.46 × − eV , for an exposure of10 years at ICAL. The results for µ − and µ + events are shown separately, in the left andright panels, respectively. The light (dark) gray bands show the ranges when the errorson other oscillation parameters and the impact of systematic uncertainties are included(excluded), as discussed in Sec. 4.2. For ﬁxed-parameter case, we use the benchmarkoscillation parameters given in Table 2.independently with Eq. 5.1. The 100 values of m thus obtained provides the expecteduncertainties on ∆ m . In Fig. 9, the gray bands show the the 90% C.L. ranges expectedto be inferred by this method in 10 years. The ﬁgure also shows the (small) deteriorationin the 90% range of ∆ m due to the incorporation of the present uncertainties in theoscillation parameters, and all the ﬁve systematic errors, as discussed in Sec. 4.2. Withsystematic uncertainties and error in other oscillation parameters, we get the expected90% C.L. allowed range for | ∆ m | from µ − events as (2.24 – 2.82) × − eV and from µ + events as (2.25 – 2.75) × − eV .An advantage of the method proposed in this section over the original proposal in[29] is that this is a one-step ﬁtting procedure as opposed to the two-step ﬁtting proceduretherein. Note that the results obtained by this method are also slightly better than the onesobtained in [29], as the shape of the valley is explicitly ﬁtted to. Moreover, we can extendthis method to measure other characteristics of the oscillation valley, such as identifyingthe signature of NSI, which we shall discuss in the next section. So far, we have discussed the characteristics of the oscillation valley in the absence of NSI( ε µτ = 0), and the measurement of | ∆ m | in this scenario. In this section, we explore thefeatures of the oscillation valley in the presence of NSI (non-zero ε µτ ). We have already– 20 – - - - - - rec - m q cos ( G e V ) r e c - m E = 0.1 tm e , - m U/D for - - - - - rec + m q cos ( G e V ) r e c + m E = 0.1 tm e , + m U/D for - - - - - rec - m q cos ( G e V ) r e c - m E = -0.1 tm e , - m U/D for - - - - - rec + m q cos ( G e V ) r e c + m E = -0.1 tm e , + m U/D for

Figure 10 : The distributions of U/D ratio in ( E rec µ , cos θ rec µ ) plane with non-zero | ε µτ | (+0.1 in upper panels and -0.1 in lower panels) as expected at ICAL in 10 years. Leftand right panels are for µ − and µ + events respectively. The white solid and dashed linescorrespond to contours with ﬁtted U/D ratio 0.4 and 0.5, respectively, obtained after ﬁttingof the oscillation valley with function F α (see Eq. 6.1). The black solid and dashed linescorrespond to log [ L rec µ /E rec µ ] = 2.2 and 3.1. We consider normal mass ordering, and thebenchmark oscillation parameters given in Table 2.observed in Sec. 2.2 that in the presence of non-zero ε µτ , the oscillation valley in theoscillogram has a curved shape that depends on the extent of ε µτ , its sign, and whetherone is observing the µ − or µ + events. Therefore, the curvature of oscillation valley isexpected to be an useful parameter for the determination of ε µτ . Fig. 10 presents the distributions of U/D ratios of µ − and µ + events separately in the( E rec µ , cos θ rec µ ) plane, for two non-zero ε µτ values ( i.e. ± . ε µτ = 0), the valley iscurved in the presence of NSI. The direction of bending of the oscillation valley is decidedby sign of ε µτ as well as whether it is µ − or µ + . An important point to note is that thenature of curvature of the oscillation valley, that we observed in the ( E ν , cos θ ν ) plane withneutrino variables (see Fig. 2) with NSI, is preserved in the reconstructed oscillation valleyin the plane of reconstructed muon observables ( E rec µ , cos θ rec µ ). This is due to the excellentenergy and angular resolutions of the ICAL detector for muons in the multi-GeV energyrange.We retrieve the information on the bending of the oscillation valley by ﬁtting with anappropriate function, which is a generalization of Eq. 5.1. We introduce an additional freeparameter α that characterizes the the curvature of the oscillation valley. The function wepropose for ﬁtting the oscillation valley in the presence of ε µτ is F α ( E rec µ , cos θ rec µ ) = Z α + N α cos (cid:18) m α cos θ rec µ E rec µ + α cos θ rec µ (cid:19) , (6.1)where Z α , N α , m α , and α are the free parameters to be determined from the ﬁtting of theU/D ratio in the ( E rec µ , cos θ rec µ ) plane. Equation 6.1 is inspired by the survival probabilityas given in Eq. 2.5. In Eq. 6.1, the parameter α enters multiplied with cos θ rec µ , where onefactor of cos θ rec µ includes the L dependence of oscillation. The other factor of cos θ rec µ is totake into account the baseline dependence of the matter potential. While this is a rathercrude approximation, it is suﬃcient to provide us a way of calibrating ε µτ , as will be seenlater in the section. In Eq. 6.1, the parameters Z α and N α ﬁx the depth of valley, while m α and α determine the orientation and curvature of oscillation valley, respectively. Therefore, m α and α together account for the combined eﬀect of the mass-squared diﬀerence ∆ m and the NSI parameter ε µτ in oscillations.As in the analysis in Sec. 5, we restrict the range of log [ L rec µ /E rec µ ] to 2.2 – 3.1, toremove the unoscillated part ( < .

2) as well as the region with signiﬁcant matter eﬀects( > . F α works well to reproduce the curved nature of the oscillationvalley with non-zero ε µτ . We see that the oscillation valley for µ − with positive (negative) ε µτ closely resembles that for µ + with negative (positive) ε µτ . This happens due to thedegeneracy in the sign of ε µτ and the sign of matter potential V CC – these appear in thecombination ε µτ V CC in the survival probability expression in Eq. 2.5. One can see thatfor µ + ( µ − ) with ε µτ = 0 . − . µ + ( µ − ) with ε µτ = 0 . − . ε µτ from oscillation valley We saw in the previous section that the shape of the oscillation valley in the muon observ-ables is well-approximated by the function in Eq. 6.1, with the two parameters α and m α .Therefore, we propose to use these two parameters for retrieving ε µτ . At ICAL, since the– 22 – - - - aD - - - - - - a m D = - . tm e = - . tm e = . tm e = . tm e NO90% C.L.

Figure 11 : The blue curve with colored circles shows the calibration of ε µτ in the(∆ α, ∆ m α ) plane, obtained from 1000-yr MC. The black stars are the values of (∆ α, ∆ m α )obtained from 100 independent data sets, each for 10 years data of ICAL, when ε µτ = 0.The gray band represents the expected 90% C.L. region in this plane for a 10-year exposureat ICAL, when ε µτ = 0. The part of the blue calibration line covered by the gray bandgives the expected 90% C.L. bounds on ε µτ . We consider normal mass ordering, and thebenchmark oscillation parameters given in Table 2.data on µ − and µ + events are available separately, we can have two independent ﬁts forthe U/D ratios in µ − and µ + events. This would give us independent determinations ofthe parameter pairs ( α − , m − α ) and ( α + , m + α ), for µ − and µ + events, respectively. One canthen get the calibration of ε µτ in the plane of α and m α , for µ − and µ + events, separately.For other experiments that do not have charge identiﬁcation capability, there would be asingle calibration curve of ε µτ with α and m α , based on the U/D ratio of muon eventswithout charge information.Building up on our experience in the analysis of the oscillation dip in the presence ofNSI (see Sec. 4.1), in this analysis of oscillation valley, we use the following observablesthat quantify the µ − vs. µ + contrast in the values of the parameters α and m α :∆ α = α − − α + , ∆ m α = m α − − m α + . (6.2)The blue line in Fig. 11 shows the calibration of ε µτ in (∆ α, ∆ m α ) plane, obtained byusing the 1000-year MC sample. It may be observed that the calibration of ε µτ is almostindependent of ∆ m α , and ε µτ is almost linearly proportional to ∆ α .– 23 –n order to estimate the extent to which ε µτ may be constrained at ICAL with anexposure of 10 years, we simulate the statistical ﬂuctuations by generating 100 independentdata sets, each corresponding to 10-year exposure and ε µτ = 0. The black dots in Fig. 11are the values of ∆ α and ∆ m α determined from each of these data sets. The ﬁgure showsthe results, with the benchmark oscillation parameters in Table 2.It is observed that with a 10-year data, the statistical ﬂuctuations lead to a strongcorrelation between the values of ∆ α and ∆ m α , and we have to employ a two-dimensionalﬁtting procedure. We ﬁt the scattered points in (∆ α , ∆ m α ) plane with a straight line.The best ﬁtted straight line is shown by the red color in Fig. 11, which corresponds tothe average of these points, and indeed passes through the calibration point obtained from1000-year MC events for ε µτ = 0. The 90% C.L. limit on ε µτ is then obtained based onthe measure of the perpendicular distance of every point from the red line. The regionwith the gray band in Fig. 11 contains 90% of the points that are closest to the red line.The expected 90% C.L. bounds on ε µτ can then be determined from the part of the bluecalibration line covered by the gray band. Any value of the pair (∆ α, ∆ m α ) that lies outsidethe gray band would be a smoking gun signature for nonzero ε µτ at 90% C.L. Figure 11indicates that if ε µτ is indeed zero, our method, as applied at the ICAL detector, wouldbe able to constrain it to the range − . < ε µτ < .

021 at 90% C.L., with an exposureof 500 kt · yr.We explore the possible dilution of our results due to the uncertainties in oscillation pa-rameters and the ﬁve systematic errors, following the same procedure described in Sec. 4.2.This yields the expected 90% C.L. bounds as − . < ε µτ < . Atmospheric neutrino experiments are suitable for exploring the non-standard interactionsof neutrinos due to the accessibility to high energies ( E ν ) and long baselines ( L ν ) throughthe Earth, for which the eﬀect of NSI-induced matter potential is large. The NSI matterpotential produced in the interactions of neutrino and matter fermions may change theoscillation pattern. The multi-GeV neutrinos, in particular, oﬀer the advantage to accessthe interplay of NSI and Earth matter eﬀects. In this paper, we have explored the impactof the NSI parameter ε µτ on the dip and valley patterns in neutrino oscillation probabili-ties, and proposed a new method for extracting information on ε µτ from the atmosphericneutrino data.Detectors like ICAL, that have good reconstruction capability for the energy, direction,and charge of muons, can reproduce the neutrino oscillation dip and valley patterns intheir reconstructed muon data, from µ − and µ + events separately. The observables chosento reproduce these patterns are the ratios of upward-going (U) and downward-going (D)events, which minimize the eﬀects of systematic uncertainties. We analyze the U/D ratioin bins of reconstructed muon observables L rec µ /E rec µ , and in the two-dimensional plane( E rec µ , cos θ rec µ ). We show that nonzero ε µτ shifts the locations of oscillation dips in thereconstructed µ − and µ + data in opposite directions in the L rec µ /E rec µ distributions. At the– 24 –ame time, nonzero ε µτ manifests itself in the curvatures of the oscillation valleys in the( E rec µ , cos θ rec µ ) plane, this bending being in opposite directions for µ − and µ + events. Thedirection of shift in the dip location and the direction of the bending of the valley alsodepend on the sign of ε µτ , as well as on the mass ordering. We present our results in thenormal mass ordering scenario, the analysis for the inverted mass ordering scenario will beexactly analogous.In the oscillation dip analysis, we introduce a new observable ∆ d , corresponding tothe diﬀerence in the locations of dips in µ − and µ + event distribution. We show thecalibration of ε µτ against ∆ d using 1000-yr MC sample at ICAL, and demonstrate that itis almost independent of the actual value of ∆ m . We incorporate the eﬀects of statisticalﬂuctuations corresponding to 10-yr simulated data, uncertainties in the neutrino oscillationparameters, and major systematic errors, by simulating multiple data sets. Including allthese features, it is possible to constrain the NSI parameter in the range − . < ε µτ < .

024 at 90% C.L. with 500 kt · yr exposure.In the oscillation valley analysis, we propose a function F α that quantiﬁes the cur-vature ( α ) and orientation ( m α ) of the oscillation valley in the ( E rec µ , cos rec µ ) plane. Thisanalysis may be performed for the µ − and µ + events separately, leading to independentmeasurements of ε µτ . However, we further ﬁnd that the diﬀerences in these parameters forthe µ − and µ + events are quantities that are sensitive to ε µτ , and show the calibration for ε µτ in the (∆ α, ∆ m α ) plane using 1000-yr MC sample at ICAL. Incorporating the eﬀectsof statistical ﬂuctuations, uncertainties in the neutrino oscillation parameters, and majorsystematic errors, the NSI parameter can be constrained in the range − . < ε µτ < . · yr exposure.A special case of the function F α , with α = 0, corresponds to a valley without anycurvature. It is found that the method of ﬁtting for this function F to the U/D ratio inthe ( E rec µ , cos rec µ ) plane leads to a more robust and precise measurement of the value of∆ m than the method proposed earlier [29].Our analysis procedure is quite straightforward and transparent, and is able to capturethe physics of the NSI eﬀects on muon observables quite eﬃciently. It builds upon the majorfeatures of the impact of ε µτ on the oscillation dip, and identiﬁes the contrast in the shiftin dip locations of µ − and µ + as the key observable that would be sensitive to NSI. Itexploits the major pattern in the complexity of the two-dimensional event distributionsin the ( E rec µ , cos rec µ ) plane by ﬁtting to a function approximately describing the curvatureof the valley. The two complementary methods of using the oscillation dips and valleysprovide a check of the robustness of the results. The shift of the dip location and thebending of the oscillation valley, and survival of these eﬀects in the reconstructed muondata, are deep physical insights obtained through the analysis in this paper.Note that muon charge identiﬁcation plays a major role in our analysis. In a detectorlike ICAL, the measurement of ε µτ is possible independently in both the µ − and µ + chan-nels. The oscillation dips and valleys in µ − and µ + event samples have diﬀerent locationsand shapes. As a result, the information available in the dip and valley would be severelydiluted in the absence of muon charge identiﬁcation. But more importantly, the quantitiessensitive to ε µτ in a robust way turn out to be the ones that reﬂect the diﬀerences in– 25 –he properties of oscillation dips and valleys in µ − and µ + events. For example, ∆ d isthe observable identiﬁed by us, whose calibration against ε µτ is almost independent of theactual value of ∆ m . For ∆ d to be measured at a detector, the muon charge identiﬁca-tion capability is necessary. The data from ICAL will thus provide the measurement of acrucial observable that is not possible at other large atmospheric neutrino experiments likeSuper-K or DeepCore/ IceCube. At future long-baseline experiment like DUNE, that willhave access to 1 – 10 GeV neutrino energy and separate data sets for neutrino and antineu-trino, our methodology to probe NSI using oscillation-dip location can also be adopted,by replacing the U/D ratio with the ratio of observed number of events and the predictednumber of events without oscillations.We expect that further exploration of the features of the oscillation dips and valleys,observable at ICAL-like experiments, would enable us to probe neutrino properties in novelways. Acknowledgments

This work is performed by the members of the INO-ICAL collaboration to suggest a newapproach to probe Non-Standard Interactions (NSI) using atmospheric neutrino experi-ments which can distinguish between neutrino and antineutrino events. We thank S. UmaSankar, Srubabati Goswami, and D. Indumathi for their useful and constructive commentson our work. A.K. would like to thank the organizers of the virtual Mini-Workshop on Neu-trino Theory at Fermilab, Chicago, USA during 21st to 23rd September, 2020 for giving anopportunity to present the main results of this work. We would like to thank the Depart-ment of Atomic Energy (DAE), Govt. of India for ﬁnancial support. S.K.A. is supportedby the DST/INSPIRE Research Grant [IFA-PH-12] from the Department of Science andTechnology (DST), India and the Young Scientist Project [INSA/SP/YSP/144/2017/1578]from the Indian National Science Academy.

References [1]

Particle Data Group

Collaboration, P. Zyla et al.,

Review of Particle Physics , PTEP (2020), no. 8 083C01.[2] L. Wolfenstein,

Neutrino Oscillations in Matter , Phys.Rev.

D17 (1978) 2369–2374.[3] J. Valle,

Resonant Oscillations of Massless Neutrinos in Matter , Phys. Lett. B (1987)432–436.[4] M. M. Guzzo, H. Nunokawa, P. C. de Holanda, and O. L. G. Peres,

On the massless ’just-so’solution to the solar neutrino problem , Phys. Rev.

D64 (2001) 097301, [ hep-ph/0012089 ].[5] P. Huber and J. W. F. Valle,

Nonstandard interactions: Atmospheric versus neutrino factoryexperiments , Phys. Lett.

B523 (2001) 151–160, [ hep-ph/0108193 ].[6] A. M. Gago, M. M. Guzzo, P. C. de Holanda, H. Nunokawa, O. L. G. Peres, V. Pleitez, andR. Zukanovich Funchal,

Global analysis of the postSNO solar neutrino data for standard andnonstandard oscillation mechanisms , Phys. Rev.

D65 (2002) 073012, [ hep-ph/0112060 ]. – 26 –

7] F. J. Escrihuela, M. Tortola, J. W. F. Valle, and O. G. Miranda,

Global constraints onmuon-neutrino non-standard interactions , Phys. Rev.

D83 (2011) 093002,[ arXiv:1103.1366 ].[8] M. Gonzalez-Garcia, M. Maltoni, and J. Salvado,

Testing matter eﬀects in propagation ofatmospheric and long-baseline neutrinos , JHEP (2011) 075, [ arXiv:1103.4365 ].[9] T. Ohlsson,

Status of non-standard neutrino interactions , Rept. Prog. Phys. (2013)044201, [ arXiv:1209.2710 ].[10] M. C. Gonzalez-Garcia and M. Maltoni, Determination of matter potential from globalanalysis of neutrino oscillation data , JHEP (2013) 152, [ arXiv:1307.3092 ].[11] O. G. Miranda and H. Nunokawa, Non standard neutrino interactions: current status andfuture prospects , New J. Phys. (2015), no. 9 095002, [ arXiv:1505.06254 ].[12] Y. Farzan and M. Tortola, Neutrino oscillations and Non-Standard Interactions , Front.inPhys. (2018) 10, [ arXiv:1710.09360 ].[13] P. S. Bhupal Dev et al., Neutrino Non-Standard Interactions: A Status Report , in

NTNWorkshop on Neutrino Non-Standard Interactions St Louis, MO, USA, May 29-31, 2019 ,2019. arXiv:1907.00991 .[14] J. Kopp, M. Lindner, T. Ota, and J. Sato,

Non-standard neutrino interactions in reactor andsuperbeam experiments , Phys. Rev.

D77 (2008) 013007, [ arXiv:0708.0152 ].[15] M. C. Gonzalez-Garcia and M. Maltoni,

Atmospheric neutrino oscillations and new physics , Phys. Rev.

D70 (2004) 033010, [ hep-ph/0404085 ].[16] I. Mocioiu and W. Wright,

Non-standard neutrino interactions in the mu–tau sector , Nucl.Phys.

B893 (2015) 376–390, [ arXiv:1410.6193 ].[17] I. Esteban, M. Gonzalez-Garcia, M. Maltoni, I. Martinez-Soler, and J. Salvado,

UpdatedConstraints on Non-Standard Interactions from Global Analysis of Oscillation Data , JHEP (2018) 180, [ arXiv:1805.04530 ].[18] J. Salvado, O. Mena, S. Palomares-Ruiz, and N. Rius, Non-standard interactions withhigh-energy atmospheric neutrinos at IceCube , JHEP (2017) 141, [ arXiv:1609.03450 ].[19] IceCube

Collaboration, M. Aartsen et al.,

Search for Nonstandard Neutrino Interactionswith IceCube DeepCore , Phys. Rev.

D97 (2018), no. 7 072009, [ arXiv:1709.07079 ].[20]

Super-Kamiokande

Collaboration, G. Mitsuka et al.,

Study of Non-Standard NeutrinoInteractions with Atmospheric Neutrino Data in Super-Kamiokande I and II , Phys.Rev.

D84 (2011) 113008, [ arXiv:1109.1889 ].[21]

ICAL

Collaboration, S. Ahmed et al.,

Physics Potential of the ICAL detector at theIndia-based Neutrino Observatory (INO) , Pramana (2017), no. 5 79, [ arXiv:1505.07380 Resonant ampliﬁcation of neutrino oscillations in matter andsolar neutrino spectroscopy , Nuovo Cim. C9 (1986) 17–26.[24] S. Mikheev and A. Y. Smirnov, Resonance Ampliﬁcation of Oscillations in Matter andSpectroscopy of Solar Neutrinos , Sov.J.Nucl.Phys. (1985) 913–917.[25] S. Choubey, A. Ghosh, T. Ohlsson, and D. Tiwari, Neutrino Physics with Non-StandardInteractions at INO , JHEP (2015) 126, [ arXiv:1507.02211 ]. – 27 –

26] A. Khatun, S. S. Chatterjee, T. Thakore, and S. Kumar Agarwalla,

Enhancing sensitivity tonon-standard neutrino interactions at INO combining muon and hadron information , Eur.Phys. J. C (2020), no. 6 533, [ arXiv:1907.02027 ].[27] P. Huber, M. Lindner, and W. Winter, Superbeams versus neutrino factories , Nucl.Phys.

B645 (2002) 3–48, [ hep-ph/0204352 ].[28] G. L. Fogli, E. Lisi, A. Marrone, D. Montanino, A. Palazzo, and A. M. Rotunno,

Solarneutrino oscillation parameters after ﬁrst KamLAND results , Phys. Rev.

D67 (2003) 073002,[ hep-ph/0212127 ].[29] A. Kumar, A. Khatun, S. K. Agarwalla, and A. Dighe,

From oscillation dip to oscillationvalley in atmospheric neutrino experiments , arXiv:2006.14529 .[30] Super-Kamiokande

Collaboration, Y. Ashie et al.,

A Measurement of atmosphericneutrino oscillation parameters by Super-Kamiokande I , Phys.Rev.

D71 (2005) 112005,[ hep-ex/0501064 ].[31]

IceCube

Collaboration, M. Aartsen et al.,

Energy Reconstruction Methods in the IceCubeNeutrino Telescope , JINST (2014) P03009, [ arXiv:1311.4767 ].[32] N. Dash, V. M. Datar, and G. Majumder, A Study on the time resolution of Glass RPC , arXiv:1410.5532 .[33] A. D. Bhatt, V. M. Datar, G. Majumder, N. K. Mondal, Pathaleswar, and B. Satyanarayana, Improvement of time measurement with the INO-ICAL resistive plate chambers , JINST (2016), no. 11 C11001.[34] A. Gaur, A. Kumar, and M. Naimuddin, Study of timing response and charge spectra of glassbased Resistive Plate Chamber detectors for INO-ICAL experiment , JINST (2017), no. 03C03081.[35] M. Sajjad Athar, M. Honda, T. Kajita, K. Kasahara, and S. Midorikawa, Atmosphericneutrino ﬂux at INO, South Pole and Pyh´asalmi , Phys.Lett.

B718 (2013) 1375–1380,[ arXiv:1210.5154 ].[36] M. Honda, M. S. Athar, T. Kajita, K. Kasahara, and S. Midorikawa,

Atmospheric neutrinoﬂux calculation using the nrlmsise-00 atmospheric model , Phys. Rev. D (Jul, 2015)023004.[37] D. Casper, The Nuance neutrino physics simulation, and the future , Nucl.Phys.Proc.Suppl. (2002) 161–170, [ hep-ph/0208030 ].[38] P. de Salas, D. Forero, S. Gariazzo, P. Mart´ınez-Mirav´e, O. Mena, C. Ternes, M. T´ortola, andJ. Valle, , arXiv:2006.11237 .[39] I. Esteban, M. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, The fate of hints:updated global analysis of three-ﬂavor neutrino oscillations , JHEP (2020) 178,[ arXiv:2007.14792 ].[40] F. Capozzi, E. Di Valentino, E. Lisi, A. Marrone, A. Melchiorri, and A. Palazzo, Globalconstraints on absolute neutrino masses and their ordering , Phys. Rev. D (2017), no. 9096014, [ arXiv:2003.08511