Constraints on neutrino decay lifetime using long-baseline charged and neutral current data

We investigate the status of a scenario involving oscillations and decay for charged and neutral current data from the MINOS and T2K experiments. We first present an analysis of charged current neutrino and anti-neutrino data from MINOS in the framework of oscillation with decay and obtain a best fit for non-zero decay parameter $\alpha_3$. The MINOS charged and neutral current data analysis results in the best fit for $|\Delta m_{32}^2| = 2.34\times 10^{-3}$~eV$^2$, $\sin^2 \theta_{23} = 0.60$ and zero decay parameter, which corresponds to the limit for standard oscillations. Our combined MINOS and T2K analysis reports a constraint at the 90\% confidence level for the neutrino decay lifetime $\tau_3/m_3>2.8 \times 10^{-12}$~s/eV. This is the best limit based only on accelerator produced neutrinos.

We can parametrize the decay by the ratio of the lifetime parameter τ i and the mass m i for each of the mass eigenstates i = 1, 2, 3. The role of invisible neutrino decay was investigated for the solar neutrino anomaly [15,32,33] and showed no evidence for the dominance of the decay scenario. From these we can constrain values of the decay parameter, τ 2 /m 2 > 8.7 × 10 −5 s/eV at 90% C.L., where τ 2 and m 2 are respectively the lifetime and the highest mass eigenstate in a two generation scenario [15].
In the visible decay scenario, we can search for ν e → ν e conversion using a pure ν e source such as the Sun [24]. The null results from the solar ν e appearance impose a constrain on the decay parameter: τ 2 /m 2 > 6.7×10 −2 s/eV from the KamLAND experiment [27].
Neutrinos produced in supernovas are interesting to investigate for the presence of decays due to the large distance traveled by the neutrino. From the observation of electron neutrinos from SN1987A [34,35] we should have a lower limit in neutrino lifetime. Otherwise we could not see any signal. For larger values of the mixing angle, such as the current LMA solution for the solar neutrino anomaly, no constraint is possible [36]. Other possibilities for neutrinos coming from a supernova include the neutrino decay catalyzed by a very dense media [25,29] where the matter effects can increase the decay rate of neutrinos. The diffuse supernova neutrinos (neutrinos coming from all past supernova explosions) [17], can provide very robust sensitivity in the range of τ/m < 10 10 s/eV [18,19].
Astrophysical neutrino sources megaparsec away can generate all neutrino flavors. Due to the long distance from the sources we are in the limit L → ∞, where all dependence on the lifetime parameter τ i fades away. However, if we have a precise determination of the ratio of flavors of these neutrinos we can discriminate the case with and without decay [37][38][39][40].
Concerning the accelerator and atmospheric neutrinos we can test the decay scenario of the third generation of neutrino mass eigenstates investigating how the τ 3 /m 3 decay parameter changes the ν µ → ν µ survival probability [20,21]. The MI-NOS experiment made a search for the decaying neutrino and constrained the lifetime to τ 3 /m 3 > 2.1 × 10 −12 s/eV at 90% C.L., using both neutral and charged current events [23]. The combined analysis of Super-Kamiokande atmospheric neutrinos with K2K and MINOS accelerator neutrinos show a 90% C.L. lower bound value of τ 3 /m 3 > 2.9 × 10 −10 s/eV [21].
This article is organized as follows, in Section 2 we discuss our neutrino decay scenario. Next, we introduce the χ 2 analysis developed for the neutrino charged and neutral current data of the MINOS experiment in Section 3. We then present our bounds for the neutrino lifetime based on a charged current analysis (Section 4) and on a combined charged and neutral current analysis (Section 5), comparing with previous limits.

