Probing Lorentz Invariance Violation with Atmospheric Neutrinos at INO-ICAL

The possibility of Lorentz Invariance Violation (LIV) may appear in unified theories, such as string theory, which allow the existence of a new space-time structure at the Planck scale ($M_p \sim 10^{19}$ GeV). This effect can be observed at low energies with a strength of $\sim 1/M_p$ using the perturbative approach. In the minimal Standard Model extension (SME) framework, the neutrino mass-induced flavor oscillation gets modified in the presence of LIV. The Iron Calorimeter (ICAL) detector at the proposed India-based Neutrino Observatory (INO) offers a unique window to probe these LIV parameters by observing atmospheric neutrinos and antineutrinos separately over a wide range of baselines in the multi-GeV energy range. In this paper, for the first time, we study in detail how the CPT-violating LIV parameters $(a_{\mu\tau}, a_{e\mu}, a_{e\tau})$ can alter muon survival probabilities and expected $\mu^-$ and $\mu^+$ event rates at ICAL. Using 500 kt$\cdot$yr exposure of ICAL, we place stringent bounds on these CPT-violating LIV parameters at 95\% C.L., which are slightly better than the present Super-Kamiokande limits. We demonstrate the advantage of incorporating hadron energy information and charge identification capability at ICAL while constraining these LIV parameters. Further, the impact of the marginalization over the oscillation parameters and choice of true values of $\sin^2\theta_{23}$ on LIV constraints is described. We also study the impact of these LIV parameters on mass ordering determination and precision measurement of atmospheric oscillation parameters.


