Cosmological Limits on the Neutrino Mass and Lifetime

At present, the strongest upper limit on $\sum m_{\nu}$, the sum of neutrino masses, is from cosmological measurements. However, this bound assumes that the neutrinos are stable on cosmological timescales, and is not valid if the neutrino lifetime is less than the age of the universe. In this paper, we explore the cosmological signals of theories in which the neutrinos decay into invisible dark radiation on timescales of order the age of the universe, and determine the bound on the sum of neutrino masses in this scenario. We focus on the case in which the neutrinos decay after becoming non-relativistic. We derive the Boltzmann equations that govern the cosmological evolution of density perturbations in the case of unstable neutrinos, and solve them numerically to determine the effects on the matter power spectrum and lensing of the cosmic microwave background. We find that the results admit a simple analytic understanding. We then use these results to perform a Monte Carlo analysis based on the current data to determine the limit on the sum of neutrino masses as a function of the neutrino lifetime. We show that in the case of decaying neutrinos, values of $\sum m_{\nu}$ as large as 0.9 eV are still allowed by the data. Our results have important implications for laboratory experiments that have been designed to detect neutrino masses, such as KATRIN and KamLAND-ZEN.


Introduction
Over the last few decades, a series of oscillation experiments have convincingly established that the neutrinos have masses, and determined their mass splittings.However, the actual values of the masses of the three neutrino species continue to remain a mystery.In particular, it is still not known whether the spectrum of neutrino masses is hierarchical, inverse hierarchical or quasi-degenerate.The question of whether the neutrino masses are Dirac or Majorana also remains unanswered.
At present, the strongest limit on the sum of neutrino masses, m ν < 0.12 eV, is from cosmological observations [1].These measurements are sensitive to the neutrino masses through the gravitational effects of the relic neutrinos left over from the Big Bang.In determining the size of this effect [2,3], (reviews with additional references may be found in [4][5][6][7]), it is assumed that the neutrinos are stable on timescales of order the age of the universe.In particular, if the neutrino lifetime is less than the age of the universe [8,9], or if the neutrinos have annihilated away into lighter states [10,11], this bound on the neutrino masses is no longer valid and must be reconsidered.In this paper, we explore the cosmological signals that arise from a general framework in which the neutrinos decay into dark radiation on timescales shorter than the age of the universe, and determine the bound on the sum of neutrino masses as a function of the neutrino lifetime in this scenario.Our focus is on the case in which neutrinos decay after becoming non-relativistic.
The case for neutrino decay is theoretically extremely well-motivated.Neutrino decay is in fact a characteristic feature of models in which neutrinos have masses.Even in the minimal extensions of the Standard Model (SM) that incorporate Dirac neutrino masses by adding right-handed neutrinos, or Majorana masses by including the non-renormalizable Weinberg operator, the heavier neutrinos undergo two-body decays at one loop into a lighter neutrino and a photon [12][13][14][15][16], (useful discussions may also be found in [17,18]).In these scenarios, the lifetime of the massive neutrino is given by τ ν ∼ 10 50 s (0.05 eV/m ν ) 5 , assuming the daughter neutrino mass is negligible.This is much longer than the age of the universe, and therefore these minimal frameworks do not give rise to observable cosmological signals from neutrino decay.However, in more general extensions of the SM that incorporate neutrino masses, the neutrino lifetime can be much shorter.In particular, this includes theories in which the generation of neutrino masses is associated with the spontaneous breaking of global symmetries in the neutrino sector [19][20][21][22][23], (see also [24,25]).In this framework, the heavier neutrinos can decay into a lighter neutrino and one of the Goldstone bosons associated with the spontaneous breaking of the global symmetry.The timescale for this process can be shorter than or comparable to the age of the universe, giving rise to cosmological signals.In general, neutrinos that are unstable on cosmological timescales remain an intriguing possibility due to the strong motivations for new physics that explains the smallness of neutrino masses.
In the past, the decaying neutrino scenario has been explored as a solution to the solar and atmospheric neutrino problems [26][27][28][29].However, the resulting predictions for the energy spectrum of the solar neutrinos and the decay lengths required for this proposal are now disfavored by the data [30][31][32].There has also been earlier work studying the impact of the decay of massive neutrinos on structure formation [33,34].However the range of parameter space that was considered is much above the current limits on the masses of the neutrinos.More recently, radiative neutrino decays have been proposed as an explanation of the 21 cm signal observed by the EDGES experiment [35].
The current limits on the neutrino lifetime are rather weak, except in the case of decays to final states involving photons.In this specific case, the absence of spectral distortions in the cosmic microwave background (CMB) places strong bounds on radiative decays from a heavier neutrino mass eigenstate to a lighter one, τ ν > ∼ 10 19 s for the larger mass splitting and τ ν > ∼ 4 × 10 21 s for the smaller one [36].There are also very strong, albeit indirect, limits on radiative neutrino decays based on the tight laboratory and astrophysical bounds on the neutrino dipole moment operators that induce this process [37][38][39][40][41].
In contrast, the decay of neutrinos into dark radiation that does not possess electromagnetic interactions is only weakly constrained by current cosmological, astrophysical, and terrestrial data.The most stringent bound on this scenario arises from CMB measurements.If neutrino decay and inverse decay processes are effective during the CMB epoch, they prevent the neutrinos from free streaming, leading to observable effects on the CMB [42][43][44].Current measurements of the CMB power spectra require neutrinos to free stream from redshifts z ≈ 8000 until recombination, z ≈ 1100 [45][46][47][48]. 1his can be used to set a lower bound on the neutrino lifetime τ ν ≥ 4 × 10 8 s (m ν /0.05 eV) 3 for SM neutrinos decaying into massless dark radiation [48].Several astrophysical observations have also been used to set limits on the neutrino lifetime.However, the resulting bounds are much weaker.The observation that the neutrinos emitted by Supernova 1987A did not decay prior to reaching the earth can be used to set a bound on the lifetime of the electron-neutrino, τ νe /m νe ≥ 5.7 × 10 5 s/eV [50].Similarly, the detection of solar neutrinos at the earth can be used to place a bound on the lifetime of the mass eigenstate ν 2 , τ ν /m ν 10 −4 s/eV [32,51,52].Limits on the neutrino lifetime can also be obtained from atmospheric neutrinos and long-baseline experiments, but the resulting constraints are even weaker (see e.g.[53][54][55][56]).Therefore, at present there is no evidence that neutrinos are stable on cosmological timescales, and that the cosmic neutrino background (CνB) has not decayed away into dark radiation.
The impact of non-vanishing neutrino masses on cosmological structure formation is well understood, (see [4,6] for useful reviews).
• Sub-eV neutrinos constitute radiation at the time of matter-radiation equality.Therefore, fluctuations about the background neutrino number density do not contribute significantly to the growth of structure until after neutrinos have become non-relativistic.Consequently, perturbations on scales that enter the horizon prior to neutrinos becoming non-relativistic evolve differently than scales that enter afterwards, thereby affecting the matter power spectrum.
• After neutrinos become non-relativistic, their overall contribution to the energy density redshifts away less slowly than that of a relativistic species of the same abundance.This results in a larger Hubble expansion, reducing the time available for structure formation.This leads to an overall suppression of large scale structure (LSS).
Then the leading effect of non-vanishing neutrino masses is to suppress the growth of structure on scales that entered the horizon prior to the neutrinos becoming non-relativistic.The extent of this suppression depends on the values of the neutrino masses.Since heavier neutrinos become nonrelativistic earlier and also contribute a greater fraction of the total energy density after becoming nonrelativistic, a larger neutrino mass results leads to more suppression of LSS.In the case of neutrinos that decay, this suppression now also depends on the neutrino lifetime.After neutrinos have decayed, their contribution to the energy density redshifts like that of massless neutrinos, resulting in a milder suppression of structure as compared to stable neutrinos of the same mass.It follows that there is a strong degeneracy between the neutrino mass and the lifetime inferred from the matter power spectrum.The cosmological upper bound on the neutrino mass is therefore lifetime-dependent, as was first discussed in [8,9].Neutrino masses also lead to observable effects on the CMB.Sub-eV neutrinos become nonrelativistic after CMB decoupling.The main "primary" effect on the CMB is through the early and late integrated-Sachs-Wolfe effects, as well as a modification of the angular diameter distance to the last scattering surface.Because of their impact on the growth of structure detailed above, neutrinos also affect the CMB through the "secondary" effect of lensing.At the precision of Planck, the effects of lensing drive the CMB constraints on the sum of neutrino masses.Since neutrino decay results in a milder suppression of structure as compared to stable neutrinos of the same mass, the bounds on m ν from CMB lensing are also lifetime dependent.We begin our analysis by deriving the Boltzmann equations that govern the cosmological evolution of density perturbations in the case of unstable neutrinos.We then appropriately modify the Boltzmann code CLASS 2 [57] to calculate the CMB and matter power spectra to accommodate this framework.We find that the results admit a simple analytic understanding.We then perform a Monte Carlo analysis based on CMB and LSS data (Planck+BAO+Pantheon+LSS) to determine the bounds on this scenario.We use the likelihood function from the Planck 2015 analysis [58]. 3 We find that when the stable neutrino assumption is relaxed, the limits on the neutrino masses from this data set become much weaker, with the bound on m ν increasing from 0.25 eV to 0.9 eV.Importantly, this shows that the cosmological bounds do not exclude the region of parameter space in which future experiments such as KATRIN [60], KamLAND-ZEN (KLZ) [61] and the Enriched Xenon Observatory (EXO) [62,63] are sensitive to the neutrino masses.
Our focus in this paper is on the decay of neutrinos to dark radiation, since this framework has a greater impact on the bound on m ν than the decay of heavier neutrinos to lighter ones.In particular, at present the cosmological limits on m ν only constrain quasi-degenerate neutrino spectra, so that decays of heavier neutrinos to lighter ones are not expected to alter the current bound significantly.In appendix A we present an example of a simple model in which the neutrinos decay into dark radiation on timescales of order the age of the universe.This model is consistent with all current cosmological, astrophysical and laboratory bounds, and represents a concrete realization of the scenario we are considering.However, we stress that the results presented in the body of the paper are not restricted to this specific model, but apply to any theory in which the neutrinos decay to dark radiation after becoming nonrelativistic.
The outline of this paper is as follows.In the next section we discuss the parameter space of the neutrino mass and lifetime, outlining the current bounds.In Sec. 3, we derive the Boltzmann equations that dictate the cosmological evolution of perturbations in the phase-space distribution of unstable neutrinos and their daughter radiation.While our focus is on the case in which the decaying particles are neutrinos, the formalism is more general and can be applied to the much larger class of models in which warm dark matter decays into dark radiation.In Sec.4.1, we numerically compute the growth of perturbations in the case of unstable neutrinos, and determine the effects on the matter power spectrum and on CMB lensing.To obtain a physical understanding, in Sec.4.2 we derive analytical expressions for these effects.In Sec. 5, we perform a Monte Carlo scan of the parameter space and derive constraints on the mass and lifetime of the neutrino from current data.Our conclusions are in Sec. 6.In appendix A, we present a realistic example of a model in which the neutrinos decay into dark radiation on timescales of order the age of the universe.The plot shows the current constraints in the m ν − Γ ν parameter space.The colored regions are excluded by current data while the white region is allowed.The orange dashed line separates the region of parameter space in which neutrinos decay while still relativistic from that in which they decay after becoming non-relativistic.Our study focuses on the region below this line, corresponding to the latter scenario.The light grey regions show current constraints on neutrino mass and lifetime coming from CMB free streaming and the bound on stable neutrinos (labelled "CMB+LSS (stable neutrino)").Our analysis excludes the blue region labelled "CMB+LSS (this work)" based on CMB and LSS data (Planck+BAO+Pantheon+LSS).The dashdotted line represents the approximate constraint obtained by simply requiring that the matter power spectrum be consistent with observations in the neighborhood of k = 0.1 h/Mpc with fixed H 0 .This is seen to provide a reasonable estimate to the constraints from all data.The vertical brown band shows the projected KATRIN sensitivity and also the current KLZ sensitivity.The vertical red line shows the projected KLZ-800 sensitivity in the case of a normal hierarchy.

