Enhancing Sensitivity to Non-Standard Neutrino Interactions at INO combining muon and hadron information

The neutral current non-standard interactions (NSI's) of neutrino with matter fermions while propagating through long distances inside the Earth matter can give rise to the extra matter potentials apart from the standard MSW potential due to the $W$-mediated interactions in matter. In this paper, we explore the impact of flavor violating neutral current NSI parameter $\varepsilon_{\mu\tau}$ in the oscillation of atmospheric neutrino and antineutrino using the 50 kt magnetized ICAL detector at INO. We find that due to non-zero $\varepsilon_{\mu\tau}$, $\nu_\mu\rightarrow\nu_\mu$ and $\bar\nu_\mu\rightarrow\bar\nu_\mu$ transition probabilities get modified substantially at higher energies and longer baselines, where vacuum oscillation dominates. We estimate the sensitivity of the ICAL detector for various choices of binning schemes and observables. The most optimistic bound on $\varepsilon_{\mu\tau}$ that we obtain is $-0.01<\varepsilon_{\mu\tau}<0.01 $ at 90$\%$ C.L. using 500 kt$\cdot$yr exposure and considering $E_\mu,\, \cos\theta_\mu,\,E'_{\rm had}$ as observables in their ranges [1, 21] GeV, [-1, 1], and [0, 25] GeV respectively. For the first time we show that the charge identification capability of the ICAL detector is crucial to set stringent constraints on $\varepsilon_{\mu\tau}$. We also show that when we marginalize over $\varepsilon_{\mu\tau}$ in fit in its range of -0.1 to 0.1, the mass hierarchy sensitivity deteriorates by 10$\%$ to 20$\%$ depending on the analysis mode, and the precision measurements of atmospheric parameters remain quite robust at the ICAL detector.


Introduction
The observation of neutrino oscillation proves that the neutrinos are massive and mix with each other, which is the first direct evidence for physics beyond the Standard Model [1]. Currently, the three-flavor lepton mixing is confirmed and neutrino oscillation has entered the era of precision [2][3][4]. To explain the tiny masses of neutrino and large leptonic mixing angles, the extensions of the Standard Model allow some interactions which are not possible in the Standard Model. These interactions are termed as non-standard interactions (NSI's). The presence of NSI's in nature can have subdominant effect on the oscillation of neutrino and antineutrino. Therefore, the phenomenological consequences of NSI's in three-flavor mixing using neutrino oscillation experiments are interesting and widely studied by many authors in Refs. .
In this paper, we study the impact of neutral current (NC) non-standard interactions of neutrino which may arise when atmospheric neutrinos travel long distances inside the Earth. While NC NSI's affect neutrinos during their propagation, there are charged-current NSI's which may modify the neutrino fluxes at the production stage and interaction crosssection at the detection level. In this work, we only focus on the NC NSI's, and do not consider NSI's at production and detection level. In most of the cases, NSI's come into the picture as a low-energy manifestation of high-energy theory involving new heavy states. For a detailed discussion on this topic, see the reviews [40,54,56,57]. Therefore, at low energies, a neutral current NSI can be described by a four-fermion dimension-six operator [58], where G F is the Fermi coupling constant, ε C f αβ is the dimensionless parameter which represents the strength of NSI relative to G F , and ν α and ν β are the neutrino fields of flavor α and β respectively. Here, f denotes the matter fermions, electron (e), up-quark (u), and down-quark (d). Here, P L = (1 − γ 5 )/2, P R = (1 + γ 5 )/2, and the subscript C = L, R expresses the chirality of the f f current. Due to the hermiticity of the interaction, we have ε C f αβ = (ε C f βα ) * . The NSI's of neutrino with matter fermions can give rise to the additional matter induced potentials apart from the standard MSW potential due to the W -mediated interactions in matter as denoted by V CC . The total relative strength of the matter induced potential generated by the NC NSI's of neutrinos with all the matter fermions (ν α +f → ν β +f ) can be written in the following fashion, The quantity N f denotes the number density of matter fermion f in the medium. For antineutrino, V f → −V f and V CC → −V CC . In general, the total matter induced potential in presence of all the possible NC non-standard interactions of neutrino with matter fermions can be written as In the present study, we focus our investigation to flavor violating NSI parameter ε µτ , that is, we only allow ε µτ to be non-zero in our analysis, and assume all other NSI parameters to be zero. We also consider ε µτ to be real entertaining both of its negative and positive values. Since the atmospheric neutrino oscillation is mainly governed by ν µ → ν τ transition, it is expected that NSI parameter ε µτ would have significant impact on this oscillation channel, which in turn can modify ν µ → ν µ survival probability by a considerable amount. We can study this effect by observing the atmospheric neutrinos at the proposed 50 kt magnetized ICAL detector. If we will not see any significant deviation from the standard µ − and µ + event spectra at ICAL, we can use this fact to place tight constraints on NSI parameter ε µτ . This is the main theme of our present study.
This article is organized in the following fashion. In Sec. 2, we briefly review the existing bounds on NSI parameter ε µτ from various neutrino oscillation experiments. We discuss the possible modification in oscillation probabilities of neutrino and antineutrino due to non-zero ε µτ in Sec. 3. In Sec. 4, we present the expected total µ − and µ + events and their distributions as a function of reconstructed E µ and cos θ µ for the following three cases: (i) ε µτ = 0 (SM), (ii) ε µτ = 0.05, and (iii) ε µτ = −0.05 using 500 kt·yr exposure of the ICAL detector. In Sec. 5, we discuss the numerical procedure and various binning schemes that we use in our analysis. We present all the results of our study in Sec. 6 where we show the following: (a) The possible improvement in the sensitivity of the ICAL detector in constraining ε µτ due to the inclusion of events with E µ in range of 11 to 21 GeV in addition to the events that belong to the E µ in range of 1 to 11 GeV. (b) How much the limit on ε µτ can be improved by considering the information on reconstructed hadron energy (E had ) as an additional observable along with reconstructed variables E µ and cos θ µ on an event-by-event basis. (c) We show the advantage of having charge identification (CID) capability in the ICAL detector in placing competitive constraint on ε µτ . (d) We present the expected limits on ε µτ considering different exposures of the 50 kt ICAL detector. (e) We also explore the possible impact of non-zero ε µτ in determining the mass hierarchy and in the precision measurement of atmospheric oscillation parameters. We provide a summary of this study in Sec. 7.

Existing Limits on NSI Parameter ε µτ
There are existing constraints on the NSI parameter ε µτ from various neutrino oscillation experiments. The Super-Kamiokande collaboration performed an analysis of the atmospheric neutrino data collected during its phase-I and -II run assuming only NSI's with d-quarks [59]. The following bounds at 90% C.L. are obtained: Since N d = N u = 3N e for an electrically neutral and isoscalar Earth matter, the above constraints as obtained in Ref. [59] are actually on the NSI parameters ε αβ /3. Therefore, the above constraints at 90% C.L. can be interpreted as Recently, the authors in Ref. [60] considered the possibility of NSI's in µ-τ sector in the one-year high-energy through-going muon data of IceCube. In their analysis, they included various systematic uncertainties on both the atmospheric neutrino flux and detector properties, which they incorporated via several nuisance parameters. They obtained the following limits − 6.0 × 10 −3 < ε µτ < 5.4 × 10 −3 at 90% credible interval (C.I.). (2. 3) The IceCube-DeepCore collaboration also searched for NSI's involving ε µτ [61]. Using their three years of atmospheric muon neutrino disappearance data, they placed the following constraint at 90% confidence level − 6.7 × 10 −3 < ε µτ < 8.1 × 10 −3 .

(2.4)
A preliminary analysis to constrain the NSI parameters in context of the ICAL detector was performed in Ref. [55]. Using an exposure of 500 kt·yr and considering only muon momentum (E µ , cos θ µ ) as observable, the authors in Ref. [55] obtained the following bound − 0.015 (−0.027) < ε µτ < 0.015 (0.027) at 90 (3σ) C.L. with NH . (2.5) In the present study, we estimate new constraints on ε µτ considering the reconstructed hadron energy (E had ) as an additional observable along with the reconstructed E µ and cos θ µ on an event-by-event basis at the ICAL detector.      In upper panels of Fig. 1, we present the oscillograms for ν µ survival channel in the plane of cos θ ν vs. E ν considering NH. Here, we draw the oscillograms for three different cases: i) ε µτ = −0.05 (left panel), ii) ε µτ = 0.0 (the SM case, middle panel), and iii) ε µτ = 0.05 (right panel). The lower panel depicts the same but forν µ →ν µ oscillation channel. To prepare Fig. 1, we take the following benchmark values of vacuum oscillation parameters in three-flavor framework: sin 2 θ 23 = 0.5, sin 2 2θ 13 = 0.1, sin 2 θ 12 = 0.3, ∆m 2 21 = 7.5 × 10 −5 eV 2 , ∆m 2 eff = 2.4 × 10 −3 eV 2 , and δ CP = 0 • . We estimate the value of ∆m 2 31 from ∆m 2 eff 1 , where ∆m 2 eff has the same magnitude for NH and IH with positive and negative signs respectively. It is evident from the upper panels of Fig. 1 that in the presence of negative (see left panel) and positive (see right panel) non-zero values of ε µτ , ν µ survival probabilities get modified substantially at higher energies and longer baselines, where vacuum oscillation dominates. We observe the similar changes in case ofν µ oscillation probabilities as well (see lower panels).
Another broad feature which has been emerging from Fig. 1 is that the ν µ → ν µ oscillation probabilities with positive (negative) ε µτ as shown in upper right (left) panel are similar to that ofν µ →ν µ transition probabilities with negative (positive) ε µτ as can be seen from lower left (right) panel. We can understand this behaviour with the help of following approximate analytical expression of ν µ → ν µ transition probability. Assuming ∆m 2 21 L/4E → 0 and θ 13 = 0, ν µ → ν µ oscillation channel in the presence of non-zero ε µτ and under the constant matter density approximation takes the form [11,21] where and Since for antineutrino, V CC → −V CC , following Eq. 3.7, we can write and Eq. 3.8 and Eq. 3.9 explain the broad features in Fig. 1 that we mention above. 1 The effective mass splitting is related to ∆m 2 31 as follows [63,64] ∆m 2 eff = ∆m 2 31 − ∆m 2 21 cos 2 θ12 − cos δCP sin θ13 sin 2θ12 tan θ23 .    The lower panels are forνµ →νµ oscillation channel. Here, in all the panels, we assume NH.
To have a better look at the changes induced by non-zero ε µτ as compared to the SM case, we give Fig. 2 where we plot the difference in ν µ → ν µ survival channel considering the cases ε µτ = 0 (the SM case) and ε µτ = −0.05 (see top left panel). In top right panel, we present the same for the cases of ε µτ = 0 (the SM case) and ε µτ = 0.05. The lower panels are for antineutrinos. In all the panels, we see a visible difference in ν µ survival channel due to the presence of non-zero ε µτ as compared to the SM case (ε µτ = 0.0) at higher baselines with cos θ ν in the range −1 to −0.8. This range of cos θ ν corresponds to the baseline in the range ∼ 12700 km to 10000 km where neutrino and antineutrino mostly travel through inner and outer part of the Earth's core 2 and have access to large Earth matter effect. Also, we see a trend that the impact of NSI's is large at higher energies where the three-flavor oscillations are suppressed because the oscillation lengths (L osc = 4πE ∆m 2 ij ) are large at higher energies.
To have a quantitative estimate of the difference in ν µ → ν µ oscillation channel due to non-zero ε µτ as compared to ε µτ = 0 case, we use Eq. 3.7 and obtain the following Eq. 3.10 clearly suggests that the impact of NSI is proportional to both NSI induced matter potential (V CC ε µτ ) and baseline (L). We see in upper left and lower right panels of Fig. 2 that around E ∼ 18 GeV and cos θ ν ∼ − 0.9, P diff µµ approaches to zero suggesting that the impact of NSI is negligible. We observe the opposite behaviour in upper right and lower left panels, where around E ∼ 18 GeV and cos θ ν ∼ − 0.9, P diff µµ attains a quite large value of − 0.8 suggesting that the influence of NSI is significant there. Here, cos θ ν = − 0.9 corresponds to L = 11500 km for which the line-averaged constant Earth matter density according to the PREM [62] profile is 6.8 g/cm 3 . Therefore, the standard line-averaged Earth matter potential 3 for 11500 km baseline is V CC ∼ 2.6 × 10 −13 eV. Considering ε µτ = 0.05, we get For E = 18 GeV, L = 11500 km, and ∆m 2 31 = 2.36×10 −3 eV 2 (this mass-squared difference is obtained from Eq. 3.2 using benchmark values of oscillation parameters), Thus, for E = 18 GeV and L = 11500 km, from Eq. 3.10, we obtain and Eq. 3.14 and Eq. 3.15 confirm the observations regarding P diff µµ and P dif µμ in Fig. 2 that we mention above. We know that due to its CID capability, ICAL has an edge to resolve the issue of neutrino mass hierarchy (sign of ∆m 2 31 ) by observing the Earth matter effect in µ − and µ + events separately [66]. Similarly, the four panels in Fig. 2 suggest that the CID capability of ICAL can provide useful information to determine the sign of NSI parameter ε µτ for a particular choice of mass hierarchy. 3 The standard neutrino matter potential due to the W -mediated interactions with the ambient electrons can be written as a function of matter density ρ as follows: where Ye ( Ne Np+Nn ) is the relative number density. For electrically, neutral and isoscalar medium, Ne = Np = Nn, and therefore, Ye = 0.5.

