Enhancing sensitivity to non-standard neutrino interactions at INO combining muon and hadron information

In this paper, we explore the impact of flavor violating neutral current non-standard interaction (NSI) parameter εμτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon _{\mu \tau }$$\end{document} in the oscillation of atmospheric neutrinos and antineutrinos separately using the 50 kt magnetized ICAL detector at INO. We find that due to non-zero εμτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon _{\mu \tau }$$\end{document}, νμ→νμ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _\mu \rightarrow \nu _\mu $$\end{document} and ν¯μ→ν¯μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bar{\nu }}_\mu \rightarrow {\bar{\nu }}_\mu $$\end{document} transition probabilities get modified substantially at higher energies and longer baselines, where vacuum oscillation dominates. We demonstrate for the first time that by adding the hadron energy information along with the muon energy and muon direction in each event, the sensitivity of ICAL to the NSI parameter εμτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon _{\mu \tau }$$\end{document} can be enhanced significantly. The most optimistic bound on εμτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon _{\mu \tau }$$\end{document} that we obtain is -0.01<εμτ<0.01\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\,\,0.01< \varepsilon _{\mu \tau } < 0.01 $$\end{document} at 90%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\%$$\end{document} C.L. using 500 kt·\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\cdot $$\end{document}yr exposure and considering Eμ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_\mu $$\end{document}, cosθμ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\cos \theta _\mu $$\end{document}, and Ehad′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E'_{\mathrm{had}}$$\end{document} as observables in their ranges of [1, 21] GeV, [− 1, 1], and [0, 25] GeV, respectively. We discuss for the first time the importance of the charge identification capability of the ICAL detector to have better constraints on εμτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon _{\mu \tau }$$\end{document}. We also study the impact of non-zero εμτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon _{\mu \tau }$$\end{document} on mass hierarchy determination and precision measurement of oscillation parameters.


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 a e-mail: amina@iopb.res.in b e-mail: sabya.s.chatterjee@durham.ac.uk c e-mail: tarak.thakore@ific.uv.es d e-mail: sanjib@iopb.res.in (corresponding author) 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,61]. Therefore, at low energies, a neutral current NSI can be described by a four-fermion dimension-six operator [62], (1) 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 ε 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 between electron neutrino and ambient electron, which takes the form √ 2G F N e (≡ 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, where, 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 Sect. 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 Sect. 3. In Sect. 4, we present the expected total μ − and μ + events and their distributions as a function of reconstructed E μ and cos θ μ for the following three cases: (1) ε μτ = 0 (SM), (2) ε μτ = 0.05, and (3) ε μτ = − 0.05 using 500 kt·yr exposure of the ICAL detector. In Sect. 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 Sect. 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-21 GeV in addition to the events that belong to the E μ in range of 1-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 study the impact of non-maximal θ 23 while deriving limits on ε μτ at ICAL. (d) We show the advantage of having charge identification (CID) capability in the ICAL detector in placing competitive constraint on ε μτ . (e) We present the expected limits on ε μτ considering different exposures of the 50 kt ICAL detector. (f) 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 Sect. 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 [63]. 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. [63] are actually on the NSI parameters ε αβ /3. Therefore, the above constraints at 90% C.L. can be interpreted as Recently, the authors in Ref. [64] considered the possibility of NSI's in μ-τ sector in the one-year high-energy throughgoing 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 The IceCube-DeepCore collaboration also searched for NSI's involving ε μτ [65]. 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 .
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. (8) 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.

ν μ → ν μ transition with non-zero ε μτ
This section is devoted to explore the effect of non-zero ε μτ in the oscillation of atmospheric neutrino and antineutrino propagating long distances through the Earth matter. For this, we numerically estimate the three-flavor oscillation probabilities including NSI parameter ε μτ and using the PREM profile [66] for the Earth matter density. The NSI parameter ε μτ modifies the evolution of neutrino in matter, which in the flavor basis takes the following form, where ε μτ is real in our analysis. 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: (1) ε μτ = − 0.05 (left panel), (2) ε μτ = 0.0 (the SM case, middle panel), and (3) ε μτ = 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 ε μτ , ν μ 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 In Eqs. 12 and 13, positive and negative signs in front of η μτ are associated with the normal and inverted hierarchy respectively. In case of maximal mixing (θ 23 = 45 • ), Eq. 11 boils down to the following simple expression [69] P ν μ →ν μ = cos 2 L m 2 Since for antineutrino, V CC → −V CC , following Eq. 15, we can write and Equations 16 and 17 explain the broad features in Fig. 1 that we mention above.
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 ν θ cos      Here, in all the panels, we assume NH ν θ cos     Fig. 2 The upper left panel shows the difference in ν μ → ν μ oscillation channel between the SM case (ε μτ = 0) and ε μτ = − 0.05. In the top right panel, the difference is due to the SM case and ε μτ = 0.05. The lower panels are forν μ →ν μ oscillation channel. Here, in all the panels, we assume NH with cos θ ν in the range − 1 to − 0.8. This range of cos θ ν corresponds to the baseline in the range ∼ 12,700 to 10,000 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 ) 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. 15 and obtain the following expression Equation18 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 [66] profile is 6.8 g/cm 3 . Therefore, the standard line-averaged Earth matter potential 3 for 11,00 km baseline is V CC ∼ 2.6 × 10 −13 eV. Considering ε μτ = 0.05, we get For E = 18 GeV, L = 11, 500 km, and m 2 31 = 2.36×10 −3 eV 2 (this mass-squared difference is obtained from Eq. 10 using benchmark values of oscillation parameters), 2 According to a simplified version of the PREM profile [66], the inner core has a radius of ∼ 1220 km with an average density of 13 g/cm 3 . For outer core, R min 1220 km and R max 3480 km with an average density of 11.3 g/cm 3 . Note that in our analysis, we consider the detailed version of the PREM. 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 Y e Ne N p +Nn is the relative number density. For electrically, neutral and isoscalar medium, N e = N p = N n , and therefore, Y e = 0.5.
Thus, for E = 18 GeV and L = 11, 500 km, from Eq. 18, we obtain (22) and Equations 22 and 23 confirm the observations regarding P diff μμ and P diff μμ 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 [70]. 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.