Parameter Space of the Unstable Neutrino
In this section we outline the constraints on the decay of neutrinos to dark radiation.As explained in the introduction, these bounds only place limits on a combination of the neutrino mass and the lifetime.Therefore, in this study we will map out the constraints and the signals in the two-dimensional parameter space spanned by the sum of neutrino masses ( m ν ) and the neutrino decay width (Γ ν ), as displayed in Fig. 1.In our analysis we make the simplifying assumption that all three neutrinos are degenerate in mass.As we shall see, the bounds on m ν are always much larger than the observed mass splittings, and so this is an excellent approximation in the relevant parameter space.We further assume that all three neutrinos have the same decay width Γ ν .Since the mixing angles in the neutrino sector are large, this is a good approximation in many simple models of decaying neutrinos if the spectrum of neutrinos is quasi-degenerate.In particular, the model presented in appendix A exhibits this feature.
There is a hard lower limit on the sum of neutrino masses from the atmospheric and solar mass splittings which constrain m ν ≥ ∆m 2 31 + ∆m 2 21 = 0.06 eV in the case of normal ordering and m ν ≥ 2 × ∆m 2 31 = 0.1 eV in the case of inverted ordering [7].Therefore, we present the parameter space starting from m ν = 0.06 eV.CMB observations can be used to obtain an upper bound on the sum of neutrino masses.The current CMB data constrains the effective number of neutrinos, N eff , during the epoch of acoustic oscillations to be 2.99 ± 0.17 [1], which is perfectly compatible with the SM value of 3.046.Then, if neutrinos are stable on CMB timescales, we can obtain an approximate upper bound on their masses by requiring that all three species of neutrinos are relativistic at recombination.This translates into an approximate limit, m ν 3T rec ≈ 0.9 eV.A more precise bound can be obtained from a fit to the CMB data.
The CMB can also be used to constrain the masses of neutrinos that decay prior to recombination.As mentioned in the introduction, CMB data requires the species that constitute N eff to be free streaming at redshifts below z ≈ 8000 until recombination, z ≈ 1100.This can be used to place limits on processes such as neutrino decays and inverse decays that prevent neutrinos from free streaming at late times.The resulting bound depends on the neutrino mass, and is given by τ ν ≥ 4 × 10 8 s (m ν /0.05 eV) 3 [48].This bound excludes the grey region at the top of Fig. 1.Naively, one might expect the CMB bounds from free streaming to rule out all theories in which the neutrino decays before recombination, independent of the neutrino mass.However, in the case of an ultrarelativistic mother particle, the decay process results in approximately collinear daughter particles moving in the same direction as the mother.Similarly the inverse decay process generally only involves collinear initial state particles, so that there is no significant disruption in the flow of energy even if the decay and inverse decay processes are efficient [45].The net constraint from CMB free streaming is therefore much weaker on the decays of light neutrinos.
As discussed in the introduction, massive neutrinos suppress the growth of matter perturbations by reducing the time available for structure formation.In the case of stable neutrinos, this has been used to set a constraint on the sum of neutrino masses, m ν ≤ 0.12 eV [1].Unstable neutrinos that decay after becoming non-relativistic also lead to a suppression in the growth of structure that now depends on the neutrino lifetime.In this paper we determine the resulting bound in the two dimensional parameter space spanned by m ν and the neutrino lifetime.Based on the Monte Carlo study presented in Sec. 5, CMB and LSS data (Planck+BAO+Pantheon+LSS) exclude the blue region labelled as "CMB+LSS (this work)" in Fig. 1.We have scanned the region between 0 ≤ log 10 Γν km/s/Mpc ≤ 5.5.In Fig. 1, we simply extrapolate the bound at log 10 Γν km/s/Mpc = 0 to Γ ν = 0, because the constraint on m ν is independent of Γ ν when Γ ν H 0 .The existing constraint on the masses of stable neutrinos from this data set forms the lower boundary of this region (labelled as "CMB+LSS (stable neutrino)").
The dash-dotted line that approximately envelopes the blue shaded region represents the con-straint obtained by simply requiring that the matter power spectrum be consistent with observations in the neighborhood of k = 0.1 h/Mpc with fixed H 0 , where the current LSS measurements have the best sensitivity.We see that it provides a good approximation to the true bound, except in the region of m ν 0.9 eV, where the CMB limits on N eff at recombination become important.The impact of neutrinos on the matter power spectrum depends slightly on the mass ordering as the individual mass eigenstates become non-relativistic at different times.However, since the current limits are only sensitive to quasi-degenerate spectra, we are justified in neglecting this effect.
The orange dashed line (Γ = H(z nr )) separates the region where neutrinos decay when nonrelativistic from the region where they decay while still relativistic.Here z nr , the approximate redshift at which neutrinos become non-relativistic, is defined implicitly from the relation 3T ν (z nr ) = m ν .This definition is based on the fact that for relativistic neutrinos at temperature T ν , the average energy per neutrino is approximately 3T ν .The Hubble scale at z nr is given by, Since our study assumes neutrinos decay after they become non-relativistic, we only present the constraints below this orange dashed line.
The currently allowed parameter space is represented by the white regions in Fig. 1.In the white region above the orange dashed line, even though neutrinos decay when still relativistic, their small mass allows them to evade the current CMB free streaming constraints.In this scenario their contribution to the energy density evolves in a manner similar to that of massless neutrinos, and so the effects on LSS are similar in the two cases.In the white region below the orange dashed line the neutrinos decay after becoming non-relativistic, but because their masses are too small or their lifetimes too short, the suppression of the matter power spectrum is too small to be detected with current data.
We see from this discussion that the unstable neutrino paradigm greatly expands the range of neutrino masses allowed by current data.This has important implications for current and future laboratory experiments designed to detect neutrino masses.Next generation tritium decay experiments such as KATRIN [60] are expected to be sensitive to values of m νe as low as 0.2 eV, corresponding to m ν of order 0.6 eV.A signal in these experiments would conflict with the current cosmological bound, m ν < 0.12 eV, for stable neutrinos.However, in the decaying neutrino paradigm, we have seen that the current cosmological upper bound on the sum of neutrino masses is relaxed, with the result that m ν as high as 0.9 eV is still allowed.Therefore, a signal at KATRIN can be accommodated if neutrinos are unstable on cosmological timescales.In Fig. 1, we display a brown vertical line m ν ≈ 0.6 eV that corresponds to the expected KATRIN sensitivity.
In the case of Majorana neutrinos, current data from neutrinoless double-beta decay experiments such as KLZ and EXO have already ruled out m ν 0.6 eV (brown vertical line) [61,63].An updated version of KLZ, the KLZ-800, is currently probing m ν as low as 0.17 eV [64] (red vertical line) in the case of the normal hierarchy and the entire parameter space for the inverted hierarchy.If this experiment were to see a signal, we cannot immediately conclude that hierarchy is inverted based on the current cosmological bound of m ν < 0.12 eV, since the decaying neutrino paradigm would still admit a normal hierarchy.