Expected Events at ICAL with non-zero ε µτ
The Monte Carlo based neutrino event generator NUANCE [67] is used to simulate the CC interactions of ν µ andν µ in the ICAL detector. To generate events in NUANCE, we give a simple geometry of the ICAL detector with 150 alternate layers of iron and glass plates in each module. We have three such modules to account for the 50 kt ICAL detector. As far as the neutrino flux is concerned in generating the neutrino events in the present study, we use the flux as predicted at Kamioka 4 [70]. To reduce the statistical fluctuation, we generate the unoscillated CC neutrino and antineutrino events considering a very high exposure of 1000 years and 50 kt ICAL. Then, we implement various oscillation probabilities using the reweighting algorithm. Next, we fold the oscillated events with detector response for muon and hadron as described in Ref. [71,72]. In the present study, we assume that the ICAL particle reconstruction algorithms can separate the hits due to the hadron shower from the hits originating from a muon track with 100% efficiency. It means that whenever a muon is reconstructed, we consider all the other hits to be part of the hadronic shower in order to calibrate the hadron energy. It also implies that the neutrino event reconstruction efficiency is same as the muon reconstruction efficiency. Finally, the reconstructed µ − and µ + events are scaled down to the exposure of 10 years for 50 kt ICAL. Now, we present the expected µ − and µ + events for 500 kt·yr exposure of the ICAL detector assuming the SM case (ε µτ = 0) and ε µτ = ± 0.05. To estimate these event rates, we use the values of oscillation parameters as considered in Sec. 3 to draw the oscillograms.