Expected events at ICAL with non-zero ε μτ
The Monte Carlo based neutrino event generator NUANCE [71] 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 [74]. 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. [75,76]. 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 4 Preliminary calculation of the expected fluxes at the INO site has been performed in Ref. [72,73]. 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. [73]) in future studies. 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 Sect. 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: (1) ε μτ = 0.05, (2) ε μτ = 0 (the SM case), and (3) ε μτ = − 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 to 21 GeV and reconstructed hadron energy from 15 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 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 . Table 1 Expected number of μ − and μ + events for 500 kt·yr exposure of the ICAL detector considering low-energy (LE) and high-energy (HE) binning schemes. We present the event rates for the following three cases: (1) ε μτ = 0.05, (2) ε μτ = 0 (the SM case), and (3) ε μτ = − 0.05. Apart from showing the total μ − event rates, we also give the estimates of individual event rates coming from ν μ → ν μ (P μμ ) disappearance channel and ν e → ν μ (P eμ ) appearance channel. For μ + events also, we separately show the contributions fromν μ →ν μ (Pμμ) disappearance channel andν e →ν μ (Pēμ) appearance channel. Here, we consider NH and assume the benchmark values of the oscillation parameters as mentioned in Sect. 4

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-25 GeV. In each panel, we compare the event spectra for three different cases: (1) ε μτ = 0.05 (blue line), (2) ε μτ = 0 (the SM case, black line), and (3) ε μτ = − 0.05 (red line). Before we discuss the impact of non-zero ε μτ , we would like to mention few general features that 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. [73]), but, due to the poor reconstruction efficiency of the ICAL detector along the horizontal direction [75], 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 nonzero ε μτ (±0.05) get increased for a wide range of cos θ μ . We see similar features in Fig. 2 in Sect. 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 μ θ cos    [11,21] GeV ∈ μ E μ θ cos  Fig. 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] GeV in right panel. In each panel, we consider three different cases: (1) ε μτ = 0.05 (blue line), (2) ε μτ = 0 (the SM case, black line), and (3) ε μτ = − 0.05 (red line). Here, we sum over E had in its entire range of 0-25 GeV and show the results for 500 kt·yr exposure and assuming NH events in Sect. 6.1. Next, we discuss the numerical technique and analysis procedure which we adopt to obtain the final results.

