Invisible neutrino decay at long-baseline neutrino oscillation experiments

We perform an updated analysis of long-baseline accelerator data in the framework of neutrino oscillations in presence of invisible neutrino decay. We analyze data from T2K, NOvA and MINOS/MINOS+ and show that the combined analysis of all experiments improves the previous bound from long-baseline data by approximately one order of magnitude.


I. INTRODUCTION
The presence of three neutrino oscillations is, at present, observed in solar, atmospheric, reactor and accelerator experiments requiring the presence of nonzero neutrino masses and mixing angles [1][2][3].A nonzero neutrino mass automatically opens the possibility that, besides oscillating, neutrinos can decay.This idea, firstly proposed as a mechanism to solve the solar neutrino problem [4,5], is deeply investigated in the literature.Standard Model neutrino decays, both through radiative and non-radiative processes, already have strong cosmological bounds thanks to the high precision measurement of the cosmic microwave background [6,7].On the other hand, processes involving beyond Standard Model (BSM) physics such as are much less constrained.Here ν i (i = 1, 2, 3) is a neutrino mass eigenstate with mass m i , while ν and X are particles in the final state.In particular, ν can correspond to one or more neutrinos and X represents one or more non-observable particles, typically identified as scalar or pseudoscalar fields.The BSM decay in Eq. (1) can be classified as i) visible decay, in which at least one neutrino in the final state is active; ii) invisible decay, in which all ν and X are non-observable particles.In this last case, the final state neutrino particles are sterile neutrinos ν s .We are interested in discussing invisible neutrino decay.In this case, a beam of (relativistic) a christoph.ternes@lngs.infn.itb giulia.pagliaroli@lngs.infn.itν i -neutrinos with lifetime τ i , is depleted due to the invisible neutrino decay by the factor where E is the neutrino energy, L is the distance between the source and the detector, and is the decay parameter.It is evident that, for a given ratio L/E, the neutrino decay is only sensitive to decay parameters with α i ≥ E/L.Limits on α i have been derived using different neutrino sources.Astrophysical neutrinos, travelling long distances, are most promising for testing neutrino decay.Indeed, for electron antineutrinos, the most favourable combination is provided by the supernova SN1987A (L = 50 kpc, E ∼ 20 MeV): the observation of electron antineutrinos in Kamiokande-II [8] and IMB [9] and Baksan [10] yields the lower limit α 1 ∼ α 2 ≤ 10 −5 eV/s [11,12].High-energy neutrinos observed by IceCube [13] are, in principle, more powerful.However we need to identify the source position or the production process to bound the α i exploiting the observed flavor ratio [14].The recent discovery of the first steadystate source of high-energy neutrinos, NGC 1068 [15], allowed to investigate this possibility.However the present event rates are too low and the uncertainties on emission model too large to allow to set a bound on neutrino decay [16].Neutrino decay has been also suggested as an option to alleviate the tension between cascade and track events in IceCube high energy data [17].
Atmospheric neutrinos or long-baseline neutrino data have been used so far to bound the ν 3 lifetime.As a leading example, the combined analysis of Super-Kamiokande data together with long-baseline data from K2K and MINOS provides a lower bound on the decay parameter τ 3 /m 3 ≥ 9.3 × 10 −11 s/eV at 99% confidence level (C.L.) [18].Whereas, the combined analysis of NOvA and T2K initial data sets showed a small preference for the decay scenario with a best-fit of τ 3 /m 3 = 3.16 × 10 −12 s/eV.However, at 3σ the combined constraint was found to be τ 3 /m 3 ≥ 1.5 × 10 −12 s/eV [19].
In this paper, we improve the bound on the ν 3 lifetime by combining larger data sets from the MINOS/MINOS+, NOvA and T2K experiments under the assumption of neutrino decay into invisible daughters.

II. THE FORMALISM
Let us assume that active neutrinos are subject to both standard mixing and invisible decays.Therefore, propagating neutrinos mix among flavor eigenstates in an oscillatory timereversible manner and disappear due to timeirreversible decay.In this hypothesis, the evolution process is described by the Hamiltonian given by where the first two terms correspond to the standard vacuum and matter terms, , where E is the neutrino energy, G F the Fermi constant and N e the electron number density.The last term in Eq. ( 4) represents the neutrino decay part In this paper we are interested in the decay of ν 3 only, hence we set α 1 = α 2 = 0. Neutrino decay implies that the sum of the neutrino oscillation probabilities might be different from one, This is a feature which will be used in the analysis of MINOS/MINOS+ data below, where we include the neutral current data set, which would be affected due to Eq. ( 8).For discussion of the neutrino oscillation probability or formulas derived in this scenario, we refer the interested reader to Refs.[20][21][22][23][24][25].

III. DATA ANALYSIS
In this section we describe the data that we use in the analyses of this paper.In our analyses we keep the solar parameters fixed at the best fit values from Ref. [1].In addition, we have checked that reactor experiments are only sensitive to values of α 3 larger than the ones considered here.Therefore in our analyses we include priors on θ 13 and ∆m 2  31 as obtained from the combined analysis of Daya Bay [26] and RENO [27] data [1], which read The determination of sin 2 θ 13 is the same for normal and inverted neutrino mass ordering (NO and IO, respectively).In the following subsections we briefly detail the analysis strategy for MINOS/MINOS+, T2K, and NOvA, the experiments under consideration in this paper.