Total Event Rates
First, we address the following question: can we see the signature of non-zero ε µτ in the total number of µ − and µ + events which will be collected at the ICAL detector over 10 years of running? To have an answer of this question, we estimate the total number of events for the following three cases: i) ε µτ = 0.05, ii) ε µτ = 0 (the SM case), and iii) ε µτ = −0.05. We present these numbers in Table 1 with NH and using 500 kt·yr exposure of the ICAL detector integrating over entire ranges of E µ , cos θ µ , and E had that we consider in our analysis. As far as the binning schemes are concerned, we use the low-energy (LE) and high-energy (HE) binning schemes 5 , and for both these binning schemes, we take the entire range of cos θ µ spanning over -1 to 1. The energy ranges for reconstructed E µ and E had are different in LE and HE binning schemes. For LE binning scheme, E µ ∈ [1,11] GeV and E had ∈ [0, 15] GeV. In case of HE binning scheme, E µ ∈ [1, 21] GeV, and E had ∈ [0, 25] GeV. When we increase the reconstructed muon energy from 11 GeV to 21 GeV and reconstructed hadron energy from 15 GeV to 25 GeV, the number of µ − and µ + events get increased by 300 and 150 respectively for 500 kt·yr exposure of the ICAL detector. Apart from showing the total µ − event rates in Table 1, we also present the estimates of 4 Preliminary calculation of the expected fluxes at the INO site has been performed in Ref. [68,69]. The visible differences between the neutrino fluxes at the Kamioka and INO sites appear at lower energies. The main reason behind this is that the horizontal components of the geo-magnetic field are different at the Kamioka (30µT) and INO (40 µT) locations. We plan to use these new fluxes estimated for the INO site (see Ref. [69]) in future studies. 5 For a detailed description of the two binning schemes that we consider in our analysis, see Sec. 5.1.

