Invisible neutrino decay in the light of NOvA and T2K data

We probe for evidence of invisible neutrino decay in the latest NOvA and T2K data. It is seen that both NOvA and T2K data sets are better fitted when one allows for invisible neutrino decay. We consider a scenario where only the third neutrino mass eigenstate $\nu_3$ is unstable and decays into invisible components. The best-fit value for the $\nu_3$ lifetime is obtained as $\tau_{3}/m_{3} = 3.16\times 10^{-12}$ s/eV from the analysis of the NOvA neutrino data and $\tau_{3}/m_{3} = 1.0\times 10^{-11}$ s/eV from the analysis of the T2K neutrino and anti-neutrino data. The combined analysis of NOvA and T2K gives $\tau_{3}/m_{3} = 5.01\times 10^{-12}$ s/eV as the best-fit lifetime. However, the statistical significance for this preference is weak with the no-decay hypothesis still allowed at close to 1.5$\sigma$ C.L. from the combined data sets, while the two experiment individually are consistent with no-decay even at the 1$\sigma$ C.L. At 3$\sigma$ C.L., the NOvA and T2K data give a lower limit on the neutrino lifetime of $\tau_{3}/m_{3}$ is $\tau_{3}/m_{3} \geq 7.22 \times 10^{-13}$ s/eV and $\tau_{3}/m_{3} \geq 1.41 \times 10^{-12}$ s/eV, respectively, while NOvA and T2K combined constrain $\tau_{3}/m_{3} \geq 1.50 \times 10^{-12}$ s/eV. We also show that in presence of decay the best-fit value in the $\sin^{2}\theta_{23}$ vs $\Delta m^{2}_{32}$ plane changes significantly and the allowed regions increase significantly towards higher $\sin^{2}\theta_{23}$.