Evolution of Perturbations in the Decay of Non-Relativistic Particles into Radiation
In this section we derive the set of Boltzmann equations describing the evolution of the phase-space density of massive particles decaying into massless daughter particles, working to first order in the perturbations.In contrast to the case of cold dark matter (CDM) decay (see, e.g., [65,66]), we cannot assume that the mother particles are at rest, but must take into account their non-trivial momentum distribution, as in the studies [67][68][69].This allows us to study the cosmological effects of a warm particle species, such as neutrinos or warm dark matter, decaying into radiation.We implement these new Boltzmann equations into the numerical code CLASS to generate the results in sections 4.1 and 5.
The phase-space distribution of a particle species in the expanding universe is a function of the position x, the comoving momentum q ≡ qn, and the comoving time τ .The evolution of this distribution is determined by the Boltzmann equation, where C[f ] is the collision term that accounts for all processes involving the species.We consider the case of a massive mother (with the subscript M for mother) of mass M decaying into N daughters (D i=1,2...N ).For the sake of simplicity, we restrict ourselves to the case where the mother particles decay after becoming non-relativistic, but nevertheless keep track of their non-trivial momentum distribution.In this regime, inverse-decay processes can be safely neglected.We also ignore any effects arising from Pauli blocking and spontaneous emission since f M,Di 1.The collision terms for the mother and daughter particles are then given by, where S ≡ (q 2 S + m 2 S a 2 ) 1/2 represents the comoving energy of the species S(≡ M, D i ) and d3 q ≡ d 3 q/(2π) 3 .From the definition of the decay width, the collision term for the mother particle can be simplified to where Γ denotes the decay width in the rest frame of the decaying particle, and the relativistic boost factor γ ≡ q 2 M + M 2 a 2 /(M a) accounts for time-dilation in the cosmic frame.To determine the evolution of inhomogeneities in our universe, we consider perturbations about the homogeneous and isotropic background phase space distribution functions, f S (q S , n, x, τ ) = f 0 S (q S , τ ) + ∆f S (q S , n, x, τ ), S = M, D i . (3.5)

Background: Zeroth Order
Treating ∆f M and fluctuations about the gravitational background as higher order perturbations, the zeroth order Boltzmann equations for f 0 M arising from Eq. (3.1) take the form, The formal solution to f 0 M (q, τ ) from the differential equations in Eq. (3.6) is given by, where τ i denotes the initial conformal time and f i (q) represents the initial momentum distribution.We will focus on the case where the mother decays after becoming non-relativistic.Using integration by parts, the exponent in Eq. (3.7) can be rewritten as, where we have used adτ = dt.It is computationally demanding to solve the integral for general a(τ ).However, the behavior of the exponential factor is rather simple: the exponential is close to 1 when τ is smaller than the mother lifetime ∼ γ/Γa, and f M no longer contributes when τ is much larger than the mother lifetime.The only time that the exponential factor exhibits a non-trivial adependence is when τ ∼ γ/Γa.Since our focus is on decays in the non-relativistic regime, so that γ(a) is slowly varying at the time of decay.Then the second term on the right-hand side of Eq. (3.8), which depends on the time derivative of γ(a), can be neglected in favor of the first term.This allows us to approximate the exponent as We have verified numerically that Eq. (3.9) is a good approximation to the full solution.Therefore, the mother distribution we use in this paper is It is worth pointing out that the mother distribution described by Eq. (3.10) is a general formula that can also be applied to the case of decaying CDM.This limiting case corresponds to the distribution f i (q M ) = δ(q M )N M i /(4πq 2 M ), where N M i represents the initial comoving number density of mother particles.Since this distribution is localized entirely at q M = 0, the boost factor γ(a) = 1.
Then Eq. (3.10) reduces to the known result for decaying cold dark matter [70][71][72].Our analysis is, however, more general, because it accounts for the fact that the contribution of warm dark matter to the background energy density scales with the redshift in a more complicated manner than a −3 .In addition, it takes into account the fact that, in general, particles with larger momenta live longer as a consequence of time dilation.
We now apply the above general formula to the decay of massive neutrinos.The SM neutrinos decoupled from the photon bath when they were ultra-relativistic.Therefore, their distribution prior to decay is of the Fermi-Dirac form.Therefore, f i = 1/(e q/T ν0 + 1), leading to The collision terms for the daughter particles are more challenging.However, we can simplify this set of equations by using the total integrated Boltzmann equations for the daughters.This is done by integrating the Boltzmann equations for the individual daughter species with respect to d3 q Di Di and adding them up.The resulting total integrated collision term for the daughter species is given by, The simplification in the last line follows from the covariant conservation of the energy-momentum tensor, where we have used Eq.(3.2), Eq. (3.4), and M /γ = M a to obtain this expression.In this work we focus on the case in which the mother neutrino decays into massless daughter particles.The relation in Eq. (3.12) can be used to express the Boltzmann equation for the daughters in terms of the total comoving energy density of the daughters E D and the comoving number density of the mother Since the daughter particles constitute massless radiation, we can rewrite the expression for the evolution of the daughter distribution in Eq.The right-hand side of the Eq.(3.14) is exactly the same as in the case of cold dark matter decay.While mother particles that have higher momentum have more energy, they also decay more slowly due to time-dilation in the cosmic frame.This perfect cancellation between relativistic energy and time-dilation is neatly encapsulated in the simplification M /γ = M a that was used in obtaining Eq. (3.12).