low-energy (LE)
high-energy (HE)  individual events coming from ν µ → ν µ (P µµ ) disappearance channel and ν e → ν µ (P eµ ) appearance channel. Also, for µ + events, we separately show the contributions originating fromν µ →ν µ (Pμμ) disappearance andν e →ν µ (Pēμ) appearance channels. Here, we see that only ∼ 2% of the total µ − events at the ICAL detector come via the appearance channel. Note that the differences in the total number of µ − and µ + events between the SM case (ε µτ = 0) and non-zero ε µτ of ±0.05 are not significant. But, later while presenting our final results, we see that the ICAL detector can place competitive constraints on ε µτ by exploiting the useful information contained in the spectral distributions of µ − and µ + events as a function of reconstructed observables E µ , cos θ µ , and E had . To establish this claim, now, we show how the expected µ − and µ + event spectra get modified in the presence of non-zero ε µτ in terms of reconstructed E µ and cos θ µ while integrating over entire range of E had .

Event Spectra
In Fig. 3, we show the distributions of µ − (upper panels) and µ + (lower panels) events as a function of reconstructed cos θ µ for three different ranges of E µ . The ranges of E µ that we consider in left, middle, and right panels are [3,4] GeV, [5,11] GeV, and [11,21] GeV respectively. Here, we integrate over E had in its entire range of 0 to 25 GeV. In each panel, we compare the event spectra for three different cases: i) ε µτ = 0.05 (blue line), ii) ε µτ = 0 (the SM case, black line), and iii) ε µτ = −0.05 (red line). Before we discuss the impact of non-zero ε µτ , we would like to mention few general features that µ θ cos    [11,21] GeV ∈ µ E µ θ cos  Figure 3: The distributions of µ − (upper panels) and µ + (lower panels) events for three different Eµ ranges: [3,4] GeV in left panel, [5,11] GeV in middle panel, and [11,21]  are emerging from various panels in Fig. 3. For both µ − (upper panels) and µ + (lower panels), the number of events get reduced as we go to higher energies. It happens because of ∼ E −2.7 ν dependence of the atmospheric neutrino flux. Though the neutrino fluxes are higher along the horizontal direction (cos θ µ around 0) as compared to the other directions (for detailed discussion, see Ref. [69]), but, due to the poor reconstruction efficiency of the ICAL detector along the horizontal direction [71], we see a suppression in µ − and µ + event rates around cos θ µ ∈ [−0.1, 0] irrespective of the choices of E µ ranges. Important to note that as we proceed towards higher E µ , the relative differences in µ − and µ + event rates between the SM case (ε µτ = 0) and non-zero ε µτ (±0.05) get increased for a wide range of cos θ µ . We see similar features in Fig. 2 in Sec. 3, where we show the differences in ν µ → ν µ oscillograms due to the SM case (ε µτ = 0) and non-zero ε µτ (±0.05). We show the improvement in the sensitivity to constrain ε µτ due to high energy events in Sec. 6.1. Next, we discuss the numerical technique and analysis procedure which we adopt to obtain the final results.

