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

We propose a new approach to explore the neutral-current non-standard neutrino interactions (NSI) in atmospheric neutrino experiments using oscillation dips and valleys in reconstructed muon observables, at a detector like ICAL that can identify the muon charge. We focus on the flavor-changing NSI parameter $\varepsilon_{\mu\tau}$, which has the maximum impact on the muon survival probability in these experiments. We show that non-zero $\varepsilon_{\mu\tau}$ shifts the oscillation dip locations in $L/E$ distributions of the up/down event ratios of reconstructed $\mu^-$ and $\mu^+$ in opposite directions. We introduce a new variable $\Delta d$ representing the difference of dip locations in $\mu^-$ and $\mu^+$, which is sensitive to the magnitude as well as the sign of $\varepsilon_{\mu\tau}$, and is independent of the value of $\Delta m^2_{32}$. We further note that the oscillation valley in the ($E$, $\cos \theta$) plane of the reconstructed muon observables bends in the presence of NSI, its curvature having opposite signs for $\mu^-$ and $\mu^+$. We demonstrate the identification of NSI with this curvature, which is feasible for detectors like ICAL having excellent muon energy and direction resolutions. We illustrate how the measurement of contrast in the curvatures of valleys in $\mu^-$ and $\mu^+$ can be used to estimate $\varepsilon_{\mu\tau}$. Using these proposed oscillation dip and valley measurements, the achievable precision on $|\varepsilon_{\mu\tau}|$ at 90% C.L. is about 2% with 500 kt$\cdot$yr exposure. The effects of statistical fluctuations, systematic errors, and uncertainties in oscillation parameters have been incorporated using multiple sets of simulated data. Our method would provide a direct and robust measurement of $\varepsilon_{\mu\tau}$ in the multi-GeV energy range.


Introduction and Motivation
The phenomenon of neutrino oscillations has now been well-established, from measurements at the solar, atmospheric, reactor as well as accelerator experiments with short and long baselines [1]. Neutrino oscillations are the consequences of mixing among different neutrino flavors and non-degenerate values of neutrino masses, with at least two neutrino masses nonzero. However, neutrinos are massless in the Standard Model (SM) of particle physics, and therefore, physics beyond the SM (BSM) is needed to accommodate nonzero neutrino masses and mixing. Many models of BSM physics suggest new non-standard interactions (NSI) of neutrinos [2], which may affect neutrino production, propagation, and detection in a given experiment. The possible impact of these NSI at neutrino oscillation experiments have been studied extensively, for example see Refs. [3][4][5][6][7][8][9][10][11][12][13]. In this paper, we propose a new method for identifying NSI at atmospheric neutrino experiments which can reconstruct the energy, direction, as well as charge of the muons produced in the detector due to chargedcurrent interactions of ν µ andν µ .
We shall focus on the neutral-current NSI, which may be described at low energies via effective four-fermion dimension-six operators as [2] L NC−NSI = −2 √ 2 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 strength of 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 flavors. The subscript C ∈ {L, R} represents the chiral projection operator P L = (1 − γ 5 )/2 or P R = (1 + γ 5 )/2. The hermiticity of the interactions demands ε f βα,C = (ε f αβ,C ) * . The effective NSI parameter relevant for the neutrino propagation through matter is (1.2) where N f is the number density of fermion f . In the approximation of a neutral and isoscalar Earth, the number densities of electrons, protons, and neutrons are identical, which implies N u ≈ N d ≈ 3N e . Thus, ε αβ ≈ ε e αβ + 3 ε u αβ + 3 ε d αβ .
(1. 3) In the presence of NSI, the modified effective Hamiltonian for neutrino propagation through matter is where U is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix that parametrizes neutrino mixing, and ∆m 2 ij ≡ m 2 i − m 2 j are the mass-squared differences. The quantity V CC ≡ √ 2G F N e is the effective matter potential due to the coherent elastic forward scattering of neutrinos with electrons inside the medium via the SM gauge boson W . For antineutrinos, V CC → −V CC , U → U * , and ε αβ → ε * αβ . In the present study, we suggest a novel approach to unravel the presence of flavorchanging neutral current NSI parameter ε µτ , via its effect on the propagation of multi-GeV atmospheric neutrinos and antineutrinos through Earth matter. We choose ε µτ as the NSI parameter to focus on, since it can significantly affect the evolution of the mixing angle θ 23 and the mass-squared difference ∆m 2 32 in matter [14], which in turn would alter the survival probabilities of atmospheric muon neutrinos and antineutrinos substantially [15,16]. In general, ε µτ can be complex, i.e. ε µτ ≡ |ε µτ |e iφµτ . However, in the disappearance channels ν µ → ν µ andν µ →ν µ that dominate in our analysis 1 , ε µτ appears only as |ε µτ | cos φ µτ at the leading order [14]. Thus, a complex phase only changes the effective value of ε µτ to a real number between −|ε µτ | and +|ε µτ | at the leading order. We take advantage of this observation, and restrict ourselves to real values of ε µτ in the range −0.1 ≤ ε µτ ≤ 0.1. From the arguments given above, this covers the whole range of complex values of ε µτ with |ε µτ | ≤ 0.1.
Based on the global neutrino data analysis in [17], where the possible contributions to NSI from only up and down quark have been included, the bound on |ε µτ | turns out to be |ε µτ | < 0.07 at 2σ C.L.. The existing bounds on ε µτ from various neutrino oscillation experiments are listed in Table 1. An important point to note is the energies of neutrinos involved 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 the Super-K experiment is more efficient in the sub-GeV energy range [20].
The proposed 50 kt magnetized Iron Calorimeter (ICAL) detector at the India-based Neutrino Observatory (INO) [21,22] would be sensitive to multi-GeV neutrinos, since it can efficiently 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 would also be in a unique position to detect any interplay between the matter effects and NSI. Another important feature is that ICAL can explore physics in neutrinos and antineutrinos separately, unlike Super-K and IceCube/DeepCore. The studies of physics potential of ICAL for detecting NSI have shown that, using the reconstructed muon momentum, it would be possible to obtain a bound of |ε µτ | < 0.015 at 90% C.L. [25] with 500 kt·yr exposure. When information on the reconstructed hadron energy in each event is also included, the expected 90% C.L. bound improves to |ε µτ | < 0.010 [26]. The results in [25,26] are obtained using a χ 2 analysis with the pull method [15,27,28] The wide range of neutrino energies and baselines available in atmospheric neutrino experiments offer an opportunity to study the features of "oscillation dip" and "oscillation valley" in the reconstructed µ − and µ + observables, as demonstrated in [29]. These features can be clearly identified in the ratios of upward-going and downward-going muon events at ICAL. If the muon neutrino disappearance is solely due to non-degenerate masses and non-zero mixing of neutrinos, then the valley in both µ − and µ + is approximately a straight line. The location of the dip, and the alignment of the valley, can be used to determine ∆m 2 31 [29]. These features may undergo major changes in the presence of NSI, and can act as 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 oscillation dip feature, we note that non-zero ε µτ shifts the oscillation dip location in opposite directions for µ − and µ + . We demonstrate that this opposite shift in dip location due to NSI 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 a whole new analysis methodology to extract the information on ε µτ using the dip locations. For this, we define a new variable exploiting the contrast between the shifts in reconstructed dip locations, which eliminates the dependence of our results on the actual value of ∆m 2 32 . For the oscillation valley feature, we notice that the valley becomes curved in the presence of 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 of the 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 the presence of non-zero ε µτ , and discuss the shifts in the dip locations as well as the bending of the oscillation valleys in the survival probabilities of ν µ andν µ . In Sec. 3, we investigate the survival of these two striking features in the reconstructed L rec µ /E rec µ distributions and in 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 the shifts of dip locations in µ − and µ + . This variable leads to the calibration of ε µτ , and is used to find the expected bound on ε µτ from a 500 kt·yr exposure of ICAL. In Sec. 5, we come up with a new procedure for determining the alignment of the oscillation valley and estimating the value of ∆m 2 32 in the absence of NSI, which we extend to the NSI analysis in Sec. 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 the valley analysis at ICAL. Finally, in Sec. 7, we summarize our findings and offer concluding remarks.