Perturbations: First Order
In the synchronous gauge, the metric perturbations can be parametrized as, where dτ = dt/a(τ ) and the indices i and j run over the three spatial coordinates, (i, j = 1, 2, 3).It is convenient to work in Fourier space, where k is conjugate to x and k is the unit vector.In Fourier space the first order terms in Eq. (3.1) for the mother particle can be collected as, where µ ≡ k • n and P l (µ) are the Legendre polynomials.
As usual, we can expand the angular dependence of the perturbations as a series in Legendre polynomials, Here X represents any of the perturbations ∆f M,D j , ∆E D or ∆N M , which are defined as in Eqs.(3.5) and (3.13).Exploiting the orthonormality of the Legendre polynomials, we arrive at a Boltzmann hierarchy of moments in which any moment is related only to its neighboring moments.The diminishing importance of the higher moments allows us to cutoff the calculation at some l = l max , where the choice of l max depends on our desired level of accuracy.We use the improved truncation scheme from Ref. [73], which has been generalized to spatial curvature in Ref. [74].
The Boltzmann hierarchy for the perturbations of the mother particle becomes, In the limit that the decay term is set to zero, these equations reduce to the standard equations for massive neutrinos in the synchronous gauge [73], as expected.
For the Boltzmann hierarchy of daughter particles, we integrate with respect to d3 q D j q D j P l (µ D j ) on both sides of Eq. (3.1) for each daughter particle and add them up.The collision term becomes Again, our focus is on the case in which the mother particle decays after becoming non-relativistic.Then, up to corrections of order q M /(M a) arising from the motion of the mother particle, the decay into daughters is isotropic, so that there is no correlation between the directions of the mother and daughter momenta (n M,D ).Given that the perturbations of the daughter particles give only a small contribution to structure formation, we can ignore this subleading correction in q M /(M a) and assume that µ M and µ D j are uncorrelated.In this case, the angular integrals over the Legendre polynomials can be performed independently, so that l This implies that in the daughter equations, only the zeroth moment of the source term from the mother particle decay (∆f M (0) ) survives in the limit of non-relativistic decay.The source term shows up in the equation for ∆f D(0) .We can therefore take l = l = 0 and simplify the collision term to get a source term similar to that in Eq. (3.12), but with f 0 M replaced by the perturbation ∆f M (0) .Therefore, the Boltzmann hierarchy for the daughter energy perturbations, ∆E D(l) , in terms of the ∆N M (l) and the metric perturbations h and η is given by, Similar equations can also be found in [65,66,68,71].Again, we neglect the source terms with ∆N M (l>0) due to the additional q M /(M a) suppressions in these terms.Other quantities such as the overdensity, perturbed pressure, energy flux/velocity-divergence, and shear stress can be calculated from these moments in the usual manner, to be fed into the perturbed Einstein field equations as detailed in [73].

Cosmological Signals of Neutrino Decay
In this section we determine the impact of decaying neutrinos on the matter power spectrum and on CMB lensing.In Sec.4.1, we solve the Boltzmann equations of the previous section numerically using CLASS, and determine the matter power spectrum and the CMB lensing potential C φφ as a function of the neutrino mass and lifetime.This allows us to establish numerically that there is indeed a degeneracy in the matter power spectrum between neutrino mass and lifetime.In Sec.4.2, we determine the matter power spectrum analytically, after making certain well-motivated approximations.
We show that the results closely reproduce those based on the numerical study, and admit a physical interpretation of the effects of decaying neutrinos.

Numerical Results
To simplify the analysis, we assume that the three neutrinos have degenerate masses and lifetimes.This extends the parameter space of the ΛCDM model to include two additional parameters; the sum of neutrino masses, m ν , and the logarithm of the decay width, log 10 Γ ν .In our analysis, we fix the following cosmological parameters to their central values from the Planck 2015 TT, TE, EE+low-P data:{ω b = 0.022032, ω cdm = 0.12038, ln(10 10 A s ) = 3.052, n s = 0.96229, τ reio = 0.0648}.The impact of neutrino masses on the matter power spectrum looks different depending on whether θ s or H 0 is kept fixed [75].This is because, to keep θ s fixed, H 0 must be adjusted within CLASS, leading to an overall shift of the matter power spectrum.While fixing H 0 is more conventional, fixing θ s gives a better reflection of the constraining effects of a combined analysis of CMB+LSS data, since CMB data pins θ s down very precisely.In the following, we will show results with either H 0 = 67.56km/s/Mpc or 100 × θ s = 1.043, explicitly stating in each case what convention is chosen.Since the galaxy power spectrum is known to trace the CDM and baryon overdensities, we focus on the power spectrum where ρcb (δρ cb ) is the average (perturbation) of the sum of CDM and baryon energy densities. 4In Fig. 2, we display the residuals of P cb (left) and the CMB lensing potential (right) with respect to the case of massless neutrinos for m ν fixed at 0.25 eV, keeping the value of H 0 fixed.We compare three different values of Γ ν and the limiting case of stable neutrinos.The curves run from top to bottom in order of decreasing Γ ν .The analytic results are shown as dashed lines in the plot, and are seen to agree reasonably well at large k or with the numerical results, shown as solid lines.These plots demonstrate that the main effect of a non-zero decay rate of neutrinos is to reduce the power suppression at large k arising from their mass.Moreover, they establish that the gravitational effects of unstable relic neutrinos can indeed give rise to observable signals in LSS, provided that the decays occur sufficiently long after the neutrinos have become non-relativistic.
Because of the effects of nonlinearities at large k (small scales) and cosmic variance at small k (large scales), current experiments are sensitive only to a narrow range of k in the neighborhood of 0.1h/Mpc.We see from Fig. 2 that in this region there are no qualitative features in P cb | z=0 or C φφ that would allow unstable neutrinos to be distinguished from stable ones.Although P cb | z=0 and C φφ are more suppressed in the stable case, as expected, this effect can be mimicked if the the neutrino masses in the unstable scenario are suitably heavier.This results in a strong parameter degeneracy between the neutrino lifetime and the sum of neutrino masses as determined from P cb | z=0 and C φφ .
In Fig. 3 we show an explicit example of the degeneracy between mass and lifetime in the values of P cb and C φφ at fixed H 0 .We consider a model with stable neutrinos of mass m ν = 0.2 eV, and a different model with unstable neutrinos of mass m ν = 0.36 eV and width Γ ν = 10 4 km/s/Mpc.In the P cb (z = 0) case, we see from the figure that the blue (stable neutrino) and purple (unstable neutrino) curves cannot be distinguished by measurements such as SDSS DR7 (used later in sec.5), whose sensitivity is shown in grey.However, we note that the lensing power spectrum can potentially help in breaking the degeneracy, because it receives its dominant contribution at higher z ≈ 3 [76].
We will explore the possibility of breaking the degeneracy by using next generation measurements at different redshifts in future work.Finally, we show in Figs. 4 the effects of neutrino masses and decay at fixed θ s on P cb , C φφ and C T T,EE .This fixes the peak locations in the CMB power spectra and only generates negligible deviations away from the massless neutrino case in C T T,EE [75].The same choices of parameters, however, do generate sizeable deviations in C φφ and P cb away from the massless neutrino case that are close to the current sensitivities.This demonstrates that as expected, for sub-eV m ν , it is the CMB- lensing and matter power spectrum measurements that provide the constraining power.Additionally, note that the change in H 0 required to keep θ s fixed leads to an overall shift of P cb .This makes the BAO in the three models out of phase and leads to small oscillations at large k on top of the power suppression.

Analytic Understanding
In this section we provide an analytic derivation of the effects of neutrino decay on CMB and LSS observables.We begin by showing how the results in the literature for the effects of massive neutrinos on the matter power spectrum (P cb (k)) and CMB lensing (C φφ ) can be reproduced analytically.We improve on the existing analytical treatment of the cosmological effects of massive neutrinos by taking into account their momentum distribution.We then build on this to derive an expression for the evolution of overdensities in scenarios with unstable neutrinos.Once neutrinos become non-relativistic, their contribution to the background energy density leads to an increase the Hubble rate, leaving less time for structure formation as compared to a universe with massless neutrinos.The net result is an overall suppression of power at small scales in the matter power spectrum.The size of this effect can be determined by studying the evolution of density perturbations.Consider δ i = δρ i /ρ i for particle species i, for a mode that is already deep inside the horizon when neutrinos become non-relativistic at z = z nr .In the matter dominated era, the Einstein equation for the density perturbation with wavenumber k can be approximated as Here φ is the metric perturbation in the conformal Newtonian gauge [73]. 5We assume baryons have already decoupled from photons.This allows us to combine the baryon contribution to the matter density with that of CDM to simplify the discussion.Since δ ν δ cb for perturbation modes that enter the horizon before z nr , we can write, where τ is the comoving time and ρtot ≡ ρcb + ρν .Inserting this expression into the Boltzmann equation for CDM perturbations yields, where the dots represent derivatives with respect to τ .Deep in the matter dominated era, neutrinos only contribute up to a few percent of the total energy density.Therefore, throughout this derivation, we work to leading order in f ν ( 1).We look for a solution of the form, where now the function h(τ ) is to be determined.Inserting this expression into Eq.(4.4) and dropping the term proportional to f 2 ν , we obtain the following differential equation for h(τ ).
Thus far we have not made any assumption about the redshift dependence of f ν .For massless or ultrarelativistic neutrinos in the matter dominated era, we have In this case we can solve for the function h(τ ) as, This leads to the following approximate solution for perturbations in the case of massless or ultrarelativistic neutrinos, (4.9) In the limit that neutrinos are non-relativistic, f ν (τ ) goes to a constant value.Then Eq. (4.6) admits a solution where h(τ ) is constant.This implies that in the case of massive neutrinos, the h-function can be approximated as The result is almost identical to Eq. (4.8) since in both cases the exponent is dominated by f mν ν (τ i ) (and f mν ν after τ > τ nr is much smaller than the expansion parameter f mν ν ) This means that the solution in the case of massive neutrinos can be approximated as, Then the ratio of the perturbations in the two cases is given by,