Decay model for neutrinos
We are going to introduce the neutrino evolution equation in which the neutrino can decay. This is made by putting an imaginary part related to the neutrino lifetime, which is the ratio α 3 ≡ m 3 /τ 3 , in the evolution equation. We are going to assume the decay of the heaviest state, ν 3 → ν s + φ, where both final products are invisible. The two-generation system is considered, which is adequate to describe the muon neutrino and anti-neutrino data from MINOS. The evolution equation is where the stateν ≡ ν µ ν τ , and E is the neutrino energy. U is the usual rotation matrix, where c 23 ≡ cos θ 23 and s 23 ≡ sin θ 23 . The same evolution equation applies for anti-neutrinos as well. From Eq. (1) we obtain the muon neutrino survival probability as P(ν µ → ν µ ) = cos 2 θ 23 + sin 2 θ 23 e − α 3 L where L is the distance traveled by the neutrinos. You can notice that a non-zero decaying parameter α 3 changes only the amplitudes: the constant amplitude (first term of the equation above), and the oscillation amplitude (second term). Both amplitudes are damped but the oscillation phase does not change. In the two-ν standard oscillation probability formula we have the symmetry cos 2 θ 23 ↔ sin 2 θ 23 , but in Eq. (3) the symmetry is broken, and then we should scan the parameter space of the variable sin 2 θ 23 in the range (0, 1). This broken symmetry will appear in our plots later.
The limiting case where the oscillations are induced only by decay, ∆m 2 32 → 0, can also be tested. In this case the probability assumes the simple form where now we have two free parameters, the mixing amplitude sin 2 θ 23 and the decay parameter α 3 . Even at this limit we also have an asymmetry between sin 2 θ 23 and cos 2 θ 23 .
In the standard neutrino oscillation scenario, the sum of the probabilities over the active states is equal to unity. Then the spectrum of neutral current (NC) events is not effected by active oscillation, which means that the expected number of NC events is the same with or without oscillations. But in the extended scenario involving sterile neutrinos the sum of probabilities, β P(ν µ → ν β ), where β = µ, τ, obviously does not sum up to 1. This is discussed in the general context of non-unitary neutrino evolution in Ref. [41].
We can compute the conversion probability for the two-ν oscillations with decay scenario, Thus, from the Eq. (6) we observe that there will be an effect on the neutral current interaction events under the oscillations with decay model.

Analysis of MINOS Charged and Neutral Current Data
We have performed a combined analysis using the published data of charged and neutral current MINOS analyses. MINOS is a long-baseline neutrino experiment [42] using two detectors and exposed to a neutrino beam produced at Fermilab. The NuMI beam line is a two-horn-focused neutrino beam that can be configured to produce muon neutrinos or anti-neutrinos. The Near Detector is located at Fermilab, around 1 km from the NuMI target and the Far Detector is 735 km far from the target.
The first data set used in our analysis comprises the charged current (CC) contained-vertex neutrino disappearance data [9] using the ν µ enhanced beam with exposure of 10.71 × 10 20 protons on target (POT), with 23 points of non-equally divided bins of energy up to 14 GeV; the second is the CC contained-vertex anti-neutrino disappearance data [9] from theν µ enhanced beam with 3.36×10 20 POT, using 12 points for energies up to 14 GeV as well; and the third is the neutral current data [43], based on 7.07 × 10 20 POT. The spectrum of NC events is described as a function of a reconstruted energy, E reco . We then use the information from a previous MINOS analysis [23] to separate the NC data into two bins: (a) events with E reco < 3.0 GeV and (b) events with 3.0 < E reco < 20.0 GeV, which have median neutrino energies of 3.1 and 7.9 GeV, respectively.
For the χ 2 calculation we use the following function in bins of energy, where N th were read off from References [9] and [43].
The total error used for both neutrino and anti-neutrino CC data sets is given by where σ data i is the total error of the data, σ stat,th i and σ syst,th i are the statistical and systematic error of the prediction, respectively. We have σ data [44]. For the NC data set we do not have the total error of the data events. We then calculate it summing in quadrature the statistical and systematic errors. The former is the square root of the number of data events, and the latter is an estimate based on the systematic error of the extracted expectation of events.
We scan the parameter region in the variables τ 3 /m 3 , sin 2 θ 23 and ∆m 2 32 and for the β parameter of our function χ 2 = χ 2 (τ 3 /m 3 , sin 2 θ 23 , ∆m 2 32 , β). The χ 2 is then marginalized over the nuisance β parameter to obtain the effective χ 2 eff as a function of the other parameters To have a guess about the range of the decay parameter we can probe, we will assume the argument of the exponential term in Eq. (3) is of the order of unity, where we can get the most sensitivity. Using the MINOS distance and the energy range of E ν = (0.5, 10) GeV, we get τ 3 /m 3 ∼ 10 −13 − 10 −11 s/eV.