A. MINOS/MINOS+
The accelerator-based neutrino oscillation experiment MINOS studied neutrinos produced at Fermilab at the NuMI beam facility and detected them at two detectors placed at 1.04 km and 735 km distance.First the neutrino beam peaked at an energy of 3 GeV, which later was increased to a peak energy of 7 GeV.These phases were called MINOS and MINOS+, respectively.In the analysis of this paper we consider the full data set of both phases, corresponding to an exposure of 10.56 × 10 20 protons-on-target (POT) in MINOS and 5.80×10 20 POT in MINOS+ [28].We use the publicly available MINOS/MINOS+ code used in the analysis for the search for light sterile neutrinos, see Ref. [28].We modify this code to account for neutrino decay instead of active-sterile neutrino oscillations.The statisti-cal analysis is performed by using the following χ 2 definition as a function of the oscillation parameter set ⃗ p: where N dat,i and N exp,i are the observed and the predicted event rates, in the energy bin i.The covariance matrix V cov.accounts for the contributions from several sources of systematic uncertainties [28].The last term, ) where X = NO, IO for normal (inverted) ordering, contains the priors discussed above, see Eq. ( 9).The first two sums are taken over detectors (near and far) and over the different channels used in MINOS/MINOS+, which include disappearance and appearance channels using charged current events, and also the neutral current data which is sensitive to the effect described in Eq. ( 8).

B. T2K and NOvA
We also consider data collected by the modern experiments T2K [29] and NOvA [30].Both experiments observe events due to interactions of neutrinos and antineutrinos.In the case of T2K an exposure at Super-Kamiokande of 1.97×10 21  POT in neutrino mode and 1.63×10 21 POT in antineutrino mode has been reached up to now.NOvA, instead, has reached 1.36×10 21 POT in neutrino mode [31] and 1.25×10 21 POT in antineutrino mode.The analysis for these experiments is performed using GLoBES [32,33].We assume Gaussian smearing for the energy reconstruction in both experiments.In addition we use bin-to-bin efficiencies, which are adjusted to reproduce the best-fit spectra reported by the T2K and NOvA collaborations.Also included are systematic uncertainties to account for uncertainties in the predicted signal and background spectra, see Refs.[29,34] and [30].In the case of T2K and NOvA the χ 2 function is given by where again N dat,i and N exp,i are the data and prediction in bin i and the first sum is taken over the different oscillation channels: ν µ → ν µ , ν µ → ν µ , ν µ → ν e and ν µ → ν e .The second term contains the penalties σ i for the systematic uncertainties δ i , while the last term is equivalent to the one in the χ 2 function for MI-NOS/MINOS+.
In addition to the individual analyses of the T2K, NOvA and MINOS/MINOS+ experiments, we will also discuss a combined analysis of all accelerator data.Note that in this case the penalty term χ 2 pen. is added only once to the overall χ 2 function.FIG.1: The bounds on τ 3 3 obtained from the analysis of NOvA (red), T2K (blue), MINOS/MINOS+ (green) data, and the result from the combined fit (black).The solid (dashed) lines correspond to the analysis using normal (inverted) neutrino mass ordering.Also shown is the bound from Daya Bay (orange), which, as explained in the text, is weaker.

