Constraining the invisible neutrino decay with KM3NeT-ORCA

Several theories of particle physics beyond the Standard Model consider that neutrinos can decay. In this work we assume that the standard mechanism of neutrino oscillations is altered by the decay of the heaviest neutrino mass state into a sterile neutrino and, depending on the model, a scalar or a Majoron. We study the sensitivity of the forthcoming KM3NeT-ORCA experiment to this scenario and find that it could improve the current bounds coming from oscillation experiments, where three-neutrino oscillations have been considered, by roughly two orders of magnitude. We also study how the presence of this neutrino decay can affect the determination of the atmospheric oscillation parameters $\sin^2\theta_{23}$ and $\Delta m_{31}^2$, as well as the sensitivity to the neutrino mass ordering.


Introduction
Over the last two decades or so, we have found overwhelming evidence for oscillating neutrinos. The oscillatory behavior can be described in terms of six parameters: the solar and atmospheric mass splittings ∆m 2 21 and ∆m 2 31 , the solar angle θ 12 , the atmospheric angle θ 23 , the reactor angle θ 13 and ules (DOMs) and 20 m horizontal spacing between the detector lines. This will result in an energy threshold of only a few GeV, enhancing its sensitivity to lower energies with respect to ANTARES, that had a threshold of around 20 GeV. Therefore, ORCA seems a very promising candidate not only to improve the current sensitivity to neutrino oscillation parameters, but also to look for signals of physics beyond the Standard Model. The presence of new physics might change the well established picture of neutrino oscillations and, hence, it is crucial to improve the precision of current measurements to look for potential deviations of the standard scenario, which would only arise at sub-leading order.
One of these new physics scenarios, able to alter the neutrino oscillation pattern, is based on the existence of unstable neutrinos. Taking, for example, the Majoron model [12][13][14][15][16], a neutrino ν i can decay into a lighter neutrino ν j and a new boson, the Majoron J, through ν i → ν j + J or ν i → ν j + J. Likewise, Dirac neutrinos could decay through ν iL → ν jR + ξ into a scalar ξ and a light right-handed neutrino ν jR . If the decay product ν j is an active neutrino, we talk about a visible neutrino decay, otherwise it is an invisible decay. In this work we will focus on the latter one. For previous studies on visible neutrino decay at current or future experiments, we refer the interested reader to Refs. [17][18][19][20][21][22].
Actually, the idea of unstable neutrinos is not new, since it was already proposed to explain the solar neutrino anomaly with a decaying mass state ν 2 [23]. However, neutrino decay alone could not explain the solar neutrino deficit and flavor oscillations were needed anyway [24]. Therefore, it is usually assumed that this process can appear at subdominant level in combination with neutrino oscillations, as studied in Refs. [25][26][27][28]. For the invisible neutrino decay, the best bound from oscillation experiments on the ν 2 lifetime comes from the combination of solar and reactor data and corresponds to approximately τ 2 /m 2 > 2 × 10 −3 s/eV at 90% C.L. [29]. See also Ref. [30] for a similar result and Ref. [31] for a recent forecast analysis on these parameters using dark matter detectors.
The decay of the third neutrino mass eigenstate, ν 3 , was also considered as an attempt to explain the atmospheric neutrino problem, but it was found that also here it can contribute only with sub-leading effects [32]. Several studies using atmospheric and long-baseline experiments have been performed in this direction. For instance, the analysis in Ref. [33] combines neutrino data from Super-Kamiokande, K2K and MINOS to obtain the limit τ 3 /m 3 > 2.9 × 10 −10 s/eV at 90% C.L. The authors of Ref. [34], on the other hand, combine long-baseline neutrino data from T2K and MINOS, obtaining the bound τ 3 /m 3 > 2.8 × 10 −12 s/eV at 90% C.L. Note, however, that these results have been derived under the two-neutrino approximation and, therefore, a full three-neutrino analysis might loosen this bound. Recently, following a three-neutrino approach, a new constraint on the neutrino decay lifetime has been calculated from the combination of T2K and NOνA data [35], giving approximately τ 3 /m 3 > 2 × 10 −12 s/eV at 90% C.L.
Prompted by the good sensitivity of the forthcoming experiment KM3NeT-ORCA to the atmospheric neutrino oscillation parameters in the GeV energy range, in this work we study whether it can provide better bounds on the invisible neutrino decay, in comparison to other neutrino oscillation experiments. We comment on current bounds from other astrophysical and cosmological probes in Sec. 5.
This letter is organized as follows. Sec. 2 provides an introduction to the unstable neutrino scenario considered in this work. In Sec. 3, we explain how the analysis and simulation of the ORCA experiment is performed. In Sec. 4, our main results are presented and discussed in detail. Finally, in Sec. 5, we summarize the most relevant outcome of our work and give some concluding words.