INTRODUCTION
The neutrino oscillation physics has entered its precision age. Among the six oscillation parameters associated with the standard neutrino oscillation physics, ∆m 2 21 and θ 12 have been measured precisely from the solar neutrino experiments [1] and KamLAND [2]. The atmospheric neutrino experiments [3] first showed evidence of ν µ − ν τ flavor transformation and gave a measurement of |∆m 2 32 | and θ 23 . The mixing angle θ 13 is also now precisely determined from reactor experiments DayaBay [4], Double CHOOZ [5] and RENO [6]. Currently the unknown quantities in the neutrino oscillation physics are the neutrino mass hierarchy/ordering i.e., whether the lightest neutrino state is ν 1 or ν 3 , the precise measurement of θ 23 and its octant and the CP-violating phase δ CP . NOvA [7,8] and T2K [9] are the presently running long baseline experiments. These are complimentary to the previous experiments and are expected to shed light on the unknown neutrino parameters.
The NOvA long-baseline experiment is in USA. The ν µ flux is generated by the NuMI beam at Fermilab. The experiment employs two identical Totally Active Scintillator Detector (TASD) which are different only in terms of their mass. The Near detector (ND) is a 100 m deep 290 ton detector while the Far Detector (FD) is on the surface with mass 14 kton. The neutrino flux is measured first at the ND located at 1 km away from the target. The * Electronic address: sandhya@hri.res.in † Electronic address: debajyotidutta@hri.res.in ‡ Electronic address: dipyamanpramanik@hri.res.in neutrino beam is next detected at the FD situated 810 km away near Ash River, Minnesota at an off-axis angle of 14.6 mrad. The 14.6 mrad off-axis location of the FD gives a narrow energy spectrum at the FD with a peak near 2 GeV, which is tuned close to the first oscillation maximum at this baseline.
T2K (Tokai to Kamioka) is a long-baseline experiment in Japan. Here the ν µ beam is produced at the J-PARC accelerator complex in Japan by impinging a 30 GeV proton beam onto a carbon target. This experiment also employs a two detector set-up that are off-axis compared to the beam direction. The neutrino produced from the beam are detected first at the ND, ND280, at 280 m from the target. The FD is the Super-Kamiokande detected with fiducial mass 22.5 kton and situated 295 km away from the source at 2.5 • off-axis from the main beam axis, giving a narrow beam peaked at around 600 MeV, which for T2K's baseline corresponds to the first oscillation maximum.
NOvA and T2K have both presented their initial results. The T2K experiment announced their first result with 1.43 × 10 20 protons on target (POT) on electron appearance in 2011 [10]. With six observed electron candidate events and 1.5 expected backgrounds, T2K gave the first direct evidence of non-zero θ 13 at 2.5σ C.L. The first results announcement on muon disappearance came a year later in 2012 [11] with the same POT of 1.43×10 20 and gave a best-fit of ∆m 2 32 = 2.65 × 10 −3 eV 2 and sin 2 2θ 23 = 0.98. The first anti-neutrino result from T2K was published in [12], where they usedν µ beam with 4.01 × 10 20 POT and obtained best fit of sin 2 θ 23 = 0.45 and ∆m 2 32 = 2.51 × 10 −3 eV 2 , very consistent with the measurements obtained using the ν µ beam. The T2K collaboration has published their results periodically since the first announcements and their data-sets and best-fit oscillation parameters have remained consistent, with θ 23 close to maximal. The latest result form T2K [13] used 7.482 × 10 20 POT for neutrino mode and 7.471 × 10 20 POT for anti-neutrino data. This currently gives the best fit ∆m 2 32 = 2.52 ± 0.08(2.51 ± 0.08) × 10 −3 eV 2 and sin 2 θ 23 = 0.55 +0.05 −0.09 (0.55 +0.05 0.08 ) for normal(inverted) ordering, using both electron appearance as well as muon disappearance data. The first result of muon-neutrino disappearance from NOvA came in 2016 [14], where they used 2.74 × 10 20 POT and got the best fit ∆m 2 32 = (2.52 +0. 20 −0.18 ) × 10 −3 eV 2 and sin 2 θ 23 = 0.43 and 0.60. The was immediately followed with first result on electron appearance data [7] with the same exposure. Next disappearance data came in 2017 [15] which used 6.05×10 20 POT and gave ∆m 2 32 = (2.67±0.11)×10 −3 eV 2 , while for sin 2 θ 23 they obtained two statistically degenerate values 0.404 +0.030 −0.022 and 0.624 +0.022 −0.030 and claimed that the NOvA data disfavours maximal mixing at 2.6σ. Results from the combined analysis of NOvA's appearance and disappearance data was presented in [8]. So while T2K prefers maximal mixing for θ 23 , the early analysis from the NOvA showed 2.6σ preference for non-maximal mixing. This tension between the two experiments led several authors to propose new physics ideas to explain the tension between the two datasets. However, the NOvA collaboration has recently done an improved re-analysis of their disappearance dataset [16]. The newer analysis mainly addresses better the energy resolution of the hadron sample leading to an improved neutrino energy resolution. They have divided the muon events into four quantiles of different resolutions from ∼ 6% to ∼ 12%, based on their hadronic energy fraction. This approach changed the measurement of θ 23 at NOvA with maximal θ 23 mixing being preferred by NOvA as well. Hence, the tension between NOvA and T2K has been resolved for now and the datasets seem to be consistent with standard three-generation flavor oscillations.
Although the data appears to be consistent with the standard expected three-generation paradigm, there can still be new physics effects present in these experiments. One such new physics scenario is neutrino decay. The active neutrino state could decay into another lighter active neutrino state and boson(s), or it could decay into a sterile fermion and boson(s). The former scenario is called visible neutrino decay [17][18][19].(since the final state fermion is active and hence "visible" to the detector) while the latter is known as invisible neutrino decay (since the final state fermion is sterile and hence "invisible" to the detector). The bosonic state(s) are assumed to be invisible in both class of models. In this paper we will consider neutrino decay into all invisible states only. There are two possibilities for such models: (i) The neutrino could decay into a Dirac fermion [20,21] ν j →ν iR + χ, whereν iR is a right-handed singlet and χ is an iso-singlet scalar carrying a lepton number. (ii) The neutrinos could decay into a Majorana fermion and a Majoron ν j → ν s +J [22,23]. To evade the constraints of the Z decay to invisible particles from LEP data, the Majoron should be dominantly singlet [24]. First idea of decaying neutrinos was proposed very early in order to explain the solar neutrino problem [25]. Later neutrino oscillations along with decay solution to the solar neutrino problem was studied in [21,[26][27][28][29][30][31][32]. These studies obtained bound on τ 2 by considering ν 2 to be unstable. The bound on τ 2 from the solar data is τ 2 /m 2 > 8.7 × 10 −7 s/eV at 99 % C.L. [31]. (See [32,33] for a more recent study). The bound on τ 2 from SN1987A are much more stringent [34]. Atmospheric and long-baseline (LBL) neutrinos give the bound on the ν 3 lifetime. To solve the atmospheric neutrino problem a pure neutrino decay was proposed in [35], however this gave a very poor fit to the data. Authors of [36,37] considered neutrino decay with mixing and claimed that it could somewhat reproduce the SK results, however, the zenith angle dependent SK data gave a poor fit [38]. In [36,38] the ∆m 2 dependent terms were averaged out. However, if the unstable neutrino state is allowed to decay into a sterile state with which it does not mix then the constraints on ∆m 2 can be relaxed. Two scenarios have been studied in this context. In [39] ∆m 2 << 10 −4 eV 2 was considered and it was claimed to fit the data better, however the analysis by the SK collaboration showed that this gives poorer fit than the only oscillation case [40]. The second scenario was considered in [41], where ∆m 2 was kept unconstrained. For SK data, it was shown that a small non-zero decay parameter and ∆m 2 ∼ 0.003 eV 2 gave a better fit to the data. The global analysis of atmospheric and MINOS data was performed in [42]. Although only oscillation gave the best-fit, but the decay plus oscillation scenario also gave a good fit. The bound obtained from the analysis is τ 3 /m 3 ≥ 2.9 × 10 −10 s/eV at the 90 % C.L. Prospects of constraining the ν 3 lifetime with atmospheric neutrino events at INO was studied in [43].
The case for visible neutrino decay in long-baseline experiments was considered for T2K in [44] and DUNE in [45]. However, the visible neutrino decay scenario is very tightly constrained from data from other experiments. Therefore, we consider in this paper only the case for invisible neutrino decay. The analysis [46] considered MI-NOS and T2K data with two generation of neutrinos and gave the bound, τ 3 /m 3 ≥ 2.8 × 10 −12 s/eV at 90 % C.L. The expected results at DUNE in the invisible neutrino decay scenario was worked out in [47].
In this paper we present the current constraints on τ 3 /m 3 from the recent data from NOvA and T2K, using the full three generation oscillation framework with matter effects. We also study the effect of decay on the measurement of θ 23 and ∆m 2 32 . As was shown in [47], there can be significant impact on the measured value of θ 23 if decay is present. We will study how presence of decay changes the best-fit values of θ 23 and ∆m 2 32 as well as the C.L. contours allowed by the T2K and NOvA data.
The paper is organised as follows. In the next section we give the theoretical description of neutrino propaga-tion when the ν 3 state decays into a sterile state with which it does not mix. Also given in section 2 are the details of our simulation framework needed for the real data analysis of T2K and NOvA. We next give our result in section 3 and finally we conclude in section 4.

INVISIBLE NEUTRINO DECAY AND SIMULATION FRAMEWORK
We assume that the ν 3 state decays into a sterile neutrino and a singlet scalar (ν 3 →ν 4 + J).We also assume that the neutrino mass eigenstate decays, therefore the mass matrix as well as the decay matrix can be simultaneously diagonalised. In this case the mixing between the flavor and the mass eigenstates can be written as, where the greek indices represent the standard flavor states i.e., e, µ, τ and the latin indices represent the mass eigenstates. U is the standard PMNS matrix. Here we assume NH i.e., m 3 > m 2 > m 1 . In presence of decay the evolution equation in matter becomes: where A = 2 √ 2G F n e E represents the matter potential due to neutrino electron scattering in matter, G F is the Fermi coupling, E is the neutrino energy and n e is the electron density. We solve Eq. (2) numerically using Runge-Kutta with constant matter density.
The simulation is done using a modified version of GLoBES, with modifications needed for real data analysis.
For the analysis of NOvA we have taken a 14 kt detector at a baseline of 812 km with constant matter density. We have taken 8.5% energy resolution for electron events and 6% resolution for muon events. The signal efficiency is chosen to be 99 % for electron events and 91 % for muon events. We normalize the number of events to match the best fit event spectra given in [8] for electrons (6.04×10 20 POT) and in [16] for muons (8.85 × 10 20 POT).
For the analysis of T2K, we have taken a 22.5 kt detector at a baseline of 295 km with constant matter density. The energy resolution is taken 8.5 %. The signal efficiency is chosen to be 51.5 % for electron events and 90 % for muon events. We normalise the event spectra to match the event spectra given in [13] which corresponds to 7.482 × 10 20 POT in neutrino mode and 7.471 × 10 20 POT for anti-neutrino mode.

RESULTS
In Fig. 1 we show the constraint on τ 3 /m 3 from the current data of NOvA and T2K, where the latest appearance as well as disappearance data sets of both experiments have been taken into consideration in the analysis. The green solid curve is obtained using NOvA data alone, the dark red dotted curve is obtained using T2K data alone, while the black dashed curve is obtained from 3 σ a combined analysis of NOvA and T2K data. The parameters θ 23 and ∆m 2 32 are marginalised in the fit over their 3σ allowed ranges and δ CP is marginalised over its full range. It can be seen from the Fig. 1, that for both experiments the no decay scenario is slightly disfavoured and the best-fit value from the fit comes for non-zero decay. NOvA disfavors the no decay scenario at 0.7σ and the best-fit value is τ 3 /m 3 = 3.16 × 10 −12 s/eV. T2K disfavors the no decay case at slightly more than 1σ and the best-fit value is 1 × 10 −11 s/eV. For the  combined analysis case, the no decay case is disfavoured at 1.5 σ and the best-fit obtained from the combination of the two data is τ 3 /m 3 = 5.01 × 10 −12 s/eV. The minimum χ 2 (χ 2 min ) for NOvA, T2K and the combined case are 10.38, 69.34 and 87.19, respectively, which are slightly less than the standard oscillation fit, for which the χ 2 min are 10.93, 70.39 and 88.65, respectively. Therefore, the invisible decay scenario we consider in this work, fits the data slightly better than the standard oscillation case. The data sets also set a lower bound on the lifetime. The 3σ lower bound on τ 3 /m 3 from NOvA data is seen to be τ 3 /m 3 ≥ 7.22 × 10 −13 s/eV, while from T2K it is τ 3 /m 3 ≥ 1.41 × 10 −12 s/eV. The 3σ combined constraint from both experiments taken together is τ 3 /m 3 ≥ 1.50 × 10 −12 s/eV. It can be clearly seen from Fig. 1, that the sum of the two ∆χ 2 is less than the ∆χ 2 for the combined analysis. This points at a synergy between the two experiments. This synergy results in an improved fit for the decay scenario compared to standard oscillation when we perform the combined analysis of the two experiments. Fig. 2 shows the allowed region in the sin 2 θ 23 vs ∆m 2 32 plane at 95 % C.L. The solid curves are for only oscillation case without decay whereas the dashed curves are for the case where ν 3 are allowed to decay. The fit is marginalised over δ CP in both cases and over τ 3 /m 3 as well for the case of decay plus oscillation. The blue curves are for NOvA, the black curves are for T2K and the red curves are for the combined analysis. For the standard case the best-fit points (sin 2 θ 23 , ∆m 2 32 ) are (0.45, 2.41 × 10 −3 eV 2 ), (0.52, 2.56 × 10 −3 eV 2 ) and (0.46, 2.51 × 10 −3 eV 2 ) for NOvA, T2K and the combined cases, respectively. On the other hand, for the case with decay and oscillation the corresponding bestfit points are: (0.48,2.39 × 10 −3 eV 2 ), (0.62, 2.62 × 10 −3 eV 2 ) and (0.48, 2.52 × 10 −3 eV 2 ) for NOvA, T2K and the combined cases, respectively. The interesting point to notice here is that for all cases, the allowed region of the parameter space increases significantly when decay is considered along with oscillation. Also note that with inclusion of decay, the best-fit shifts towards higher values of sin 2 θ 23 . This behaviour is very similar to what was seen in [46] in the context of MINOS and T2K and in [47] in the context of DUNE. The shift of θ 23 to higher values in presence of decay can be understood in terms of the survival probability given in [46] in the two-generation approximation neglecting matter effect, In Eq. (3), there is an exponential suppression due to neutrino decay in both the oscillatory as well as the nonoscillatory term. Therefore for a given θ 23 , the survival probability for the decay case will be less than the standard oscillation case. Hence, when decay is considered in the fit, the value of sin 2 θ 23 increases in order to reproduce the same probability obtained for the standard case. Fig. 3 gives the allowed region in the sin 2 θ 23 vs τ 3 /m 3 plane at 95 % C.L. with ∆m 2 32 marginalised over its current 3σ range and δ CP marginalised over full range. The blue solid, black dashed and red dashed-dotted curves are for NOvA, T2K and the combined analysis, respectively. The blue cross gives the best-fit for NOvA (0.47, 3.16 × 10 −12 s/eV), the black plus gives the best-fit for T2K (0.61, 5.011 × 10 −12 s/eV) and the red star gives the best-fit for the combined case (0.49, 5.011 × 10 −12 s/eV). Again as in Fig. 2, the best-fit is seen to be for finite τ 3 /m 3 . Fig. 4 gives the ∆χ 2 vs sin 2 θ 23 for the standard case and the decay plus oscillation case. The fit is marginalised over ∆m 2 32 and δ CP for the standard case and over ∆m 2 32 , τ 3 /m 3 and δ CP for the decay plus oscillation case. The left panel is for NOvA, the middle panel is for T2K and the right panel is for the combined analysis. In all the panels the blue solid curves are for the standard oscillation case while the red dashed curves represent the decay with oscillation case. For T2K our standard oscillation best-fit sin 2 θ 23 = 0.52 matches very well with the best-fit obtained by the T2K collaboration [13]. For NOvA on the other hand, our best-fit for standard oscillation comes at sin 2 θ 23 = 0.45 while the NOvA in the higher octant. The reason for this mild mis-match could be because our experimental simulation is based on GLoBES which cannot include all systematics in a rigorous manner. However, the figure shows that the χ 2 differences between the minima of the two octant is less than 1, and hence the best-fit obtained by us for the higher octant is not so different from the correct value obtained from the collaboration. Fig. 5 gives the allowed region in the ∆m 2 32 vs τ 3 /m 3 plane at 95 % C.L. with θ 23 marginalised over current 3σ allowed region and δ CP marginalised over full range. The blue curve is for NOvA, the black is for T2K and the red is for the combined case. The blue cross gives the best-fit for NOvA (2.36 × 10 −3 , 3.16 × 10 −12 ), the black plus for the T2K (2.61 × 10 −3 , 5.011 × 10 −12 ) and the red star is for the combined best-fit (2.51 × 10 −3 , 5.011 × 10 −12 ). The above best-fit values are in units of (eV 2 , s/eV). Fig. 6 gives the ∆χ 2 vs ∆m 2 32 with θ 23 and δ CP marginalised for the standard and θ 23 , δ CP and τ 3 /m 3 marginalised for the case of decay plus oscillation. Just as in Fig. 4, the left panel is for NOvA, the middle is for T2K and the right panel is for the combined case. The blue solid curves are for the standard case whereas the red dashed curves are for the decay plus oscillation case. It can be seen from the figure that for NOvA the allowed range of ∆m 2 32 increases on both sides while for T2K the best-fit ∆m 2 32 shifts towards higher values and as a result the allowed range shifts towards the right.

SUMMARY & CONCLUSION
We analysed the recent NOvA and T2K data for invisible decay of the neutrino. We considered a framework where the ν 3 is unstable and it decays into some lighter sterile state with which the active neutrinos do not mix. We used modified GLoBES for simulating the T2K and NOvA experiments taking into account the experimental exposure, systematic uncertainties and resolution functions for the current data sets. We considered the latest disappearance and appearance data given by the T2K collaboration which corresponds to 7.482 × 10 20 POT for neutrino and 7.471×10 20 POT for anti-neutrinos [13]. For NOvA we consider the disappearance data and simulation framework announced by the collaboration in [16] and appearance data presented in [8]. To obtained the oscillation probabilities we considered the full three generation framework with matter effects and decay and solved the evolution equation for the neutrinos numerically. We performed a two-pronged analysis. On one hand we looked at how compatible the current T2K and NOvA data are with neutrino decay. On the other hand we checked how the best-fit values and allowed ranges on ∆m 2 32 and sin 2 θ 23 change when decay is included in the fit.
We found that both T2K and NOvA data give better fit when neutrino decay in included. The best-fit lifetime corresponding to the T2K data is τ 3 /m 3 = 1.0 × 10 −11 s/eV. This can be compared with the best-fit lifetime of τ 3 /m 3 = 1.6 × 10 −12 s/eV obtained in [46] using the older T2K data [48] which corresponded to 6.57 × 10 20 POT for the neutrino mode and did not have the antineutrino data. The main difference between the analysis performed in [46] and this work is the use of the anti-neutrino data, more exposure in the neutrino data and the inclusion of the electron appearance data. Our best-fit can also be compared to the best-fit lifetime of τ 3 /m 3 = 1.2 × 10 −12 s/eV obtained in [46] from the combined fit of T2K and MINOS data. The best-fit ν 3 lifetime from NOvA data corresponds to τ 3 /m 3 = 3.16 × 10 −12 s/eV, which is about an order of magnitude smaller than the T2K best-fit. The combined fit of T2K and NOvA returns a best-fit τ 3 /m 3 = 5.01 × 10 −12 s/eV. The datasets also put a lower bound on the ν 3 lifetime. The 3σ bound put by T2K, NOvA and T2K+NOvA are τ 3 /m 3 ≥ 1.41 × 10 −12 s/eV, τ 3 /m 3 ≥ 7.22 × 10 −13 s/eV and τ 3 /m 3 ≥ 1.50 × 10 −12 s/eV, respectively.
We also studied the effect of decay on the measurement of the standard parameters θ 23 and ∆m 2 32 . We found that if we include decay in our fit, the best-fit values for sin 2 θ 23 and ∆m 2 32 change significantly. The best-fit (sin 2 θ 23 , ∆m 2 32 ) obtained for the standard oscillation case from analysis of NOvA, T2K and both experiments combined are (0.45, 2.41 × 10 −3 eV 2 ), (0.52, 2.56×10 −3 eV 2 ) and (0.46, 2.51×10 −3 eV 2 ) respectively. On including decay in the fit, the corresponding best-fit points become (0.48,2.39 × 10 −3 eV 2 ) for NOvA, (0.62, 2.62 × 10 −3 eV 2 ) for T2K and (0.48, 2.52 × 10 −3 eV 2 ) for the NOvA and T2K combined. The best-fit sin 2 θ 23 is seen to be shifting to higher values. We also give the 95 % C.L. contours in the two-parameter space and the ∆χ 2 vs sin 2 θ 23 and ∆χ 2 vs ∆m 2 32 plots from which 1, 2 and 3σ ranges of these parameters can be read for both hypothesis, with and without decay. Decay is seen to shift the allowed range of sin 2 θ 23 significantly to higher values, thereby extending the allowed ranges in the higher octant. The reason for this behavior was discussed. The allowed range of ∆m 2 32 is also seen to change with inclusion of delay, albeit very mildly.
In conclusion, both T2K and NOvA, and in particular NOvA, seem to favor neutrino decay. Even though this conclusion is not statistically significant yet, it will be interesting to see the results from the forthcoming next-generation long baseline experiments like DUNE and T2HK. Invisible neutrino decay also results in shifting θ 23 to higher values and this would be again an interesting phenomenon to study at the next-generation experiments.