Oscillation Dip and Valley in the Presence of NSI
In the limit of θ 13 → 0, and the approximation of one mass scale dominance scenario [∆m 2 21 L/(4E) ∆m 2 32 L/(4E)] and constant matter density, the survival probability of ν µ when traveling a distance L ν is given by [15] and where β ≡ sgn(∆m 2 32 ). That is, β = +1 for normal mass ordering (NO), and β = −1 for inverted mass ordering (IO).
In the limit of maximal mixing (θ 23 = 45 • ), Eq. 2.1 reduces to the following simple expression [16]: We further find that the correction due to the deviation of θ 23 from maximality, χ ≡ θ 23 − π/4, is second order in the small parameter χ, and hence can be neglected. Here, one notes that the NSI parameter ε µτ primarily modifies the wavelength of neutrino oscillations. It also comes multiplied with V CC , which increases at higher baselines inside the Earth's matter. Thus, the modification in ν µ survival probabilities due to NSI varies with the baseline L ν (or cos θ ν ).
2.1 Effect of ε µτ 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 fixed values of cos θ ν (i.e. cos θ ν = −0.4, −0.8), and for three values of the NSI parameter ε µτ (i.e. ε µτ = +0.1, 0.0, −0.1). The other benchmark values of the oscillation parameters are given in Table 2.  The benchmark values of oscillation parameters that we use in the analysis. Normal mass ordering corresponds to m 1 < m 2 < m 3 . Figure 1 also indicates the sensitivity ranges for the atmospheric neutrino experiments Super-K, ICAL, and IceCube. Note that these ranges are different for cos θ ν = −0.4 and −0.8, which correspond to L ν around 5100 km and 10200 km, respectively While calculating these L ν /E ν ranges, the energy ranges chosen are those for which the detectors perform very well: we use the E ν range of 100 MeV -5 GeV for Super-K [30], 1-25 GeV for ICAL [21], and 100 GeV -10 PeV for IceCube [31]. Note that these ranges are only indicative. The following observations may be made from the figure.
• For log 10 [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 difference do not develop for such small value of L ν /E ν , i.e. ∆m 2 32 L ν /E ν 1, the disappearance of ν µ is possible due to the L ν ε µτ V CC term, which can be high for large baselines (for core passing neutrinos, V CC is large). For example, log 10 [L ν /E ν ] = 1 may correspond to L ν = 5000 km and E ν = 500 GeV, so that ∆m 2 32 L ν /E ν ≈ 0.025. However, for this baseline, the average density of the Earth is ρ ≈ 3.9 g/cc, and hence for ε µτ = 0.1, we have L ν ε µτ V CC ≈ 0.93, which takes the oscillation probability away from unity. The effect of NSI at such small L ν /E ν is energy-independent, and can be seen at detectors like IceCube [18] due to its better performance at high energy.
• As we go to higher value of log 10 [L ν /E ν ] (> 2), the term containing neutrino mass splitting becomes comparable to the Lε µτ V CC term in Eq. 2.5, and the competition  Figure 1: The survival probabilities of ν µ andν µ as functions of log 10 (L ν /E ν ) in left and right panels, respectively. The black lines correspond to ε µτ = 0 (SI), whereas red and blue lines are with ε µτ = 0.1, and −0.1, respectively. The neutrino direction taken in upper (lower) panels is cos θ ν = −0.4 (−0.8). The horizontal bands shown here are the indicative 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 due to the mass splitting are the dominating contribution, the change in oscillation wavelength due to non-zero ε µτ results in a shift of the oscillation minima towards left or right side (lower or higher values of L ν /E ν , respectively), depending on the amplitude of ε µτ and its sign. The direction of shift in the dip location depends on whether it is neutrino or antineutrino (since the sign of V CC for them are opposite), and on the neutrino mass ordering. This effect on the shift in the dip location is discussed in next section in detail. This region is relevant for Super-K, INO, and most of the long-baseline experiments.
We now discuss the modification of the oscillation dip due to non-zero ε µτ . First, using the approximate expression of ν µ survival probability from Eq. 2.5, we obtain the value of L ν /E ν at the first oscillation minimum (or the dip): The above expression may be written by expressing V CC in terms of the line-averaged matter density 2 ρ and taking care of units, . (2.8) Here a positive sign in the denominator corresponds to neutrinos, whereas a negative sign is for antineutrinos. This approximation is useful for understanding the shift of dip position with 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 lower value of L ν /E ν than the SI. On the other hand, with negative ε µτ , the oscillation dip would occur at a higher value of L ν /E ν . These two features can be seen clearly in the left panels of Fig. 1. For antineutrinos and the same mass ordering, since the matter potential has the opposite sign, the shift of oscillation dip is in the opposite direction as compared that for neutrinos, given the same ε µτ . Expanding Eq. 2.8 to first order in ε µτ , one can write (2.9) where β ≡ sgn(∆m 2 32 ), as defined 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 neutrinos and antineutrinos, as well as for the two mass orderings. We have checked that for E ν < 25 GeV and typical line-averaged density of the Earth, indeed ε µτ < 0.1 is small enough for the above approximation to hold.
2.2 Effect of ε µτ on the oscillation valley in the (E ν , cos θ ν ) plane In Fig. 2, we show the oscillograms of survival probabilities for ν µ andν µ , in the plane of neutrino energy and cosine of neutrino zenith angle (cos θ ν ), with NO as the neutrino mass ordering. We show only upward-going neutrinos (cos θ ν < 0), since for downwardgoing neutrinos (cos θ ν > 0), the baseline L ν is very small 3 and neutrino oscillations do 2 We can write VCC approximately as a function of matter density ρ VCC ≈ ±7.6 × Ye × ρ 10 14 g/cm 3 eV . (2.7) Here, Ye = Ne Np+Nn is the relative electron number density. In an electrically neutral and isoscalar medium, Ye = 0.5. The positive (negative) sign is for neutrino (antineutrino). 3 The zenith angle θν is related to the baseline Lν by 10) Figure 2: The oscillogram for the survival probabilities in (E ν , cos θ ν ) plane for neutrino and antineutrino in left and right panels, respectively, with normal mass ordering and the benchmark 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 neutrinos are 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 affect the νµ oscillation probabilities because R h d.
not develop, making P µµ ≈ 1.0. The differences observed among the upper, middle, and lower panels are due to different values of ε µτ . A major impact of NSI is observed on the feature corresponding to the first oscillation minima, seen in the figure as a broad blue/black diagonal band. We refer to this broad band as the oscillation valley [29]. In the SI case, the oscillation valley appears like a triangle with straight edges, while in the presence of NSI, its edges seem to acquire a curvature. The sign of this curvature is opposite for neutrinos and antineutrinos, as can be seen by comparing the left and right plots of upper and lower panels of Fig. 2.
The modification in the oscillation valley due to NSI may be explained from the following relation between E ν and cos θ ν at the first oscillation minima. We rewrite Eq. 2.6 with L ν ≈ 2R| cos θ ν |: Taking care of units and expressing V CC in terms of line-averaged constant matter density ρ, the above expression may be written as Here, the negative sign in the denominator corresponds to neutrinos, whereas the positive sign is for antineutrinos. Putting ε µτ = 0 in Eq. 2.12 gives the condition of the first oscillation minima to be E ν = |(5.08/π) · ∆m 2 32 · R[km] · cos θ ν |, and this relation clearly shows 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 ε µτ > 0 in Eq. 2.12, then for neutrinos, 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 oscillation valley bends in the opposite direction, for both neutrinos and antineutrinos, as shown in the bottom panels of Fig. 2. In the inverse mass ordering scenario, the bending of the oscillation valley will be in the direction opposite to that in the normal mass ordering scenario.
From Eqs. 2.9 and 2.12, the shift in the oscillation minima and the bending in the oscillation valley are in opposite directions for normal and inverted mass ordering. The effect of ε µτ therefore depends crucially on the mass ordering. In this paper, we present our analysis in the scenario where the mass ordering is known to be NO. The analysis for IO may be performed in an exactly analogous manner. Our results are presented with the exposure corresponding to 10 years of ICAL data-taking. It is expected that the mass ordering will be determined from the neutrino experiments (including ICAL) by this time.