.13)
This ratio can be expressed in terms of the scale factor as, where ρν (a) ≡ ρν,mν (a) − ρν, mν (a) represents the difference in the neutrino energy between the two scenarios.If all the neutrinos are stable and become non-relativistic instantly at a i , ρν (a)/ρ tot (a) = ρν,mν /ρ tot is a constant, and Eq.(4.14) recovers the well-known result for the ratio of perturbations in the massive and massless neutrino scenarios, We can improve on this estimate by incorporating a more precise expression for the neutrino energy in Eq. (4.14), Here q = a p ν denotes the neutrino's conformal momentum, and f (q) = e q/T ν0 + 1 −1 represents the momentum distribution of neutrinos.ρν (a) exhibits non-trivial redshift dependence since the neutrino energy goes from being radiation-like to being matter-like.In fig. 5, we show the evolution of the ratio in Eq. (4.14) as a function of redshift (black dashed curves) for two different values of the neutrino mass.We start our approximation from a i = 2 × 10 −3 to make sure we are deep inside the matter dominated era so that the assumptions leading to Eq. (4.14) are justified.We stress, however, that the result is quite insensitive to order one changes in a i .As we can see, Eq. (4.14) is a good approximation to the full numerical results (black solid curves), and describes the evolution of the δ cb ratio from the relativistic to the non-relativistic regime much better than the approximation based on Eq. (4.15) (black dotted curves).Using this, we can estimate the ratio of the power spectrum between the two scenarios, (4.17) The density perturbation grows much slower in the cosmological constant dominant era, and we take the final scale factor to be at a f = 0.7 for a good approximation to the power spectrum ratio today 6 .We now turn our attention to the effects of massive neutrinos on CMB lensing.The difference in the density perturbation δρ cb between the massive and massless neutrino scenarios results in a change in the gravity perturbation φ.The photons are therefore deflected differently in the CMB lensing process.The correlation function of the lensing potential, C φφ ∼ φφ , parameterizes the size of angular deflection of CMB photons.The ratio of C φφ in the massive neutrino case to that in the massless case can be approximated using Limber's formula [77,78] Here τ * ≈ 2.8 × 10 2 Mpc is the conformal time at last scattering, while τ f ≈ 1.4 × 10 4 Mpc is the conformal time today.The value of τ f differs a bit between the massive and massless neutrino scenarios, since the contribution of neutrinos to the total energy density is different in the two cases.However, since the neutrino mass only results in a significant difference in the contributions to the background energy in the short period of time between the neutrinos becoming non-relativistic and the universe becoming dominated by the cosmological constant, the difference in χ * between the two scenarios can be neglected.Then, the difference between C φφ in the two cases primarily arises from differences in the evolution of φ.
According to the Einstein Eq. (4.2), the ratio of φ between the two scenarios for large modes at a given value of the scale factor is, q q j B i J I o S f 0 g l 6 d R + f Z e X P e Z 6 0 Z Z z 5 z i H 7 B + f g G Z 9 e S c Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D F / 0 r T P c x A W R y 9 p A 5 b B J 9 q q j B i J I o S f 0 g l 6 d R + f Z e X P e Z 6 0 Z Z z 5 z i H 7 B + f g G Z 9 e S c Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D F / 0 r T P c x A W R y 9 p A 5 b B J 9 q q j B i J I o S f 0 g l 6 d R + f Z e X P e Z 6 0 Z Z z 5 z i H 7 B + f g G Z 9 e S c Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D F / 0 r T P c x A W R y 9 p A 5 b B J 9 < l a t e x i t s h a 1 _ b a s e 6 4 = " D F / 0 r T P c x A W R y 9 p A 5 b B J 9 L P d Y n M = " > A A A B 9 X i c b V B N S w M x E M 3 W r 1 q / q h 6 9 B I v g q e y 2 B b 0 q q j B i J I o S f 0 g l 6 d R + f Z e X P e Z 6 0 Z Z z 5 z i H 7 B + f g G Z 9 e S c Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D F / 0 r T P c x A W R y 9 p A 5 b B J 9 L P d Y n M = " > A A A B 9 X i c b V B N S w M x E M 3 W r 1 q / q h 6 9 B I v g q e y 2 B b 0 q q j B i J I o S f 0 g l 6 d R + f Z e X P e Z 6 0 Z Z z 5 z i H 7 B + f g G Z 9 e S c Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D F / 0 r T P c x A W R y 9 p A 5 b B J 9 L P d Y n M = " > A A A B 9 X i c b V B N S w M x E M 3 W r 1 q / q h 6 9 B I v g q e y 2 B b 0 W g p X P H I J f s D 6 + A U q C l x s = < / l a t e x i t > Dashed: our approx.Dotted: naive approx.
. Evolution of the ratio of the CDM+baryon density perturbation with respect to the case of a massless neutrino, δ mν cb /δ mν cb .The results are shown for the case of a single massive neutrino with m ν = 60 meV.All the solid curves are obtained from numerical calculations using the modified CLASS code discussed in Sec. 3. The black curve is for the stable neutrino scenario, and the blue (orange) curve is for the neutrino with decay rate Γ ν = 10 4 (10 3 ) km /s/Mpc.The dashed curves represent the first approximations to the solid curves, based on the derivation in Eq. (4.14).The dotted curves are based on the approximation method in Eq. (4.15), where we assume a i to be the value when 80% of neutrinos have their momenta lower than m ν and a f = a dec .As we see, Eq. (4.14) provides a much better approximation to the full numerical result.
Since C φφ receives its dominant contribution close to z ≈ 3 [76], we can estimate the ratio of the C φφ as, C φφ Based on a very similar analysis, we can predict the suppression of P cb (k) and C φφ for large k and in the unstable neutrino case.We consider a scenario with a single massive neutrino species that becomes non-relativistic after last scattering and decays into dark radiation.After the decay, the energy density of the daughter particles redshifts more quickly than that of a stable neutrino of the same mass as the mother.We work in the instantaneous decay approximation and assume that all neutrinos decay at the same time, corresponding to a scale factor a dec , which is implicitly determined by the equation, The difference in energy density ρν between an unstable neutrino and a massless neutrino evolves in a more complicated way than in the case of a stable neutrino.The instantaneous decay approximation allows us to separate the evolution into two parts.On timescales shorter than the proper lifetime of the neutrino, the difference in energy density follows the equation, In the instantaneous decay approximation, the energy density in non-relativistic neutrinos is immediately transferred into radiation energy at a dec .It subsequently redshifts with an extra (a dec /a) factor as compared to a non-relativistic neutrino, so that The ratio of CDM density perturbations in the case of unstable neutrinos can be obtained by inserting the energy density ratios in Eqs.(4.22) and (4.23) into Eq.(4.14).Then the ratios of P (k) and C φφ in the limit of large k and can be obtained from Eqs. (4.17) and (4.20)In Fig. 5, we show the ratio of δ cb from the numerical calculation described in Sec. 3 for both the decaying (blue and orange) and stable (black) neutrinos.The plots are for a single massive neutrino with m ν = 60 meV (upper) and 80 meV (lower), and a decay rate Γ ν = 10 4 (10 3 ) km /s/Mpc for the blue (orange) curves.In this scenario, more than 80% of the neutrinos have momenta p ν < m ν after a > 0.012 (a > 0.0096) for m ν = 60 (80) meV neutrino.It is at this point, when most of the neutrinos have become non-relativistic, that the major suppression of δ cb begins.During this period the δ cb ratio drops with the power described in Eq. (4.15) (grey line).The blue (orange) dotted lines give the value of the δ cb -suppression if the later contributions of daughter particles to the energy density shown in Eq. (4.23) are ignored.As we see, this underestimates the suppression of δ cb , showing that the contributions of daughter particles to the energy density cannot be neglected.It is clear from the figures that Eqs.(4.22) and (4.23) provide a good description of the δ cb evolution in unstable neutrino scenarios (dashed blue and orange), both before and after neutrino decay.This shows that the effects of neutrino decay on the evolution of δ cb on these length scales primarily arise from the contributions of the unstable neutrinos and their daughter particles to the background energy density, and not from their perturbations.

Current Limits on the Neutrino Mass and Lifetime from Monte Carlo Analysis
In this section we perform a Monte Carlo analysis to determine the current bounds on the neutrino mass and lifetime.