Simulation Method
In the present study, we produce all the results with low-energy (LE) and high-energy (HE) binning schemes. Table 2 shows the detailed information about the LE binning scheme for the three reconstructed observables E µ , cos θ µ , and E had . Table 3 portrays the same for the HE binning scheme. In case of LE binning scheme, the range of E µ is [1,11] GeV with total 10 bins each having a width of 1 GeV. In case of HE binning scheme, we extend the range of E µ up to 21 GeV by adding two additional bins in the range of 11 to 21 GeV, where each bin has a width of 5 GeV. As far as reconstructed E had is concerned, in case of LE (HE) binning scheme, the considered range is 0 to 15 GeV (0 to 25 GeV). We can see from Table 2 and Table 3 that the first three bins for E had are same for both the binning schemes, whereas the last bin extend from 4 to 15 GeV (4 to 25 GeV) for LE (HE) binning scheme. For both these binning schemes, we consider the entire range of cos θ µ from -1 to 1. For upward going events, that is cos θ µ ∈ [-1, 0], we consider 10 uniform bins each having width of 0.1. For downward going events, that is cos θ µ ∈ [0, 1], we consider 5 uniform bins each having width of 0.2. Important to note that the downward going events do not have enough path length to oscillate, but, these events play an important role to increase the overall statistics and to minimize the effect of normalization uncertainties in atmospheric neutrino fluxes. Here, we would like to mention that we have not optimized these binning schemes to obtain the best sensitivities, but we ensure that there are sufficient statistics in most of the bins.

Numerical Analysis
In our analysis, the χ 2 function gives us the median sensitivity of the experiment in the frequentist approach [73]. We use the following Poissonian χ 2 − for µ − events in our statistical analysis considering E µ , cos θ µ , and E had as observables (the so-called "3D" analysis  as considered in [74]): In the above equations, N data ijk and N theory ijk denote the observed and expected number of µ − events in a given [E µ , cos θ µ , E had ] bin. In case of LE (HE) binning scheme, N Eµ = 10 (12), N cos θµ = 15, and N E had = 4. In Eq. 5.2, N 0 ijk represents the number of expected events without systematic uncertainties. Following Ref. [75], we consider five systematic errors in our analysis: 20% flux normalization error, 10% error in cross-section, 5% tilt error, 5% zenith angle error, and 5% overall systematics. We incorporate these systematic uncertainties in our simulation using the well known "pull" method [21,76,77]. In Eq. 5.1 and Eq. 5.2, the quantities ζ l denote the "pulls" due to the systematic uncertainties.
When we produce results with only E µ and cos θ µ as observables and do not use the information on hadron energy E had (the so-called "2D" analysis as considered in Ref. [75]), the Poissonian χ 2 − for µ − events takes the form For both the "2D" and "3D" analyses, the χ 2 + for µ + events is determined following the same technique described above. We add the individual contributions from µ − and µ + events to estimate the total χ 2 for both the "2D" and "3D" schemes: In our analysis, we simulate the prospective data considering the following benchmark values of the oscillation parameters: sin 2 θ 23 = 0.5, sin 2 2θ 13 = 0.1, sin 2 θ 12 = 0.3, ∆m 2 21 = 7.5 × 10 −5 eV 2 , and |∆m 2 eff | = 2.4 × 10 −3 eV 2 . To estimate the value of ∆m 2 31 from ∆m 2 eff , we use the Eq. 3.2, where ∆m 2 eff has the same magnitude for NH and IH with +ve and -ve signs respectively. In the fit, we first minimize χ 2 ICAL (see Eq. 5.5) with respect to the "pull" variables ζ l , and then marginalize over the oscillation parameters sin 2 θ 23 in the range [0.36, 0.66], |∆m 2 eff | in the range [2.1, 2.6] × 10 −3 eV 2 , and over both the choices of mass hierarchy, NH and IH, while keeping θ 12 , ∆m 2 21 , sin 2 2θ 13 fixed at their benchmark values. We consider δ CP = 0 • throughout our analysis.