Impact of NSI on the Event Distribution at ICAL
While the effect of NSI on the neutrino and antineutrino survival probabilities was discussed in the last section, it is important to confirm whether the features present at the probability level can survive in the observables at a detector, and if they can be reconstructed. Here is where the response of the detector plays a crucial role. Neutrinos cannot be directly detected in an experiment, however the charged leptons produced from their chargedcurrent interactions in the detector contains information about their energy, direction, and flavor, which can be recovered depending on the nature of the detector.
The upcoming ICAL detector at the India-based Neutrino Observatory (INO) will be composed of 151 layers of magnetized iron plates; the 50 kt iron acting as the target, and around 30,000 glass resistive plate chambers as detection elements. The magnetic field of strength 1.5 Tesla and the time resolution of less than 1 ns [32][33][34] will enable ICAL to distinguish µ − from µ + events in the multi-GeV energy range. In this study, we generate the charged-current interactions of ν µ andν µ , similar to Ref. [29], using neutrino and antineutrino fluxes calculated at the Theni site by Honda et.al. [35,36] and the neutrino generator NUANCE [37]. The reconstructed energy (E rec µ ) and zenith angle (θ rec µ ) of muons produced in neutrino and antineutrino interactions will be used in the analysis. The quantity L rec µ associated with the reconstructed muon direction cos θ rec µ is defined (similar to Eq. 2.10) as Note that L rec µ is simply a proxy observable, and there is no need to associate it directly with the distance traveled by the incoming neutrino.    3 presents the log 10 [L rec µ /E rec µ ] 4 distribution of the upward-going µ − and µ + events (U, cos θ rec µ < 0) with SI (ε µτ = 0) and with SI + NSI (ε µτ = ±0.1) expected with 10 yr exposure at ICAL. We include the statistical fluctuation in the number of events by simulating 100 independent data sets, and calculating the mean and root-mean-square deviation We consider normal mass ordering, and the benchmark oscillation parameters given in Table 2. in each bin. For log 10 [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 10 [L rec µ /E rec µ ] in the range 0-4. The downward-going events (D, cos θ rec µ > 0) are not affected significantly by oscillations or by additional NSI interactions. The log 10 [L rec µ /E rec µ ]-distribution of µ − and µ + events is identical to the one shown in the upper panel of Fig. [4] in Ref. [29], and we do not repeat it here. Fig. 3 clearly shows that the log 10 [L rec µ /E rec µ ] range of 2.4-3.0 would be important for NSI studies, since non-zero ε µτ has the largest effect in this range. With positive ε µτ , the number of µ − events is lower as compared to that of SI case, whereas the number of µ + events is higher. If ε µτ is negative, then the modifications of the number of µ − and µ + events are the other way around. As a consequence, if a detector is not able to distinguish between the µ − and µ + events, the difference between SI and NSI would get diluted substantially. The charge identification capability of the magnetized ICAL detector would be crucial in exploiting this observation.

Distributions in the (E rec
µ , cos θ rec µ ) plane In order to study the effect 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 16 non-uniform E rec µ bins in the range 1 -25 GeV, and the whole range of −1 ≤ cos θ rec µ ≤ 1 is divided into 20 uniform bins. The reconstructed E rec µ and cos θ rec µ are binned with the same binning scheme as used in [29], and is shown in Table 4. For demonstrating event distributions, we scale the 1000-year MC sample to an exposure of 10 years, and the 1 -1 -1 -1 1 1   Figure 4: The (E rec µ , cos θ rec µ ) distributions of difference of events between including nonzero NSI and only SI (ε µτ = 0) expected with 500 kt·yr of ICAL. Left and right panels are for µ − and µ + events, respectively, whereas upper and lower panels correspond to µτ = 0.1 and −0.1, respectively. We consider normal mass ordering, and the benchmark oscillation parameters given in Table 2.
difference in the number of events with SI + NSI (|ε µτ | = 0.1) and the number of events with SI, for µ − and µ + events are shown in Fig. 4 Table 4: The binning scheme adopted for E rec µ and cos θ rec µ of µ − and µ + events.
From this figure, we can make the following observations.
• The events in the bins with with E rec µ > 5 GeV and cos θ rec µ < −0.2 are particularly useful for NSI searches at ICAL. This is true for both µ − and µ + events.
• There is almost no difference in the number of events between SI+NSI and SI at E rec µ < 5 GeV. The events with reconstructed muon energy less than 5 GeV do not seem to be sensitive to NSI due to small NSI-induced matter effects: L ν ε µτ V CC ∆m 2 /(2E ν ) at small energies.
• As we go to higher energies, the mass-squared difference-induced neutrino oscillations die down, while the NSI-induced matter effect term L ν ε µτ V CC increases. Therefore, large effects of NSI are observed at high energies.
• The effect of NSI is also larger at higher baselines, since in this case, neutrinos pass through inner and outer cores that have very high matter densities (around 11 g/cm 3 and 13 g/cm 3 , respectively), and hence higher values of L ν ε µτ V CC .
• In many of the (E rec µ , cos θ rec µ ) bins, the difference in the number of events in µ − and µ + events has the opposite sign. Therefore, muon charge information is a crucial ingredient in the search for NSI. We have also seen this feature in the log 10 [L rec µ /E rec µ ] distributions of µ − and µ + events in Sec. 3.1.

Identifying NSI through the Shift in Oscillation Dip Location
To study the effect 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 event ratio U/D is defined for cos θ rec µ < 0, as [29] U/D(E rec µ , cos θ rec where the numerator corresponds to the number of upward-going events and the denominator is the number of downward-going events, with energy E rec µ and direction | cos θ rec µ |. We associate the U/D ratio with cos θ rec µ < 0 as well as the corresponding L rec µ of upward-going events. 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 ν µ → ν µ oscillation channel [29]. One advantage of the use of the U/D ratio would be to nullify the effects of systematic uncertainties like flux normalization, cross sections, overall detection efficiency, and energy dependent tilt error, as can be seen later. Note that the up-down symmetry of the detector geometry and detector response play an important role in this.  Figure 5: The log 10 (L rec µ /E rec µ ) distributions of ratio of upward-going and downwardgoing µ − (left panel) and µ + (right panel) events using 10 years data of ICAL. The black lines correspond to SI (ε µτ = 0), whereas red and blue lines are for ε µτ = 0.1 and -0.1, respectively. The statistical fluctuations shown by shaded boxes are the root-mean-square deviation of 100 independent distributions of U/D ratio, each for 10 years, whereas the mean of these distributions are shown by solid lines. We consider normal mass ordering, and the benchmark oscillation parameters given in Table 2. modification in the U/D ratio due to the presence of non-zero ε µτ (±0.1), particularly in the location of the oscillation dip. The dip shifts towards left or right (i.e. to lower or higher values of L rec µ /E rec µ , respectively) from that of the SI case with non-zero ε µτ . This shift for µ − and µ + events is observed to be in opposite directions, as is expected from the 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, the oscillation dip shifts to higher values of L rec µ /E rec µ . For ε µτ < 0, this shift is in the opposite directions. Moreover, as the oscillation dip shifts towards the higher value of L rec µ /E rec µ , the dip becomes shallower due to the effect of rapid oscillation at high L ν /E ν . This effect is visible in both, µ − and µ + . It can be easily argued that the modification in oscillation dips of µ − and µ + due to non-zero ε µτ is the reflection of how the survival probabilities of neutrino and antineutrino changes in the presence of non-zero ε µτ , as expected from Eq. 2.6 in Sec. 2.1.
Note that for a neutrino detector that is blind to the charge of muons, the dip itself would get diluted since different average inelasticities of neutrino and antineutrino events at these energies would lead to slightly different dip locations in µ − and µ + distributions. Moreover, the shift of oscillation dip towards opposite directions in µ − and µ + events due to non-zero ε µτ would further contribute to the dilution of NSI effects on the dip location. On the other hand, a detector like ICAL that has the charge identification capability, not only 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. The calibration of ε µτ with difference in dip location of µ − and µ + . Note that for −0.1 < ε µτ < −0.07, the dip location is simply the minimum U/D ratio since the fitting is not stable for these values. We consider normal mass ordering, and the benchmark oscillation parameters given in Table 2.

A novel variable ∆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 location in either of these distributions may be used to estimate ε µτ , independently. We have already developed a dip-identification algorithm [29], where the dip location is not simply the lowest value of U/D ratio, but takes into account information from the surrounding bins that have a U/D ratio below a threshold value. The algorithm essentially identifies the cluster of contiguous bins that have a U/D ratio lower than any surrounding bins, and fits these bins with a quadratic function whose minimum corresponds to the dip. The value of log 10 [L rec µ /E rec µ ] corresponding to this minimum is termed as the dip location. We denote the dip location for µ − and µ + events as d − and d + , respectively.
We 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 curves of ε µτ with d − (left panel) and d + (right panel), for three different values of ∆m 2 32 . The straight line nature of the calibration curves can be understood qualitatively using Eq. 2.6. For a majority of events (E rec µ < 25 GeV), the ε µτ V CC E ν term is much smaller than ∆m 2 32 term, and therefore the linear expansion as shown in Eq. 2.9 is valid. In Eq. 2.9, the dominating ∆m 2 32 dependence of the dip location is through the first term. We can get rid of this dependence by introducing a new variable which is expected to be proportional to ε µτ , and can be used to calibrate for ε µτ . As may be seen in the bottom panel of Fig. 6, this indeed removes the dominant ∆m 2 32 -dependence in the calibration of ε µτ . Note that some ∆m 2 32 -dependence still survives, which may be seen in the slightly different slopes of the calibration curves for different ∆m 2 32 values. However, this dependence is clearly negligible, as it appears in Eq. 2.9 as a multiplicative correction to ε µτ . We have thus obtained a calibration for ε µτ that is almost independent of ∆m 2 32 . Note that ε µτ = 0 does not imply ∆d = 0, since different inelasticities and matter effects in the neutrino and antineutrino channels mean that the dip locations in these two locations are not identical even in the absence of NSI.

Constraints on ε µτ from the measurement of ∆d
We saw in the last section that one can infer the value of ε µτ in the observed data from the calibration curve between ∆d and ε µτ , independent of the actual value of ∆m 2 32 . We obtain the calibration curve using the 1000-year MC-simulated sample, and show it with the solid blue line in the panels of Fig. 7.
For estimating the expected constraints on ε µτ , we simulate 100 independent sets of data, each for 10 years exposure of ICAL with the benchmark oscillation parameters as given in Table 2, and ε µτ = 0. In Fig. 7, the range indicated along the x-axis by solid lines is the expected 90% C.L. range of the measured value of ∆d, while the range indicated along the y-axis by the dark gray band is the expected 90% C. L. limit on ε µτ . From the figure, we find that one can expect to constrain ε µτ to −0.024 < ε µτ < 0.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 2 32 , we expect that even if the actual value of ∆m 2 32 were not exactly known, the results will not change. We have confirmed this by generating 100 independent sets of 10-yr exposure, where the ∆m 2 32 values used as inputs follow a Gaussian distribution ∆m 2 32 = (2.46 ± 0.03) × 10 −3 eV 2 , in accordance to the range allowed by the global fit [38][39][40] of available neutrino data.
We also checked the impact of uncertainties in the values of the other oscillation parameters. We first simulated 100 statistically independent unoscillated data sets. Then for each of these data sets, we take 20 random choices of oscillation parameters, according to the Gaussian distributions   Table 2. sin 2 2θ 12 = 0.855 ± 0.020 , sin 2 2θ 13 = 0.0875 ± 0.0026 , sin 2 θ 23 = 0.50 ± 0.03 , (4.4) guided by the present global fit [41]. We keep δ CP = 0, since its effect on ν µ survival probability is known to be highly suppressed in the multi-GeV energy range [14]. This procedure effectively enables us to consider the variation of our results over 2000 different combinations of oscillation parameters, to take into account the effect of their uncertainties.
We do not observe any significant dilution of results due to the variation over oscillation parameters. This is expected, since (i) solar oscillation parameters ∆m 2 21 and θ 12 do not contribute significantly to the ν µ survival probability in the multi-GeV range, (ii) the reactor mixing angle θ 13 is already measured to a high precision, (iii) the mixing angle θ 23 does not affect 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 2 32 uncertainty, as shown in Sec. 4.1.
In addition, we take into account the five major systematic uncertainties in the neutrino fluxes and cross sections that are used in the standard ICAL analyses [21]. These five uncertainties are (i) 20% in overall flux normalization, (ii) 10% in cross sections, (iii) 5% in the energy dependence, (iv) 5% in the zenith angle dependence, and (iii) 5% in overall systematics. For each of the 2000 simulated data sets, we modify the number of events in each (E rec µ , cos θ rec µ ) bin as 5) where N (0) is the theoretically predicted number of events, and E 0 = 2 GeV. Here (δ 1 , δ 2 , δ 3 , δ 4 , δ 5 ) is an ordered set of random numbers, generated separately for each simulated data set, with the Gaussian distributions centered around zero and the 1σ widths given by (20%, 10%, 5%, 5%, 5%). The normalization uncertainties and energy tilt uncertainty are expected to be canceled in the U/D ratio, while the zenith angle distribution uncertainty affects the upward-going and downward-going events differently, and hence would be expected to affect the U/D ratio. We explicitly check this and indeed find this to be true. When the oscillation parameter uncertainties and all the five systematic uncertainties mentioned above are included, it is found that the method based on the shift in the dip locations may constrain ε µτ to −0.025 < ε µτ < 0.024 at 90% C.L.. The results denoting the effects of these uncertainties have been shown in Fig. 7.

A New Method for Determining ∆m 2 32 in the Absence of NSI
The reconstruction of the oscillation valley using the distributions of the upward-going and downward-going muons in (E rec µ , cos θ rec µ ) plane is discussed in detail in Ref. [29]. It has been shown that the alignment of the valley can provide a measurement of ∆m 2 32 . In this paper, we use an alternative procedure for getting the alignment of the oscillation valley, which is more robust and more well-motivated as compared to what was used in Ref. [29], and determine the expected precision that the ICAL detector can achieve on the atmospheric oscillation parameter ∆m 2 32 , using this method. Fig. 8 shows the distributions of the ratio of upward-going and downward going events at ICAL in the (E rec µ , cos θ rec µ ) plane for 10 years of ICAL data. The left and right panels show the U/D ratios for µ − and µ + events, respectively. For the sake of demonstration (in order to reduce the statistical fluctuations in the figure and emphasize the physics), we show the binwise average of 100 independent data sets, each with 10 years exposure of ICAL. The blue bands in Fig. 8 correspond to the reconstructed oscillation valley. Note that a good reconstruction of the oscillation valley would be an evidence for the fidelity of the ICAL detector.
To get the alignment of the oscillation valley, we fit the U/D distribution for µ − or µ + independently with a functional form where Z 0 , N 0 , and m 0 are the independent parameters to be determined from the fitting of U/D distributions. The parameter Z 0 quantifies the minimum depth of the fitted U/D ratio, N 0 is the normalization constant, whereas m 0 is the slope of the oscillation valley.  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 data sets, each for 10 years exposure at ICAL, with ε µτ = 0. The white solid and dashed lines correspond to the contours with fitted U/D ratio 0.4 and 0.5, respectively, obtained after the fitting of oscillation valley with the function F 0 as given in Eq. 5.1. The black solid and dashed lines correspond to log 10 [L rec µ /E rec µ ] = 2.2 and 3.0. We consider normal mass ordering, 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 0 , that resembles Eq. 2.5, would be suitable for fitting the oscillation valley in the (E rec µ , cos θ rec µ ) plane. Here, the slope m 0 can be used directly to calibrate ∆m 2 32 . The parameters N 0 and Z 0 also contain information about the atmospheric mixing parameters, however, they cannot be connected easily to the determinations of these parameters. Figure 8 also shows the contour lines for the fitted U/D ratio of 0.4 and 0.5. The contour lines for µ + are seen to reproduce the data better than those for µ − . The reason behind this is the matter effects, which appear in µ − data since the mass ordering is assumed to be NO here. Note that the fitting function F 0 is designed for only the vacuum oscillation, and we focus only on the region where vacuum oscillation is expected to dominate. In order to be safe from matter effects, we remove the events with log 10 [L rec µ /E rec µ ] > 3.0, where the matter effects are significant. We also discard the events with log 10 [L rec µ /E rec µ ] < 2.2, since these would have barely oscillated. The latter cut also gets rid of most of the events near horizon, where the distance traveled by the neutrinos has large errors. In order to minimize the data from the bins with clearly large fluctuations, we use another cut on the maximum value of the U/D ratio. This cut is taken to be U/D < 0.9 for both µ − and µ + events.
We calibrate for ∆m 2 32 by fitting the U/D ratio of 1000-year MC data with input  For fixed-parameter case, we use the benchmark oscillation parameters given in Table 2. independently with Eq. 5.1. The 100 values of m 0 thus obtained provides the expected uncertainties on ∆m 2 32 . In Fig. 9, the gray bands show the the 90% C.L. ranges expected to be inferred by this method in 10 years. The figure also shows the (small) deterioration in the 90% range of ∆m 2 32 due to the incorporation of the present uncertainties in the oscillation parameters, and all the five systematic errors, as discussed in Sec. 4.2. With systematic uncertainties and error in other oscillation parameters, we get the expected 90% C.L. allowed range for |∆m 2 32 | from µ − events as (2.24 -2.82)×10 −3 eV 2 and from µ + events as (2.25 -2.75)×10 −3 eV 2 .
An advantage of the method proposed in this section over the original proposal in [29] is that this is a one-step fitting procedure as opposed to the two-step fitting procedure therein. Note that the results obtained by this method are also slightly better than the ones obtained in [29], as the shape of the valley is explicitly fitted to. Moreover, we can extend this method to measure other characteristics of the oscillation valley, such as identifying the signature of NSI, which we shall discuss in the next section.

Identifying NSI through Oscillation Valley
So far, we have discussed the characteristics of the oscillation valley in the absence of NSI (ε µτ = 0), and the measurement of |∆m 2 32 | in this scenario. In this section, we explore the features of the oscillation valley in the presence of NSI (non-zero ε µτ ). We have already  µ /E rec µ ] = 2.2 and 3.1. We consider normal mass ordering, and the benchmark oscillation parameters given in Table 2. observed in Sec. 2.2 that in the presence of non-zero ε µτ , the oscillation valley in the oscillogram has a curved shape that depends on the extent of ε µτ , its sign, and whether one is observing the µ − or µ + events. Therefore, the curvature of oscillation valley is expected to be an useful parameter for the determination of ε µτ .
6.1 Curvature of oscillation valley as a signature of NSI 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.1). For the sake of demonstration (in order to reduce the statistical fluctuations in the figure and emphasize the physics), we show the binwise average of 100 independent data sets, each with 10 years exposure of ICAL. In contrast to the linear nature of oscillation valley with SI (ε µτ = 0), the valley is curved in the presence of NSI. The direction of bending of the oscillation valley is decided by sign of ε µτ as well as whether it is µ − or µ + . An important point to note is that the nature of curvature of the oscillation valley, that we observed in the (E ν , cos θ ν ) plane with neutrino variables (see Fig. 2) with NSI, is preserved in the reconstructed oscillation valley in the plane of reconstructed muon observables (E rec µ , cos θ rec µ ). This is due to the excellent energy and angular resolutions of the ICAL detector for muons in the multi-GeV energy range.
We retrieve the information on the bending of the oscillation valley by fitting with an appropriate function, which is a generalization of Eq. 5.1. We introduce an additional free parameter α that characterizes the the curvature of the oscillation valley. The function we propose for fitting the oscillation valley in the presence of ε µτ is where Z α , N α , m α , and α are the free parameters to be determined from the fitting of the U/D ratio in the (E rec µ , cos θ rec µ ) plane. Equation 6.1 is inspired by the survival probability as given in Eq. 2.5. In Eq. 6.1, the parameter α enters multiplied with cos 2 θ rec µ , where one factor of cos θ rec µ includes the L dependence of oscillation. The other factor of cos θ rec µ is to take into account the baseline dependence of the matter potential. While this is a rather crude approximation, it is sufficient to provide us a way of calibrating ε µτ , as will be seen later in the section. In Eq. 6.1, the parameters Z α and N α fix the depth of valley, while m α and α determine the orientation and curvature of oscillation valley, respectively. Therefore, m α and α together account for the combined effect of the mass-squared difference ∆m 2 32 and the NSI parameter ε µτ in oscillations.
As in the analysis in Sec. 5, we restrict the range of log 10 [L rec µ /E rec µ ] to 2.2 -3.1, to remove the unoscillated part (< 2.2) as well as the region with significant matter effects (> 3.1). The cut on the U/D value is applied at 0.9. In Fig. 10, the contours of the fitted U/D ratio of 0.4 and 0.5, shown with white solid and dashed lines, respectively, clearly indicate that the function F α works well to reproduce the curved nature of the oscillation valley with non-zero ε µτ . We see that the oscillation valley for µ − with positive (negative) ε µτ closely resembles that for µ + with negative (positive) ε µτ . This happens due to the degeneracy in the sign of ε µτ and the sign of matter potential V CC -these appear in the combination ε µτ V CC in the survival probability expression in Eq. 2.5. One can see that for µ + (µ − ) with ε µτ = 0.1 (−0.1), the contour with the fitted U/D ratio of 0.4 is absent. This is because the valley is shallower in these cases. This feature corresponds to a similar one observed in Fig. 5, where the oscillation dip is shallower for µ + (µ − ) with ε µτ = 0.1 (−0.1).

