Sensitivity to neutrino decay with atmospheric neutrinos at INO

Sensitivity of the magnetised Iron CALorimeter (ICAL) detector at the proposed India-based Neutrino Observatory (INO) to invisible decay of the mass eigenstate $\nu_3$ using atmospheric neutrinos is explored. A full three-generation analysis including earth matter effects is performed in a framework with both decay and oscillations. The wide energy range and baselines offered by atmospheric neutrinos are shown to be excellent for constraining the $\nu_3$ lifetime. We find that with an exposure of 500 kton-yr the ICAL atmospheric experiment could constrain the $\nu_3$ lifetime to $\tau_3/m_3>1.51\times10^{-10}$ s/eV at the 90\% C.L. This is two orders of magnitude tighter than the bound from MINOS. The effect of invisible decay on the precision measurement of $\theta_{23}$ and $|\Delta{m^2_{32}}|$ is also studied.


I. INTRODUCTION
Neutrino oscillation experiments spanning various energy ranges and baselines have helped in establishing the fact that neutrinos oscillate from one flavor to another. Most of the neutrino oscillation parameters have been pinned down and are now known rather precisely. 1 The main open questions remaining in neutrino oscillation physics are neutrino mass hierarchy, octant of the mixing angle θ 23 and the value of the CP phase δ CP . Several experiments are running or are being planned in order to answer the above-mentioned questions. The leading experimental proposals for the future include the long-baseline experiments DUNE [3] and T2HK [4], reactor experiments JUNO [5] and RENO50 [6], and atmospheric neutrino experiments PINGU [7], ORCA [8] and ICAL [9,10]. It is expected that the neutrino oscillation probabilities would change in the presence of new physics. This could be used to constrain new physics scenarios at neutrino oscillation experiments. At the same time, a given new physics scenario could also interfere with the measurement of the standard neutrino oscillation parameters and hence pose a challenge to the proposed experiments, unless ways are found to cancel out their effects through synergistic measurements at multiple experiments. One such new physics scenario is the decay of neutrino during its flight from the source to the detector.
While there is no observational evidence in support for unstable neutrinos, since they are massive, its not unlikely that they would decay. Radiative decays of neutrinos are very severely constrained by cosmological data. Since the measured neutrino masses suggest that the neutrinos would radiatively decay in the microwave energy range, the most stringent bounds are provided by cosmic microwave background data [11], making radiative decay of neutrinos totally uninteresting for neutrino oscillation experiments. However, there still remains the possibility that neutrinos could decay into a lighter fermion state and a beyond standard model boson. The Majoron model [12][13][14] for instance allows the following decay modes for Majorana neutrinos: ν i → ν j +J or ν i →ν j +J, where ν j andν j are lighter neutrino and anti-neutrino states and J is a Majoron. The Majoron in principle could belong to either a singlet or a triplet representation of the standard model gauge group. But the triplet model is severely constrained [13,14]and hence J must predominantly be an electroweak singlet. If the final state fermion is a lighter active neutrino, the decay is called visible decay. On the other hand, if the final state fermion is a sterile state with no standard model interaction, then the decay scenario is termed invisible decay. Even for Dirac neutrinos in extensions of the standard models one could write down terms in the Lagrangian coupling neutrinos with a light scalar boson and light right-handed neutrinos allowing the decay mode ν i →ν iR + χ, whereν iR is a right-handed singlet neutrino and χ is an iso-singlet scalar carrying lepton number +2 [15,16]. In this paper, we will work in a scenario where the final state particles remain invisible to the detector.
The lifetime of ν 2 (and ν 1 ) is constrained by the solar neutrino experiments. Neutrino decay as a solution to the solar neutrino deficit problem was suggested in [17], however, now we know that neutrino decay alone cannot explain this deficit. Attempts to constrain the neutrino lifetime by considering neutrino decay as a subdominant effect along with the leading LMA-MSW solution was done in [18][19][20][21][22][23][24][25][26]. Most of these studies considered the invisible decay scenario. Since U e3 is small, the ν e state mostly resides in the ν 2 and ν 1 states and hence all of these studies worked in the two-generation framework. Bounds on the lifetime of ν 2 was obtained from a global analysis of solar neutrino data in [24] where the impact of the Sudbury Neutrino Observatory neutral current data was highlighted. It was shown that the bound on ν 2 lifetime was τ 2 /m 2 > 8.7 × 10 −5 s/eV at 99% C.L. for a 3 parameter fit. This bound was revisited in [25] (see also [26]) where the authors obtained the 95% C.L. limit τ 2 /m 2 > 7 × 10 −4 s/eV for both normal and inverted mass hierarchy and τ 1 /m 1 > 4 × 10 −3 s/eV for inverted mass hierarchy. These results are very consistent with the earlier analysis of [24] where the 95% C.L. limit for a one parameter fit is seen to be τ 2 /m 2 > 4.4 × 10 −4 s/eV. The corresponding constraints from SN1987A are stronger [27].
Limits on the lifetime of ν 3 come from the atmospheric and long-baseline neutrino experiments. Like in the case of solar neutrinos, any fit with neutrino decay alone [28,29] is unable to explain the atmospheric neutrino zenith angle data. A lot of work has gone into considering decay along with oscillations. The analyses can be broadly classified into two categories depending on the model used. If one considers decay of ν 3 to a state with which it oscillates, then the bounds coming from K-decays [30] restrict the corresponding mass squared difference between them to ∆m 2 > 0.1 eV 2 [31]. However, if the state to which ν 3 decays is a sterile state then the ∆m 2 driving the leading oscillations of ν µ is unconstrained. The former case is that of decay to active neutrinos and was studied in the context of atmospheric neutrinos in [31,32] and no good fit was found. The latter is the invisible decay scenario to sterile neutrinos and was analysed against the atmospheric neutrino data in [33][34][35][36]. The invisible decay case can be again classified into two. In one case we can make the assumption that ∆m 2 ≪ 10 −4 eV 2 , causing it to drop out of the oscillation probability. The authors of [34] argued that this could explain the Super-Kamiokande atmospheric neutrino data, however, the Super-Kamiokande collaboration itself reported [35] that this scenario was not supported by their data. The other case of invisible decay is when ∆m 2 is left free in the fit to be determined by the data. This case was first proposed by some of us in [33]. The results of [33] were updated in [36] where the authors obtained the limit τ 3 /m 3 > 2.9 × 10 −10 s/eV for invisible decay at the 90% C.L. from a combined analysis of Super-Kamiokande atmospheric and MINOS data. More recently, the analysis of oscillation plus invisible decay scenario with unconstrained ∆m 2 was performed in [37] in the context of MINOS and T2K data and gave a bound τ 3 /m 3 > 2.8 × 10 −12 s/eV at 90% C.L. The constraint for the visible decay scenario using the MINOS and T2K charged as well as neutral current data was performed in [38]. The bounds on neutrino lifetime could be improved considerably by observations at IceCube using cosmological baselines [39][40][41][42].
All the above mentioned papers which considered neutrino decay alongside oscillations performed their analysis in the framework of two-generations and did not take earth matter effects into account. Recently a three-generation analysis including earth matter effect and decay in the context of the Deep Underground Neutrino Experiment (DUNE) was performed in [43] for visible decays and [44] for invisible decays. It was shown that DUNE could improve the bound on τ 3 /m 3 for the invisible decay case by at least an order of magnitude compared to the current limits from MINOS and T2K. In this work, we consider invisible neutrino decay within a three-generation oscillation framework in the context of atmospheric neutrinos and include earth matter effects. Atmospheric neutrinos span many orders of magnitude in energy and baseline. Since the effect of neutrino decay increases for lower energies and longer baselines, atmospheric neutrino experiments are expected to give a tighter bound on τ 3 /m 3 than the proposed long-baseline experiments. We will study the sensitivity of the atmospheric neutrinos at INO to neutrino decay.
The India-based Neutrino Observatory (INO) is a proposed underground laboratory in India, which plans to house a 50 kton magnetised Iron CALorimeter (ICAL). The detector will be mainly sensitive to muon type neutrinos, which are detected through the observation of a muon track and the accompanying hadron shower in a charged current interaction. The detector response to muons [45][46][47][48] and hadrons [49][50][51][52] have been performed via the Geant4based [53][54][55] detector simulation code for ICAL. This detector owing to its magnetisation can distinguish between neutrino and anti-neutrino events which makes it an excellent detector to determine the neutrino mass hierarchy [9,[56][57][58][59][60]. ICAL will also perform precision measurements of |∆m 2 32 | and the mixing angle θ 23 [9,58,59,[61][62][63][64]. In addition, there are a variety of new physics scenarios which could be constrained and/or discovered at ICAL. Some of the new physics scenarios studied by the INO collaboration include, CPT violation [65], dark matter [66], non-standard neutrino interactions [67] and sterile neutrino oscillations [68]. In this work we will study in detail the sensitivity of ICAL to invisible neutrino decay using the full physics analysis simulation framework of ICAL. We will also study the effect of invisible neutrino decay on the precision measurement of |∆m 2 32 | and the mixing angle θ 23 .
The paper is organised as follows. The scenario of invisible decay plus oscillations for three-generation mixing and oscillations in earth matter are discussed in Section II. The simulation of events and χ 2 analysis are explained in Section III. In Section IV we present our results for the sensitivity to the decay parameter τ 3 /m 3 . The effects of the presence of decay on the precision measurements of sin 2 θ 23 and |∆m 2 32 | are discussed in Sections V A and V B respectively. The exclusion contours are presented in Section V C. Conclusions are presented in Section VI.

II. INVISIBLE DECAY AND OSCILLATIONS IN THE PRESENCE OF MATTER
In this section we consider the oscillations and decay of ν 3 in the presence of matter. Let the state ν 3 decay invisibly via ν 3 → ν s + J, where J is a pseudo-scalar and ν s is a sterile neutrino. Since ν s does not mix with the three active neutrinos, the mixing matrix U in vacuum [69][70][71] is given by : where c ij = cos θ ij , s ij = sin θ ij ; θ ij are the mixing angles and δ is the CP violating phase. The mass of ν s is such that when the hierarchy is normal, m s < m 1 < m 2 < m 3 . Since ν s does not mix with the active neutrinos, the propagation equation is not affected by this. The effect of decay is included in the three-flavor evolution equation in the presence of earth matter as follows : where E is the neutrino energy, α 3 = m 3 /τ 3 is the decay constant in units of eV 2 , m 3 is the mass of ν 3 and τ 3 its rest frame life time. Since the term α 3 appears in the propagation equation along with ∆m 2 31 , it has to be in units of eV 2 . The conversion factor used here is 1 eV/s = 6.58 × 10 −16 eV 2 . The matter potential is where, G F is the Fermi constant and n e is the electron number density in matter and ρ is the matter density. For anti-neutrinos, both the sign of A cc and the phase δ in Eq. (2) are reversed.

A. Effect of the decay term
The decay term is of the form of exp (−αL/E). No decay corresponds to α = 0 and the exponential term as 1 whereas complete decay will be when the exponential term tends to 0. The effect of the decay parameter α for various L/E values can be understood from The ranges of exp (−αL/E) values for various values of α accessible for the specified range of L/E for a given baseline is shown in Table I. For a given L, a broader range of E will improve the sensitivity to α; on the other hand for a given E range the sensitivity to α will increase if longer baselines are available. In principle any experiment which spans over a wide range of L/E will have a better sensitivity to decay; with larger L/Es being sensitive to smaller values of α and vice versa. Atmospheric neutrino oscillation experiments fulfill this exact requirement. If we consider the neutrino energy range of 0.5-25 GeV, atmospheric neutrinos will span the L/E range of [0.6, 25484] (km/GeV) which includes all possible baselines from 15 km to the earth's diameter. The INO ICAL detector becomes relevant in this context. Since ICAL can detect neutrinos in the range 0.5-25 GeV [59] and since it is an atmospheric neutrino experiment, it will be sensitive to a wide range of α values. As seen from Fig. 1 ICAL should give a sensitivity to α = 10 −6 eV 2 also. The sensitivity to low α values come from the low energy part of the spectrum, while the higher energy parts of the spectrum will help us rule out larger values of α. First let us consider the effect of α 3 alone for a given θ 23 . The plots for α 3 = 0 correspond to the oscillation only case and as the value of α 3 increases the effect of decay becomes prominent which can be seen from the figure. In general the effect of decay is seen to be more for the lower energy neutrinos. For the decay constant α 3 = 10 −4 eV 2 , the effect of decay increases and the neutrino probabilities show significant depletion as compared to the no decay case for neutrino energies up to ∼ 15 GeV. The presence of decay reduces the oscillation amplitude near maxima and elevates it near minima. As α 3 increases to 10 −3 eV 2 , the survival probability of the neutrino and anti-neutrinos show a difference over the entire energy range considered. We also note that the effect of decay is mainly to damp out the oscillatory behavior in the probabilities. For the large decay case the oscillatory behavior is seen to be largely washed out. From Fig. 2 it can be seen that, the relative change in the oscillation probability due to decay is more forP µµ than P µµ whereas the relative change in P eµ is more compared to that inP eµ . Hence the contribution to the α 3 sensitivity χ 2 will be more from anti-neutrino events in the former case and neutrino events in the latter case.
However since P µµ andP µµ are the dominant channels at ICAL, the major contribution to α 3 sensitivity is expected to come from anti-neutrino events in the present study. Now let us look at the effect of θ 23 alone for a given α 3 value. The effect of θ 23 is also to vary the oscillation amplitude. In general, P µµ andP µµ decrease with increase in θ 23 . However beyond 13 GeV, for α 3 = 0 and 10 −4 eV 2 , θ 23 = 45 • gives the lowest probability compared to those for 39 • and 51 • , though the relative variation is much less. From the plots in the lower panels of Fig. 2 we see that P eµ andP eµ increase with θ 23 , the increase in P eµ is larger than that inP eµ . For all values of θ 23 , P eµ andP eµ decrease.
Since both α 3 and θ 23 affect the oscillation amplitudes, when combined in the following way, similar probabilities can be obtained. The combination of θ 23 in the first octant + a larger (smaller) value of α 3 will give a probability similar to that with θ 23 in second octant + a smaller (larger) value of α 3 for P µµ andP µµ (P eµ andP eµ ). Since the event spectrum is dominated by P µµ andP µµ events, this combined effect will affect the sensitivity/discovery potential to/of α 3 and the precision measurement on θ 23 ,which is discussed in Section V.

III. DETAILS OF NUMERICAL SIMULATIONS
ICAL will be a 50 kton magnetised iron detector which is optimised for the detection of atmospheric ν µ andν µ .Both ν µ (ν µ ) and ν e (ν e ) fluxes can contribute to the ν µ (ν µ ) events observed at ICAL. Hence the number of events detected by ICAL will be : where n d is the number of nucleon targets in the detector, σ µ is the differential neutrino interaction cross section in terms of the energy and direction of the muon produced, Φ µ and Φ e are the ν µ and ν e fluxes and P m αβ is the oscillation probability of ν α → ν β in matter and in presence of decay. A sample of 1000 years of unoscillated neutrino events are generated using NUANCE-3.5 neutrino generator [74], in which the Honda 3D atmospheric neutrino fluxes [73] along with neutrino-nucleus cross-sections and a simplified ICAL detector geometry are incorporated. Each event is oscillated by multiplying with the relevant oscillation probability including decay and oscillations in Earth matter assuming PREM density profile [72]. The probabilities are obtained by solving the propagation equation in matter in presence of decay. The events are then smeared according to the resolutions and efficiencies obtained from [45,46]. These two steps are done on an event by event basis for the entire 1000 year sample. Both "data" and theory are generated via this method, "data" with the central values of the parameters as described in Table II and theory by varying them in their respective 3σ ranges. Afterwards the oscillated samples of 1000 years of events, both "data" and theory are scaled down to the required number of years, 10 for our current analysis. This is done to reduce the effect of Monte-Carlo fluctuations on sensitivity studies.
In the current analysis, the efficiencies and resolutions of muons in the central region of the detector [45,46] have been used over the entire detector. These resolutions and efficiencies have been obtained by the INO collaboration via detailed detector simulations using a GEANT4-based simulation toolkit for ICAL. The central region of the ICAL detector [45,46] has the best efficiencies and resolutions for muons, the few-GeV muons in ICAL have a momentum resolution of ∼ 10% and direction resolution of ∼ 1 • on the average. Their relative charge identification efficiencies is about ∼ 99%. However, ICAL has two more regions namely the peripheral [47,48] and side regions depending on the magnitude and strength of the magnetic field. The peripheral region which has lesser reconstruction efficiencies but only slightly worse resolutions compared to the central region, constitutes 50% of the detector. Hence, in a realistic scenario where the efficiencies and resolutions in different regions are taken appropriately, the results obtained with 10 years of running of 50 kton of ICAL will only be obtained by increasing the run time to 11.3 years, as mentioned in [59].
Since the charged current ν µ (ν µ ) interactions have µ − (µ + ) in the final state along with the hadron shower, and since ICAL is capable of measuring the energy of the hadron shower, we include in our analysis the data on those as well. It was reported in [49] from ICAL simulations that hadrons in ICAL have energy resolutions of 85% at 1 GeV and 36% at 15 GeV and the events are smeared accordingly before including them in the final 3D-binned analysis which includes muons binned in observed energy and direction and hadrons binned in energy. There are 15 bins in E obs µ between (0.5 − 25) GeV, 21 bins in cos θ obs µ between (−1, +1)and 4 bins in E ′obs had between (0 − 15) GeV, thus giving 1260 bins. More details of the binning scheme and the numerical simulations can be found in Ref. [59].
The true values and the 3σ ranges of the oscillation parameters used to generate the probabilities are given in Table II. Since ICAL is not directly sensitive to δ CP , it is taken as 0 • in this analysis and kept fixed. The 1-2 oscillation parameters ∆m 2 21 and sin 2 θ 12 are also kept fixed throughout our analysis. For the remaining parameters two types of analyses are performed, -one with fixed parameter and the other with marginalisation. In the former all parameters are kept fixed while in the latter, the parameters other than the one for which the sensitivity study is done are marginalised in their respective 3σ ranges shown in Table II.

Parameter
True value Marginalization range θ 13 8.  To statistically analyse the data, we define the following χ 2 function Here i, j, k sum over muon energy, muon angle and hadron energy bins respectively. The number of predicted (theory) events with systematic errors in each bin are given by The number of theory events without systematic errors in a bin is given by T 0± ijk and the observed events ("data") per bin are given by D ± ijk . It should be noted that both D ± ijk and T 0± ijk are obtained from the scaled NUANCE neutrino events as mentioned earlier.
The following values are taken for the systematic uncertainties [75,76]: π 1 = 20% flux normalisation error, π 2 = 10% cross section error, π 3 = 5% tilt error, π 4 = 5% zenith angle error, π 5 = 5% overall systematics and π 6 = 2.5% on Φ νµ /Φν µ ratio. These are included in the analysis via pull method. The "tilt" error is incorporated as follows. The event spectrum with the predicted values of atmospheric neutrino fluxes is calculated and then shifted according to the relation : where E 0 is 2 GeV, and δ is the 1σ systematic tilt error (5%). Flux error is included as the difference Φ δ (E) − Φ 0 (E). A prior of 8% at 1σ is added to sin 2 2θ 13 . This is the only prior in this calculation. No prior is imposed at all on the quantities whose sensitivities are to be studied, i.e on α 3 , θ 23 and |∆m 2 32 |. The contribution from prior to the χ 2 is : where, σ(sin 2 2θ 13 ) = 0.08 × sin 2 2θ true 13 . Hence, the final χ 2 for ICAL will be : where χ 2 is given by Eq. (6).

IV. SENSITIVITY OF ICAL TO α 3
The results of the sensitivity studies of ICAL to α 3 are presented in this section. We first show how the number of oscillated events change with decay as a function of zenith angle and energy. Then we proceed further to discuss the sensitivity as well as the discovery potential of ICAL to neutrino decay and the bound on α 3 from our analysis.
a. Effect of decay on the number of oscillated events: In Fig. 3, we show the zenith angle distribution of the ν µ andν µ events for different values of the decay constant α 3 . The four panels are for four different energy bins. The convention used in these plots is such that cos θ obs µ = [0, 1] indicates the up-coming neutrinos. It can be seen from the figure that both ν µ andν µ events deplete with an increase in the value of α 3 . We also note that the effect of decay is more prominent in the lower energy bins. With increase in energy, there is no significant effect of decay on the number of events if the decay parameter is less than 10 −4 eV 2 as can be seen from the lower panels.    Table II. It should be noted that the y-axes are not the same. Only up-coming events (oscillated) are shown here.
b. Sensitivity to the decay parameter α 3 : In this section, first the study of the sensitivity of ICAL to α 3 is presented with 500 kton-yr exposure of the detector taking normal hierarchy (NH) as the true hierarchy. To that end, we simulate the prospective "data" for no decay and fit it with a theory of oscillation plus decay. The corresponding χ 2 is shown as a function of α 3 (test) in the left panel of Fig. 4.  The blue dashed curve is obtained for a fixed parameter fit while the blue solid one corresponds to the sensitivity when the χ 2 is marginalised over all oscillation parameters as described in Section III. A comparison of the solid and dashed curves gives us an idea of the impact of marginalisation over the oscillation parameters on the sensitivity of the experiment to decay. From Fig. 4 it can be seen that with marginalisation of the oscillation parameters, the sensitivity decreases as expected. The right panel shows the sensitivity to decay in terms of τ 3 /m 3 in s/eV. The expected sensitivity of ICAL to α 3 are shown in Table III. The corresponding values of τ 3 /m 3 in units of s/eV are also given. Note that by sensitivity limit we mean the value of α 3 (τ 3 /m 3 ) upto which ICAL can rule out neutrino decay.  The lower bound on τ 3 /m 3 for the invisible decay scenario from MINOS data was shown to be τ 3 /m 3 > 2.8 × 10 −12 (s/eV) at 90% C.L. This corresponds to an upper limit α 3 < 2.35 × 10 −4 eV 2 . Table III shows that ICAL is expected to tighten these bounds by two orders of magnitude with just charged current ν µ andν µ events. At 90% C.L, ICAL with marginalisation is expected to give a lower bound of τ 3 /m 3 > 1.51 × 10 −10 (s/eV) which corresponds to α 3 < 4.36 × 10 −6 eV 2 .

) (test)
The expected sensitivity with fixed parameters as well as marginalisation for true IH are shown in Fig. 5. At 90% C.L, the upper bound on α 3 are α 3 < 2.78 × 10 −6 eV 2 with fixed parameters and α 3 < 5.82 × 10 −6 eV 2 with marginalisation. These are only slightly worse than the sensitivities obtained with true NH. In terms of τ 3 /m 3 , these limits translate as the lower limits τ 3 /m 3 > 2.42 × 10 −10 s/eV and τ 3 /m 3 > 1.14 × 10 −10 s/eV for the fixed parameter and marginalised cases, respectively. The expected sensitivity to α 3 at different C.L. with true IH is summarised in Table IV.  The analysis discussed above gives us the sensitivity to α 3 when we fit a "data" with no decay with a theory which has decay. On the other hand, if neutrinos indeed decay into sterile components, and if the decay rate is large enough to be observed in ICAL, we will be able to discover neutrino decay at this experiment. Therefore, we next estimate how much the decay rate needs to be in order for ICAL to make this discovery. For this analysis, we simulate the "data" with different values of α 3 and fit it with a theory with no decay. The analysis was done for 500 kton-yr exposure of ICAL for fixed parameters as well as with marginalisation of the undisplayed parameters over their respective 3σ ranges. The results are shown in Fig. 6 by the red-dashed curve for the fixed parameter case and the red-solid line for the marginalized case. However, we find that for the discovery potential, the marginalization has no effect and gives the same result as the fixed parameter case. We find that ICAL will be able to discover neutrino decay at the 90% C.L. if α 3 > 2.5 × 10 −6 eV 2 . We also plot the sensitivity curves, blue dashed (solid) lines for the fixed parameter (marginalized) case, in this figure for a comparison between the 'sensitivity" and "discovery" potential of α 3 . We can see that the "sensitivity" and "discovery" limits of ICAL are very similar for fixed parameter analysis. However for the marginalised case the "discovery potential" is significantly higher than the "sensitivity" limit and is same as the fixed parameter case. The reason why the expected "sensitivity" limit worsens due to marginalisation while the expected "discovery" limit does not can be understood as follows. For the "sensitivity" analysis we generate the data for no decay and θ 23 maximal and fit it with a theory where α 3 = 0. Since the effect of decay is to reduce the number of events and suppress the event spectrum for fixed parameter there will be a difference between the data and the theory giving a higher χ 2 . For the marginalized case, this can be compensated to some extent by suitably changing the value of θ 23 from maximal and thereby reducing sin 2 2θ 23 , the leading term that controls the amplitude of oscillations in the case of muon neutrino survival probability. This can be seen in Fig. 7. In this figure the solid line denotes the "data" generated with α 3 = 0 i.e no decay and θ 23 = 45 • while the dashed (dotted) lines show the theory events for a non-zero α 3 and θ 23 = 45 • (38.65 • ). We can see that the lower value of θ 23 compensates for the depletion due to decay and can give a lower χ 2 . As a result the expected sensitivity drops when the sensitivity χ 2 is marginalised over θ 23 .     GeV (left) with α 3 = 0 eV 2 in "data" and α 3 = 1.316 × 10 −5 eV 2 in theory; (right) with α 3 = 1.316 × 10 −5 eV 2 in "data" and α 3 = 0 eV 2 in theory for the marginalised case. "D" represents data and "T" represents theory events. The blue histograms are for ν µ and the red ones are forν µ events. The theory events are generated with marginalisation of parameters except α 3 in their respective 3σ ranges.
On the other hand, for the expected "discovery" limit case we generate the data for nonzero α 3 and maximal mixing and fit it with a theory with no decay. In this case, the data has events lower than the theory due to decay. This can be seen from the second panel of Fig. 7 where the blue (red) solid line denotes the data events for muon neutrinos (antineutrinos). However, unlike the "sensitivity" case, here one cannot change θ 23 to reduce the event spectrum any further to compensate for the difference between data and theory since maximal mixing already corresponds to maximal suppression of the muon neutrino survival probability, the leading oscillation channel for atmospheric neutrinos. As a result, the fit continues to keep θ 23 at its maximal value and marginalisation fails to lower the χ 2 any further. This can also be seen from Fig. 7 where the dotted line shows the theory events obtained after marginalizationand this is higher than the data events and same as the fixed parameter case. We next look at the impact of neutrino decay on the precision measurement of the mixing angle θ 23 and the mass squared difference |∆m 2 32 | at ICAL. A comparison of the precision measurement in the presence and absence of decay is presented. In the no decay case both "data" and theory are generated without the decay parameter and in the case with decay both "data" and theory are generated with non-zero values of α 3 . For all results presented in this section, the value α 3 = 1 × 10 −5 eV 2 is used to generate the "data". In the fixed parameter analysis this is kept fixed in theory and for the marginalised case, the range over which α 3 is marginalised is taken to be α 3 = [0, 2.35 × 10 −4 ] eV 2 which corresponds to the 90% CL bound given by the MINOS analysis. The other parameters are kept fixed at their true values as shown in Table II for the fixed parameter analyses and varied in the 3σ ranges as shown in the same table for the marginalized case. The 1σ precision on a parameter λ is defined as : where λ max-2σ and λ min-2σ are the maximum and minimum allowed values of λ at 2σ and λ true is the true choice.
A. Precision on sin 2 θ 23 in the presence of decay : The sensitivity to sin 2 θ 23 in the presence and absence of ν 3 decay is shown in Fig. 8. The left panel shows the fixed parameter results whereas the right panel shows the results for the marginalised case. For the fixed parameter case, in the absence of decay, the 1σ precision on sin 2 θ 23 is ∼ 8.9%. In presence of decay the 1σ precision is ∼8.6% which is similar to the no decay case. However, it is important to note that even though the percentage precision is same, the allowed parameter space is shifted to the right when there is decay, as compared to the no decay case. The minimum and maximum values of sin 2 θ 23 at 2σs in the presence and absence of decay are shown in Table V.  In order to understand the shift of parameter space, we show in Fig. 9, the number of oscillated decay there are differences between the number of events for various θ 23 values. For the case of no decay this difference is less as compared to the case where decay is present. Comparing the figures on the left and right panels one also observes that, the difference between the number of events for θ 23 = 39 • and 45 • is more in presence of decay and the curve for θ 23 = 45 • is closer to 52 • . Now, in obtaining the precision plot the data is generated with true θ 23 = 45 • and in theory the θ 23 is kept fixed. For θ 23 in the lower octant the difference of the number of events with that for 45 • being more in presence of decay, the χ 2 for a θ 23 in the lower octant will be higher as compared to the no decay case. On the other hand, for θ 23 in the higher octant the difference in the number of events with θ 23 = 45 • being less in presence of decay, one gets a lower χ 2 as compared to the no decay case. This explains why the precision curve shifts towards higher θ 23 values.
For the marginalised case, in the presence of decay the overall precision becomes worse compared to the no decay case. The 1σ precision when decay is present is ∼10.85% whereas for no decay it is ∼8.9%. This can be explained as follows. In the marginalised case, for only oscillation we are trying to fit the "data" generated with θ 23 = 45 • , varying the other parameters in theory. In this case the θ 13 can be adjusted to give a slightly lower χ 2 . In presence of decay we generate the "data" for a particular non-zero α 3 and θ 23 = 45 • . But now in theory we vary α 3 as well as the other parameters. For θ 23 in the lower octant, the theory events will be higher than the "data" events as can be seen by comparing the events in the second panel of Fig. 9. However, in this case the α 3 can be increased to give a better fit and a lower χ 2 . On the other hand for θ 23 in the higher octant, the data events are higher than the theory events and α 3 can be decreased in theory to match the data better and give a lower χ 2 . This explains the widening of the χ 2 vs θ 23 curve in presence of decay. Note that this is more for the higher octant because the difference of the events for θ 23 = 45 • and say 52 • is less as compared to θ 23 in the lower octant, say 39 • . This gives a lower χ 2 thus allowing more θ 23 values in the higher octant.    : The precision on the magnitude of the mass square difference |∆m 2 32 | in the presence and absence of invisible decay of ν 3 is presented in Fig. 10. NH is taken as the true hierarchy. The relative 1σ precision on |∆m 2 32 | with oscillations only and with decay is ∼2.5% for the fixed parameter case. When marginalisation is done this becomes ∼2.6% for both the cases. Thus it can be seen that the presence of decay does not affect the precision on |∆m 2 32 | much. This is because decay mainly affects the amplitude of the oscillations and not the phase which is determined by |∆m 2 32 |. The minimum and maximum values of sin 2 θ 23 at 2σs in the presence and absence of decay are shown in Table VI.

VI. SUMMARY AND DISCUSSIONS
The expected sensitivity of ICAL to the decay lifetime of the mass eigenstate ν 3 , when it decays via the invisible decay mode was presented. The analysis was performed in the threegeneration neutrino oscillation framework including decay as well as earth matter effects. The decay was parameterised in terms of α 3 = m 3 /τ 3 , where, m 3 is the mass and τ 3 the lifetime at rest of the mass eigenstate ν 3 . With 500 kton-yr of exposure, ICAL is expected to constrain the invisible decay rate to α 3 < 4.36 × 10 −6 eV 2 at 90% C.L., which is two orders of magnitude tighter than the bound obtained in [37] for MINOS. In [37] both charged current (CC) and neutral current (NC) events were considered where as in our study only atmospheric CC ν µ andν µ events were used. For invisible neutrino decay, the NC background will be less. Hence the sensitivity to α 3 is expected to improve.
The effect of decay on the 2-3 oscillation parameters was also studied. Since the amplitude of oscillations is affected most by the presence of decay, it was found that decay affected the precision measurement of sin 2 θ 23 . For 500 kton-yrs of exposure assuming NH as the true hierarchy, the 1σ precision on sin 2 θ 23 was found to worsen to 10.85% when α 3 = 10 −5 eV 2 was assumed. This is worse as compared to the 8.87% obtained with oscillation only hypothesis. In the case of |∆m 2 32 | the 1σ precision without decay is 2.5% whereas the inclusion of invisible decay does not affect it at all. The effect of α 3 on the sensitivity to neutrino mass hierarchy and octant of θ 23 will be studied elsewhere [77].
It is also noteworthy that the sensitivity to smaller α 3 comes mainly from the lower energy bins below 2 GeV. Hence, if we can improve the efficiencies and resolutions of the detector, especially for muons in the lower energy region, we will be able to put a better limit on α 3 . Reduction of the energy threshold for the detection of low energy neutrinos in future will also help probing phenomena like decay with increased precision. This is important since the atmospheric neutrino flux peaks at lower energies and by being able to detect and analyse more events we will further improve our sensitivities to all parameters including α 3 .