Expected Bounds on NSI parameter ε µτ
We quantify the statistical significance of the analysis to constrain the NSI parameter ε µτ in the following fashion Here, χ 2 ICAL (SM) and χ 2 ICAL (SM + ε µτ ) are calculated by fitting the prospective data with zero (the SM case) and non-zero value of NSI parameter ε µτ respectively. In our analysis procedure, statistical fluctuations are suppressed, and therefore, χ 2 ICAL (SM) ≈ 0. Let us first identify the regions in cos θ µ and E µ plane which give significant contributions towards ∆χ 2 ICAL−NSI . In Fig. 4, we show the distribution 6 of ∆χ 2 − from µ − events in the reconstructed [cos θ µ -E µ ] plane using 500 kt·yr exposure of the ICAL detector and assuming NH. In all the panels of  Fig. 5, we show the sensitivity in the plane of reconstructed cos θ µ and E µ for the "2D" analysis, where we do not use any information on hadrons. But, in right panels of these figures, we portray the sensitivity in the plane of reconstructed cos θ µ and E µ for the "3D" case, where the events are further divided into four sub-bins depending on the reconstructed hadron energy for LE (see Table 2) and HE binning schemes (see Table 3). 6 In Fig. 4, we do not consider the constant contributions in χ 2 coming from the term which involves five pull parameters ζ 2 l in Eq. 5.1 and Eq. 5.3. Also, we do not marginalize over the oscillation parameters in the fit to produce these figures. We adopt the same strategy for Fig. 5 as well. Note that we show our final results considering full pull contributions and marginalizing over the oscillation parameters in the fit as mentioned in previous section. µ θ cos    The common features which are emerging from all the panels in Fig. 4 and Fig. 5 are that most of the sensitivity towards the NSI parameter ε µτ stems from higher energies and longer baselines where the matter effect term 2 √ 2G F N e E becomes sizeable. We observe similar trends in Fig 2 where we plot the differences in ν µ → ν µ oscillation probabilities for the cases ε µτ = 0 and ε µτ = ±0.05. The event spectra as shown in Fig. 3 also confirm this fact. Fig. 4 and Fig. 5 clearly demonstrate while going from LE to HE binning scheme that the sensitivity towards the NSI parameter ε µτ get enhanced due to the increment in the range of E µ from 11 GeV to 21 GeV and for extending the fourth E had bin from 15 GeV to 25 GeV. We can also observe from these figures that with the addition of hadron energy information, the area in the E µ -cos θ µ plane which contributes significantly to ∆χ 2 ± increases, consequently enhancing the net ∆χ 2 ± for both LE and HE binning schemes. Here, we would like to mention that the increase in χ 2 ± is not just due to the information contained in E had , but also due to the valuable information coming from the correlation between E had and muon momentum (E µ , cos θ µ ).
In Fig. 6, we show the sensitivity of the ICAL detector to constrain ε µτ using an exposure of 500 kt·yr and assuming NH as the true mass hierarchy. We obtain these results µ θ cos  (GeV (GeV after performing marginalization over θ 23 , ∆m 2 eff , and both the choices of mass hierarchy as discussed in Sec. 5.2. In the left (right) panel, the results are shown for the LE (HE) binning scheme. In each panel, the red solid line shows the sensitivity for the "3D" case where we consider E µ , cos θ µ , and E had as observables. The black dashed line in each panel portrays the sensitivity for the "2D" case considering E µ and cos θ µ as observables. We see considerable improvement in the sensitivity for both the LE and HE binning schemes when we add E had along with E µ and cos θ µ as observables. We see significant gain in the sensitivity when we increase the E µ range from 11 GeV to 21 GeV and extend the fourth E had bin from 15 GeV to 25 GeV. It is evident from both the panels in Fig. 6 that for the [HE, 3D] case, we obtain the best sensitivity towards the NSI parameter ε µτ , whereas the [LE, 2D] mode gives the most conservative limits.
The 3σ (90%) confidence level bounds on the flavor violating NSI parameter ε µτ obtained using 500 kt·yr exposure of the ICAL are listed in Table 4. The results are shown for true NH (3rd column) and true IH (4th column). For the [HE, 3D] case, we expect the best limit of −0.01 < ε µτ < 0.01 at 90% C.L. using 500 kt·yr exposure of the ICAL detector and irrespective of the choices of true mass hierarchy. For the [LE, 2D] mode, we obtain
the most conservative limit of −0.03 < ε µτ < 0.034 at 90% confidence level assuming NH as true choice. So far we have considered sin 2 2θ 13 = 0.1 as our benchmark choice both in data and theory. If we consider the current best fit value of sin 2 2θ 13 = 0.085 [2][3][4], we have checked that our results will remain almost unaltered. For an instance, if we consider ε µτ = 0.02 in the fit, then we obtain ∆χ 2 = 9.49 assuming sin 2 2θ 13 = 0.1 for the [HE,3D] (fit) mode (see the red curve in the right panel of Fig. 6). Under the same condition, if we take sin 2 2θ 13 = 0.085, then the ∆χ 2 changes to 9.38.

Advantage of having Charge Identification Capability
The ICAL detector is expected to have a uniform magnetic field of strength around 1.5 Tesla over the entire detector. It will enable the ICAL detector to identify the µ − and µ + events separately by observing the bending of their tracks in the opposite directions in the presence of the magnetic field. We label this feature of ICAL as the charge identification capability. In Ref. [71], it has been demonstrated that the ICAL detector will have a very good CID efficiency over a wide range of reconstructed E µ and cos θ µ . In this work, we estimate for the first time the gain in the sensitivity that ICAL may have in constraining the NSI parameter ε µτ due to its CID capability. In each panel of Fig. 7, we show the expected sensitivity of ICAL in constraining ε µτ with (red solid line) and without (black dashed line) CID capability using 500 kt·yr exposure and assuming NH. While preparing these plots, we keep the oscillation parameters fixed in the fit and depict the result for the 2D: E µ , cos θ µ (3D: E µ , cos θ µ , E had ) mode in the left (right) panel assuming the HE binning scheme. It is apparent from Fig. 7 that the CID capability of ICAL in distinguishing µ − and µ + events plays an important role to make it sensitive to the NSI parameter ε µτ like the mass hierarchy measurements [66,74]. In the following, we quote the 90% confidence level limits on ε µτ that the ICAL detector can place with and without CID capabilities for [HE, 2D] and [HE, 3D] modes. with CID : −0.01 < ε µτ < 0.011 at 90% C.L. , without CID : −0.018 < ε µτ < 0.025 at 90% C.L. (6.3) The limits on ε µτ mentioned in Eq. 6.2 and Eq. 6.3 clearly demonstrate the improvement that the ICAL detector can have in constraining the NSI parameter ε µτ due its CID capability.
6.3 Limits on ε µτ for various exposures Fig. 8 shows the 3σ limit on ε µτ as a function of run-time 7 for 50 kt ICAL. The left (right) panel is for LE (HE) binning scheme. In each panel, the black and red lines depict the results for 2D (E µ , cos θ µ ) and 3D (E µ , cos θ µ , E had ) modes respectively. In Fig. 8, various sensitivity curves are drawn keeping all the oscillation parameters fixed in the fit and assuming NH. Here, we give the results only for positive values of ε µτ . We have checked that the results look similar if we consider negative values of ε µτ as well. If we take total 250 kt·yr exposure (50 kt ICAL with a run-time of 5 years), the expected bound is |ε µτ | 0.28 at 3σ C.L. assuming NH in [HE, 3D] binning scheme. It suggests that ICAL can place competitive constraints on ε µτ even for less exposure.

Impact of non-zero ε µτ on Mass Hierarchy Determination
This section is devoted to study how the flavor violating NSI parameter ε µτ affects the mass hierarchy measurement which is the prime goal of the ICAL detector. We quantify the performance ICAL to rule out the wrong hierarchy by adopting the following χ 2 expression: Here, we obtain χ 2 ICAL (true MH) and χ 2 ICAL (false MH) by performing the fit to the prospective data assuming true and false mass hierarchy respectively. Since the statistical fluctuations are suppressed in our analysis, χ 2 ICAL (true MH) ≈ 0. First, we estimate the sensitivity of the ICAL detector to determine the neutrino mass hierarchy by adopting the procedure as outlined in Ref. [74] for the standard case, which we denote as "∆χ 2 ICAL−MH (SM)" in the third column of Table 5. Next, to estimate the mass hierarchy sensitivity in the presence of non-zero ε µτ , we adopt the following strategy. We generate the data with a given mass   Table 5: The mass hierarchy sensitivity of the ICAL detector using 500 kt·yr exposure. For the "SM" case (third column), we do not consider εµτ in data and in fit. For the "SM + εµτ " case (fourth column), we introduce εµτ in the fit and marginalize over it in the range of [-0.1, 0.1] along with oscillation parameters θ23 and ∆m 2 eff . Last column shows how much the mass hierarchy sensitivity deteriorates in presence of εµτ as compared to the SM case. We present our results for various choices of binning schemes and observables assuming both true NH and true IH.
hierarchy assuming ε µτ = 0. Then, while fitting the prospective data with the opposite hierarchy, we introduce ε µτ in the fit and marginalize over it in the range of -0.1 to 0.1 along with the oscillation parameters θ 23 and ∆m 2 eff in their allowed ranges as mentioned in Sec. 5. We label this result as "∆χ 2 ICAL−MH (SM + ε µτ )" in the fourth column of Table 5. We show our results for various choices of binning schemes and observables assuming both true NH and true IH. We consider 500 kt·yr exposure of the ICAL detector. We can see from Table 5 that depending on the choice of true mass hierarchy and the analysis mode, the mass hierarchy sensitivity of ICAL gets reduced by 10% to 20% due to the presence of non-zero ε µτ in the fit.
Next, we study the impact of non-zero ε µτ in the precision measurement of atmospheric parameters in the following fashion. We again generate the prospective data considering the true values of sin 2 θ 23 and |∆m 2 32 | as mentioned above. Then, while estimating the allowed regions in sin 2 θ 23 -|∆m 2 32 | (test) plane, we introduce ε µτ in the fit and marginalize over it in the range of [-0.1, 0.1]. We present these results for the "SM + ε µτ " case at 90% C.L. (2 d.o.f.) with the help of dashed lines in Fig. 9 for various analysis modes. We do not see any appreciable change in the contours when we introduce ε µτ in the fit and vary in its ±10% range. It suggests that the precision measurement of atmospheric oscillation parameters at the ICAL detector is quite robust even if we marginalize over ε µτ in the fit. Similar results were obtained by the Super-Kamiokande Collaboration in Ref. [59], where they studied the impact of NSI's in ν µ -ν τ sector using their Phase I and Phase II atmospheric data.

Summary and Concluding Remarks
In this paper, we explore the possibility of lepton flavor violating neutral current nonstandard interactions (NSI's) of atmospheric neutrino and antineutrino while they travel long distances inside the Earth matter before reaching to the ICAL detector. During the propagation of these neutrinos, we allow an extra interaction vertex with ν µ as the incoming particle and ν τ as the outgoing one and vice versa. With such an interaction vertex, the neutral current non-standard interaction of neutrino with matter fermions gives rise to a new matter potential whose relative strength as compared to the standard matter potential (V CC ) is denoted by ε µτ .
We show that the ICAL detector would be able to place tight constraints on the NSI parameter ε µτ considering reconstructed hadron energy and muon momentum as observables. We find that with E µ ∈ [1,11] GeV and with [E µ , cos θ µ ] as observables, the expected limit on ε µτ at 90% C.L. is −0.03 < ε µτ < 0.03. If we increase the muon energy range from 11 to 21 GeV (E µ ∈ [1, 21] GeV) and consider the reconstructed hadron energy (E had ) as an extra observable on top of the four momenta of muon (E µ , cos θ µ ), we find a significant improvement in the limit which is −0.01 < ε µτ < 0.01 at 90% C.L. using 500 kt·yr exposure of the ICAL detector. We observe that the charge identification capability of the ICAL detector plays an important role to obtain these tight constraints on ε µτ as mentioned above.
Assuming 1 to 21 GeV reconstructed muon energy range and considering E µ , cos θ µ , and E had as observables, we find that the mass hierarchy sensitivity at the ICAL detector deteriorates by ∼10% if we introduce the NSI parameter ε µτ in the fit and marginalize over it in the range of -0.1 to 0.1 along with other standard oscillation parameters. On the other hand, the precision measurement of atmospheric oscillation parameters at the ICAL detector is quite robust even if we marginalize over the NSI parameter ε µτ in fit in the range -0.1 to 0.1.