Constraints on ε µτ from oscillation valley
We saw in the previous section that the shape of the oscillation valley in the muon observables 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  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 exposure at ICAL, when ε µτ = 0. The part of the blue calibration line covered by the gray band gives the expected 90% C.L. bounds on ε µτ . We consider normal mass ordering, and the benchmark oscillation parameters given in Table 2. data on µ − and µ + events are available separately, we can have two independent fits for the U/D ratios in µ − and µ + events. This would give us independent determinations of the parameter pairs (α − , m − α ) and (α + , m + α ), for µ − and µ + events, respectively. One can then get the calibration of ε µτ in the plane of α and m α , for µ − and µ+ events, separately. For other experiments that do not have charge identification capability, there would be a single calibration curve of ε µτ with α and m α , based on the U/D ratio of muon events without charge information.
Building up on our experience in the analysis of the oscillation dip in the presence of NSI (see Sec. 4.1), in this analysis of oscillation valley, we use the following observables that quantify the µ − vs. µ + contrast in the values of the parameters α and m α : The blue line in Fig. 11 shows the calibration of ε µτ in (∆α, ∆m α ) plane, obtained by using the 1000-year MC sample. It may be observed that the calibration of ε µτ is almost independent of ∆m α , and ε µτ is almost linearly proportional to ∆α.
In order to estimate the extent to which ε µτ may be constrained at ICAL with an exposure of 10 years, we simulate the statistical fluctuations by generating 100 independent data sets, each corresponding to 10-year exposure and ε µτ = 0. The black dots in Fig. 11 are the values of ∆α and ∆m α determined from each of these data sets. The figure shows the results, with the benchmark oscillation parameters in Table 2.
It is observed that with a 10-year data, the statistical fluctuations lead to a strong correlation between the values of ∆α and ∆m α , and we have to employ a two-dimensional fitting procedure. We fit the scattered points in (∆α, ∆m α ) plane with a straight line. The best fitted straight line is shown by the red color in Fig. 11, which corresponds to the average of these points, and indeed passes through the calibration point obtained from 1000-year MC events for ε µτ = 0. The 90% C.L. limit on ε µτ is then obtained based on the measure of the perpendicular distance of every point from the red line. The region with 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 blue calibration line covered by the gray band. Any value of the pair (∆α, ∆m α ) that lies outside the gray band would be a smoking gun signature for nonzero ε µτ at 90% C.L. Figure 11 indicates that if ε µτ is indeed zero, our method, as applied at the ICAL detector, would be able to constrain it to the range −0.022 < ε µτ < 0.021 at 90% C.L., with an exposure of 500 kt·yr.
We explore the possible dilution of our results due to the uncertainties in oscillation parameters and the five systematic errors, following the same procedure described in Sec. 4.2. This yields the expected 90% C.L. bounds as −0.024 < ε µτ < 0.020, thus keeping them almost unchanged.