Invisible neutrino decay
Here we discuss how the presence of an invisible neutrino decay would affect the calculation of the neutrino oscillation probability. Besides the three light known neutrinos, we consider the presence of a fourth sterile neutrino, ν 4 . Along this work, we will assume the decay of the heaviest mass eigenstate (ν 3 in NO) to this new neutrino state where J is a pseudo-scalar singlet, or Majoron. We assume that there is no mixing among the three active neutrinos and the sterile one, so it can not oscillate back into an active state. Therefore, the neutrino mixing matrix in vacuum is given by the standard three-family mixing matrix U , where the Greek index α = e, µ, τ indicates the flavor eigenstates and the Latin index k = 1, 2, 3, the mass eigenstates. Because of the absence of active-sterile mixing, the propagation of the active states is not affected by the presence of ν 4 . For the mass spectrum we assume normal mass ordering for the active states and a fourth state ν 4 = ν s , lighter than the decaying one m 4 < m 3 . Note that this is not in tension with the results of Refs. [36][37][38], since we assume that the sterile state cannot oscillate back into active states.
To take into account the neutrino decay in the evolution process, we modify the neutrino Hamiltonian, including a decay constant α 3 = m 3 /τ 3 , where m 3 is the heaviest neutrino mass and τ 3 is its rest-frame lifetime. Hence, the full neutrino Hamiltonian can be written as where the first two terms correspond to the standard vacuum and matter terms, namely with V = 2E √ 2G F N e . E is the neutrino energy, G F the Fermi constant and N e the electron number density. Finally, the last term in Eq. (3) represents the neutrino decay part Then, effectively, the only change to the standard oscillation picture is a shift in the 33 entry of the Hamiltonian in the mass basis, from ∆m 2 31 to ∆m 2 31 − iα 3 . Note, however, that this implies that the sum of the neutrino oscillation probabilities might be different from one, Therefore, if the heaviest neutrino mass eigenstate ν 3 decays, apart from a changed oscillatory pattern, we could also have missing neutrinos. In order to show the effect of the decay, we present in Fig. 1 the difference in the survival probability (left panel) of atmospheric muon neutrinos with and without decay, ∆P αβ = P decay αβ − P standard αβ , for a value of α 3 = 10 −5 eV 2 . The right panel shows the analogous result for the electron neutrino appearance probability. Note that this value of the decay constant is rather large and we have chosen it for illustrative purposes only. As it is seen in the figures, the main effect of the decay is concentrated in the region of neutrino energies close to the resonance (∼ 3-8 GeV) for values of the zenith angle at which matter effects are more relevant. This results in a softening of the oscillation pattern, as better illustrated in Fig. 2. There, we show the muon neutrino survival probability P µµ as well as the electron neutrino appearance probability P µe for a particular arrival direction, corresponding to cos θ Z = −0.82. The solid lines in the plot correspond to the standard case, without neutrino decay, while the dashed, dashed-dotted and dotted lines have been obtained assuming values of the decay constant α 3 equal to 10 −5 , 10 −4 and 10 −3 eV 2 , respectively. In this figure, one can see how the presence of the invisible neutrino decay leads to a reduction of the oscillatory behavior, that becomes almost suppressed for larger values of α 3 , with special impact for the values of cos θ Z and neutrino energies close to the matter effect resonance at around 6-7 GeV. An interesting result from this suppression is the increment of the survival probability at some regions in the plane cos θ Z -E, which seems contradictory with the idea of decaying neutrinos. This comes from the softening of flavor oscillations in the presence of the ν 3 decay, since it washes out oscillations controlled by ∆m 2 31 while keeping untouched those driven by ∆m 2 21 . As a consequence, an almost averaged oscillation pattern appears, which enhances the survival muon neutrino probability at its minima. However, Eq. (6) holds and the sum of neutrino probabilities below unity tells us that they are actually decaying.

Numerical analysis
In this section, we present the numerical procedure followed to simulate the neutrino signal in ORCA. First, we explain how to calculate the number of events and the χ 2 functions. Next, we provide some details about how to handle the systematic uncertainties in the data analysis.