The Data and Analysis Pipeline
Our analysis makes use of various combinations of the following datasets.
• CMB: We include Planck 2015 CMB high-TT, TE, and EE and low-TEB power spectra [79], as well as the lensing reconstruction power spectrum [80].
• BAO: We use measurements of the volume distance from 6dFGS at z = 0.106 [81] and the MGS galaxy sample of SDSS at z = 0.15 [82].We include the anisotropic measurements from the CMASS and LOWZ galaxy samples from the BOSS DR12 at z = 0.38, 0.51, and 0.61 [83].
• Growth Function: The BOSS DR12 measurements also include measurements of the growth function f , defined by where σ (vd) 8 measures the smoothed density-velocity correlation, analogous to σ 8 ≡ σ (dd) 8 that measures the smoothed density-density correlation.
• Pantheon: we use the Pantheon supernovae dataset [84], which includes measurements of the luminosity distance of 1048 SNe Ia in the redshift range 0.01 < z < 2.3.
• LSS: We use the measurement of the halo power spectrum from the Luminous Red Galaxies SDSS-DR7 [85] 7 and the tomographic weak lensing power spectrum by KiDS [86].
Our baseline analysis makes use of Planck+BAO+Growth Function+Pantheon data (i.e.data that relies on background cosmology or perturbations in the linear regime mostly).We then add LSS information to gauge the constraining power of such surveys.
Using the public code MONTEPYTHON-V38 [87,88], we run Monte Carlo Markov chain analyses using the Metropolis-Hastings algorithm assuming flat priors on all parameters.Our ΛCDM parameters are, {ω cdm , ω b , θ s , ln(10 10 A s ), n s , τ reio } , to which we add the sum of neutrino masses m ν and the logarithm of the neutrino lifetime Log 10 Γ ν .
In our analysis we assume 3 degenerate, unstable neutrino species that decay into dark radiation.
Although not detailed for brevity, there are many nuisance parameters that we analyze together with these cosmological parameters.To that end, we employ a Cholesky decomposition to handle the large  number of nuisance parameters [89], and use the default priors that are provided by MONTEPYTHON-V3.

Current Limits on the Neutrino Mass and Lifetime
In order to perform meaningful comparisons and to check the accuracy of our modified version of CLASS, we begin by running the case of stable neutrinos.Our baseline constraint on the neutrino mass, obtained with Planck+BAO+Growth Function+Pantheon, is m ν < 0.28 eV (95% C.L.).This is in good agreement with the result reported in [58].The inclusion of SDSS DR7 and KiDS improves the constraint by ∼ 10%, bringing the limit down to m ν < 0.25 eV (95% C.L.).This constraint when LSS data is included is also in good agreement with what is reported in Ref. [90].
In Fig. 6 we show the 1D and 2D marginalized posterior distribution of m ν and log 10 Γ ν for both datasets, cutting the parameter space between small decay rate log 10 Γ ν /(km/s/Mpc) ∈ [0, 3] (left panel) and large decay rate log 10 Γ ν /(km/s/Mpc) ∈ [3, 5.5] (right panel) to accelerate convergence.Strikingly, once the neutrino lifetime is let free to vary, the constraint on m ν is driven by our prior on log 10 Γ ν .We recall that this was chosen in order to ensure that neutrinos decay while non-relativistic.Interestingly, the constraint stays quite stable for log 10 Γ ν /(km/s/Mpc) < 2.5, but relaxes to m ν < 0.9 eV (with Planck+BAO+Growth Function+Pantheon) for higher values of the decay rate.We note that the limit only marginally improves with the addition of current LSS data, especially at high decay rates (right panel), for which the improvement is below numerical noise.
Our study allows us to obtain a bound on the sum of neutrino masses as a function of the neutrino lifetime.We see that m ν can be as large as 0.90 eV for neutrinos that decay close to recombination.However, given our restricted prior enforcing non-relativistic decays, our analysis does not set a true upper bound on the neutrino mass.In order to derive the true upper bound we would need to correctly incorporate relativistic decays, taking into account inverse decay processes.We refer to Ref. [45,48] for a discussion of that regime, and defer to future work a reanalysis of that region of parameter space in light of the latest Planck results.

Conclusions
The fact that the couplings of neutrinos to the other SM particles are so weak makes it extremely difficult to study their properties.Even though it has been over six decades since neutrinos were first directly observed in the laboratory, several of their fundamental properties, including their masses and lifetimes, remain to be determined.However, neutrinos are also among the most abundant particles in the universe, and their gravitational pull has effects on cosmological observables.The universe is therefore an excellent laboratory for studying the detailed properties of neutrinos.
In this paper, we have explored the cosmological signals arising from the theoretically wellmotivated scenario in which neutrinos decay into invisible dark radiation on timescales less than the age of the universe.We have studied the effects of neutrino decay on the evolution of density perturbations, both analytically and numerically, and used the results to generalize the bound on the sum of neutrino mass to the case when the lifetime of the neutrino is less than the age of the universe.We have shown that the existing mass bound from CMB and LSS measurements, which assumes that neutrinos are stable, gets weakened if neutrinos decay, so that values of m ν as large as 0.9 eV are still allowed by the data.This provides strong motivation to continue the current efforts to measure the neutrino masses directly in the lab, in spite of the limited reach of these experiments.Our analytical results show that the signals of neutrino decay in LSS and CMB-lensing primarily arise from the contributions of neutrinos and their daughters to the overall energy density, and are quite insensitive to their contributions to the fluctuations about the background.Although the bounds we obtain based on the existing data do not set independent constraints on the neutrino mass and lifetime, next generation measurements of the matter power spectrum at different redshifts will provide useful information that may help in breaking this degeneracy.We will explore this in the future work [91].
where cos θ is exactly the same as in Eq. (A.4).Now the neutrino masses are given by, Assuming that the Goldstone bosons from Σ are heavier than the massive neutrinos due to some external source of explicit breaking, the dominant decay modes of the massive neutrinos are to a massless sterile neutrino and either φ or φ .Following the discussion above, the total neutrino decay width is given by where f and f are as defined after Eq. (A.8).One characteristic feature of this model is that the widths of the neutrinos scale as the cube of their masses, Γ ν i /Γ ν j = m 3 i /m 3 j .In the case of quasidegenerate neutrinos, m 1 ≈ m 2 ≈ m 3 , it is clear that all neutrinos have almost the same total width.Assuming f = f , we find that the total width is of order H 0 for f ∼ 10 5 GeV and neutrino masses of order 0.1 eV, The parameter space of this model is constrained by astrophysical, cosmological and laboratory data.These limits are very similar to those on conventional Majoron models, and can be expressed in terms of bounds on the decay constants f and f .In the case of massless Goldstone bosons, the bounds from cosmology and astrophysics are the most severe.A strong cosmological constraint arises from requiring consistency with the observation that the cosmic neutrinos are free streaming at temperatures below an eV [45][46][47][48].Neutrino-neutrino scattering mediated by Goldstone boson exchange can prevent the neutrinos from free streaming, impacting the heights and locations of the CMB peaks.This translates into constraints on f and f of order 100 keV [92].A stronger although somewhat model-dependent constraint, f, f 100 MeV, may be obtained by requiring that the Goldstone bosons and right-handed neutrinos not contribute significantly to the energy density in radiation at the time of Big Bang nucleosynthesis (BBN), or during the CMB epoch.
The strongest astrophysical bounds arise from the effects of Goldstone bosons on supernovae.The large chemical potential for electron neutrinos inside the supernova means that these particles can now decay into final states containing a Goldstone boson and a right-handed neutrino.This has the effect of deleptonizing the core, preventing the explosion from taking place.In addition, the free streaming of Goldstone bosons out of the supernovae core can lead to overly rapid energy loss.The resulting constraints are at the level of f, f 100 keV [93][94][95][96][97].There are also bounds on the couplings of neutrinos to Goldstone bosons from laboratory experiments, such as neutrinoless double beta decay [98,99], meson decays [93,100], charged lepton decays [101] and tritium decay [102].These constraints arise from corrections to the energy spectrum of the visible final states due to Goldstone boson emission.However, in all these cases, the limits are weaker than astrophysical and cosmological bounds on massless Goldstone bosons.Clearly, our benchmark values of f, f ∼ 10 5 GeV are easily consistent with all current bounds.