Introduction and motivation
The Lorentz symmetry has been preserved in the fundamental theories of physics, such as the general theory of relativity and the quantum field theory. However, this symmetry could be challenged at the Planck scale physics (M p ∼ 10 19 GeV), where the unification of gravity with the gauge fields of the Standard Model (SM) of particle physics is expected. Quantum loop gravity [1][2][3][4][5] and String theory [6][7][8][9][10][11] attempt such unification by allowing small perturbation of Lorentz symmetry breaking, so-called Lorentz Invariance Violation (LIV).
neutrino [78,[81][82][83][84][85][86][87][88][89][90]. ICAL provides a unique avenue to probe a plethora of BSM scenarios [49,53,[91][92][93][94][95][96][97][98]. In the context of searching the signature of CPT violation at ICAL, a dedicated work has been performed in Ref. [99] by introducing the new CPT-violating parameters in the effective Hamiltonian as suggested by the authors in Ref. [100]. Since ICAL can measure the atmospheric mass-squared differences using atmospheric neutrinos and antineutrinos separately, any possible differences in their measurements would be a smoking gun signature of the CPT-violation, which may also indicate the possible violation of Lorentz Invariance. This possibility has been explored by the INO Collaboration in Refs. [101,102]. In the present work, for the first time, we explore in detail the impact of CPT-violating LIV parameters in the SME framework using the atmospheric neutrinos at the ICAL detector. The sensitivity of ICAL towards the presence of CPT-violating LIV parameters is estimated using the ICAL detector with 500 kt·yr exposure, and stringent constraints are placed. We will also demonstrate the impact of the presence of LIV on the measurement of neutrino mass ordering and oscillation parameters.
In section 2, we describe the theoretical background of LIV in the context of the neutrino oscillations. The experimental attempts to probe LIV are summarized in section 3. The impacts of LIV on neutrino oscillograms are presented in section 4. In section 5, we elaborate the method to simulate neutrino events at the ICAL detector and modifications in the event distributions due to the presence of LIV. The method to perform the statistical analysis is explained in section 6 which is followed by the results in section 7 where we constrain the CPT-violating LIV parameters and show the impact of LIV on the measurement of mass ordering and precision measurement. Finally, we summarize our findings and conclude in section 8. In appendix A, we describe the properties of gauge invariant LIV parameters. The effect of LIV parameters on the appearance channel is discussed in appendix B. The appendix C presents some additional plots identifying effective regions in the plane of energy and direction of reconstructed muons while constraining LIV parameters.

Lorentz and CPT violation in neutrino oscillations
The study of atmospheric neutrinos provides an avenue to probe neutrino oscillations as an interferometer with a length scale of the diameter of Earth, which may also get affected by the presence of LIV. In order to understand the effect of LIV on neutrino oscillation probabilities, we first look at the modified Hamiltonian in the presence of LIV keeping only the terms which are gauge-invariant and renormalizable under the minimal SME framework. In the ultra-relativistic limit, the effective Hamiltonian H eff ij describing the propagation of left-handed neutrino in vacuum can be written in the following fashion: where, i and j are the neutrino flavor indices, whereas p µ and E are the four momenta and energy of neutrino, respectively. In Eq. 2.1, the first two terms correspond to the standard kinematics 1 , whereas the last two terms denote the LIV interactions governed 1 Here, m 2 ij = U · diag(m 2 1 , m 2 2 , m 2 3 ) · U † , where U is the unitary neutrino mixing matrix known as the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [21,103,104]. by the parameters a µ L (CPT-violating LIV parameters) and c µν L (CPT-conserving LIV parameters). While going from neutrino to antineutrino, the CPT-violating LIV parameters change their sign, whereas the sign of the CPT-conserving parameters remain unchanged (see the discussion in detail in appendix A).
So far, we consider the neutrino propagation in vacuum, but in reality, the neutrinos propagate through Earth's matter. The electron neutrinos undergo W -mediated chargedcurrent interactions with the ambient electrons. This coherent and forward-scattering of ν e with matter electrons gives rise to the so-called "Wolfenstein matter potential" [29,31,105]. Note that the Wolfenstein matter potential doesn't get affected by the presence of non-zero LIV parameters. In the present work, we focus on the isotropic 2 component (µ = 0, ν = 0) of the LIV parameters and denote a 0 L ≡ a, c 00 L ≡ c, p 0 ≡ E. here, U is the standard 3×3 neutrino mixing matrix called PMNS matrix, ∆m 2 ab (≡ m 2 a − m 2 b ) are the mass-squared splittings of the neutrino mass eigenstates, and √ 2G F N e is matter potential where G F is Fermi weak coupling constant, and N e is the electron number density inside the Earth's matter. The contributions to H eff from LIV are given by the terms containing CPT-violating parameters a αβ and CPT-conserving parameters 3 c αβ . Note that the effective Hamiltonian H eff in Eq. 2.3 is written for neutrino, if we go to antineutrino then, U → U * , √ 2G F N e → − √ 2G F N e , a αβ → −a * αβ , and c αβ → c * αβ , where the negative sign for CPT-violating LIV parameter a αβ appears as a multiplicative factor due to the construction of LIV formalism as discussed in appendix A. However, the sign of CPT-conserving LIV parameter c αβ does not change when we go to antineutrino (see appendix A).
In this work, we focus our attention on CPT-violating LIV parameters (a αβ ), whereas the CPT-conserving LIV parameters (c αβ ) will be studied separately in another work. As far as the CPT-violating LIV parameters are concerned, in the present work, we only focus on the off-diagonal parameters: a µτ ≡ |a µτ |e iφµτ , a eµ ≡ |a eµ |e iφeµ , and a eτ ≡ |a eτ |e iφeτ . Our analysis is mainly dominated by ν µ → ν µ andν µ →ν µ disappearance oscillation 2 One of the choices for the isotropic frame is the Sun-centered celestial-equatorial frame [106]. In this frame, we ignore the anisotropy generated due to the Earth's boost, which is suppressed by a factor of ∼ 10 −4 . 3 In Eq. 2.3, there is a multiplicative factor 4/3 in the terms containing c αβ . This factor arises due to the choice of a frame which is rotational invariant. See Eq. A.25 and the related discussion in appendix A. channels where the off-diagonal LIV parameter a µτ appears only in the form of |a µτ | cos φ µτ at the leading order (one can see a similar discussion in the context of NSI in Ref. [107]). Therefore, a complex phase only modifies the effective value of a µτ to a real number between −|a µτ | to |a µτ | at the leading order. We exploit this observation by considering only real values of a µτ in the range of −|a µτ | to |a µτ | for the present work. From these arguments, we can say that our method effectively covers the entire range of complex values of a µτ which is equivalent to the variation of φ µτ in the full range of −π to π. As far as the off-diagonal LIV parameters a eµ and a eτ are concerned, they appear at the subleading order in ν µ → ν µ survival channel and are always θ 13 -suppressed. The imaginary parts associated with the phases are further suppressed in the survival channel compared to the real parts and have a negligible impact ( 1%). Thus, the real values of a eµ and a eτ cover the entire complex parameter space. Note that the CPT-violating LIV parameters (a αβ ) appear in an analogous fashion to that of neutral-current (NC) NSI parameters (ε αβ ) in the effective Hamiltonian (see Eq. 2.3). One can obtain the following relation between the NC-NSI and CPT-violating LIV parameters [23,27] Here, we would like to mention that though these two new physics scenarios show a correspondence as given in the above equation, their physics origins are completely different. Further, the NC-NSI which appears during neutrino propagation is driven by the nonstandard matter effects and it does not exist in vacuum. On the other hand, the LIV is an intrinsic effect that can be present even in vacuum. Note that for long-baseline experiments such as DUNE [108], where a line-averaged constant Earth matter density may be a valid approximation, this one-to-one correspondence between CPT-violating LIV and NC-NSI parameters as shown in Eq. 2.4 may hold. But in the present paper, our focus is on atmospheric neutrino experiments where we deal with a wide range of baselines spanning from 15 km to 12757 km. For each of these baselines, the scaling between CPT-violating LIV and NC-NSI parameters is expected to be different. Moreover, for a large fraction of these baselines (5721 km to 12757 km), neutrinos pass through the inner mantle or even core, where the density changes significantly following the PREM [109] profile with a sharp jump in the density at the core-mantle boundary. Therefore, for these trajectories, the line-averaged constant Earth matter density approximation is no longer valid and a simple scaling between CPT-violating LIV and NC-NSI parameters does not work.
While showing the approximate analytical expressions of neutrino oscillation probabilities in the presence of non-zero CPT-violating LIV parameters a αβ , we use the existing expressions in the literature for NC-NSI [53,107] by replacing the NC-NSI parameters ε αβ with the CPT-violating LIV parameters a αβ following Eq. 2.4. Note that to study the impact of CPT-violating LIV parameters on the atmospheric neutrinos at INO-ICAL, we indeed perform a full-fledged numerical simulation from scratch considering three-flavor neutrino oscillation probabilities in the presence of the PREM profile of Earth and obtain all the results that we show in the present paper.
In this work, we focus on the survival channel (ν µ → ν µ ) and appearance channel (ν e → ν µ ) because these two channels have significant contribution to the neutrino events at the ICAL detector. The survival channel probability P (ν µ → ν µ ) is dominantly affected by LIV parameter a µτ whereas the effects of a eµ and a eτ are subdominant. Using approximations of one mass scale dominance (OMSD) [∆m 2 21 L/(4E) ∆m 2 32 L/(4E)] in the limit of θ 13 → 0 with constant matter density, the survival probability for ν µ can be expressed as [53,110], where and where, we have replaced ε µτ with a µτ using Eq. 2.4. For a given value of CPT-violating parameter a αβ , the parameter ω = +1 for neutrino and ω = −1 for antineutrino due to the construction of LIV formalism as explained in appendix A. Here, β ≡ sgn(∆m 2 32 ) which is the sign of atmospheric mass-squared difference. We have β = +1 for normal ordering (NO, m 1 < m 2 < m 3 ) and β = −1 for inverted ordering (IO, m 3 < m 1 < m 2 ). In the limit of maximal mixing (θ 23 = 45 • ), the Eq. 2.5 boils down to [111], As far as the effects of LIV parameters a eµ and a eτ on survival channel (ν µ → ν µ ) are concerned, it appears in subleading order terms and non-trivial to express analytically. The effects of a eµ and a eτ with similar strengths become dominant in appearance channel (ν e → ν µ ) because it occurs at leading order terms. We provide the expression for P (ν e → ν µ ) in the presence of both LIV parameters a eµ and a eτ at-a-time in appendix B. The effect of LIV parameter a µτ on P (ν e → ν µ ) is not significant, hence, it is not present in Eq. B.4. Note that these approximate expressions for oscillation probabilities are just for understanding purposes. However, in this work, we use numerically calculated full threeflavor oscillation probabilities with LIV in the presence of Earth's matter with the PREM profile [109].

A brief history of the search for LIV in neutrino oscillations
In this section, we discuss how the search for LIV parameters in neutrino oscillations progressed. In 2004, a general formalism for LIV and CPT-violation was developed in the neutrino oscillation sector, where possible definitive signals are predicted in the light of prevailing and proposed neutrino experiments [112]. After that, a few sample and global models are proposed to illustrate the key physical effects of LIV by considering both with and without neutrino mass, along that, a few generalized models with operators of arbitrary dimension are discussed to accommodate the signals observed in various neutrino experiments [113][114][115][116][117][118][119][120][121][122]. Now, we summarize the experimental attempts to probe the LIV: • A proposal of possible LIV signal was made in the context of the results obtained from the short-baseline experiments, and in 2005, the first experimental attempt to probe Lorentz and CPT violation was made by the liquid scintillator neutrino detector (LSND) experiment [123] to accommodate the oscillation excess. But null results of LIV in LSND signal put bounds on (a L ) and (E × c L ) of the order of 10 −19 GeV, which is in an expected scale of suppression.
• In 2008, a search for the sidereal modulation in the MINOS near detector neutrino data was performed [124]. No significant evidence was found, resulting in the bounds on the sidereal components of LIV parameters to the order of 10 −20 GeV. Consequently, in 2010, the same search was performed with the MINOS far detector neutrino data (Run I, II, and III) [125]. No signature of sidereal effect was found, and the bounds are further improved to the order of 10 −23 GeV.
• In 2010, the IceCube experiment [126] had reported the results of the search for a Lorentz-violating sidereal signal with the atmospheric neutrinos. No directiondependent variation was found, and the constraints improved by a factor of 3 for (a L ) and an order of 3 for (c L ).
• In 2011, the result from the MiniBooNE Collaboration [127] was consistent with MINOS near detector bounds.
• In 2012, the MINOS Collaboration searched the sidereal variation with muon type antineutrino data at the near detector [128], and the result was consistent with their earlier work with neutrino data [124].
• In the later part of 2012, the Double Chooz Collaboration [129] showed their results of the search for LIV with reactor antineutrinos. Data indicated no sidereal variations, so their work set the first limits on 14 LIV coefficients associated with e-τ sector and set two competitive limits associated with e-µ sectors.
• The presence of LIV allows the neutrino and antineutrino mixings, hence to see such a signal, in early 2013, MINOS has performed a search analysis [130]. There was no evidence of such a signal, resulting in the limit of the appropriate 66 LIV coefficients (H and g) 4 , which govern the neutrino-antineutrino mixing.
• Further search for neutrino-antineutrino mixing was performed using the disappearance of reactor antineutrinos in the Double Chooz experiment [131] and set limits on another 15 LIV coefficients (H and g). A similar study in the context of solar neutrinos was performed by the authors in Ref. [132].
• In 2015, the Super-Kamiokande (Super-K) Collaboration published their result to search Lorentz Invariance with atmospheric neutrinos with a large range of baselines and wide range of energies [22]. No evidence of Lorentz violation was observed, so limits were set on the renormalizable isotropic LIV coefficients in the e-µ, µ-τ , and e-τ sectors, with an improvement of seven orders of magnitude.
• In 2017, the T2K Collaboration searched LIV using sidereal time dependence of neutrino flavor transitions with the T2K on-axis near detector [133]. The results were consistent with the limits put by other short-baseline experiments.
• In 2018, a simultaneous measurement of LIV coefficients was done using the Daya Bay reactor antineutrino experiment [134]. The bounds on the appropriate parameters were agreed to the suppression of the Planck mass scale (M p ).
• The IceCube experiment analyzed the LIV hypothesis using the high-energy part of atmospheric neutrinos in the year 2018 [25]. They perform the analysis using a perturbative approach in an effective two-flavor neutrino oscillation framework and obtained the limits on individual LIV coefficients (a L and c L ) with mass dimensions ranging from 3 to 8. These are the most stringent limits on Lorentz violation (in the µ-τ sector) set by any physical experiment.
The current constraints on all the relevant LIV/CPT-violating parameters are nicely summarized in Ref. [106]. In a recent review [135], the authors have discussed the phenomenological effects of CPT and Lorentz symmetry violations in the context of particle and astroparticle physics.

Effects of Lorentz Invariance Violation on oscillograms
In this section, we describe how does the neutrino oscillation probabilities get affected due to the presence of LIV. The atmospheric neutrinos are produced as a result of the interaction of primary cosmic rays with the nuclei of the atmosphere and mostly consist of electron and muon flavors. The atmospheric neutrinos with the multi-GeV energy range have access to a wide range of baselines starting from about 15 km to 12757 km. The upward-going neutrinos with longer baselines experience the Earth's matter effect, which results in the modification of oscillation probabilities. The effect of LIV on neutrino oscillations can be explored at several L/E values available for atmospheric neutrinos.
In the present analysis, we numerically calculate the neutrino oscillation probability in the three-flavor paradigm with Earth's matter effect using the PREM profile [109]. We use the benchmark values of oscillation parameters given in Table 1. The value of ∆m 2 31 is obtained from the effective mass-squared difference 5 ∆m 2 eff . To consider NO as mass 5 The effective mass-squared difference is defined in terms of ∆m 2 31 as follows [140,141]: ∆m 2 eff = ∆m 2 31 − ∆m 2 21 (cos 2 θ12 − cos δCP sin θ13 sin 2θ12 tan θ23).   ordering, we use positive value of ∆m 2 eff whereas for IO, ∆m 2 eff is taken as negative with the same magnitude.
The neutrino events at ICAL are contributed by both disappearance (ν µ → ν µ ) channel as well as appearance (ν e → ν µ ) channel. Since, more than 98% events at ICAL are contributed by disappearance (ν µ → ν µ ) channel, we discuss the impact of LIV on P(ν µ → ν µ ) in this section. In Fig. 1, we show the impact of non-zero value of CPT-violating LIV parameter a µτ on ν µ → ν µ survival probability oscillograms in the plane of energy (E ν ) and direction 6 (cos θ ν ) of neutrino. We consider three different choices of a µτ which are −1.0 × 10 −23 GeV, 0, and 1.0 × 10 −23 GeV as shown in the left, middle and right columns, respectively, in Fig. 1. The panels in top and bottom rows correspond to neutrinos and antineutrinos, respectively.
In Fig. 1, the effect of LIV parameter a µτ can be observed at energies above 10 GeV and baselines with cos θ ν < −0.6. This region is dominated by vacuum oscillations. The dark blue diagonal band, which is termed as oscillation valley in Ref. [53,89], bends in the presence of LIV parameter a µτ . The direction of bending depends on the sign of a µτ . For example, the oscillation valley for neutrino bends in the downward direction for negative a µτ in the top left panel, whereas for positive a µτ , the bending is in the upward direction in the top right panel. For a given value of a µτ , the oscillation valley curves in the opposite directions for neutrinos and antineutrinos. In summary, we can infer that the bending of oscillation valley for neutrino with positive a µτ (top right panel) is similar to that for antineutrino with negative a µτ (bottom left panel) and vice versa.
We now discuss how the modification of the oscillation valley due to non-zero a µτ can be understood analytically. The oscillation valley corresponds to the first oscillation minimum in ν µ survival probability in Eq. 2.9 which results in the following relation between E ν and cos θ ν [53], where, we assume L ν ≈ 2R| cos θ ν |. For the case of standard interactions (SI) with a µτ = 0, Eq. 4.3 results in a linear relation between E ν and cos θ ν which is manifested as oscillation valley with straight line for SI case in the middle panels of Fig. 1.
In the case of NO (β = +1), if a µτ is positive then for neutrino (ω = +1), E ν increases compared to that for SI case for a given cos θ ν as observed in the top right panel of Fig. 1. If we consider negative value of a µτ then for neutrino, E ν decreases as seen in the top left panel of Fig. 1. As far as antineutrino is concerned, where ω = −1, positive (negative) value of a µτ results in decrease (increase) in E ν as shown in the bottom right (left) panel of Fig. 1. Since, β and ω appear as product, the effect of positive (negative) a µτ with neutrino is same as that of negative (positive) a µτ with antineutrino. If we consider the case of IO (β = −1), then all these effects are reversed with the opposite curvatures of oscillation valley.
In Fig. 2, we present the impact of CPT-violating LIV parameter a eµ on ν µ → ν µ survival probability oscillograms in the plane of neutrino energy (E ν ) and direction (cos θ ν ). We assume three different values of a eµ which are −2.5 × 10 −23 GeV, 0, and 2.5 × 10 −23 GeV as shown in the left, middle and right columns, respectively, in Fig. 2. In the top and bottom panels, we portray plots for neutrinos and antineutrinos, respectively. Unlike 6 The relation between neutrino zenith angle θν and neutrino baseline Lν is given as where, R, h, and d represent the radius of Earth, the average production height for neutrinos, and the depth of the detector, respectively. In the present study, we consider R = 6371 km, h = 15 km, and d = 0 km.

Figure 2:
The νµ survival probability oscillograms in (Eν , cos θν ) plane for neutrino (top panels) and antineutrino (bottom panels) with Earth's matter effect considering the PREM profile. We use NO and benchmark oscillation parameters given in Table 1. We consider the values of aeµ = −2.5 × 10 −23 GeV, 0, and 2.5 × 10 −23 GeV in the left, middle and right panels, respectively.
the case of a µτ , here, we do not observe any bending in the oscillation valley due to the presence of non-zero a eµ . However, a distortion in the oscillation valley can be observed due to a eµ . The regions with matter effect also get modified. For neutrino, the effect due to positive a eµ is significantly more than that due to negative a eµ . We also observe that the effect of a eµ for antineutrino is not significant. We would like to highlight that the effect of a eµ appears at the subleading order in the expression of P (ν µ → ν µ ), which are non-trivial to express analytically. We can analyze the same by studying the analytical expression for the appearance channel P (ν e → ν µ ) as given in Eq. B.4 in appendix B, where the effect of a eµ appears in the leading order terms. The effect of a eµ is dominantly contributed by the fifth term in Eq. B.4, which has the following form where (∆m 2 31 − a CC ) factor in the denominator causes the matter-driven resonance effect for the case of neutrino with NO. Since this term has positive sign for the case of neutrino (ω = +1), the positive value of a eµ increases P (ν e → ν µ ). This can be translated as the decrease in P (ν µ → ν µ ) around the region of matter effect as shown in the top right panel antineutrino (bottom panels) with Earth's matter effect considering the PREM profile. We use NO and benchmark oscillation parameters given in Table 1. We consider the values of aeτ = −2.5 × 10 −23 GeV, 0, and 2.5 × 10 −23 GeV in the left, middle and right panels, respectively.
of Fig. 2 because we have, where P (ν µ → ν e ) = P (ν e → ν µ ) for δ CP = 0. Note that ν µ → ν τ oscillation channel also gets affected due to matter effects in certain ranges of energies and baselines (see Ref. [142]). But as far as the impact of LIV parameter a eµ is concerned, it appears only at the subleading order in ν µ → ν τ oscillation channel. For negative value of a eµ , P (ν e → ν µ ) decreases and P (ν µ → ν µ ) increases, hence we observe dilution of matter effect in the top left panel of Fig. 2. We can see that the effect of a eµ is larger in the case of neutrino than antineutrino. This happens because a CC becomes negative for antineutrino and matterdriven resonance condition is not fulfilled for NO due to which the above-mentioned term does not contribute significantly. Similar to LIV parameter a eµ , now we depict the impact of a eτ on ν µ → ν µ survival probability oscillograms in the plane of neutrino energy (E ν ) and direction (cos θ ν ) in Fig. 3. We take three different values of a eτ which are −2.5×10 −23 GeV, 0, and 2.5×10 −23 GeV as shown in the left, middle and right columns, respectively, in Fig. 3. The features present in the oscillograms due to a eτ in Fig. 3 are analogous to that due to a eµ in Fig. 2. Here also, we can observe the distortions in the oscillation valley as well as the matter effect regions. The effect of a eτ is more for neutrinos than that for antineutrino. If we focus on neutrinos, we see that the effect of positive a eτ (top right panel) is larger than that of negative a eτ (top left panel). Following the similar arguments as for a eµ in the previous paragraph, these features can be explained using the seventh term in Eq. B.4. These asymmetric effects for CPT-violating LIV parameter a eµ and a eτ can also be observed in our result section where we constrain them using 500 kt·yr exposure at the ICAL detector.
Till now, we have described the impact of CPT-violating LIV parameters at the probability level, where we observe some interesting features. In the next section, we explain the procedure to simulate neutrino events at the ICAL detector and explore how the event distribution modifies in the presence of LIV.

Event generation at ICAL
The proposed ICAL detector at INO [78] consists of a 50 kt magnetized iron stack of size 48 m × 16 m × 14.5 m which would be able to detect atmospheric neutrinos and antineutrinos separately in the multi-GeV range of energies over a wide range of baselines. The vertically stacked 151 layers of iron of thickness 5.6 cm act as the target material for neutrino interactions. The charged-current interactions of neutrinos and antineutrinos with iron nuclei produce charged muons, µ − and µ + respectively. The 2 m × 2 m Resistive Plate Chambers (RPCs) [143][144][145] sandwiched between iron layers act as active detector elements. The charged particles deposit energies while passing through RPCs, which result in the production of electronic signals. The perpendicularly arranged pickup strips of RPC give the X and Y coordinates of the hit, whereas the Z coordinate is obtained from the layer number of RPC.
The muons are minimum ionizing particles in the multi-GeV range of energies, and hence, they pass through many layers producing hits in each of those layers in the form of tracks. The ICAL can efficiently reconstruct muon energies and directions in the reconstructed muon energy range of 1 to 25 GeV [80]. Thanks to the magnetic field of 1.5 Tesla [79], ICAL will be able to distinguish µ − and µ + from the opposite curvatures of their tracks. This charge identification (CID) capability of ICAL helps us to separately identify neutrinos (ν µ ) and antineutrinos (ν µ ). Due to ns time resolution of RPCs [146][147][148], ICAL can efficiently distinguish the upward-going and downward-going muon events. The resonance scattering and deep-inelastic scattering of neutrinos with iron nuclei result in the production of hadrons which produce showers inside the detector.
In the present work, the neutrino events are simulated using NUANCE Monte Carlo (MC) neutrino event generator [149] using ICAL-geometry as a target. While simulating data, we use atmospheric neutrino flux at the INO site [150,151] in the Theni district of Tamil Nadu, India, with a rock coverage of 1 km in all directions. Due to the rock coverage of about 1 km (3800 m water equivalent), the downward-going cosmic muon background gets reduced by ∼ 10 6 [152]. The application of fiducial volume cut further vetoes these events. We incorporate the effect of solar modulation on atmospheric neutrinos flux by considering solar maxima for half exposure and solar minima for another half. Since we estimate the median sensitivities in our analysis, we simulate unoscillated neutrino events at ICAL for a large exposure of 1000-year MC to minimize the statistical fluctuation. The neutrino oscillation with LIV in the three-flavor framework in the presence of matter with the PREM profile [109] is incorporated using reweighting algorithm following Refs. [82,83,85].
The detector response to muons and hadrons is simulated using GEANT4 simulations of the ICAL detector as described in Refs [80,153]. The track-like events associated with muons are fitted with the Kalman filter technique to obtain information about the energy, direction, and charge of the reconstructed muons [154]. Although, the configuration of ICAL is mainly optimized to measure the four-momenta of muon, it can also measure the hadron energy in each event. In the resonance and deep-inelastic scattering processes, a significant fraction of neutrino energy is carried away by the hadrons and it is defined as E had = E ν − E µ . The hadron energy deposited in a shower is estimated using the total number of hits that are not part of the muon track. Since we can obtain information on the hadron energy apart from precisely measuring the four-momenta of muon on an event-by-event basis, we consider the binning in three separate observables: muon direction cos θ µ , muon energy E µ , and hadron energy E had . The use of both E µ and E had as separate observables, allows us to indirectly probe the incoming neutrino energy. Note that if we add these two observables to reconstruct the neutrino energy, we may not be able to take the advantage of precise measurement of muon energy. This additional information on hadron energy is very important to improve the constraints on the LIV parameters, as we show in our result section. The reconstructed µ − and µ + events at ICAL detector with observables cos θ rec µ , E rec µ , and E rec had are obtained after folding with detector properties using the migration matrices [80,153] provided by ICAL collaboration following the procedure mentioned in Refs [82,83,85]. These reconstructed µ − and µ + events at 50 kt ICAL detector are scaled from 1000-yr MC to 10-yr MC for analysis. Now, we present the reconstructed µ − and µ + events expected for 500 kt·yr exposure at ICAL for the case of SI and SI with LIV.

Total event rates
First of all, we would like to address the question: can we observe the signal of nonzero CPT-violating LIV parameter (a µτ , a eµ , and a eτ ) in the total number of µ − and µ + events at 50 kt ICAL detector for 10-year exposure? To answer this question, we estimate total reconstructed µ − and µ + events at ICAL for 500 kt·yr exposure for the cases of SI and SI with non-zero CPT-violating LIV parameters one at-a-time considering a µτ = ± 1.0 × 10 −23 GeV, a eµ = ± 2.5 × 10 −23 GeV, and a eτ = ± 2.5 × 10 −23 GeV. We present these numbers in Table 2 assuming NO as true mass ordering while using the values of benchmark oscillation parameters given in Table 1. While calculating total number of events, we integrate over reconstructed muon energy E rec µ in the range 1 to 25 GeV, reconstructed muon direction cos θ rec µ in the range -1 to 1, and reconstructed hadron energy E rec had in the range 0 to 25 GeV. In Table 2, we present total µ − events along with the estimates of individual events contributed from ν µ → ν µ disappearance channel and ν e → ν µ appearance channel. Similarly, we also show total µ + events originating fromν µ →ν µ disappearance andν e →ν µ appearance channels. Here, we observe that about 2% of total µ − events at ICAL for SI  are contributed by the appearance channel. The tau lepton may get produced during the interaction of ν τ inside the detector via ν µ → ν τ and ν e → ν τ channels. The number of muon events produced during tau decay inside the detector is small (about 2% of the total upward going muon events produced from ν µ interactions [155]), and the energies of these muon events are lower and mostly below the 1 GeV energy threshold of the ICAL detector.
In the present work, we do not consider this small contribution. We would like to point out that the difference between total µ − (µ + ) events for the cases of SI and SI with non-zero CPT-violating LIV parameters (a µτ , a eµ , and a eτ ) one ata-time is not large. But while presenting our results, we demonstrate that ICAL can place competitive bounds on LIV parameters using the information contained in the distributions of reconstructed µ − and µ + events as a function of E rec µ , cos θ rec µ , and E rec had . To validate this claim, now we show the impact of non-zero LIV parameters one at-a-time on the spectral distributions of reconstructed µ − and µ + events as a function of E rec µ , and cos θ rec µ while integrating over E rec had .

Reconstructed event distributions
In the present analysis, we use the binning scheme for reconstructed observables E rec µ , cos θ rec µ , and E rec had as shown in Table 3. We have total 13 bins for E rec µ in the range of 1 to 25 GeV, 15 bin for cos θ rec µ in the range -1 to 1, and 4 bins for E rec had in the range 0 to 25 GeV. Although no significant oscillations are developed for downward-going events, these events are included in our analysis because they help in increasing the overall statistics and minimizing the normalization errors in the atmospheric neutrino events. In this way, we also take care of those upward-going (near horizon) neutrino events, which result in the downward-going reconstructed muon events due to angular smearing during the interaction  of neutrinos as well as reconstruction. We have considered the same binning scheme for reconstructed µ − as well as µ + events.
In Fig. 4, we present the event distributions of reconstructed µ − events as a function of cos θ rec µ in the range -1 to 0 for three different ranges of reconstructed muon energies. For E rec µ , we consider the energy range of [3,5] GeV, [5,11] GeV, and [11,25] GeV as shown in the top, middle, and bottom panels, respectively. We integrate over E rec had in the range of 0 to 25 GeV. The error bars show the statistical uncertainties, which are obtained by taking the square root of the number of events. The left, middle and right panels show the impact of a µτ , a eµ , and a eτ , respectively, considering only one parameter at-a-time and remaining parameters are assumed to be zero. In the same fashion, we demonstrate the impact of a µτ , a eµ , and a eτ on reconstructed µ + events in Fig. 5.
Before we discuss the impact of non-zero LIV parameters a µτ , a eµ , and a eτ , we would like to mention few general features of the event spectra in Figs. 4 and 5. The number of reconstructed µ − and µ + events decrease as we go towards higher energies. This happens because of power law where the atmospheric neutrino flux has an energy dependence of E −2.7 ν . Though the neutrino and antineutrino fluxes along horizontal direction (cos θ ν ∼ 0) are higher than the vertical direction (cos θ ν ∼ ±1), we observe lesser reconstructed µ − and µ + events around cos θ rec µ ∈ [−0.1, 0] due to the poor reconstruction efficiency of the ICAL detector for horizontal direction irrespective of the choice of E rec µ range. Note that the number of reconstructed µ + events are smaller than the reconstructed µ − events due to the lower interaction cross-section for antineutrinos.
In the atmospheric sector, the impact of a µτ is largest among all the CPT-violating LIV parameters as shown in the left panels of Fig. 4. In left panels, we compare three cases corresponding to a µτ = −1.0 × 10 −23 GeV (red curves), 0 (black curves), and 1.0 × 10 −23 GeV (blue curves). The impact of a µτ is significant at higher energies of E rec µ ∈ [5, 11] GeV and [11,25] GeV. We also observe that the polarity of difference between SI and SI with non-zero a µτ changes as move from cos θ rec µ of -1 to 0. Therefore, the total rate shows a diluted effect of non-zero a µτ and the good directional resolution of the ICAL detector plays an important role while constraining LIV parameter a µτ . The impact of the LIV Figure 4: The distributions of reconstructed µ − events at the ICAL detector for 500 kt·yr exposure for three different range of E rec µ : [3,5] GeV (top panels), [5,11] GeV (middle panels) and [11,25] GeV (bottom panels) where y-ranges are different among these rows. Here, error bars show the statistical uncertainties.
Note that, although the events are binned in (E rec µ , cos θ rec µ , E rec had ) binning scheme given in Table 3, the reconstructed events in the hadron energy bins are integrated in the range 0 to 25 GeV in these plots. The left, middle and right panels show the effect of aµτ , aeµ, and aeτ , respectively one parameter at-a-time. In the left panels, the red, black, and blue curves represent aµτ = −1.0 × 10 −23 GeV, 0, and 1.0 × 10 −23 GeV, respectively. In the middle panels, the red, black, and blue curves refer to aeµ = −2.5 × 10 −23 GeV, 0, and 2.5 × 10 −23 GeV, respectively. Similarly, in the right panels, the red, black, and blue curves stand for aeτ = −2.5 × 10 −23 GeV, 0, and 2.5 × 10 −23 GeV, respectively. We assume sin 2 θ23 = 0.5 and NO as true mass ordering. The values of other oscillation parameters are taken from Table 1. parameter a eµ is shown in the panels of the middle column where the red, black, and blue curves refers to a eµ = −2.5 × 10 −23 GeV, 0, and 2.5 × 10 −23 GeV, respectively. We can observe that the difference between SI and SI with non-zero a eµ is significant for E rec µ ∈ [5,  [11,25] GeV ∈ + µ rec E Figure 5: The distributions of reconstructed µ + events at the ICAL detector for 500 kt·yr exposure for three different range of E rec µ : [3,5] GeV (top panels), [5,11] GeV (middle panels) and [11,25] GeV (bottom panels) where y-ranges are different among these rows. Here, error bars show the statistical uncertainties.
Note that, although the events are binned in (E rec µ , cos θ rec µ , E rec had ) binning scheme given in Table 3, the reconstructed events in the hadron energy bins are integrated in the range 0 to 25 GeV in these plots. The left, middle and right panels show the effect of aµτ , aeµ, and aeτ , respectively one parameter at-a-time. In the left panels, the red, black, and blue curves represent aµτ = −1.0 × 10 −23 GeV, 0, and 1.0 × 10 −23 GeV, respectively. In the middle panels, the red, black, and blue curves refer to aeµ = −2.5 × 10 −23 GeV, 0, and 2.5 × 10 −23 GeV, respectively. Similarly, in the right panels, the red, black, and blue curves stand for aeτ = −2.5 × 10 −23 GeV, 0, and 2.5 × 10 −23 GeV, respectively. We assume sin 2 θ23 = 0.5 and NO as true mass ordering. The values of other oscillation parameters are taken from Table 1. 11] GeV. Similarly, we show the effect of non-zero a eτ on the distribution of reconstructed µ − events in the right panels of Fig. 4 where the red, black, and blue curves stand for a eτ = −2.5 × 10 −23 GeV, 0, and 2.5 × 10 −23 GeV, respectively. The difference between SI and SI with non-zero a eτ is largest for E rec µ ∈ [5,11] GeV, but we also see significant difference for E rec µ ∈ [3, 5] GeV. As far the case of µ + is concerned, we observe that the impact of a µτ is significant in the higher muon energy range of [5,11] GeV and [11,25] GeV as shown in the left panels in Fig. 5. By comparing the left panels of Figs. 4 and 5, we observe that for a given value of a µτ , the effect is in the opposite direction for reconstructed µ − and µ + events. For example, the reconstructed µ − events increase for the negative value of a µτ and the red curve is higher than black curve in the bottom left panel of Fig. 4 for E rec µ ∈ [11,25] GeV whereas the corresponding curve is lower in the case of µ + as shown in Fig. 5. The addition of reconstructed µ − and µ + events will dilute this feature. This observation indicates that the charge identification efficiency of the ICAL detector will play an important role while constraining LIV parameter a µτ . The effect of other LIV parameters a eµ and a eτ on the event spectra of reconstructed µ + events is significantly lower than that for the case of µ − . This happens because the effect of a eµ and a eτ is driven by matter resonance term, which occurs for the case of neutrino for NO but not antineutrino as explained in the section 4. Next, we discuss the numerical method and analysis procedure that we use to obtain our final results.

Simulation method
To obtain the median sensitivity of the ICAL atmospheric neutrino experiment in the frequentist approach [156], we define the following Poissonian χ 2 − [157] for reconstructed µ − events in E rec µ , cos θ rec µ , and E rec had observables (the so-called "3D" analysis as considered in [85]): In these equations, N data ijk and N theory ijk stand for the expected and observed number of reconstructed µ − events in a given (E rec µ , cos θ rec µ , E rec had ) bin, respectively. We used the binning scheme mentioned in Table 3 where N E rec µ = 13, N cos θ rec µ = 15, and N E rec had = 4. In Eq. 6.2, N 0 ijk corresponds to the number of expected events without considering systematic uncertainties. In this analysis, we consider five systematic uncertainties following Refs. [82,83,158]: 20% error in flux normalization, 10% error in cross section, 5% energy dependent tilt error in flux, 5% zenith angle dependent tilt error in flux, and 5% overall systematics. We use the well-known method of pulls [110,159,160] to incorporate these systematic uncertainties in our simulation. The pulls due to the systematic uncertainties are denoted by ξ l in Eqs. 6.1 and 6.2.
When we do not incorporate the hadron energy information E rec had in our analysis and use only E rec µ , and cos θ rec µ (the so-called "2D" analysis as considered in Ref. [82]), the Poissonian χ 2 − for µ − events boils down to In Eq. 6.3, N data jk and N theory jk correspond to the observed and expected number of reconstructed µ − events in a given (E rec µ , cos θ rec µ ) bin, respectively. The expected number of events without systematic uncertainties are represented by N 0 jk in Eq. 6.4. In case of binning scheme mentioned in Table 3, N E rec µ = 13, and N cos θ rec µ = 15. Following the same procedure, as described above, the χ 2 + for reconstructed µ + events is determined for both the "2D" and "3D" analyses. The total χ 2 is estimated by adding the contribution from µ − and µ + : In this analysis, we simulate the prospective data using the benchmark values of oscillation parameters given in Table 1. We use Eq. 4.1 to estimate the value of ∆m 2 31 from ∆m 2 eff where ∆m 2 eff has the same magnitude for NO and IO with positive and negative signs, respectively. To incorporate the systematic uncertainties, the χ 2 ICAL is minimized with respect to pull parameters ξ l in the fit. After that the marginalization is performed over atmospheric mixing angle sin 2 θ 23 in the range (0.36, 0.66), atmospheric mass square difference |∆m 2 eff | in the range (2.1, 2.6) ×10 −3 eV 2 , and over both choices of mass ordering, NO and IO, in the theory. We do not consider any priors for sin 2 θ 23 and |∆m 2 eff | in our analysis. While performing fit, we keep the solar oscillation parameters sin 2 2θ 12 and ∆m 2 21 fixed at their values as mentioned in Table 1. Since the reactor mixing angle is already very well measured [136][137][138][139], we use a fixed value of sin 2 2θ 13 = 0.0875 both in data and theory. Throughout our analysis, we take δ CP = 0 both in data and theory.

Effective regions in (E rec
µ , cos θ rec µ ) plane to constrain LIV parameters For simulating the prospective data for statistical analysis, we assume the oscillations with standard interactions only where LIV parameters are taken as zero. The statistical significance of the analysis for constraining the LIV parameters (a µτ , a eµ , a eτ ) is quantified in the following fashion: Here, χ 2 ICAL (SI) and χ 2 ICAL (SI + a αβ ) are calculated by fitting the prospective data with only SI case (no LIV) and with SI in the presence of non-zero LIV parameter a αβ , respectively. We consider only one LIV parameter at-a-time where a αβ can be a µτ , a eµ , or a eτ . Since, the statistical fluctuations are suppressed to estimate the median sensitivity, we have χ 2 ICAL (SI) ∼ 0. As we will demonstrate in the later part of the result section that ICAL being an atmospheric neutrino experiment dominated by ν µ → ν µ survival channel places tightest constraint on a µτ among all the LIV parameters. Here, we discuss what regions in (E rec µ , cos θ rec µ ) plane contribute significantly towards the sensitivity of ICAL for a µτ . A similar discussion for other two off-diagonal LIV parameters a eµ and a eτ is provided in the appendix C. The sensitivity of ICAL towards the LIV parameter a µτ stems from both µ − and µ + events for both the choices of mass ordering. While considering a µτ = 1.0 × 10 −23 GeV (−1.0 × 10 −23 GeV) in theory and a µτ = 0 in data, the contribution of µ − and µ + towards fixed-parameter 7 ∆χ 2 ICAL-LIV for 500 kt·yr exposure at ICAL for NO is 43.5 (43.9) and 25.0 (23.8), respectively.
Before we present our final results, let us first identify the effective regions in (E rec µ , cos θ rec µ ) plane which contribute significantly towards ∆χ 2 ICAL-LIV . In Fig. 6, we show the distribution of fixed-parameter ∆χ 2 − ( ∆χ 2 + ) without pull penalty term 8 from reconstructed µ − (µ + ) events in (E rec µ , cos θ rec µ ) plane for 500 kt·yr exposure at ICAL assuming NO at true mass ordering. For demonstration purpose, we have added the ∆χ 2 contribution from all E rec had bins for each (E rec µ , cos θ rec µ ) bin while using binning scheme mentioned in Table 3. In the top (bottom) panels in Fig. 6, we take non-zero LIV parameter a µτ = 1.0 × 10 −23 GeV (−1.0 × 10 −23 GeV) in theory while considering a µτ = 0 in the prospective data. The left and the right panels show the distribution of ∆χ 2 − and ∆χ 2 + , respectively. In all the panels, we observe that a significant contribution is received from bins with 7 GeV < E rec µ < 17 GeV and cos θ rec µ < −0.4. These are the regions around the oscillation valley as observed in Fig. 1.

Advantage of hadron energy information in constraining LIV parameters
In Fig. 7, we constrain the LIV parameters a µτ , a eµ , and a eτ using 500 kt·yr exposure at the ICAL detector and describe the advantage of incorporating the hadron energy information. In the prospective data, we assume the case of SI where the values of all LIV parameters are considered to be zero, whereas in theory, we consider SI + LIV case where the value of the LIV parameters a µτ , a eµ , and a eτ are varied one-at-a-time. In the left, middle and right panels, we constrain the LIV parameters a µτ , a eµ , and a eτ , respectively. In all the panel, the red curves represent the 2D case where we use only E rec µ and cos θ rec µ variables without considering any hadron energy information. For reconstructed variables of E rec µ and cos θ rec µ , we use binning scheme mentioned in Table 3, where the events are integrated 7 Please note that in the fixed-parameter scenario, we minimize only over the systematic uncertainties, but we keep the oscillation parameters fixed in both theory and data at their benchmark values as mentioned in Table 1. 8 We do not include pull penalty term 5 l=1 ξ 2 l (see Eq. 6.1) while calculating ∆χ 2 − and ∆χ 2 + to explore contributions from each bin in the plane of E rec µ and cos θ rec µ for µ − and µ + events, respectively. over E rec had bins in the hadron energy range of 0 to 25 GeV. On the other hand, the black curves represent the 3D case considering reconstructed variables of E rec µ , cos θ rec µ , and E rec had where we use the hadron energy information available at the ICAL detector. For 3D case, we used the binning scheme given in Table 3 where E rec had is divided into 4 bins in the hadron energy range of 0 to 25 GeV.
We can observe in Fig. 7 that the incorporation of hadron energy information results in the improved constraints for all three cases of LIV parameters a µτ , a eµ , and a eτ . The ICAL detector places the tightest constraint on LIV parameter a µτ among all the three off-diagonal LIV parameters. For the case of a eµ and a eτ , we observe that the bounds are asymmetric with stronger constraints for the positive values than that for the negative values. These observations are consistent with the effect of a eµ and a eτ on ν µ survival probability oscillograms in Figs. 2 and 3, respectively. In section 4, we use Eq. B.4 to explain the reason behind stronger effects of a eµ and a eτ for positive values. In the next [GeV] (fit) 23 10 × [GeV] (fit) 23 10 ×  Figure 7: The sensitivities to constrain the LIV parameters aµτ , aeµ, and aeτ using 500 kt·yr exposure at the ICAL detector as shown in the left, middle, and right panels, respectively. In each panel, the red lines represent the 2D analysis using reconstructed observables (E rec µ , cos θ rec µ ) whereas the black lines refers to the 3D analysis using reconstructed observables (E rec µ , cos θ rec µ , E rec had ). We have used the benchmark value of oscillation parameters given in Table 1. In theory, we have marginalized over oscillation parameters sin 2 θ23, |∆m 2 eff |, and both choices of mass orderings. [GeV] (fit) 23 10 × [GeV] (fit) 23 10 × [GeV] (fit) 23 10 ×  the ICAL detector as shown in the left, middle, and right panels, respectively. In each panel, the black lines represent the case when the charge identification capability of the ICAL detector is used while estimating the sensitivities, whereas the red lines refer to the case when the charge identification capability of ICAL is absent. We have used the benchmark value of oscillation parameters given in Table 1. In theory, we have marginalized over oscillation parameters sin 2 θ23, |∆m 2 eff |, and both choices of mass orderings.
section, we explore the impact of CID capability of the ICAL detector in constraining the LIV parameters.

Advantage of charge identification capability to constrain LIV parameters
The presence of a 1.5 T magnetic field enables the ICAL detector to distinguish µ − and µ + events, which leads to the capability of the ICAL detector to separately identify their parent particles neutrinos and antineutrinos, respectively. In Fig. 8, we discuss the advantage of charge identification capability of the ICAL detector while constraining the CPT-violating LIV parameters a µτ , a eµ , and a eτ using 500 kt·yr exposure. In each panel of Fig. 8, black curves show the sensitivity of the ICAL detector where CID capability is exploited, and two separate sets of bins for reconstructed µ − and µ + events are used. On the other hand, the red curves represent the case where the CID capability of the ICAL detector is not utilized, and a single set of combined bins for both reconstructed µ − and µ + events is considered. We can observe in Fig. 8 that the incorporation of CID capability of the ICAL detector significantly improves the sensitivity to constrain the LIV parameters a µτ , a eµ , and a eτ . In section 4, we discussed that for a given value of a µτ , the oscillation valley bends in the opposite directions for neutrinos and antineutrinos. The combined binning of both reconstructed µ − and µ + events in a single set of bins will lead to a dilution of these features, which are caused due to the presence of non-zero LIV parameter. Thus, it is important to separately identify reconstructed µ − and µ + events to preserve the information about LIV. All these observations validate the fact that the charge identification capability of the ICAL detector will play a crucial role in constraining LIV parameters.
In the subsections 7.2 and 7.3, we discussed the advantage of incorporating hadron energy information and presence of CID while constraining CPT-violating LIV parameters a µτ , a eµ , and a eτ . Now, Table 4 nicely summarizes the findings from these two studies in tabular form, and at the same time, we compare the performance of ICAL with other experiments. In Table 4, we present the constraints on LIV parameters a µτ , a eµ , and a eτ at 95% C.L. using 500 kt·yr exposure at the ICAL detector. By comparing the results shown in the second row with respect to that in the first row, we can infer that the incorporation of hadron energy information improves the bounds on all of these LIV parameters. The improvement due to the presence of CID capability can be observed by comparing the second row with respect to the third row. While doing comparison with existing constraints, we have mentioned the results from Super-K analysis [22] for LIV parameters a µτ , a eµ , and a eτ at 95% C.L. with real and imaginary parts separately. The IceCube collaboration [25] has performed analysis only for a µτ . We show constraints for the real and imaginary parts of a µτ separately at 99% C.L. as provided by the IceCube collaboration.

Impact of marginalization on constraining LIV parameters
In the previous sections, we have performed complete marginalization over the uncertain neutrino oscillation parameters in their 3σ ranges given by the present global fit study, but we expect in the coming decades, the precisions of neutrino oscillation parameters are going to improve because of the currently running and upcoming experiments. To understand the impact of uncertainties in oscillation parameters on our results, in this section, we marginalize over oscillation parameters while constraining the CPT-violating LIV parameters a µτ , a eµ , and a eτ using the ICAL detector with 500 kt·yr exposure and compared this with the fixed-parameter scenario as shown in Fig. 9. The red curves in Fig. 9 represent the fixed-parameter case where the oscillation parameters are kept fixed 9 In Super-K analysis [22], the authors have scanned the parameters on a logarithmic scale where the only positive values of parameters are considered assuming 10 −28 GeV to be the minimum value which is equivalent to the case of no LIV. 10 To compare with our results, we mention the dimension-three results for aµτ as given by IceCube [25]. Super-K 9 (95% C.L.) [22] Re(a µτ ) < 0.65 Re(a eµ ) < 1.8 Re(a eτ ) < 4.1 Im(a µτ ) < 0.51 Im(a eµ ) < 1.8 Im(a eτ ) < 2.8 IceCube 10 (99% C.L.) [25] |Re(a µτ )| < 0.29 --|Im(a µτ )| < 0.29 Table 4: Constraints on the CPT-violating LIV parameters aµτ , aeµ, and aeτ at 95% C.L. using 500 kt·yr exposure at the ICAL detector. The second row shows results using 2D variables (E rec µ , cos θ rec µ ) with CID, whereas the third row present the results for the case of 3D variables (E rec µ , cos θ rec µ , E rec had ) with CID. The results in the third row use 3D variables but do not use the CID capability of the ICAL detector. We have used the benchmark value of oscillation parameters given in Table 1. Note that we have marginalized over oscillation parameters sin 2 θ23, |∆m 2 eff |, and both choices of mass orderings in theory. For comparison, the last two rows mention the existing constraints on LIV parameters obtained from Super-K and IceCube experiments. All these results consider only the time component assuming isotropic nature.
[GeV] (fit) 23 10 ×  [GeV] (fit) 23 10 ×   Table 1. in theory. Please note that in the fixed-parameter scenario, we do minimize the systematic errors using the pull method. The black curves show the case where we marginalize over [GeV] (fit) 23 10 × [GeV] (fit) 23 10 × [GeV] (fit) 23 10 ×  Figure 10: The sensitivities to constrain the LIV parameters aµτ , aeµ, and aeτ using 500 kt·yr exposure at the ICAL detector as shown in the left, middle, and right panels, respectively. In each panel, the red, black, and green curves represent the cases when the true value of sin 2 θ23 is taken as 0.4, 0.5, and 0.6, respectively. We have used the benchmark value of other oscillation parameters given in Table 1. In theory, we have marginalized over oscillation parameters sin 2 θ23, |∆m 2 eff |, and both choices of mass orderings.
oscillation parameters sin 2 θ 23 , |∆m 2 eff |, and both choices of mass orderings in theory as explained in section 6.
We can observe in the left panel in Fig. 9 that the marginalization over oscillation parameters do not affect the constraints by ICAL on LIV parameter a µτ by a large amount. The constraints on LIV parameters a eµ and a eτ show some deterioration after marginalization over oscillation parameters. The largest impact of uncertainties of oscillation parameters can be observed for the case of a eτ as shown in the right panel of Fig. 9. Thus, we can infer that the more precise determination of the oscillation parameters in the future will improve the constraints on LIV parameters a eµ and a eτ using 500 kt·yr exposure at ICAL.

Impact of true values of sin 2 θ 23 on constraining LIV parameters
In the results presented so far, we have assumed sin 2 θ 23 = 0.5 which corresponds to the maximal mixing , i.e. θ 23 = 45 • . The current global fit of neutrino data indicates that θ 23 may be non-maximal where θ 23 can be found in lower octant for θ 23 < 45 • or higher octant for θ 23 > 45 • . At present, θ 23 is the most uncertain oscillation parameter apart from δ CP . Thus, it is an important question to ask how the constraints on CPT-violating LIV parameters a µτ , a eµ , and a eτ change, if θ 23 (true) is found to be non-maximal in nature. To answer this question, we present Fig. 10 where we show the impact of non-maximal θ 23 on constraining CPT-violating LIV parameters a µτ , a eµ , and a eτ as shown in the left, middle, and right panels, respectively. In each panel, we demonstrate constraints for three different true values of sin 2 θ 23 = 0.4 (red curves), 0.5 (black curves), and 0.6 (blue curves).
The left panel of Fig. 10 shows that the constraints on LIV parameter a µτ deteriorates for non-maximal values of θ 23 (true) where sin 2 θ 23 (true) = 0.4 and 0.6. We can also observe that the constraints on a µτ are the same for sin 2 θ 23 (true) = 0.4 and 0.6 which is consistent with the Eq. 2.5 where the term containing a µτ in P (ν µ → ν µ ) is proportional to sin 2 2θ 23 . The middle and right panels in Fig. 10 illustrate that the constraints on both a eµ , and a eτ improves for sin 2 θ 23 (true) = 0.6 and deteriorates for sin 2 Table 5: The sensitivity of the ICAL detector to determine mass ordering with 500 kt·yr exposure.
For the SI case (first row), we do not consider LIV in data and theory. For the cases of LIV parameters, we introduce aµτ (second row), aeµ (third row), and aeτ (fourth row) in the theory one-at-a-time and marginalize over them along with oscillation parameters sin 2 θ23 and ∆m 2 eff while assuming SI case in data. In the third and fifth columns, we show how much the mass ordering sensitivity deteriorates due to the presence of LIV parameters compared to the SI case. We present our results assuming NO (IO) in data as given in the second and third (fourth and fifth) columns. We have used the benchmark value of oscillation parameters given in Table 1.
can be explained using Eq. B.4 for appearance channel which can also be translated to the effect on ν µ survival channel following the similar arguments as in section 4. The fifth term in Eq. B.4 with significant contribution to the effect of a eµ is proportional to sin 3 θ 23 . In a similar fashion, the seventh term in Eq. B.4 having dominant effect for a eτ is proportional to sin 2 θ 23 cos θ 23 . From all these observations in Fig. 10, we can conclude that if θ 23 is found to be lying in the higher octant in nature then the sensitivity of ICAL for constraining a eµ , and a eτ will enhance whereas it will reduce for a µτ .

Impact of non-zero LIV parameters on mass ordering determination
In this section, we are going to study the impact of non-zero LIV parameters a µτ , a eµ , and a eτ on the sensitivity of the ICAL detector to determine the neutrino mass ordering. For statistical analysis, we simulate the prospective data assuming a given mass ordering. The sensitivity of ICAL to rule out the wrong mass ordering is calculated in the following fashion: where, we calculate χ 2 ICAL (true MO) and χ 2 ICAL (false MO) by fitting the prospective data assuming true and false mass ordering, respectively. Since, the statistical fluctuations are suppressed while calculating the median sensitivity, we have χ 2 ICAL (true MO) ∼ 0. We calculate the sensitivity of ICAL to measure the neutrino mass ordering with 500 kt·yr exposure using neutrino flux at the INO site following the procedure mentioned in Ref. [85]. The sensitivity towards the neutrino mass ordering is found to be 7.55 (7.48) for true NO (IO) as mentioned in the first row of Table 5.
In order to estimate the impact of non-zero LIV parameters, we generate the prospective data with a given mass ordering assuming no LIV where a µτ = 0, a eµ = 0, and a eτ = 0. Then, we fit the prospective data with opposite mass ordering assuming non- zero LIV parameter a µτ in the theory and perform marginalization over a µτ in the range 11 [−0.23, 0.23] × 10 −23 GeV along with the oscillation parameters sin 2 θ 23 and ∆m 2 eff in the ranges as mentioned in section 6. The mass-ordering sensitivity with 500 kt·yr exposure at ICAL for the cases of a µτ is mentioned in the second row of Table 5. Similarly, the third (fourth) row in Table 5 shows the sensitivity of ICAL towards neutrino mass ordering in the presence of a eµ (a eτ ) while perfoming marginalization over the range of [−2.0, 1.5] × 10 −23 GeV ([−2.8, 1.6] × 10 −23 GeV). In the third and fifth columns, we show the deterioration in the sensitivity for determining mass ordering due to the presence of non-zero LIV parameters with respect to the SI case. We observe in Table 5 that depending upon true mass ordering, the results deteriorate by 15 to 50 % due to the presence of non-zero LIV parameters when considered one parameter at-a-time. 2 32 |sin 2 θ 23 ) plane with non-zero LIV parameters Now, we explore the impact of CPT-violating LIV parameters a µτ , a eµ , and a eτ on the precision measurement of atmospheric oscillation parameters sin 2 θ 23 and |∆m 2 32 | using the ICAL detector with 500 kt·yr exposure. First of all, we simulate the prospective data assuming sin 2 θ 23 (true) = 0.5 and |∆m 2 32 |(true) = 2.46 × 10 −3 eV 2 for case of SI with no LIV. The statistical significance for this is quantified using the following expression:

Allowed regions in (|∆m
where, χ 2 ICAL (sin 2 θ 23 , |∆m 2 32 |) is estimated by fitting the prospective data with theory for a given value of sin 2 θ 23 , and |∆m 2 32 |. χ 2 0 is the minimum value of χ 2 ICAL (sin 2 θ 23 , |∆m 2 32 |) in the allowed range of oscillation parameters sin 2 θ 23 and |∆m 2 32 |. Here, χ 2 0 ∼ 0 because the statistical fluctuations are suppressed. We estimate the allowed regions in (|∆m 2 32 |sin 2 θ 23 ) plane for the case of SI at 95% C.L. (2 d.o.f.) as shown by black line in Fig. 11. Now, we discuss the impact of non-zero LIV parameters a µτ , a eµ , and a eτ one-at-atime on the precision measurement of atmospheric oscillation parameters. We simulate the prospective data considering the true values of sin 2 θ 23 and |∆m 2 32 | as mentioned above without LIV. Then, we estimate allowed regions in (|∆m 2 32 | -sin 2 θ 23 ) plane with non-zero LIV parameters a µτ , a eµ , and a eτ one-at-a-time in the fit. Further, we marginalize over a µτ in the range [−0.23, 0.23] × 10 −23 GeV, a eµ in the range [−2.0, 1.5] × 10 −23 GeV, and a eτ in the range [−2.8, 1.6] × 10 −23 GeV one-at-a-time in the fit (see footnote 11). We show these results for the cases of non-zero LIV parameters a µτ , a eµ , and a eτ at 95 % C.L. (2 d.o.f.) using dotted red, dashed green, and dash-dotted blue lines, respectively, in Fig. 11. We do not observe a significant change in the allowed regions at 95% C.L. when we introduce a µτ in fit and marginalize over it. It indicates that the precision measurement of atmospheric oscillation parameters using the ICAL detector is quite robust against the presence of a µτ in the fit. However, the presence of non-zero LIV parameter a eµ in fit deteriorate the allowed region in the direction of sin 2 θ 23 . As far as the presence of LIV parameter a eτ in fit is concerned, it deteriorates allowed region in the direction of |∆m 2 32 |.

Correlation between off-diagonal LIV parameters
So far, we have studied the neutrino oscillations in the presence of only one CPT-violating LIV parameter at-a-time. In this section, we explore the correlation between different CPT-violating LIV parameters. To begin with, we generate prospective data with 500 kt·yr exposure at the ICAL detector considering no LIV with NO as true mass ordering and benchmark value of oscillation parameters mentioned in Table 1. While constraining two LIV parameters at-a-time, we quantify the statistical significance using the following expression: where, χ 2 ICAL (SI) and χ 2 ICAL (SI + (a µτ , a eµ )) are calculated by fitting the prospective data with only SI case (no LIV) and with SI + LIV (a µτ , a eµ ) case, respectively. Since, the statistical fluctuations are suppressed, we have χ 2 ICAL (SI) ∼ 0. Then, we estimate the allowed regions in the plane of CPT-violating LIV parameters (a µτ , a eµ ) at 95% C.L. (2 d.o.f.) as shown in the left panel of Fig. 12. The true choices of LIV parameters a µτ = 0, a eµ = 0, and a eτ = 0 in data are depicted as red dots in Fig. 12. In theory, we vary two LIV parameters at-a-time with (without) marginalization over the oscillation parameters [GeV] (fit) 23 10 × [GeV] (fit) 23 10 × τ µ [GeV] (fit) 23 10 × [GeV] (fit) 23 10 × τ µ [GeV] (fit) 23 10 ×  Table 1. In theory, we vary two LIV parameters at-a-time with marginalization over the oscillation parameters sin 2 θ23, ∆m 2 eff , and both choices of mass ordering as shown by black curves. The blue curves represent the fixed-parameter scenario where we do not marginalize over any oscillation parameter in theory, but we do marginalize over systematic uncertainties. sin 2 θ 23 , ∆m 2 eff , and both choices of mass ordering as shown by black (blue) curve. Similarly, we present the allowed regions in the plane of LIV parameters (a µτ , a eτ ) and (a eτ , a eµ ) in the middle and left panels, respectively, in Fig. 12.

Summary and concluding remarks
In this paper, for the first time, we explore the Lorentz Invariance Violation (LIV) through the mass-induced flavor oscillations of atmospheric neutrino in the multi-GeV range of energies over a wide range of baselines using the proposed ICAL detector at INO. In the present analysis, we focus our attention on the isotropic CPT-violating LIV parameters a µτ , a eµ , and a eτ where we consider their real values only with both positive and negative signs.
We perform a detailed analysis by demonstrating the effect of non-zero CPT-violating LIV parameters (a µτ , a eµ , and a eτ ) one-at-a-time on probability oscillograms of the disappearance (ν µ → ν µ ) channel, which has a dominant contribution (more than 98%) to the reconstructed muon events at the ICAL detector. We observe that the vacuum oscillation valley bends due to the presence of non-zero LIV parameter a µτ where the curvature of valley depends on the sign of a µτ . For a given a µτ , the curvatures of bendings are opposite for neutrino and antineutrino. This fact indicates that the charged identification capability of the ICAL detector is a crucial property while probing LIV parameter a µτ . If we turn our attention to the other off-diagonal LIV parameters a eµ and a eτ , we observe that they do not result in any bending of oscillation valley, but they do disturb the vacuum oscillation valley as well as matter effect regions.
Next, we present the impact of LIV on the reconstructed event distributions at the ICAL detector for 500 kt·yr exposure. We observe that a µτ affects the event distributions by the largest amount among all off-diagonal CPT-violating LIV parameters. Since the effects of LIV parameters depend on energy and direction, the good energy and direction resolution of the ICAL detector is going to play an important role.
After describing the modification in reconstructed event distributions due to LIV, we calculate statistical significance for the presence of LIV using the ICAL detector at the χ 2 level. In order to calculate ∆χ 2 , we simulate the prospective data at ICAL with 500 kt·yr exposure for SI case with no LIV and perform fitting assuming SI+LIV case in theory where we consider the CPT-violating LIV parameters a µτ , a eµ , and a eτ one-at-a-time. First of all, we identify the energy and direction bins that contribute significantly to ∆χ 2 . Further, we constrain the LIV parameters a µτ , a eµ , and a eτ at 95% C.L. using 500 kt·yr exposure at the ICAL detector.
We show the advantage of incorporating hadron energy information and the CID capability of ICAL while constraining LIV parameters. We also demonstrate that the constraints on a µτ are robust against the uncertainties in oscillation parameters where the estimated bounds on a eµ and a eτ are expected to improve with more precise determination of oscillation parameters. As far as the non-maximal value of θ 23 is concerned, it deteriorates the constraints on a µτ whereas constraints on a eµ and a eτ improve (worsen) if θ 23 is found to be lying in the higher (lower) octant.
The aim of the ICAL detector at INO is to determine the mass ordering and the precise determination of atmospheric oscillation parameters θ 23 and |∆m 2 32 |. In the present analysis, we also explore how the sensitivity to measure mass ordering gets affected due to the presence of non-zero LIV parameters considered one-at-a-time. We observe that depending upon true mass ordering, the results get deteriorated by about 15 to 50% due to the presence of LIV. If we look at the effect of LIV on precision measurement of atmospheric oscillation parameters θ 23 and |∆m 2 32 |, we find that the precision measurement is quite robust against the marginalization over a µτ in the fit. However, the presence of a eµ (a eτ ) worsens the allowed region in the direction of sin 2 θ 23 (|∆m 2 32 |). The above-mentioned analyses consider one LIV parameter at-a-time. In order to explore the correlation among various off-diagonal CPT-violating LIV parameters, we considered two LIV parameter at-a-time and constrain (2 d.o.f.) them using 500 kt·yr exposure at the ICAL detector. We have obtained allowed regions in three different planes of (a µτ , a eµ ), (a µτ , a eτ ) and (a eτ , a eµ ) while marginalizing over oscillation parameters.
Before we conclude, we would like to highlight the fact that the present analysis has been performed in the multi-GeV range of energies which is complementary to energy ranges used in the analyses of Super-K and IceCube. We would also like to emphasize that the LIV bounds obtained in the present analysis using the ICAL detector are quite competitive compared to the existing bounds from Super-K and IceCube due to the charge identification capability of ICAL as well as the good resolutions in reconstructed muon energy and direction. In the absence of CID, the constraints on LIV parameters deteriorate significantly due to the dilution of LIV features while combining µ − and µ + events. This unique capability of CID enables ICAL to probe LIV separately in neutrino and antineutrino modes which may not be possible at any other existing or planned neutrino detector based on water, ice, or argon.

A Some properties of gauge invariant LIV parameters
In this section, we describe the properties of CPT-violating and conserving LIV parameters. Let us consider a set of Dirac spinors representing the neutrino field, ψ = {ν e , ν µ , ν τ } and corresponding charge conjugated field, ψ c = {ν c e , ν c µ , ν c τ }. Here, the charge conjugated field transforms as ψ c = Cψ T , and C is the charge conjugation operator. Now, lets define a Langragian density L = L o −L , where L o represents the weak interactions in the Standard Model and L represents the new physics interactions induced due to the violation of Lorentz symmetry. In order to address the new physics interactions, we are adopting the spontaneous Lorentz violation mechanism that is proposed in the string theory. At the Planck scale (M p ), spontaneous breaking of the Lorentz symmetry can be expected when a tensor field of non-perturbative vacuum in the proposed model acquires the non-zero vacuum expectation values (VEVs). These VEVs effectively can act as a fixed background to a given observer's frame of reference. The interactions induced due to this background can have the boost dependency which breaks the Lorentz symmetry 12 .
At low energy, the induction strength of Lorentz violation is expected to be suppressed by an order of 1/M p [19,20,112,161]. In the low-energy effective field theory, the LIV interaction can have many coupling coefficients which effectively maintain the power-counting renormalizability as follows [9]: where, λ stands for a dimensionless coupling constant, M p denotes the Planck mass scale, T are the non-zero tensor VEVs, Γ represent some gamma-matrix structure, and k is an integer power where the dominant terms with k ≤ 1 are renormalizable. Considering only the gauge-invariant terms, the CPT-violating and CPT-conserving parameters at low energies are obtained by choosing the mass dimension to be zero and one, respectively. For k = 0, T ∼ m 2

Mp
and for k = 1, T ∼ m, where m is the mass of fermionic field, ψ. The CPT-violating and CPT-conserving terms appear in the following fashion [17]: Here, the quantities a µ , b µ , c µν , and d µν are real because of their mode of origin via spontaneous symmetry breaking mechanism followed by the hermitian nature of the underlying theory [20]. Note that the running of these Lorentz-violating couplings between the Planck scale and the low-energy scale is discussed in Ref. [162]. Now, by interchanging the space-time indices µ ↔ ν, we get for the CPT-conserving parameters Note, L CPT−conserving must be invariant under the interchange of space-time indices (µ ↔ ν), which implies that c µν and d µν are the antisymmetric tensors, i.e., Thus, L is defined as [27,121,135], where, a µ and b µ are the CPT-violating LIV parameters, whereas c µν and d µν are the CPT-conserving LIV parameters. The 1/2 factor in Eq. A.10 appears due to the cannonical renormalization of the neutrino field.
A single field χ consisting of the ψ and ψ c fields can be defined as The Lagrangian L is defined in terms of field χ as [112] where, 14) The hermiticity demands that Γ µ = γ 0 Γ † µ γ 0 , and M = γ 0 M † γ 0 . It is expected that at low energies, the LIV observables don't have the significant effects on the SM weak interactions. Hence, the active neutrino can be treated as left-handed and the antineutrino as righthanded. The structure of the equation of motion should remain unaltered under the charge conjugation operation, which is denoted by the operator, C = iγ 2 γ 0 . Under the charge conjugation operation, Γ µ = -C(Γ µ ) T C −1 and M = CM T C −1 , which further implies Now, expanding A.15a, Similarly, expanding A.15b, Following the same procedure for A.15c and A.15d, we obtain, The experiments to probe LIV parameters can be categorized on the phenomena of i) Coherent, ii) Interferometric, or iii) Extreme effects. For the mass-induced neutrino oscillations, the LIV parameters have interference effects. Thus, the individual CPT-violating parameters a µ and b µ , as well as the CPT-conserving parameters c µν and d µν , are not the observable quantities. These interfering parameters can be redefined as the observable quantities decomposed along with the polarization of active left-handed neutrinos and right-handed antineutrinos. Therefore, we define where, indices are supressed for simplicity. The relation between a L and a R can be evaluated using A.16e and A.17f, respectively. The expectation value of a L is χ|a L |χ . Hence, we obtain Similarly using A.18a and A.18b, In simple terms, we can conclude that when we transform from neutrino to antineutrino, the CPT-violating parameters flip their signs, whereas CPT-conserving parameters do not. This is an important property of LIV parameters for the case of neutrino oscillations.
Another crucial property of the CPT-conserving LIV parameters can be explored by choosing a rotational invariant coordinate with isotropic condition. Now, using the properties of Eq. A.7 and Eq. A.8, we obtain T r c µν L g µν = 0, where g µν is the metric tensor, that implies c 00 L = 1,3 i c ii L . Applying the condition for isotropic symmetry, we got B Effect of a eµ and a eτ on appearance channel P (ν e → ν µ ) The CPT-violating LIV parameters a eµ and a eτ affects the appearance channel (ν e → ν µ ) significantly and have contribution at leading orders. The authors in Ref. [107] have presented the approximate analytic expression for P (ν µ → ν e ) in the presence of neutralcurrent NSI parameters ε eµ and ε eτ which occur during neutrino propagation. We use the Eq. [33] in Ref. [107] to obtain approximate expression for P (ν e → ν µ ) in the presence of the LIV parameters a eµ and a eτ where we replace ε αβ with a αβ with the help of Eq. 2.4, where, a CC = 2 √ 2G F N e E ν . We assume δ CP = 0 and focus on real a αβ by taking the LIV phase φ αβ = 0, and π which ensures that both positive and negative values of a αβ are considered and we have, a eµ ≡ |a eµ | cos(φ eµ ) (B.2) a eτ ≡ |a eτ | cos(φ eτ ). (B. 3) The resulting expression for P (ν e → ν µ ) in the presence of the LIV parameters a eµ and a eτ is given as, where, s ij = sin θ ij , c ij = cos θ ij , s 2×ij = sin 2θ ij , c 2×ij = cos 2θ ij . For CPT-violating parameter a αβ , the parameter ω = +1 for the case of neutrino whereas for antineutrino ω = −1. The effective mixing angle θ 13 in matter is defined as, In Equation B.4, the first three terms are driven by standard interactions where the first term has the dominant effect. The remaining terms contain the effect of CPT violating-LIV parameters a eµ and a eτ . Now let us try to understand the effect of a eµ on appearance channel P (ν e → ν µ ) from the approximate expression in Eq. B.4. We observe that the effect of a eµ is dominantly contributed by the fifth term, + 16ωa eµ E νs13 s 3 where, (∆m 2 31 −a CC ) factor in the denominator causes the matter-driven resonance effect for the case of neutrino. Since this term has positive sign for the case of neutrino (ω = +1), Figure 13: The distribution of fixed-parameter ∆χ 2 − (∆χ 2 + ) in the plane of (E rec µ , cos θ rec µ ) without pull penalty term using 500 kt·yr exposure of the ICAL detector as shown in the left (right) panels. Note that ∆χ 2 − and ∆χ 2 + are presented in the units of GeV −1 sr −1 where we have divided them by 2π × bin area. In data, aeµ = 0 with NO (true) using the benchmark value of oscillation parameters given in Table 1. In theory, aeµ = −2.5 × 10 −23 GeV and 2.5 × 10 −23 GeV in the top and bottom panels, respectively. Here, in the fixed-parameter scenario, we do not marginalize over any oscillation parameter in theory, but we do marginalize over systematic uncertainties. the positive (negative) value of a eµ increases (decreases) P (ν e → ν µ ). We can see that the effect of a eµ is larger in the case of neutrino than for antineutrino. This happens because a CC becomes negative for antineutrino, and matter-driven resonance condition is not fulfilled, due to which the above-mentioned term does not contribute significantly. A similar effect is also observed for the case of a eτ which can be explained using the seventh term in Eq. B.4. C Effective regions in (E rec µ , cos θ rec µ ) plane to constrain a eµ and a eτ In this section, we identify the regions in (E rec µ , cos θ rec µ ) plane which contribute significantly towards the sensitivity of ICAL for the CPT-violating LIV parameters a eµ and a eτ with 500 kt·yr exposure. Note that here, we consider one LIV parameter at-a-time. In Fig. 13, we Figure 14: The distribution of fixed-parameter ∆χ 2 − (∆χ 2 + ) in the plane of (E rec µ , cos θ rec µ ) without pull penalty term using 500 kt·yr exposure of the ICAL detector as shown in the left (right) panels. Note that ∆χ 2 − and ∆χ 2 + are presented in the units of GeV −1 sr −1 where we have divided them by 2π × bin area. In data, aeτ = 0 with NO (true) using the benchmark value of oscillation parameters given in Table 1. In theory, aeτ = −2.5 × 10 −23 GeV and 2.5 × 10 −23 GeV in the top and bottom panels, respectively. Here, in the fixed-parameter scenario, we do not marginalize over any oscillation parameter in theory, but we do marginalize over systematic uncertainties.
show the distribution of sensitivity of ICAL for a eµ in terms of fixed-parameter ∆χ 2 − (∆χ 2 + ) without pull penalty term contributed from reconstructed µ − (µ + ) events in (E rec µ , cos θ rec µ ) plane following the same procedure as described in section 7.1 for a µτ . While estimating sensitivity of ICAL for a eµ in the top (bottom) panels in Fig. 13, we assume non-zero value of a eµ = −2.5 × 10 −23 GeV (+2.5 × 10 −23 GeV) in theory while taking a eµ = 0 in data for the case of SI. We keep the oscillation parameters fixed in theory and data at their benchmark values as mentioned in Table 1. We consider 500 kt·yr exposure at ICAL with NO as true mass ordering. We use the binning scheme mentioned in Table 3 where we add contribution from all E rec had bins for each (E rec µ , cos θ rec µ ) bins. The left and right panels show the distribution of ∆χ 2 − and ∆χ 2 + , respectively. We can observe in Fig. 13 that the contribution towards the sensitivity of ICAL for a eµ mainly stems from µ − as shown in left panels, whereas the contribution from µ + is negligible as shown in the right panels. This happens because the term with a significant contribution of a eµ has a matter-driven resonance effect for the case of neutrino as described in section 4 while explaining the oscillograms in Fig. 2. We also observe in the left panels in Fig. 13 that ∆χ 2 − is significantly larger for the positive value (top left panel) of a eµ compared to that for negative value (bottom left panel). This observation is also consistent with the asymmetric effect of a eµ on ν µ survival oscillograms shown in Fig. 2, and the reason behind this effect is explained in detail in section 4. We would like to highlight the observation that the sensitivity of ICAL for a eµ is significantly contributed by the reconstruction muon energy and direction bins corresponding to the region of matter effect, which is not surprising because the impact of a eµ is being driven by the resonance effect with the matter as explained in section 4.
Similar to the case of a eµ , we present the distribution of sensitivity of ICAL for a eτ in terms of fixed-parameter ∆χ 2 − (∆χ 2 + ) with contribution from reconstructed µ − (µ + ) events in (E rec µ , cos θ rec µ ) plane in Fig. 14. To estimate the sensitivity of ICAL for a eτ in the top (bottom) panels in Fig. 14, we take non-zero value of a eτ = −2.5 × 10 −23 GeV (+2.5 × 10 −23 GeV) in theory while keeping a eτ = 0 in data for the case of SI. Analogous to the case of a eµ , we observe that the sensitivity for a eτ is more in the case of µ − (left panels) than that for µ + (right panels). The asymmetric effect of a eτ leads to higher sensitivity for positive a eτ as shown in the bottom right panel of Fig. 14. These features can be explained in the same fashion as we did for a eµ in the previous paragraph. The reasons behind these effects of a eτ on ν µ survival oscillograms are described in section 4.