Simulation of the neutrino signal in ORCA
The simulation of the neutrino signal expected in ORCA requires the knowledge of the atmospheric neutrino flux, that we take from [39] and further modulate with flavor oscillations. The conversion probability, from the neutrino creation point in the atmosphere to the detector after traversing  the Earth, is numerically calculated considering three-neutrino oscillations in matter. In order to do so, we discretize the neutrino path and consider the matter density at each point using the PREM profile [40]. The sensitivity tests presented in this work are performed simulating the number of events detected by ORCA in a binned area of the parameters cos θ z,rec and log 10 (E ν,rec /GeV), where θ z,rec is the reconstructed zenith angle 1 and E ν,rec is the reconstructed neutrino energy. Following the indications in Ref. [11], we divide the reconstructed parameter ranges, cos θ z,rec ∈ [−1, 0] and ε rec = log 10 (E ν,rec /GeV), in 20 bins each, where E ν,rec ∈ [1, 21] GeV. We have also considered larger values for the maximum reconstructed neutrino energy, but the results are essentially unchanged.
Given the incapability of ORCA to distinguish neutrinos from antineutrinos, both contributions are summed in each bin. Our analysis, however, makes a distinction between the two different topologies produced in a neutrino interaction with the sea water molecules, namely track-like or showerlike 2 . A track is an elongated signal of deposited energy in ORCA's photomultipliers (PMTs), which is mostly produced when a ν µ (or its antiparticle) interacts through charged-current (CC) interactions, producing a muon that travels a long distance before losing all its energy. The same topology can be produced in a ν τ CC interaction if the generated tau decays into a muon, which produces the track. Therefore, this constitutes an unavoidable background for muon neutrinos detected through CC. On the other hand, in all other cases (ν e CC, ν τ CC with the tau not decaying into a muon and all flavor neutral-current interactions) the neutrino energy is quickly lost into an electromagnetic cascade, a hadronic cascade or both, depending on the interaction, giving rise to a shower-like topology. In current large-volume neutrino telescopes, like ANTARES or IceCube, the directionality of an event producing a shower is measured with a large uncertainty. However, thanks to the multi-PMT characteristic of ORCA's DOMs, the directionality of a showerlike event will be known with a precision of a few degrees, where the exact accuracy depends on the incoming neutrino energy (less than 10 degrees for energies larger than 5 GeV) [11].
The expected number of events per bin, N cα ij , corresponding to a given interaction channel, c, and a neutrino flavor, α, contributing to a given topology, is calculated from the convolution of the neutrino flux at the detector with the corresponding cross section, detector resolutions and detector effective mass, where t is the total time of data acquisition, m p the proton mass, ε true = log 10 (E ν,true /GeV). M cα eff , R cα and φ cα det are the detector effective mass, the detector resolution and the number of neutrinos per second at the detector for the corresponding interaction channel and neutrino flavor, respectively. The indices i and j refer to the ith bin in reconstructed zenith angle θ z,rec and the jth bin in reconstructed energy ε rec . The number of ν α per second at the detector for the interaction channel c is given by where σ cα is the cross section for ν α in the interaction channel c [42,43], φ 0 β is the atmospheric neutrino flux [39] and P ν β →να is the probability for a ν β to oscillate into a ν α when arriving at ORCA. To simulate the atmospheric neutrino signal in ORCA, we have fixed the values of the neutrino oscillation parameters to their best-fit values found in [1,2] and summarized in Tab. 1. Note that, to establish ORCA's sensitivity to the invisible neutrino decay, we set the decay parameter α 3 to zero in the simulated data. The detector resolutions in Eq. 7, R cα , include both the zenith angle and energy resolutions R cα (θ z,rec , E ν,rec , θ z,true , E ν,true ) = r cα Eν (E ν,rec , E ν,true ) r cα θz (θ z,rec , θ z,true , E ν,true ), (9) where the dimensions of the individual resolutions are given by r cα Eν = GeV −1 and r cα θz = rad −1 . For simplicity, here we have neglected the dependence of the energy resolution on the arrival direction, that is actually very small [11].  Finally, with all these ingredients, one can calculate the total number of events for a given topology (track-like or shower-like) expected in ORCA. The total event number is obtained by convolving each interaction channel with the corresponding particle identification performance: T cα pid for tracks and S cα pid for showers. Each of them represents the probability that ORCA identifies an event produced from a ν α via an interaction channel c with the given topology. Thus, the number of events inside the ij-bin identified with a topology T is with T cα pid = T cα pid , S cα pid . The detector-dependent quantities (M cα eff , R cα and T cα pid ) corresponding to ORCA have been obtained by fitting the information provided in [11] for a configuration of 9 m spacing between DOMs in a line, chosen as the final experimental setup by the KM3NeT collaboration.