Summary and Concluding remarks
Atmospheric neutrino experiments are suitable for exploring the non-standard interactions of neutrinos due to the accessibility to high energies (E ν ) and long baselines (L ν ) through the Earth, for which the effect of NSI-induced matter potential is large. The NSI matter potential produced in the interactions of neutrino and matter fermions may change the oscillation pattern. The multi-GeV neutrinos, in particular, offer the advantage to access the interplay of NSI and Earth matter effects. In this paper, we have explored the impact of the NSI parameter ε µτ on the dip and valley patterns in neutrino oscillation probabilities, and proposed a new method for extracting information on ε µτ from the atmospheric neutrino 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 in their reconstructed muon data, from µ − and µ + events separately. The observables chosen to reproduce these patterns are the ratios of upward-going (U) and downward-going (D) events, which minimize the effects of systematic uncertainties. We analyze the U/D ratio in 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 the reconstructed µ − and µ + data in opposite directions in the L rec µ /E rec µ distributions. At the same 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. The direction of shift in the dip location and the direction of the bending of the valley also depend on the sign of ε µτ , as well as on the mass ordering. We present our results in the normal mass ordering scenario, the analysis for the inverted mass ordering scenario will be exactly analogous.
In the oscillation dip analysis, we introduce a new observable ∆d, corresponding to the difference in the locations of dips in µ − and µ + event distribution. We show the calibration of ε µτ against ∆d using 1000-yr MC sample at ICAL, and demonstrate that it is almost independent of the actual value of ∆m 2 32 . We incorporate the effects of statistical fluctuations corresponding to 10-yr simulated data, uncertainties in the neutrino oscillation parameters, and major systematic errors, by simulating multiple data sets. Including all these features, it is possible to constrain the NSI parameter in the range −0.025 < ε µτ < 0.024 at 90% C.L. with 500 kt·yr exposure.
In the oscillation valley analysis, we propose a function F α that quantifies the curvature (α) and orientation (m α ) of the oscillation valley in the (E rec µ , cos rec µ ) plane. This analysis may be performed for the µ − and µ + events separately, leading to independent measurements of ε µτ . However, we further find that the differences in these parameters for the µ − 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 effects of statistical fluctuations, uncertainties in the neutrino oscillation parameters, and major systematic errors, the NSI parameter can be constrained in the range −0.022 < ε µτ < 0.021 at 90% C.L. with 500 kt·yr exposure.
A special case of the function F α , with α = 0, corresponds to a valley without any curvature. It is found that the method of fitting for this function F 0 to the U/D ratio in the (E rec µ , cos rec µ ) plane leads to a more robust and precise measurement of the value of ∆m 2 32 than the method proposed earlier [29]. Our analysis procedure is quite straightforward and transparent, and is able to capture the physics of the NSI effects on muon observables quite efficiently. It builds upon the major features of the impact of ε µτ on the oscillation dip, and identifies the contrast in the shift in dip locations of µ − and µ + as the key observable that would be sensitive to NSI. It exploits the major pattern in the complexity of the two-dimensional event distributions in the (E rec µ , cos rec µ ) plane by fitting to a function approximately describing the curvature of the valley. The two complementary methods of using the oscillation dips and valleys provide a check of the robustness of the results. The shift of the dip location and the bending of the oscillation valley, and survival of these effects in the reconstructed muon data, are deep physical insights obtained through the analysis in this paper.
Note that muon charge identification plays a major role in our analysis. In a detector like ICAL, the measurement of ε µτ is possible independently in both the µ − and µ + channels. The oscillation dips and valleys in µ − and µ + event samples have different locations and shapes. As a result, the information available in the dip and valley would be severely diluted in the absence of muon charge identification. But more importantly, the quantities sensitive to ε µτ in a robust way turn out to be the ones that reflect the differences in the properties of oscillation dips and valleys in µ − and µ + events. For example, ∆d is the observable identified by us, whose calibration against ε µτ is almost independent of the actual value of ∆m 2 32 . For ∆d to be measured at a detector, the muon charge identification capability is necessary. The data from ICAL will thus provide the measurement of a crucial observable that is not possible at other large atmospheric neutrino experiments like Super-K or DeepCore/ IceCube. At future long-baseline experiment like DUNE, that will have access to 1 -10 GeV neutrino energy and separate data sets for neutrino and antineutrino, 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 predicted number 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 novel ways.