Results of the Charged Current Analysis
To test the correctness of our χ 2 analysis, we first consider the standard oscillation scenario, corresponding to the limit α 3 → 0 in Eq. (3). The results shown here are from neutrino and antineutrino charged current data obtained from Ref. [9]. However, we also used our calculation on the neutrino charged current data from Ref. [45]. For both data sets we obtain good concordance with the MINOS allowed region and best fit.
In Table 1 we show the consistency between the best fit parameters obtained by our standard oscillation analysis (accelerator data only) and MINOS [9] (which includes accelerator and atmospheric data). We obtain that the best fit of the oscillation parameters is for non-maximal mixing angle, sin 2 2θ 23 = 0.92 and for |∆m 2 32 | = 2.38 × 10 −3 eV 2 . The χ 2 of the best fit point is 19.80 for 31 degrees of freedom, which corresponds to a G.O.F. (Goodness of Fit) of 94.0%. The ∆χ 2 for the no-oscillation hypothesis compared to the standard oscillation scenario, standard oscillation is equal to 304.9, which means that we can exclude the nooscillation hypothesis with remarkable precision. The allowed regions in the (sin 2 2θ 23 -∆m 2 32 ) plane for neutrinos and anti-neutrinos under the standard oscillations model is shown in Fig. 1. The anti-neutrino data is compatible with the neutrino data but does not contribute significantly to improve the region of parameters in the combined CC analysis due to its small statistics compared to the neutrino one. Our results show that we are able to reproduce the allowed regions of square mass difference and mixing angle of the MINOS analyses.
Having confirmed the correctness of the analysis, we use it to test the oscillation with decay scenario. The procedure is the same as in the standard oscillation case, but now we have three free parameters, ∆m 2 32 , sin 2 θ 23 and τ 3 /m 3 . For the combined (neutrino and anti-neutrino) CC analysis under the oscillation with decay framework we find a finite value for the best fit of   Table 1), The effect of a finite value of τ 3 /m 3 results in a slightly lower value of ∆m 2 32 than the one obtained for standard oscillation. The χ 2 obtained is 19.28 for 30 degrees of freedom, which corresponds to a G.O.F. of 93.4%. The pure decay hypothesis is also considered in our CC analysis resulting in a χ 2 of 51.18 for 31 degrees of freedom, excluding this hypothesis at the 5.6 standard deviation level.
The projections of ∆χ 2 for each parameter are obtained minimizing the function χ 2 = χ 2 (τ 3 /m 3 , sin 2 θ 23 , ∆m 2 32 , β) accordingly to the parameters. We then use the function ∆χ 2 ≡ χ 2 − χ 2 min , where χ 2 min is the global minimum value for the oscillation with decay scenario. Next, we present the projections of ∆χ 2 as a function of τ 3 /m 3 , ∆m 2 32 and sin 2 θ 23 . The one-dimensional projection for the τ 3 /m 3 parameter, shown in Fig. 2 (dashed curve), allows us to get a lower bound on τ 3 /m 3 or equivalently an upper bound on α 3 . The allowed values for the neutrino decay lifetime is τ 3 /m 3 > 2.0 × 10 −12 s/eV (10) at the 90% C.L. (Table 1). We can compare the value we obtain for the oscillation with decay model to the one from MI-NOS [23], which used both charged and neutral current data in their analysis, but with lower statistics than the extracted data used here. We see that our combined CC result does not improve the MINOS lower limit which is τ 3 /m 3 > 2.1 × 10 −12 s/eV. An interesting behavior appears in Fig. 3 (left), where we compare the allowed values for ∆m 2 32 with and without decay. When we include a finite value for the τ 3 /m 3 parameter, smaller values of ∆m 2 32 are allowed for a certain confidence level. Then, under this scenario, we can explain the muon disappearance signal seen by MINOS as due partially to oscillation and partially to decay, since the limit ∆m 2 32 → 0 in Eq. (3) corresponds to the pure decay model.
The broken symmetry, cos 2 θ 23 ↔ sin 2 θ 23 mentioned after Eq. (3), is manifested in Fig. 3 (right), where the curve is broader for larger values of sin 2 θ 23 . This effect can be understood investigating the coefficient on the first term of Eq. (3). For sin 2 θ 23 > cos 2 θ 23 the decay term is more relevant, and for the opposite case the decay contribution is suppressed.
We also present two-dimensional projections of the allowed three-dimensional region after normalization with respect to the undisplayed parameter. Figure 4 shows the 90% C.L. allowed region in the (τ 3 /m 3 -∆m 2 32 ) and (τ 3 /m 3 -sin 2 θ 23 ) planes for the oscillation with decay scenario and our best fit point (CC data only). We can observe that for smaller values of τ 3 /m 3 , when the effects of the decay are larger, the contours allow smaller values of ∆m 2 32 and higher values of sin 2 θ 23 (the same behaviour shown in Fig. 3).
The 90% C.L. allowed regions for the oscillation parameters, ∆m 2 32 and sin 2 θ 23 , is presented in Fig. 5 both with and without decay. In agreement with the information shown in the discussion of Fig. 3 and Fig. 4 we can see that the decay allows a region of smaller values of ∆m 2 32 and larger values of sin 2 θ 23 than the standard oscillation model.
The best fit points are also shown in Fig. 5 for both models (with and without decay). The standard oscillation analysis obviously shows two possible values of θ 23 due to the symmetry cos 2 θ 23 ↔ sin 2 θ 23 . Since the oscillation with decay model does not manifest such a symmetry we find that the best fit point for this scenario is in the θ 23 > 45 • octant, with a value of sin 2 θ 23 = 0.63.