The role of systematics
We estimate ORCA's sensitivity to the invisible neutrino decay defining a χ 2 function in terms of the most relevant parameters: the decay constant α 3 and the atmospheric parameters, ∆m 2 31 and sin 2 θ 23 . With this χ 2 , we fit the simulated signal in ORCA, obtained as explained in the previous subsection. Besides the details commented there, one should also consider the presence of systematic uncertainties that may affect the simulation of the experiment. For instance, the uncertainties on the determination of the atmospheric neutrino flux or the limited knowledge of the detector response will certainly modify the calculation of the expected number of events in a given experiment. These systematic uncertainties are usually included in the Systematic Expectation value (µ k ) numerical analysis with some nuisance parameters i , implemented in the definition of the χ 2 function Here i (j) indicates the ith (jth) bin in azimuth angle (energy), N dat ij is the simulated number of events in this bin (in analogy with the observed one for a running experiment) and N ij (sin 2 θ 23 , ∆m 2 31 , α 3 ; ) is the event number for the oscillation parameters sin 2 θ 23 and ∆m 2 31 , the decay constant α 3 and the nuisance parameters = ( 1 , 2 , . . .). These parameters are fitted after the minimization of the χ 2 function, that includes a pull factor penalizing large deviations from their corresponding expectation values, µ k , compared to the associated errors, σ k . The systematics used in our analysis are listed in Tab. 2, together with their corresponding expectation values and errors. As explained above, we consider systematic uncertainties related to the detector functionality and to the theoretical predictions of the atmospheric neutrino fluxes. The first four entries in Tab. 2 are related to the atmospheric neutrino flux as, for instance, a global normalization uncertainty, N . We also consider γ, the spectral index (with a pivot point in neutrino energy at 20 GeV), and two systematics related to the composition of the atmospheric flux, namely the fraction of muon neutrinos to electron neutrinos φ νµ /φ νe in the original fluxes: f µe , and another one regarding the fraction of neutrinos to antineutrinos φ ν /φν: f ν/ν . These last two modify the flux composition, such that φ →φ. In combination with the other systematics, the flux is changed as The rest of the systematics (f r,shower , f r,track and f E ) are associated to detectorrelated effects that might affect, respectively, the resolution of shower and track events, as well as the reconstructed energy. The first two simply modify the detector response function in Eq. (9), discussed in the last subsection, while f E modifies the incoming (or true) neutrino energy. This systematic uncertainty affects the event number calculation in a more complicated way, since it replaces E true by E true (1+f E ) everywhere in the simulation and, therefore, modifies directly the oscillation probabilities. Note that all the systematics named with an f account for small deviations from their corresponding central value. Therefore, their effect can be accounted for by modifying the corresponding variable X to which they affect, such that X → X(1 + f X ).

Results and discussion
In this section, we present the results of our analysis, assuming that the ORCA experiment is running for three or ten years. First, we present the bounds it could put on the invisible decay of ν 3 from the observation of atmospheric neutrino oscillations. Next, we show how the invisible neutrino decay can affect the determination of the atmospheric neutrino parameters when marginalizing over α 3 . The systematic errors considered in this analysis were summarized in Sec. 3.2. Since current global data prefer normal neutrino mass ordering [6,44], here we will focus only on this case.