Binning scheme in E μ , cos θ μ , E had plane
In the present study, we produce all the results with lowenergy (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-15 GeV (0-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 Table 2 The low-energy (LE) binning scheme adopted for the reconstructed observables E μ , cos θ μ , and E had for each muon polarity. The last column shows the total number of bins taken for each observable 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 [77]. 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 [78]): with N theory i jk In the above equations, N data i jk and N theory i jk 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. 25, N 0 i jk represents the number of expected events without systematic uncertainties. Following Ref. [79], we consider five systematic errors in our analysis: 20% flux normalization error, 10% error in crosssection, 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,80,81]. In Eqs. 24 and 25, 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. [79]), the Poissonian χ 2 − for μ − events takes the form with N theory jk In Eq. 26, N data jk and N theory jk indicate the observed and expected number of μ − events in a given [E μ , cos θ μ ] bin. In Eq. 27, N 0 jk stands for the number of expected events without systematic errors. In case of LE (HE) binning scheme, N E μ = 10 (12) and N cos θ μ = 15.
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. 10, 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. 28) 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 χ 2 ICAL−NSI = χ 2 ICAL SM + ε μτ − χ 2 ICAL (SM) . (29) Here, χ 2 ICAL (SM) and χ 2 ICAL SM + ε μτ are calculated by fitting the prospective data with zero (the SM case) and nonzero value of NSI parameter ε μτ respectively. In our analysis procedure, statistical fluctuations are suppressed, and therefore,    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 . We show the distribution of χ 2 + from μ + events in the plane of reconstructed cos θ μ and E μ for these four different cases in Fig. 5 considering ε μτ = 0.05 in the fit. In left panels of Figs. 4 and 5, we show the sensitivity in the plane of reconstructed cos θ μ and E μ for the "2D" 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. 24 and Eq. 26. 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. 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).
The common features which are emerging from all the panels in Figs. 4 and 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. Figs. 4 and 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 to 21 GeV and for extending the fourth E had bin from 15 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, con-μ θ cos  (GeV (GeV and [E μ , cos θ μ , E had ] respectively. In all the panels, we use 500 kt·yr exposure and assume NH in both data and theory sequently 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 after performing marginalization over θ 23 , m 2 eff , and both the choices of mass hierarchy as discussed in Sect. 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 to 21 GeV and extend the fourth E had bin from 15 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] mode (see the red curve in the right panel of Fig. 6).

NH (true) yr
• 500 kt HE Fig. 6 The sensitivity of the ICAL detector to set bounds on the NSI parameter ε μτ using 500 kt·yr exposure and assuming NH. Left (right) panel is with LE (HE) binning scheme. In each panel, the red solid line shows the sensitivity for the "3D" 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.
These results are obtained after performing marginalization over θ 23 , m 2 eff , and both choices of mass hierarchy Table 4 The expected bound on ε μτ for four different choices of binning schemes and observables at 3σ and 90% C.L. obtained using 500 kt·yr exposure of the ICAL detector. We give results for the both choices of true mass hierarchy.
To obtain these constraints, we marginalize over θ 23 , m 2 eff , and both the choices of mass hierarchy in the fit

Constraints on ε μτ for non-maximal θ 23
Global fits of world neutrino data point towards non-maximal θ 23 giving rise to two degenerate solutions: one lies in the lower octant (LO) where sin 2 θ 23 < 0.5, and the other belongs to the higher octant (HO) where sin 2 θ 23 > 0.5 [2][3][4]. Therefore, it would be quite interesting to see how the expected constraints on the NSI parameter ε μτ at ICAL would be affected if θ 23 turns out to be non-maximal in Nature. To study this, we consider the [HE, 3D] binning scheme for which we obtain the best limit on ε μτ in the previous section assuming sin 2 θ 23 (true) = 0.5. In Fig. 7, we exhibit the performance of the ICAL detector to constrain ε μτ by simulating the prospective data with sin 2 θ 23 (true) = 0.4 as a benchmark choice in the LO (see the blue dashed line). We also depict the expected constraints on ε μτ at ICAL assuming sin 2 θ 23 (true) = 0.6 as a benchmark choice in the HO (see the brown dotted line). We compare these results with sin 2 θ 23 (true) = 0.5 case (see the red solid line). To generate all these results, we consider an exposure of 500 kt·yr, assume NH as the true mass hierarchy, and perform marginalization over θ 23 , m 2 eff , and both the choices of mass hierarchy as discussed in Sect. 5.2. We observe from Fig. 7 that the expected limits on ε μτ for non-maximal θ 23 values deteriorate slightly as compared to the maximal mixing choice. Also, we notice that the bounds are similar for the octant-symmetric true values of sin 2 θ 23 in the LO (0.4) and HO (0.6) since the survival probabilities of muon neutrino and antineutrino in the presence of non-zero ε μτ largely depend on sin 2θ 23 instead of sin θ 23 (see Eqs. 11, 12, and 13 ). The ICAL detector is expected to have a uniform magnetic field of strength around 1.5 T 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. [75], 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. 8, 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. 8 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 [70,78]. In the following, we quote the 90% confidence level limits on ε μτ that the ICAL detector can place with and without CID capabilities The limits on ε μτ mentioned in Eqs. 30 and 31 clearly demonstrate the improvement that the ICAL detector can have in constraining the NSI parameter ε μτ due its CID capability.  Figure 9 shows the 3σ limit on ε μτ as a function of runtime 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. 9, 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. [78] 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 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 Sect. 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-20% due to the presence of non-zero ε μτ in the fit.

Precision measurement of atmospheric parameters with non-zero ε μτ
Next, we turn our attention to the precise measurement of atmospheric oscillation parameters sin 2 θ 23 and | m 2 32 | using 500 kt·yr exposure of the ICAL detector. We quantify this performance indicator using the following expression: ICAL−PM sin 2 θ 23 , | m 2 32 | 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.  | plane for 500 kt·yr exposure of the ICAL detector assuming NH. The brown dot represents the true choices of sin 2 θ 23 and | m 2 32 |. The solid lines show the results for the "SM" case, where we do not consider ε μτ in data and in fit. The dashed lines portray the results when we introduce ε μτ in the fit and marginalize over its ±10% range. For other details, see text where χ 2 0 is the minimum value of χ 2 ICAL in the allowed parameter range. Since we suppress the statistical fluctuations, we have χ 2 0 ≈ 0. First, considering sin 2 θ 23 (true) = 0.5 and | m 2 32 | (true) = 2.4 × 10 −3 eV 2 , we estimate the allowed regions in sin 2 θ 23 -| m 2 32 | (test) plane in the absence of ε μτ at 90% C.L. (2 d.o.f.). We show these results for the "SM" case using solid lines in Fig. 10 for various analysis modes. For the [HE, 3D] case, we achieve the best precision for the atmospheric parameters, and for the [LE, 2D] case, we have the most conservative results.
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. 10 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. [63], 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 non-standard 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 exhibit for the first time that by adding the hadron energy information along with the reconstructed muon energy and muon direction in each event, the sensitivity of ICAL to the NSI parameter ε μτ can be enhanced significantly. We find that considering reconstructed E μ (in the range of 1-11 GeV) and cos θ μ (in the range of -1 to 1) as observables, the expected limit on ε μτ at 90% C.L. is − 0.03 < ε μτ < 0.03. If we increase the muon energy range further 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 also demonstrate for the first time that the charge identification capability of the ICAL detector plays an important role to obtain these tight constraints on ε μτ as mentioned above.
Assuming 1-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.
Before we conclude, we would like to emphasize that though the expected limit on ε μτ at ICAL is weaker than IceCube and DeepCore bounds, but, these results are complementary to each other. IceCube extracts information on the NSI parameters using neutrinos whose energies are above 300 GeV. At these high energies (which are relevant for Ice-Cube), oscillation probabilities are completely dominated by NSIs and there is hardly any impact of vacuum oscillations. As far as DeepCore is concerned, they use neutrinos in their analysis whose energies are above 10 GeV or so, and at these energies, NSIs start to dominate the flavor transition. The proposed ICAL detector is very efficient in the neutrino energy range of 1-10 GeV and it is going to provide a complementary information on the NSI parameter ε μτ in this energy range which is not accessible by IceCube and DeepCore. Moreover, at these energies (where ICAL is very effective), both vacuum oscillations and matter driven NSIs affect the flavor transitions. Therefore, it is quite important to probe these NSI parameters at different energies using different detectors, since they provide complementary information on these parameters.
Another important feature of ICAL is that due to the presence of magnetic field, it provides an opportunity to probe the NSI parameters in neutrino (by observing μ − events) and antineutrino (by observing μ + events) modes separately. Therefore, ICAL provides an unique platform to test some fundamental symmetries of Nature such as Charge-Parity-Time (CPT) symmetry, whose tiny violation may give rise to different mass, mixing, and NSI parameters for neutrinos and antineutrinos. This study may not be possible using any other existing or planned water-, ice-, scintillator-, or argon-based detectors.