Results of the Combined Charged and Neutral Current Analysis
In this section we present the results obtained for the combined charged and neutral current data. We expect an effect on the neutrino lifetime due to the inclusion of the neutral current data in the analysis as described in Eq. 6. First, we calculate the ratio for each bin of reconstructed energy extracted from the NC data, where N data , N bg , and N th are the number of data, background and expected events, respectively. The ratios we obtain are in good agreement with the ones from MINOS [43].
Under the oscillation with decay model the best fit point for the combined CC and NC analysis can be found in Table 1. The |∆m 2 32 | and sin 2 θ 23 values are consistent with the ones for standard oscillations. We obtain a neutrino lifetime τ 3 /m 3 → ∞, which corresponds to the case of no decay. If we compare this result to the best fit of the CC analysis, the decay effect becomes less relevant due to the inclusion of NC data, implying that the best fit value of |∆m 2 32 | increases (from 2.30 × 10 −3 eV 2 to 2.35 × 10 −3 eV 2 ).
The χ 2 of the best fit for the combined CC and NC analysis is 20.52 for 31 degrees of freedom resulting in a G.O.F. of 92.4%. The pure decay model is also used in our analysis, being excluded at the 8.8 standard deviation level. Our best fit value of τ 3 /m 3 for the pure decay scenario is close to the one from MINOS [23].
The one-dimensional projection for the τ 3 /m 3 parameter in the combined CC and NC analysis is shown in Fig. 2 (solid curve). We find a lower limit of τ 3 /m 3 > 5.7 × 10 −12 s/eV (12) at the 90% C.L., which improves the MINOS [23] limit by a factor of 2.7.
The two-dimensional projections in Fig. 4 show the effect of including the NC data into our analysis. The contours at the 90% C.L. become more symmetric when compared to the contours for CC data only (see also  Figure 6: The extracted charged current neutrino (left) and anti-neutrino (right) spectra of MINOS shown together with the curves for the following hypotheses: (i) no-oscillation, (ii) standard oscillation at the best fit parameters, and (iii) oscillation with decay at the best fit for ∆m 2 32 and sin 2 θ 23 and the 90% C.L. value of τ 3 /m 3 for the CC data only (τ 3 /m 3 = 2.0 × 10 −12 s/eV) and for the combined CC with NC data (τ 3 /m 3 = 5.7 × 10 −12 s/eV). The hatched areas show the background from neutral current events.
We can understand the NC effect by observing that this data set does not suggest evidence for decay (τ 3 /m 3 → ∞), differently from the CC analysis where we find a finite value for τ 3 /m 3 . Thus, since the decay effect becomes less relevant, the symmetry between cos 2 θ 23 and sin 2 θ 23 manifests itself again, resulting in a contour of the oscillation with decay analysis under the combined CC and NC data close to the contour of the standard oscillation scenario.
The extracted MINOS spectrum data for CC neutrinos and anti-neutrinos are shown in Fig. 6. In both plots we also show the curves for the no-oscillation hypothesis, the best fit of standard oscillations, and the background from NC events. In addition, we also show the curves for the oscillation with decay scenario at the best fit for ∆m 2 32 and sin 2 θ 23 and using the 90% C.L. value of the τ 3 /m 3 parameter for the CC data only (τ 3 /m 3 = 2.0 × 10 −12 s/eV) and for the combined CC with NC data (τ 3 /m 3 = 5.7 × 10 −12 s/eV). The reason for using the 90% C.L. values of τ 3 /m 3 is that the curves for the best fit parameters show no significant difference if compared to the standard oscillations curve. We can observe that the curve using the higher value of τ 3 /m 3 , which corresponds to smaller effect of decay, is closer to the one for standard oscillations. We also observe that the inclusion of the decay scenario into the oscillation model worsens the fit for both neutrino and anti-neutrino analyses.

Conclusions
The present scenario of standard neutrino oscillation observed by different experiments shows a strong case for nonzero neutrino masses and mixing. A direct consequence of non-zero neutrino masses is that the neutrino can decay, but usually with much longer lifetimes for decays like ν → 3ν or ν → ν + γ. The only expectation to observe sizeable effects from neutrino decays is from ν → ν s + φ, where the decay products are sterile states.
We use the more updated charged and neutral current data of the MINOS experiment for accelerator produced neutrinos and anti-neutrinos to study the impact of a non-zero decay parameter, α 3 = m 3 /τ 3 , in the ν µ → ν µ channel. Using the charged current data only we have found that the best fit scenario is for a finite τ 3 /m 3 parameter, corresponding to a neutrino lifetime τ 3 /m 3 = 1.2 × 10 −11 s/eV.
Based on the combined charged and neutral current analysis the best fit for the oscillation with decay model indicates a zero velue for the α 3 parameter, that is equivalent to the standard oscillation model. Using the combined analysis we could improve the previous constraint on the allowed neutrino decay lifetime τ 3 /m 3 > 5.7 × 10 −12 s/eV at the 90% C.L., which is the best limit based only on accelerator produced neutrino data.
We have found that the 90% C.L. range for |∆m 2 32 | is (1.95 − 2.54) × 10 −3 eV 2 and for the mixing angle is sin 2 θ 23 = (0.32 − 0.75) in the oscillation with decay scenario. We could understand the effect of the non-zero decaying parameter for the CC data analysis, implying smaller values of ∆m 2 32 and larger values of sin 2 θ 23 than the ones in standard oscillation scenario. The inclusion of NC data, which is compatible to standard oscillations, constrains the decay parameter and makes the decay scenario more constrained than in the CC analysis only.

Acknowledgments
We would like to thank C. Castromonte for helping on computational tasks and M. Goodman for valuable discussion and careful reading of the manuscript. A.L.G.G. was supported by