Estimated sensitivity to the invisible decay
Our results for the sensitivity to the decay constant α 3 are presented in Fig. 3. We have followed two approaches to estimate ORCA's sensitivity to the invisible decay. First, we have fixed the oscillation parameters, ∆m 2 31 and sin 2 θ 23 , to their best fit point and we have calculated the χ 2 function following Eq. (11). The results of this analysis are reported with dashed lines in Fig. 3, where we have considered 3 and 10 years of running time. Then, we have followed a more realistic procedure, minimizing the χ 2 function defined in Eq.(11) over the oscillation parameters for each value of α 3 , The minimization is performed over all possible values of ∆m 2 31 and sin 2 θ 23 , where no priors on the parameters have been used. Our findings for this case correspond to the solid lines plotted in Fig. 3. Thanks to the expected good resolution of ORCA to ∆m 2 31 and θ 23 , the prior knowledge on the oscillation parameters affects the sensitivity to neutrino decay only for very large values of α 3 , as can be seen by comparing the solid and dashed lines in the two panels of Fig. 3. This worsening of the sensitivity comes mostly from a jump in the preferred value of θ 23 from the second octant (corresponding to the best-fit value of θ 23 , used to simulate ORCA's data) to the first one. We have tested that the change of behavior in the sensitivity curves actually disappears when a value of sin 2 θ 23 < 0.5 is chosen to simulate the fake data. The reason of this feature is somewhat subtle, but it is related to the fact that larger modifications of the oscillation probabilities due to neutrino decay are expected close to the matter effect resonance, as explained in Sec. 2. In this region, there is a large influence coming from the neutrino appearance probabilities, P νe,µ→νµ,e , mainly driven by a term proportional to sin 2 θ 23 (contrary to the muon neutrino disappearance probability, dominated by sin 2 2θ 23 ) [11]. Because a larger decay essentially implies less track-like events (ν µ contains more ν 3 than ν e does), and the θ 23 dependence makes the muon disappearance probability almost flat around sin 2 θ 23 = 0.5 when compared to the appearance probabilities, this reduction of track-like signatures due to neutrino decay can be compensated with a lower rate of flavor oscillations in the ν µ → ν e channel, achieved switching the fitted atmospheric mixing angle from the upper to the lower octant. From Fig. 3, one can read off the expected limits on the invisible neutrino decay from ORCA, summarized in Tab. 3 for 3 and 10 years of running time, too. Using the relation α 3 = m 3 /τ 3 , it is straightforward to convert the bounds on α 3 to limits on τ 3 /m 3 , in order to compare with the existing limits in the literature, discussed in Sec. 1. We find that, within the threeneutrino framework, ORCA could improve the current bounds on α 3 from oscillation experiments by approximately two orders of magnitude.   Table 3: Expected 90% C.L. limits from ORCA after 3 and 10 years of running time.

Effect of decay on standard oscillation parameters
The existence of invisible neutrino decay could affect the determination of the standard oscillation parameters in ORCA. Here we perform an analysis with the χ 2 function in Eq. (11) marginalized with respect to the decay constant α 3 , for each pair of values (sin 2 θ 23 , ∆m 2 31 ), where α 3 is varied freely in the fit. The result of our simulation is presented in Fig. 4, where one can see that neutrino decay does not affect the sensitivity of ORCA to the mass splitting, and only worsens very slightly the sensitivity to the atmospheric angle.
We have also tested the implications of neutrino decay in the ability of ORCA to determine the neutrino mass ordering. To do so, we have simulated data for a given true ordering (TO) of neutrino masses and different assumed true values of the decay constant, α true 3 : N dat ij (sin 2 θ BF 23 , ∆m 2 BF 31 , α true 3 , TO), where we have fixed the oscillation parameters to their best fit value in Tab. 1. Then, we evaluate the χ 2 function in Eq. (11) assuming the other (wrong) mass ordering (WO) and marginalize over the two oscillation parameters as well as the decay constant α 3 , with N ij (sin 2 θ 23 , ∆m 2 31 , α 3 , WO). Our results, presented in Fig. 5, show ORCA's sensitivity to the ordering of the neutrino mass spectrum when the true ordering is assumed to be normal (blue lines) and inverted (magenta lines). In the limit of very small α 3 , we recover the standard ORCA sensitivity to the mass ordering, as expected. From the figure, we can also appreciate a reduction of ORCA's mass ordering discrimination power for α 3 ∼ 10 −4 eV 2 , when true normal mass ordering is assumed. For values of α 3 larger than 10 −4 eV 2 , a general increase on the neutrino mass ordering sensitivity is obtained. Note, however, that these values of the decay constant are already excluded by current oscillation limits. On the other hand, this enhanced sensitivity to the mass ordering arises from the different oscillatory patterns between NO and IO induced by the neutrino decay, rather than by matter effects and neutrino flux differences. In any case, in the region of interest where the experiment is most sensitive, corresponding to α 3 ∼ 10 −6 eV 2 , the sensitivity to the neutrino mass ordering is robust and does not show a large dependence on the decay. Notice that the robustness of ORCA's sensitivity to the neutrino mass ordering is mostly due to the broad number of baselines and energies accessible to atmospheric neutrinos. In consequence, this result can not be automatically extended to the case of accelerator-based long-baseline oscillation experiments, where the fixed baselines and narrow energy spectrum might hinder the ability of disentangling the effect of flavor oscillations from that of neutrino decay. Dedicated analysis will be required to investigate the impact of neutrino decay on the mass ordering sensitivity of these experiments.