Figure 1 .
Figure1.The plot shows the current constraints in the m ν − Γ ν parameter space.The colored regions are excluded by current data while the white region is allowed.The orange dashed line separates the region of parameter space in which neutrinos decay while still relativistic from that in which they decay after becoming non-relativistic.Our study focuses on the region below this line, corresponding to the latter scenario.The light grey regions show current constraints on neutrino mass and lifetime coming from CMB free streaming and the bound on stable neutrinos (labelled "CMB+LSS (stable neutrino)").Our analysis excludes the blue region labelled "CMB+LSS (this work)" based on CMB and LSS data (Planck+BAO+Pantheon+LSS).The dashdotted line represents the approximate constraint obtained by simply requiring that the matter power spectrum be consistent with observations in the neighborhood of k = 0.1 h/Mpc with fixed H 0 .This is seen to provide a reasonable estimate to the constraints from all data.The vertical brown band shows the projected KATRIN sensitivity and also the current KLZ sensitivity.The vertical red line shows the projected KLZ-800 sensitivity in the case of a normal hierarchy.
(3.1)  in terms of the background daughter energy density ρD ≡ 4πa −4 E 0 D and the background mother number density nM ≡ 4πa −3 N 0 M , where E 0 D and N 0 M are defined as in Eq. (3.13) after expanding out f M and f D i as in Eq. (3.5), ∂ ρD ∂τ + 4aH ρD = aΓM nM .(3.14)

Figure 2 .)Figure 3 .
Figure2.Plots of the fractional difference in the CDM+Baryon power spectrum P cb (left) and CMB-lensing potential C φφ (right) for various decaying (and stable) massive neutrino scenarios with respect to the case of massless neutrinos.The solid lines show the results from numerical simulations of the decaying neutrino scenario for three values of the decay width, Γ ν = 10 4.0, 3.5, 3.0 (km/s/Mpc) (top to bottom), and also the stable neutrino scenario, holding m ν = 0.25 eV and H 0 fixed.The dashed lines represent the corresponding analytic estimates from Sec. 4.2.