IV. RESULTS AND DISCUSSION
In this section we discuss the results obtained in this paper.When performing the analyses we always consider either normal or inverted neutrino mass ordering.When plotting the results they are always plotted with respect to the best fit value obtained in each ordering separately.The main result is presented in Fig. 1.Apart from long-baseline bounds we also show the bound from an analysis of Daya Bay data.The analysis procedure is the same as in Ref. [35], but with the decoherence probability replaced with the one including neutrino decay.As anticipated in Sec.III, reactor experiments are sensitive to smaller (larger) values of τ 3 /m 3 (α 3 ).This justifies to include the priors of Eq. ( 9) in our analyses of accelerator data.The different colors correspond to the analyses of the individual experiments, while solid (dashed) lines are obtained using normal (inverted) neutrino mass ordering.As can be seen, the χ 2 profiles are very similar for both orderings.The only major difference appears in the analysis of T2K data, where the bound for IO is weaker than the bound for NO.This happens, because in the case of inverted neutrino mass ordering T2K data prefer quite small values of |∆m 2  31 |, see Ref. [34].Therefore, there is a tension between T2K data (in inverted ordering) and the reactor prior imposed in our analysis.This does not happen in normal ordering, where there is better agreement.Therefore the analysis in IO is penalized resulting in a weaker bound.Note that in the case of NOvA we find a best fit value with finite lifetime.As can be seen from the profile, this best fit is statistically not very significant, and it disappears after combining with the other experiments.
In Fig. 2 we show the contours in the τ 3 /m 3 − sin 2 θ 23 -plane at 90% confidence level for two degrees of freedom for the different experiments (colored lines) and for the combined analysis (shaded region).From this figure we see that there is a correlation between τ 3 /m 3 and sin 2 θ 23 .Note, in particular in the contours for MINOS/MINOS+, that values of sin 2 θ 23 away from maximal mixing result in weaker bounds on τ 3 /m 3 .This is the reason why the combined analysis of all data provides a better bound than MINOS/MINOS+ data alone.The inclusion of NOvA and T2K data pushes the allowed region in sin 2 θ 23 closer to maximal mixing 1 which results in a stronger bound on τ 3 /m 3 .This feature can also be seen by observing how the best fit point from MINOS/MINOS+, indicated by the green arrow in Fig. 2, is outside of the allowed parameter space from T2K and NOvA.Therefore, upon combining all data we are left with a more restrictive parameter space for τ 3 /m 3 .
We summarize the bounds that we obtain in this paper in Tab.I. Our analyses of NOvA and T2K data update the analysis performed in Ref. [19], while the MINOS analysis updates Ref. [36].The individual bounds from NOvA, T2K and MINOS/MINOS+ data improve the previous bounds all by a factor of 3-4.From our combined analysis we find that long-baseline accelerator data provide a bound of τ 3 /m 3 ≥ 2.4 × 10 −11 s/eV at 90% C.L. (13) for both neutrino mass orderings.Therefore, we find that our global analysis of long-baseline accelerator data improves the previous bound from Ref. [19] (τ 3 /m 3 ≥ 2.4×10 −12 s/eV at 90% C.L.) by one order of magnitude.Note that, unlike in 1 At this point it should be noted that the results of this paper are not in tension with the ones in Ref. [1].The reason why in the present analysis MINOS/MINOS+ data disfavor maximal mixing, while it is allowed in Ref. [1], is the inclusion of the prior in Eq. ( 9).In Refs.[1,28], ∆m 2 31 is also fitted with the data resulting in a different behavior.The same argument applies to the analyses of T2K and NOvA data.We always make sure that our analyses reproduce the official results when making the same assumptions as in the original references.our case, no prior on ∆m 2  31 was used in Ref. [19].However, we have verified that the bound obtained in this analysis, Eq. ( 13), does not depend on this prior.The same bound is obtained when allowing ∆m 2  31 to vary freely in the analysis.This is expected since the determination of ∆m 2  31 from long-baseline accelerator data is more precise than the one from reactor data.

Experiment
As a final remark, let us mention that invisible neutrino decay does not have any impact on the preferred neutrino mass ordering.The preference obtained from the combined analysis is the same as the one obtained in the standard analysis.In addition, we also found that invisible neutrino decay neither helps to alleviate the disagreement in the measurement of δ in T2K and NOvA.

V. CONCLUSIONS
We have performed an analysis of longbaseline accelerator data to place a new bound on invisible neutrino decay.We find that the bound is driven by MINOS/MINOS+ data, and then further improved by the inclusion of NOvA and T2K data.This improvement happens due to the partial breaking of a correlation between sin 2 θ 23 and τ 3 /m 3 in the MINOS/MINOS+ analysis.Basically, the inclusion of NOvA and T2K data drives sin 2 θ 23 closer towards maximal mixing, which improves the bound on τ 3 /m 3 .
The individual up-to-date bounds on τ 3 /m 3 obtained in this analysis improve the former bounds from Refs.[19,36] by a factor of 3-4.A stronger bound than the one obtained here has been reported in Ref. [18] using an older MINOS dataset and also data from Super-Kamiokande.However, it should be noted that the bound in this reference shows non-Gaussian behavior due to the presence of a local minimum.Not taking into account this bound, Eq. ( 13) is the strongest bound on the neutrino lifetime obtained from neutrino oscillation data so far2 .
This laboratory bound will be further improved at future facilities.It has been shown in Refs.[22,39,40] that DUNE and T2HK can improve our bound by a factor of about 2. A bound of similar strength is also expected at the JUNO experiment [21,41].The future facility ESSνSB would only produce a bound of similar strength as the current one [42].Even stronger bounds can be expected from future atmospheric experiments, as shown in Refs.[43][44][45].
Finally, we must remark that terrestrial experiments are not the only way to obtain bounds on invisible neutrino decay.Constraints on the neutrino lifetime have also been derived from astrophysical and cosmological observations, see e.g.Refs.[46][47][48][49][50][51][52].These are significantly more stringent, up to a level [51] of τ /m ν > 8 × 10 6 s/eV, than the bounds obtained here or that can be expected from any other current or future neutrino oscillation experiment.

FIG. 2 :
FIG.2: Contours in the τ 3 /m 3 − sin 2 θ 23 -plane at 90% confidence level for two degrees of freedom for normal (upper panel) and inverted (lower panel) neutrino mass ordering.We show the contours from the analysis of T2K (blue), NOvA (red), MINOS/MINOS+ (green) data and from the combination (shaded region).The arrows and the dot indicate the best fit values found in each analysis.

TABLE I :
The bounds at 90% C.L. on τ 3 /m 3 obtained in the analyses of this paper.