Other future prospects for invisible neutrino decay
The sensitivity to the invisible neutrino decay at different future experiments has been recently estimated at several studies. Here we will summarize the most relevant results we have found in the literature.
Ref. [45] analyzes the sensitivity of the future long-baseline neutrino experiment DUNE. Considering a run of 5 years in neutrino mode plus other 5 years in antineutrino mode, it is shown that a bound of τ 3 /m 3 > 4.50 × 10 −11 s/eV could be obtained at 90% C.L. for the case of normal mass ordering.  The authors of Ref. [46], on the other hand, study the sensitivity to neutrino decay focusing on reactor experiments. They find that a bound of τ 3 /m 3 > 9.1 × 10 −11 s/eV at 90% C.L. could be obtained after 5 years of run time in JUNO. Note, however, that in this work the authors performed their analysis marginalizing over 1σ ranges of neutrino oscillation parameters and, therefore, a marginalization over the 3σ ranges could worsen this bound.
Concerning atmospheric neutrinos, the sensitivity to the invisible neutrino decay at the proposed India-based Neutrino Observatory (INO) has been analyzed in Ref. [47]. After 10 years of data taking, it is found that INO could put a limit of τ 3 /m 3 > 1.51 × 10 −10 s/eV at 90% C.L.
Therefore, assuming that ORCA will run for ten years, we find that, besides improving the current bounds on the invisible neutrino decay by approximately two orders of magnitude, it will be, at least, as competitive as DUNE, JUNO and INO, if not a bit better.

Conclusions and final remarks
In this letter we have performed an analysis of the ORCA experiment in the context of invisible neutrino decay. We find that the data on atmospheric neutrinos obtained in ORCA in the GeV energy range could improve the bounds on the decay constant from current oscillation experiments, when including three-neutrino oscillations, by roughly a factor of 100. After ten years, ORCA could constrain the parameter τ 3 /m 3 > 2.5 × 10 −10 s/eV at 90% C.L. This means that, in the context of future neutrino oscillation experiments, ORCA will better constrain the invisible neutrino decay in comparison to other experiments such as DUNE, JUNO or INO, although not very significantly.
We also show that the decay of the heaviest neutrino does not affect the determination of the atmospheric oscillation parameters at ORCA. In particular, the determination of the neutrino mass ordering is very robust against the presence of invisible neutrino decay.
To conclude this work, we must remark that terrestrial experiments are not the only way to obtain bounds on invisible neutrino decay. Complementary constraints on the neutrino lifetime have been derived from astrophysical and cosmological observations which, in many cases, are significantly more stringent than the bounds presented in this work or from other current or future oscillation experiments.
Concerning cosmology, the standard Big Bang model predicts the existence of a relic background of neutrinos that, after the weak decoupling process at MeV temperatures, free stream during later epochs of the expanding universe. However, if non-standard neutrino interactions exist, depending on their strength, neutrinos could behave like a relativistic fluid, instead of free-streaming particles with anisotropic stress. These interactions could be tested via their imprint on cosmological observables (see e.g. [48][49][50][51][52][53][54][55]), such as the anisotropies of the cosmic microwave background. In particular, several cosmological analyses have shown that very stringent bounds apply to the invisible neutrino decay [49,[51][52][53]. For instance, in [53] a cosmological limit on the neutrino lifetime was found at the level of log 10 [(τ /m)/(s/eV)] 11, several orders of magnitude above the values we are considering. However, it is expected that these bounds could be relaxed if only one of the neutrino states decays (and the others free stream like in the standard case) or other relativistic non-interacting particles are present, such as thermal axions.
A future detection of the neutrino mass with cosmological precision data, at the level guaranteed by flavor oscillations, would lead to a huge improvement on the existing limits on the neutrino lifetime by many orders of magnitude. This conclusion applies to the invisible process, and to any kind of neutrino decay, because it is independent of the neutrino decay products [56]. Likewise, lifetimes up to a similar level would be tested if the diffuse neutrino flux emitted by all past supernovae is measured in the near future [57][58][59].