TTFigure 4 .
Figure 4.The fractional differences in the CMB-lensing potential C φφ (top left), CDM+Baryon power spectrum P cb (top right), C T T (bottom left), and C EE (bottom right) for an unstable (purple) and a stable (blue) neutrino scenario with respect to the case of massless neutrinos (black) at fixed θ s .The grey regions show the 1σ uncertainties from Planck and SDSS DR7 respectively.
8 L X p / B / U t t y P c / 1 z r Y L + w e j O L J g B a y C D e C B X b A P j s E p q A I C b s E 9 e A R P z p 3 z 4 D w 7 L 5 + t G W c 0 s w x + w H n 9 A F 3 v n E A = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " U 8 6 0 H I Z m 8 8 L X p / B / U t t y P c / 1 z r Y L + w e j O L J g B a y C D e C B X b A P j s E p q A I C b s E 9 e A R P z p 3 z 4 D w 7 L 5 + t G W c 0 s w x + w H n 9 A F 3 v n E A = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " U 8 6 0 H I Z m 8 8 L X p / B / U t t y P c / 1 z r Y L + w e j O L J g B a y C D e C B X b A P j s E p q A I C b s E 9 e A R P z p 3 z 4 D w 7 L 5 + t G W c 0 s w x + w H n 9 A F 3 v n E A = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " U 8 6 0 H I Z m 8 8 L X p / B / U t t y P c / 1 z r Y L + w e j O L J g B a y C D e C B X b A P j s E p q A I C b s E 9 e A R P z p 3 z 4 D w 7 L 5 + t G W c 0 s w x + w H n 9 A F 3 v n E A = < / l a t e x i t > m ⌫ = 60 meV < l a t e x i t s h a 1 _ b a s e 6 4 = " m b 2 7 F C g j f 5 + o 0 g l 6 d R + f Z e X P e Z 6 0 Z Z z 5 z i H 7 B + f g G Z 9 e S c Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D F / 0 r T P c x A W R y 9 p A 5 b B J 9 L P d Y n M = " > A A A B 9 X i c b V B N S w M x E M 3 W r 1 q / q h 6 9 B I v g q e y 2 B b 0 I R Q / 2 W M H W Q n c t s 2 n a h i b Z J c k q Z e n / 8 O J B E a / + F 2 / + G 9 M v U O u D g c d 7 M 8 z M C 2 P O t H H d L y e z s r q 2 v p H d z G 1 t 7 + z u 5 f c P m j p K F K E N E v F I t U L Q l D N J G 4 Y Z T l u x o i B C T u / C 4 d X E v 3 u g S r N I 3 p p R T A M B f c l 6 j I C x 0 r 1 / D U J A J / V l M r 6 o d f I F t + h O g S 0 p l y u l M v Y W y o I U 0 B z 1 T v 7 T 7 0Y k E V Q a w k H r t u f G J k h B G U Y 4 H e f 8 R N M Y y B D 6 t G 2 p B E F 1 k E 6 v H u M T q 3 R x L 1 K 2 p M F T 9 e d E C k L r k Q h t p w A z 0 H + 9 i f i f 1 0 5 M 7 z x I m Y w T Q y W Z L e o l H J s I T y L A X a Y o M X x k C R D F 7 K 2 Y D E A B M T a o n A 1 h 6 e V l 0 i w V P a / o 3 V Q K 1 c t 5 H F l 0 h I 7 R K f L Q G a q i G q q j B i J I o S f 0 g l 6 d R + f Z e X P e Z 6 0 Z Z z 5 z i H 7 B + f g G Z 9 e S c Q = = < / l a t e x i t > ⌫ = H < l a t e x i t s h a 1 _ b a s e 6 4 = " W u y K y D A E w J 5 U b z 7 0 2 X Z + T L B 2 g v w = " > A A A B 9 X i c d V D L S g M x F M 3 U V 6 2 v q k s 3 w S K 4 G j K 2 t e 1 C K L q w y w r 2 A Z 1 a M m l a Q 5 P M k G S U M v Q / 3 L h Q x K 3 / 4 s 6 / M X 0 I K n r g w u G c e 7 n 3 n i D i T B u E P p z U 0 v L K 6 l p 6 P b O x u b W 9 k 9 3 d a + o w V o Q 2 S M h D 1 Q 6 w p p x J 2 j D M c N q O F M U i 4 L Q V j C 6 m f u u O K s 1 C e W 3 G E e 0 K P J R s w A g 2 V r r x L 7 E Q u J f 4 M p 6 c 1 X r Z H H Ir h X y x V I b I z a N i P l + x B H m V I j q F n o t m y I E F 6 r 3 s u 9 8 P S S y o N I R j r T s e i k w 3 w c o w w u k k 4 8 e a R p i M 8 J B 2 L J V Y U N 1 N Z l d P 4 J F V + n A Q K l v S w J n 6 f S L B Q u u x C G y n w O Z W / / a m 4 l 9 e J z a D c j d h M o o N l W S + a B B z a E I 4 j Q D 2 m a L E 8 L E l m C h m b 4 X k F i t M j A 0 q Y 0 P 4 + h T + T 5 o n r u e 5 3 l Uh V z 1 f x J E G B + A Q H A M P l E A V 1 E A d N A A B C j y A J / D s 3 D u P z o v z O m 9 N O Y u Z f f A D z t s n x O i S s Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " W u y K y D A E w J 5 U b z 7 0 2 X Z + T L B 2 g v w = " > A A A B 9 X i c d V D L S g M x F M 3 U V 6 2 v q k s 3 w S K 4 G j K 2 t e 1 C K L q w y w r 2 A Z 1 a M m l a Q 5 P M k G S U M v Q / 3 L h Q x K 3 / 4 s 6 / M X 0 I K n r g w u G c e 7 n 3 n i D i T B u E P p z U 0 v L K 6 l p 6 P b O x u b W 9 k 9 3 d a + o w V o Q 2 S M h D 1 Q 6 w p p x J 2 j D M c N q O F M U i 4 L Q V j C 6 m f u u O K s 1 C e W 3 G E e 0 K P J R s w A g 2 V r r x L 7 E Q u J f 4 M p 6 c 1 X r Z H H I r h X y x V I b I z a N i P l + x B H m V I j q F n o t m y I E F 6 r 3 s u 9 8 P S S y o N I R j r T s e i k w 3 w c o w w u k k 4 8 e a R p i M 8 J B 2 L J V Y U N 1 N Z l d P 4 J F V + n A Q K l v S w J n 6 f S L B Q u u x C G y n w O Z W / / a m 4 l 9 e J z a D c j d h M o o N l W S + a B B z a E I 4 j Q D 2 m a L E 8 L E l m C h m b 4 X k F i t M j A 0 q Y 0 P 4 + h T + T 5 o n r u e 5 3 l U h V z 1 f x J E G B + A Q H A M P l E A V 1 E A d N A A B C j y A J / D s 3 D u P z o v z O m 9 N O Y u Z f f A D z t s n x O i S s Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " W u y K y D A E w J 5 U b z 7 0 2 X Z + T L B 2 g v w = " > A A A B 9 X i c d V D L S g M x F M 3 U V 6 2 v q k s 3 w S K 4 G j K 2 t e 1 C K L q w y w r 2 A Z 1 a M m l a Q 5 P M k G S U M v Q / 3 L h Q x K 3 / 4 s 6 / M X 0 I K n r g w u G c e 7 n 3 n i D i T B u E P p z U 0 v L K 6 l p 6 P b O x u b W 9 k 9 3 d a + o w V o Q 2 S M h D 1 Q 6 w p p x J 2 j D M c N q O F M U i 4 L Q V j C 6 m f u u O K s 1 C e W 3 G E e 0 K P J R s w A g 2 V r r x L 7 E Q u J f 4 M p 6 c 1 X r Z H H Ir h X y x V I b I z a N i P l + x B H m V I j q F n o t m y I E F 6 r 3 s u 9 8 P S S y o N I R j r T s e i k w 3 w c o w w u k k 4 8 e a R p i M 8 J B 2 L J V Y U N 1 N Z l d P 4 J F V + n A Q K l v S w J n 6 f S L B Q u u x C G y n w O Z W / / a m 4 l 9 e J z a D c j d h M o o N l W S + a B B z a E I 4 j Q D 2 m a L E 8 L E l m C h m b 4 X k F i t M j A 0 q Y 0 P 4 + h T + T 5 o n r u e 5 3 l U h V z 1 f x J E G B + A Q H A M P l E A V 1 E A d N A A B C j y A J / D s 3 D u P z o v z O m 9 N O Y u Z f f A D z t s n x O i S s Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " W u y K y D A E w J 5 U b z 7 0 2 X Z + T L B 2 g v w = " > A A A B 9 X i c d V D L S g M x F M 3 U V 6 2 v q k s 3 w S K 4 G j K 2 t e 1 C K L q w y w r 2 A Z 1 a M m l a Q 5 P M k G S U M v Q / 3 L h Q x K 3 / 4 s 6 / M X 0 I K n r g w u G c e 7 n 3 n i D i T B u E P p z U 0 v L K 6 l p 6 P b O x u b W 9 k 9 3 d a + o w V o Q 2 S M h D 1 Q 6 w p p x J 2 j D M c N q O F M U i 4 L Q V j C 6 m f u u O K s 1 C e W 3 G E e 0 K P J R s w A g 2 V r r x L 7 E Q u J f 4 M p 6 c 1 X r Z H H I r h X y x V I b I z a N i P l + x B H m V I j q F n o t m y I E F 6 r 3 s u 9 8 P S S y o N I R j r T s e i k w 3 w c o w w u k k 4 8 e a R p i M 8 J B 2 L J V Y U N 1 N Z l d P 4 J F V + n A Q K l v S w J n 6 f S L B Q u u x C G y n w O Z W / / a m 4 l 9 e J z a D c j d h M o o N l W S + a B B z a E I 4 j Q D 2 m a L E 8 L E l m C h m b 4 X k F i t M j A 0 q Y 0 P 4 + h T + T 5 o n r u e 5 3 l U h V z 1 f x J E G B + A Q H A M P l E A V 1 E A d N A A B C j y A J / D s 3 D u P z o v z O m 9 N O Y u Z f f A D z t s n x O i S s Q = = < / l a t e x i t > Solid: full numerical s t a b l e n e u t r i n o m ⌫ > p ⌫ (> 80%) < l a t e x i t s h a 1 _ b a s e 6 4 = " M HZ I 5 f U K a O 7 7 L r F c L E / 1 1 J m j u z Q = " > A A A C B X i c b V D L S g M x F M 3 U V 6 2 v U Z e 6 C J Z C B S k z I t h V K b p x W c E + o D M M m T R t Q 5 P M k G S E M n T j x l 9 x 4 0 I R t / 6 D O / / G T D s L r Z 4 Q 7 u G c e 0 n u C W N G l X a c L 6 u w s r q 2 v l H c L G 1 t 7 + z u 2 f s H H R U l E p M 2 j l g k e y F S h F F B 2 p p q R n q x J I i H j H T D y X X m d + + J V D Q S d 3 o a E 5 + j k a B D i p E 2 U m A f 8 y D 1 R D J r x I v q n W W n 2 q g 7 X u U 0 s M t O z Z k D / i V u T s o g R y u w P 7 1 B h B N O h M Y M K d V 3 n V j 7 K Z K a Y k Z m J S 9 R J E Z 4 g k a k b 6 h A n C g / n W 8 x g x W j D O A w k u Y K D e f q z 4 k U c a W m P D S d H O m x W v Y y 8 T + v n + h h 3 U + p i B N N B F 4 8 N E w Y 1 B H M I o E D K g n W b G o I w p K a v 0 I 8 R h J h b Y I r m R D c 5 Z X / k s 5 5 z X V r 7 u 1 F u X m V x 1 E E R + A E V I E L L k E T 3 I A W a A M M H s A T e A G v 1 q P 1 b L 1 Z 7 4 v W g p X P H I J f s D 6 + A U q C l x s = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " M H Z I 5 f U K a O 7 7 L r F c L E / 1 1 J m j u z Q = " > A A A C B X i c b V D L S g M x F M 3 U V 6 2 v U Z e 6 C J Z C B S k z I t h V K b p x W c E + o D M M m T R t Q 5 P M k G S E M n T j x l 9 x 4 0 I R t / 6 D O / / G T D s L r Z 4 Q 7 u G c e 0 n u C W N G l X a c L 6 u w s r q 2 v l H c L G 1 t 7 + z u 2 f s H H R U l E p M 2 j l g k e y F S h F F B 2 p p q R n q x J I i H j H T D y X X m d + + J V D Q S d 3 o a E 5 + j k a B D i p E 2 U m A f 8 y D 1 R D J r x I v q n W W n 2 q g 7 X u U 0 s M t O z Z k D / i V u T s o g R y u w P 7 1 B h B N O h M Y M K d V 3 n V j 7 K Z K a Y k Z m J S 9R J E Z 4 g k a k b 6 h A n C g / n W 8 x g x W j D O A w k u Y K D e f q z 4 k U c a W m P D S d H O m x W v Y y 8 T + v n + h h 3 U + p i B N N B F 4 8 N E w Y 1 B H M I o E D K g n W b G o I w p K a v 0 I 8 R h J h b Y I r m R D c 5 Z X / k s 5 5 z X V r 7 u 1 F u X m V x 1 E E R + A E V I E L L k E T 3 I A W a A M M H s A T e A G v 1 q P 1 b L 1 Z 7 4 v W g p X P H I J f s D 6 + A U q C l x s = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " M HZ I 5 f U K a O 7 7 L r F c L E / 1 1 J m j u z Q = " > A A A C B X i c b V D L S g M x F M 3 U V 6 2 v U Z e 6 C J Z C B S k z I t h V K b p x W c E + o D M M m T R t Q 5 P M k G S E M n T j x l 9 x 4 0 I R t / 6 D O / / G T D s L r Z 4 Q 7 u G c e 0 n u C W N G l X a c L 6 u w s r q 2 v l H c L G 1 t 7 + z u 2 f s H H R U l E p M 2 j l g k e y F S h F F B 2 p p q R n q x J I i H j H T D y X X m d + + J V D Q S d 3 o a E 5 + j k a B D i p E 2 U m A f 8 y D 1 R D J r x I v q n W W n 2 q g 7 X u U 0 s M t O z Z k D / i V u T s o g R y u w P 7 1 B h B N O h M Y M K d V 3 n V j 7 K Z K a Y k Z m J S 9 R J E Z4 g k a k b 6 h A n C g / n W 8 x g x W j D O A w k u Y K D e f q z 4 k U c a W m P D S d H O m x W v Y y 8 T + v n + h h 3 U + p i B N N B F 4 8 N E w Y 1 B H M I o E D K g n W b G o I w p K a v 0 I 8 R h J h b Y I r m R D c 5 Z X / k s 5 5 z X V r 7 u 1 F u X m V x 1 E E R + A E V I E L L k E T 3 I A W a A M M H s A T e A G v 1 q P 1 b L 1 Z 7 4 v W g p X P H I J f s D 6 + A U q C l x s = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " M H Z I 5 f U K a O 7 7 L r F c L E / 1 1 J m j u z Q = " > A A A C B X i c b V D L S g M x F M 3 U V 6 2 v U Z e 6 C J Z C B S k z I t h V K b p x W c E + o D M M m T R t Q 5 P M k G S E M n T j x l 9 x 4 0 I R t / 6 D O / / G T D s L r Z 4 Q 7 u G c e 0 n u C W N G l X a c L 6 u w s r q 2 v l H c L G 1 t 7 + z u 2 f s H H R U l E p M 2 j l g k e y F S h F F B 2 p p q R n q x J I i H j H T D y X X m d + + J V D Q S d 3 o a E 5 + j k a B D i p E 2 U m A f 8 y D 1 R D