Improved cosmological constraints on the neutrino mass and lifetime

We present cosmological constraints on the sum of neutrino masses as a function of the neutrino lifetime, in a framework in which neutrinos decay into dark radiation after becoming non-relativistic. We find that in this regime the cosmic microwave background (CMB), baryonic acoustic oscillations (BAO) and (uncalibrated) luminosity distance to supernovae from the Pantheon catalog constrain the sum of neutrino masses $\sum m_\nu$ to obey $\sum m_\nu<0.42$ eV at (95$\%$ C.L.). While the the bound has improved significantly as compared to the limits on the same scenario from Planck 2015, it still represents a significant relaxation of the constraints as compared to the stable neutrino case. We show that most of the improvement can be traced to the more precise measurements of low-$\ell$ polarization data in Planck 2018, which leads to tighter constraints on $\tau_{\rm reio}$ (and thereby on $A_s$), breaking the degeneracy arising from the effect of (large) neutrino masses on the amplitude of the CMB power spectrum.


Introduction
Even though neutrinos were first detected more than six decades ago, they remain among the most mysterious particles in nature, with many of their fundamental properties still to be determined. In particular, although oscillation experiments have provided convincing evidence that neutrinos have non-vanishing masses, these measurements are only sensitive to the mass-squared splittings and consequently the spectrum of neutrino masses remains unknown. The lifetimes of the neutrinos are also poorly constrained, especially in comparison to the other particles in the Standard Model (SM). The determination of the masses and the lifetimes of these mysterious particles remain some of the most important open problems in fundamental physics.
The fact that cosmic neutrinos are among the most abundant particles in the universe, contributing significantly to the total energy density at early times, provides an opportunity to measure their properties. In particular, the evolution of the cosmological density fluctuations depends on m ν , the sum of neutrino masses. This translates into characteristic effects on the cosmic microwave background (CMB) and large-scale structure (LSS) [1,2] (for reviews see [3][4][5][6]), that are large enough to allow the sum of neutrino masses to be determined in the near future. This determination is based on the observation that massive neutrinos contribute differently to cosmological observables than either massless neutrinos or cold dark matter (CDM). At early times, while still relativistic, massive neutrinos contribute to the energy density in radiation, just as in the case of massless neutrinos. However, after neutrinos become non-relativistic, their energy density redshifts as matter and therefore contributes more to the expansion rate than massless neutrinos, which would continue to redshift as radiation. As a result, over a given redshift span, the higher expansion rate reduces the time available for the growth of matter density perturbations. However, since massive neutrinos retain pressure until late times, their contribution to the density perturbations on scales below their free streaming lengths is too small to compensate for the shorter structure formation time. Therefore, if neutrinos become non-relativistic after recombination, the net effect of non-vanishing neutrino masses is a suppression of the matter power spectrum and the CMB lensing potential. Based on this, current observations are able to place a bound on the sum of neutrino masses, m ν 0.12 eV [7]. It is important to note that this result assumes that neutrinos are stable on timescales of order the age of the universe. In scenarios in which the neutrinos decay [8,9], or annihilate away into lighter species [10,11] on timescales shorter than the age of the universe, this bound is no longer valid and must be reconsidered.
Cosmological observations can also be used to place limits on the neutrino lifetime. In the case of neutrinos that decay to final states containing photons, the bounds on spectral distortions in the cosmic microwave background (CMB) can be translated into limits on the neutrino lifetime, τ ν 10 19 s for the larger mass splitting and τ ν 4 × 10 21 s for the smaller one [12]. In the case of decays to invisible final states, the limits are much weaker. For neutrinos that decay while still relativistic, the decay and inverse decay processes can prevent neutrinos from free streaming. Measurements of the CMB power spectra set a lower bound on the neutrino lifetime, τ ν ≥ 4 × 10 6 s (m ν /0.05eV) 5 , in the case of decay into dark radiation [13] (for earlier work see [14][15][16][17]). In the case of non-relativistic neutrino decays into dark radiation, the energy density of the decay products redshifts faster than that of stable massive neutrinos. Unstable neutrinos therefore have less of an effect on structure formation than stable neutrinos of the same mass. Consequently, cosmological observables depend both on the masses of the neutrinos and their lifetimes, and heavier values of m ν may still be allowed by the data provided the neutrino lifetime is short enough. In Ref. [18], Planck 2015 and LSS data were used to place constraints on the neutrino mass as a function of the lifetime, and found that values of m ν as large as 0.9 eV were still allowed by the data. Future LSS measurements at higher redshifts may be able to break the degeneracy between the neutrino mass and lifetime and measure these parameters independently [19]. It is worth noting that there are also bounds on the neutrino lifetime from Supernova 1987A [20], solar neutrinos [21][22][23][24], astrophysical neutrinos measured at IceCube [25][26][27][28][29][30], atmospheric neutrinos and long baseline experiments [31][32][33][34]. However,these constraints are in general much weaker than the limits from cosmology.
In this paper we revisit the scenario in which neutrinos decay into dark radiation after becoming non-relativistic and obtain updated limits based on the newer data from Planck 2018. In order to take advantage of the greater precision of the new data, the analysis we perform is also more accurate. We find that, under the assumption that neutrinos decay after becoming non-relativistic, the neutrino mass bound from Planck 2018 data (in combination with BOSS baryon acoustic oscillation (BAO) data and Pantheon SN1a data) is relaxed to m ν 0.42 eV (95% C.L.). 1 While this represents a remarkable relaxation of the constraints as compared to the case of stable neutrinos, we note that it is much stronger than the limit derived from Planck 2015 data for the same decaying neutrino scenario, m ν 0.9 eV at (95% C.L.). We show that the improvement of the bound arises primarily from the more precise low-polarization data from Planck 2018, which allows an improved determination of the optical depth to reionization τ reio , thereby breaking the correlation with m ν that appears (for relatively high neutrino masses) through the impact of neutrinos on the overall height of the acoustic peaks (i.e. the "early integrated Sachs-Wolfe effect") [4].
Besides using up-to-date cosmological data, we also improve the analysis from Ref. [18] by incorporating higher order corrections due to neutrino decays into the Boltzmann equations that describe the evolution of Universe's energy and metric fluctuations. Recently, Ref. [13] provided a complete set of Boltzmann equations for the neutrino decay, but did not conduct Markov Chain Monte Carlo (MCMC) runs necessary to calculate updated neutrino bounds. In this work, we derive Boltzmann equations exactly valid in the absence of 'inverse-decays' and quantum statistics. For the numerical implementation, we follow a consistent T dec /m ν expansion, where T dec is the temperature at the time of the decay, so that the analysis is under control when neutrinos decay after become non-relativistic.
This paper is organized as follows. In section 2, we present a summary of constraints on the parameter space of decaying neutrinos. In section 3, we derive the set of Boltzmann equations to describe neutrino decay that are valid in the non-relativistic regime and compare our improved analysis to past work. In section 4, we present a MCMC analysis of the decaying neutrino scenario against up-to-date cosmological data. Finally, we conclude in section 5.

Parameter space of decaying neutrinos
In this section we outline the constraints on the mass and lifetime of neutrinos decaying into dark radiation. As explained in the introduction, current cosmological observables only place limits on a combination of the sum of neutrino masses and their lifetime. Therefore, in this study we will map out the constraints in the two-dimensional parameter space spanned by the sum of neutrino masses ( m ν ) and the neutrino decay width (Γ ν ), as shown in Fig. 1. In our analysis we assume that all three neutrinos are degenerate in mass. This is a good approximation because the current bounds on m ν are larger than the observed mass splittings (see Fig. 1). 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. While this is a simple parameterization of neutrino decays, our bounds can easily be applied to specific models, as done in great details in Ref. [36].
The CMB can be used to constrain the masses and decay widths of neutrinos that decay prior to recombination.When neutrinos decay while still relativistic, decay and inverse decay can prevent neutrinos from free-streaming. If this happens before recombination, it can alter the well-known 'neutrino drag' effect that manifests as a phase-shift at high-'s in the CMB power spectrum [37][38][39][40]. The plot shows the current constraints on decaying neutrinos in the m ν − Γ ν parameter space. The colored regions are excluded by current data while the white region is allowed. The orange dashed line represents Γ ν = H(a nr ). Our study focuses on the region below this line, meaning decay happens after neutrinos have become non-relativistic. The grey region shows current constraints on neutrino mass and lifetime coming from the requirement that neutrinos are free streaming close to recombination [13]. The light grey region indicates that this bound may not be applicable when neutrino mass is larger than the temperature of recombination: m ν > T * ∼ 0.2 eV [13]. Our analysis excludes the red (blue) region labelled "Planck 2015"("Planck 2018") based on the data (Planck+BAO+Pantheon). The vertical brown line shows the projected KATRIN sensitivity. Therefore, CMB data can place a constraint on the decay width of neutrinos. The resulting bound depends on neutrino masses, and was recently updated in Ref. [13], τ ν ≥ 4 × 10 6 s (m ν /0.05 eV) 5 . This bound excludes the grey region at the top of Fig. 1.
In addition, based on the analysis in this paper, part of the 'late-decay' parameter space can also be excluded based on the gravitational impacts of massive neutrinos on the CMB and LSS. Through the Monte Carlo study presented in section 4, the blue (red) shaded region in Fig. 1 is excluded by the data combination Planck 2018(2015)+BAO+Pantheon. 2 The orange dashed line in the figure (Γ ν = H(a nr )) separates the region where neutrinos decay when non-relativistic from the region where they decay while still relativistic. Here a nr corresponds to the approximate scale factor at 2 Note that in our analysis we scanned the region between 0 ≤ log 10 Γν km/s/Mpc ≤ 6. In Fig. 1, we have extrapolated the bound at log 10 Γν km/s/Mpc = 0 to Γν = 0, because the constraint on mν is independent of Γν when Γν H0.
the time that neutrinos transition to non-relativistic, and is defined as 3T ν (a nr ) = m ν . This simple 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 a nr is given by, where T ν0 is the present neutrino temperature. Since our study focuses on the decay of neutrinos after they become non-relativistic, we only present constraints below the orange dashed line. Our analysis shows that m ν as large as 0.42 eV is still allowed by the data.
Our results have important implications for current and future laboratory experiments designed to detect neutrino masses. Next generation tritium decay experiments such as KATRIN [41] are expected to be sensitive to values of m νe as low as 0.2 eV, corresponding to m ν of order 0.6 eV. Naively, a signal in these experiments would conflict with the current cosmological bound for stable neutrinos, m ν < 0.12 eV. However, since the unstable neutrino paradigm greatly expands the range of neutrino masses allowed by current cosmological data, it is interesting to explore whether this scenario can accommodate a potential signal at KATRIN. In Fig. 1, we display a brown vertical line m ν = 0.6 eV that corresponds to the expected KATRIN sensitivity. We see that this value of m ν is too large to be accommodated in the non-relativistic decay regime, where our analysis is valid. However, our result, in combination with those from the 'relativistic decay' scenario studied in Ref. [13], leaves open the interesting possibility that neutrinos decaying with a decay width between log 10 Γν km/s/Mpc ∼ 5.5 − 9 could reconcile cosmological observations with a potential detection at KATRIN, thereby opening a large discovery potential for laboratory experiments. To confirm this conjecture, more work needs to be done to cover the 'intermediate' decay regime (i.e. where neutrinos are neither fully relativistic nor fully non-relativistic). We leave this for future work.
In recent years, a number of studies have attempted to constrain the neutrino mass ordering, showing that under the assumption of stable neutrinos, the inverted ordering is now disfavored by constraints from joint analysis of cosmological and oscillation data [42][43][44][45][46][47] (see also Refs. [48][49][50][51] for a different take) as well as from Ly-α observations [52]. However, these arguments are centered on the fact that these analysis lead to a constraint on m ν at odds with the lower bound on the sum of neutrino masses in the case of inverted ordering, m ν 0.1 eV. Our result suggests that these constraints are strongly dependent on the assumption of neutrino stability over cosmological timescales, and therefore that the inverted ordering is not robustly excluded. It would be very interesting to extend our analysis to the inclusion of Ly-α data to confirm this conclusion.

Boltzmann equations for massive neutrinos decaying into radiation
In this section, we revisit the set of Boltzmann equations describing the evolution of the phase space distribution (PSD) of massive particles decaying into daughter radiation. In our analysis, we assume the decay happens after the neutrinos have become non-relativistic so that the contribution from inverse decay processes can be safely neglected.

Derivation of the equations
We denote the phase space distribution of each species as f (q,n, x, τ ), which is a function of the comoving momentum qn, coordinates x and conformal time τ . The general time evolution of f is controlled by the Boltzmann equations, where C[f ] is the collision term that includes all the processes involving the species. This phase space distribution has the leading order contributionf (q, τ ) that only depends on q and τ , while perturbations are encoded in ∆f (q,n, x, τ ), Treating ∆f fluctuations about the homogeneous background as higher order perturbations, the zeroth order Boltzmann equations forf take the form In this work, our focus is on the case in which neutrinos decay after turning non-relativistic. In this scenario, we can neglect the effects of inverse decay processes and quantum statistics. The collision term for the neutrino and its daughters are respectively given by [18] Hered 3 q ≡ d 3 q/(2π) 3 , ≡ q 2 + a 2 m 2 represents the comoving energy and a is the scale factor. The label i(j) denotes the ith(jth) daughter. In the case of two body decays to massless daughters, the amplitude squared |M| 2 is simply related to the rest-frame decay width of the neutrino as |M| 2 = 16πΓ ν m ν . From the collision terms above, the background evolution for decaying neutrinos is given by where Γ ν is the neutrino decay width and γ is the Lorentz boost factor, The formal solution tof ν (q, τ ) from the differential equation Eq. (3.6) is where τ ini denotes the initial conformal time andf ini (q) represents the initial momentum distribution, which we take to be of the Fermi-Dirac form,f ini = 1/(e q/T ν0 + 1). The Boltzmann equations for the individual daughter particles do not have a simple form, especially when the daughters consist of more than two species. However, since the daughter particles are taken to be massless in this study the total background density of daughter radiation can be defined as wheref Di is the background phase space distribution of the ith daughter particle. With the definition in Eq. (3.9), regardless of the number of daughter particles and their spins, the Boltzmann equation for the total background daughter densityρ D has the simple form . We now turn to the Boltzmann equations describing the perturbations of the phase space distribution of decaying neutrinos and their decay products. We work in the synchronous gauge for which the metric perturbations can be parametrized as [53] (3.11) In Fourier space, H ij is given by where k is conjugate to x and h and η are the two independent scalar metric perturbations. To obtain the Boltzmann hierarchy, we expand the angular dependence of the perturbations as a series in Legendre polynomials, where P represents the th Legendre polynomial. The Boltzmann hierarchy for the perturbations of the decaying massive neutrinos ∆f ν( ) read [18] ∆ḟ Here ε ν = q 2 + a 2 m 2 ν indicates the comoving energy of the neutrinos. To study the perturbations of the daughter radiation, we focus on the case of two-body decay. In this case, the phase space distributions of the two massless particles are basically identical and they can be considered effectively as one species with a single f D . We can therefore define multipoles F D( ) as in Ref. [54], where ρ c is the critical density of Universe today. The Boltzmann hierarchy of the F D( ) can be written as,Ḟ where r D ≡ a 4ρ D /ρ c . The terms C appearing in Eq. (3.19) arise from the integrated daughter collision term in Eq. (3.5) expanded in terms of Legendre polynomials. The expression for C is given by, (3.20) The overall factor of two in the equation above arises because we are adding the collision integrals of the two massless daughters, which are of the same form. In this expression dΩ k represents the differential solid angle along the directionk, while q 1,2 are the momenta of daughter particles. Thed 3 q 2 integral can be easily evaluated using the delta function corresponding to momentum conservation. In order to perform the integral over dΩ k , we notice that the direction ofk enters only via P (q 1 ·k) and ∆f ν (q,q ·k). Now, using the Legendre expansion of ∆f ν (q,q ·k) in Eq. (3.13) and employing the identity we can evaluate the dΩ k integral to obtain Now, notice that the direction of the neutrino momentum only enters the integrand via the angle between the neutrino momentum q and the daughter momentum q 1 , defined as cos θ 1 ≡q ·q 1 . The energy conserving delta function can be expressed in terms of this angle as where The energy conservation restricts the daughter momentum to a range of values (q + 1 , q − 1 ). The edges of this range occur when the extreme values, cos θ * 1 = ±1, are reached. For these values, . (3.25) After integrating over the delta function corresponding to energy conservation, this reduces to the simpler form, Eq. (3.26) may also be obtained by taking the appropriate limit of the more general expression in Ref. [13]. The same Boltzmann hierarchy has been derived in the context of warm matter decaying into dark radiation [55]. Performing the integral over q 1 , we can obtain the following expressions for the first few C 's, Here the functions g 2 (q, ε ν ) and g 3 (q, ε ν ) are given by, Given the complicated integrals in Eq. (3.26), it is technically challenging to keep track of all the collision terms in the Boltzmann hierarchy. Instead, we choose to keep just the first few C 's for ≤ max . The idea behind this approach is that C is of O((T dec /m ν ) ) around the time of decay.
Therefore, for non-relativistic decay (T dec /m ν 1), it is self-consistent to set C > max = 0 because those terms only have negligible effect on physical observables. To understand the scaling of C , we first note that the integral over q in Eq. (3.27) receives most of its support from the region around q ∼ T ν0 because ∆f ν( ) inherits features of the Fermi-Dirac distribution fromf ini = 1/(e q/T ν0 + 1). Deep in the non-relativistic region, q ν and T ν0 am ν . In this regime, we can employ a Taylor expansion for the functions g 2 and g 3 in powers of q/ ν to obtain, Inserting Eq. (3.29) above into Eq. (3.27), it is straightforward to see that C ∝ (T ν0 /am ν ) . Moreover, if we assume decay happens deep in the non-relativistic region, we will get C ∝ (T dec /m ν ) when decay happens, where T dec = T ν0 /a dec . Therefore, C is suppressed by powers of T dec /m ν 1 for higher . To further justify this argument, we show in section 3.3 that setting max = 2 or max = 3 makes negligible difference to cosmological observables (see Fig. 4). Therefore, we only keep C ≤3 and set C >3 = 0 in our numerical study for simplicity. Physically, the expansion in the small parameter T dec /m ν corresponds to perturbing about the ultra-nonrelativistic limit in which the momentum of the mother particle has completely redshifted away, so that it has come to rest in the cosmic frame. Energy and momentum conservation is respected order by order in this expansion. The earlier work [18] approximated the Boltzmann hierarchy for daughter radiation (Eq. 3.19) by just keeping C 0 and setting all the C ≥1 = 0. It is clear from the above discussion that this is a consistent approximation to zeroth order in an expansion in the small parameter T dec /m ν . The authors in Ref. [13] argued that the Boltzmann hierarchy for daughter radiation in Ref. [18] does not reproduce the standard decaying CDM scenario and does not respect momentum conservation. Both criticisms can be addressed by considering the term C 1 . Since C 1 begins at O(T dec /m ν ), we see that the Boltzmann hierarchy in Ref. [18] does in fact reproduce the decaying CDM scenario and respects momentum conservation up to O(T dec /m ν ) corrections, consistent with the approximation. In this limit, the momenta of the daughter particles arise entirely from the rest mass of the mother. In practice, since the contributions of neutrinos to the density perturbations are small, we will see that the higher order terms do not significantly affect the constraints derived in Ref. [18] with Planck 2015 data.

Signatures of the non-relativistic neutrino decay on the CMB spectra
To make this work fully self-contained, we briefly summarize the impact of the non-relativistic invisible neutrino decays on the CMB spectra, following the discussion in Ref. [18]. In Fig. 2, we display the residuals in the CMB (lensed) TT, EE and lensing power spectra, for the sum of neutrino masses m ν = 0.6 eV and several decay widths Log 10 (Γ ν /km/s/Mpc) = 0, 2, 4, 6. In all cases, the ΛCDM parameters are set to their best-fit values from Planck 2018, that is, {100θ s = 1.04089, ω cdm = 0.1198, ω b = 0.02233, n s = 0.9652, ln(10 10 A s ) = 3.043, τ reio = 0.0540}. Our reference ΛCDM model makes use of the same parameters and assumes standard massless neutrinos.
For the value of the mass considered ( m ν = 0.6 eV) and at fixed angular size of the sound horizon θ s , neutrino masses primarily impact the lensing spectrum. Indeed, as they reduce power below the free-streaming scale, they produce a significant matter power suppression at small scales, which leads to a ∼ 20% reduction in the C φφ at large (blue curve in Fig. 2). Consequently, this power suppression decreases the smoothing in the high-part of the TT and EE spectra, which can be seen as 'wiggles' in the corresponding plots.
In addition, stable neutrinos dilute like non-relativistic matter at late times (ρ ν ∼ a −3 ), which increases the value of Ω m . As we impose the closure relation Ω m + Ω Λ = 1 at late-times, this is compensated for by a decrease in Ω Λ (later beginning of Λ-domination), and thus a reduction in the Late Integrated Sachs-Wolfe effect (LISW), leaving a signature in the low-TT spectrum. Furthermore, the modified expansion history H(z) changes quantities integrated along z, such as τ reio , which affects the multipoles at ∼ 10 in the EE spectrum.
When a non-negligible Γ ν is considered (orange, green and red curves in Fig. 2), one can see that the aforementioned effects typically become less prominent for earlier decays. This is particularly true for the high-part of the lensing spectrum (and consequently the smoothing at high-in TT and EE) since decay of neutrinos reduce their impact on structure formation. The reduction of the effect in the low-part of the TT and EE spectra is not entirely monotonic, as intermediate values of Γ ν can induce additional time variation in the gravitational potentials (thereby affecting the LISW effect), as well as time variations in H(z) (thereby affecting τ reio ). As a result, the ΛCDM limit is reached not only for small values of m ν , but also for high values of Γ ν . This will be reflected in the MCMC analysis in section 4, which shows a large positive correlation between both parameters. It is precisely this degeneracy which relaxes the neutrino mass bounds.

Consistency of the implementation of Boltzmann equations
We begin by comparing the approximation used in Ref. [18] for the background energy density of decaying massive neutrinos to the more accurate results obtained by evaluating the integral in Eq. (3.8) numerically. In Ref. [18], the phase space distribution of neutrinos in Eq. (3.8) is approximated through the following analytic formula, (3.30) Figure 3. Redshift evolution of the quantity (ρ ν +ρ D )/ρ ur (whereρ ur denotes the energy density of stable massless neutrinos), which should be equal to 1 in the limit of relativistic decays. We consider a very small value of the neutrino mass sum, m ν = 0.06 eV, and several values for the decay width, Log 10 (Γ ν /km/s/Mpc). "approx. PSD" refers to the approximated phase space distribution in Eq. (3.30) while "Full PSD" refers to the exact solutions of Eq. (3.8).
As argued in Ref. [18], this approximation is valid under the assumption that the decay happens deep in the non-relativistic regime. To see the difference between the approximation and the full result, we plot the ratio r ≡ (ρ ν +ρ D )/ρ ur in Fig. 3, for several values of the decay width Γ ν and a fixed value of the total neutrino mass m ν = 0.06 eV. Hereρ ur denotes the energy density of stable massless neutrinos. If neutrinos decay while relativistic, this ratio always gives r 1. However, if the decay happens when the neutrinos are already non-relativistic (ρ ν ∼ a −3 ), then the ratio evolves from r 1 to r ∼ a, and will eventually reach a plateau once all the neutrinos have decayed. From Fig. 3, we can see that the approximate formula in Eq. (3.30) gradually improves as we go to smaller decay widths (that is, going deeper into the regime of non-relativistic decays), as expected. The error in the case of neutrinos decaying right around the time of the non-relativistic transition (Log 10 (Γ ν /[km/s/Mpc]) 4 for 0.06 eV) is around 25%. Nevertheless, as we argue below, the impact on observables is much smaller given that neutrinos only contribute a small fraction of the total energy density for masses considered in this work. Not surprisingly, the approximate formula fails in the relativistic regime, leading to r < 1 at late-times. Therefore future work focusing on this regime should make use of the exact formula.
In Figs. 4 and 5, we show the effects of various approximations in dealing with decaying neutrinos (at the background and perturbation level) on the CMB TT, EE and lensing spectra. We compare the impact of using either the approximated or the exact PSD of neutrinos discussed above, as well as the impact of only keeping C ≤ max in the Boltzmann hierarchy of daughter particles in Eq. (3.19), where we vary max from zero to three. We show the residuals of these approximations with respect to the 'optimal' case (i.e. including all terms up to max = 3 and the exact background PSD) for a fixed value of the neutrino mass ( m ν = 0.6 eV) and two different decay widths (Log 10 (Γ ν /[km/s/Mpc]) = 5.5 in Fig. 4 and Log 10 (Γ ν /[km/s/Mpc]) = 4 in Fig. 5). Fig. 4 corresponds to decays happening around the time of the non-relativistic transition, T dec /m ν 0.3, where the effects of the approximations are expected to be largest. Fig. 5 on the other hand refers to decays happening deep in the non-relativistic regime, T dec /m ν 0.03. We also show the Planck 2018 1-σ error bars, as well as the (binned) cosmic variance.
For decays close to the non-relativistic transition T dec /m ν 0.3 shown in Fig. 4, we find that the biggest improvement in the CMB TT spectrum occurs when including C ≤1 (i.e., the contribution from the decaying neutrino bulk velocity) in the Boltzmann hierarchy of daughter radiation, which impacts the integrated Sachs-Wolfe (ISW) effect at multipoles 100. On the other hand, the approximate background distribution of neutrinos does not have a significant effect. For the CMB EE spectrum shown in the same figure, which is not sourced by the ISW effect, the impact of the approximate background distribution of neutrinos is comparable to the effect of the approximate perturbed hierarchy. Nevertheless, one can see that for max ≥ 2, additional contributions to the daughter hierarchy have negligible impacts, which justifies our choice of cutting the collision term C contribution at max = 3. Finally for the CMB lensing spectrum, the effects due to the approximate treatment of the background PSD dominate over the ones due to including higher order terms in the Boltzmann hierarchy of the dark radiation. This is expected given that the matter power spectrum suppression scales approximately withρ ν /ρ m [2,4] whereρ m is the total matter density, while neutrino perturbations are very small well below the free-streaming scale, so that their detailed dynamics is not as important as on larger scales.
The impact of the various approximations in the case of decays deep in the non-relativistic regime T dec /m ν 0.03, displayed in Fig. 5, is much less visible. In that case, one can therefore safely neglect C >0 and consider the approximate PSD, as done in Ref. [18]. In this section we perform a numerical scan over the parameter space to obtain updated limits on the neutrino mass and lifetime. We perform comprehensive MCMC analyses with the MontePython-v3 3 [56,57] code interfaced with our modified version of CLASS. We fit the decaying neutrino model to a combination of the following data-sets: • The Planck 2018 high-TT, TE, EE + low-data TT, EE + lensing data [7]. We will also compare these results with the use of Planck 2015 data to disentangle the effects of our improved formalism and that of the new data.
We adopt wide flat priors on the following six ΛCDM parameters: {ω b , ω cdm , H 0 , n s , A s , τ reio }. We assume three degenerate neutrinos decaying into massless radiation and consider flat priors on m ν /eV and Log 10 (Γ ν /[km/s/Mpc]). To accelerate convergence, we split the parameter space between Log 10 (Γ ν /[km/s/Mpc]) ∈ [0.1, 2.5] and Log 10 (Γ ν /[km/s/Mpc]) ∈ [2.5, 6.5]. In both cases we take wide priors on m ν ∈ [0.06, 1.5] eV. We assume our MCMC chains to be converged when the Gelman-Rubin criterion R − 1 < 0.05 [64]. In our baseline analysis, we do not apply any specific cut to the parameter space, even if neutrinos decay in the relativistic regime (this occurs for low m ν and high Γ ν ). In appendix A, we investigate the impact of imposing a prior that excludes the parameter space corresponding to relativistic decay from our analysis and show that the limit at 95% on m ν agrees within a few percent.

Main results: updated limit on the neutrino mass and lifetime
The results of our analyses are presented in Figs. 6. For very late decays, Log 10 (Γ ν /[km/s/Mpc]) 2.5, no relaxation of the constraints on m ν /eV is visible, in agreement with what was found in Ref. [18]. The impact of the new Planck data is visible as a significantly improved bound on the sum of neutrino mass, namely we find m ν < 0.127 eV (95%C.L.), an improvement of about ∼ 35% over 2015 data, in good agreement with Ref. [7]. For Log 10 (Γ ν /[km/s/Mpc]) 2.5, one can see that the bound relaxes as expected, although not as much with Planck 2018 data as for Planck 2015 data.
Taking the intersect of the non-relativistic decay line as our 2σ limit, we find that Planck 2018 allows neutrinos with masses up to m ν = 0.42 eV. In appendix A, we present an alternative analysis that directly imposes the non-relativistic decay criterion as a prior while performing the scan. Marginalizing over all parameters we find the same result, m ν < 0.42 eV (95% C.L.). The excellent agreement between these different analyses leads us to conclude with confidence that, within the regime of non-relativistic decay, values of m ν as large as 0.4 eV are still allowed by the data. This bound is significantly stronger than the limit from Planck 2015 data, for which m ν ∼ 0.9 eV was still allowed in the non-relativistic decay scenario.
Our result also has implications for laboratory searches. For m ν = 0.6 eV, the smallest mass scale that the KATRIN experiment is designed to probe, Planck 2018 data requires decay rate Γ ν 10 5.5 km/s/Mpc, a constraint roughly one order of magnitude stronger than from Planck 2015 data. However, this value of the decay rate is now slightly beyond the regime of validity of our work 4 , indicating that, in the event of a neutrino mass discovery at KATRIN, a more involved analysis including inverse-decays would be necessary to confirm that the decay scenario can reconcile laboratory and cosmological measurements.

Comparison with former results and the impact of Planck 2018 data
Comparing with the constraints presented in Ref. [18] for Planck 2015, we find that, while the impact of our improved treatment is clearly visible in the CMB power spectra (and will be relevant for future experiments), it has only a marginal impact on the constraints, and our bounds are in very good agreement with those derived in Ref. [18], which only included the leading order term in the daughter radiation hierarchy 5 . The bulk of the improvement is due to the newest Planck 2018 data and can be understood as follows. As shown in Fig. 2, for the masses we consider, the main effect is an almost scale independent suppression of CMB lensing spectrum. This suppression can be compensated for by increasing the primordial amplitude A s or by adjusting the matter density ω cdm (see Ref. [65] for 5 Let us note that the implementation of the BAO/f σ8 DR12 likelihood used in Ref. [18] within the MontePython code had an issue that led to constraints on mν that were somewhat milder than the true bounds. MontePython has since then been corrected, leading to an improvement on the constraints on the stable/long-lived (Γν < 10 3 km/s/Mpc) case by about 20%. However, we have verified that this bug had no impact in the short-lived case (Γν > 10 2.5 km/s/Mpc). a discussion of the correlation between { m ν , A s , τ reio , ω cdm }). Due to the well-known degeneracy between A s and e −2τ reio , Planck 2015 data, which was limited in polarization, were unable to place a tight constraint on A s , and thus the constraining power on the sum of neutrino mass and lifetime was limited. The precise measurements of low-polarization from Planck 2018 leads to constraints on τ reio that are tighter by a factor of two than those from Planck 2015. As a result, parameters degenerate with τ reio such as A s are now much better constrained. Consequently, the constraints on the sum of neutrino mass and lifetime have significantly improved with Planck 2018 data. To confirm this simple argument, we perform another MCMC run with Planck 2015 data and a tight gaussian prior on τ reio = 0.0540 ± 0.0074, chosen to match the optical depth to reionization reconstructed from Planck 2018. Given that the constraints on m ν are independent of Γ ν below Γ ν 10 3 , and the scaling above Γ ν 10 5.5 is monotonic, we focus on the parameter space Log 10 (Γ ν /[km/s/Mpc]) ∈ [3, 5.5] to accelerate convergence. Our results are presented in Fig. 7, where one can see that this simple prescription leads to constraints that are very similar to those from the full Planck 2018 data. We attribute the remaining differences to the additional constraining power of Planck 2018 data on the parameters ω cdm and ω b , which are mildly correlated with m ν (see Fig. 6, top panel). Note that our constraints are a factor of two weaker than those advocated in Ref. [35], which performed a 'model-independent' reconstruction of the neutrino mass as a function of redshift, but neglects the decay products. As we show here, including details about the daughter radiation is necessary to accurately compute the effect of neutrino decays even in the non-relativistic regime. Finally, as discussed in Refs. [19,65], a combination of CMB data with future tomographic measurements of the power spectrum by DESI [66] or Euclid [67], and an improved determination of the optical depth to reionization by 21-cm observations with SKA [68,69], could greatly increase the sensitivity of cosmological probes to neutrino masses and lifetimes.

Conclusions
Cosmological observations are known to set the strongest constraints on the sum of neutrino masses. Yet, the existing mass bound from CMB and LSS measurements, which assumes that neutrinos are stable, is significantly weakened if neutrinos decay. In this work, we provide up-to-date limits on the lifetime of massive neutrinos that decay into dark radiation after becoming non-relativistic, from a combination of CMB, BAO, growth factor measurements, and Pantheon SN1a data.
Compared to the earlier analysis [18], we have incorporated higher-order corrections up to O((T dec /m ν ) 3 ) when solving the dark radiation perturbations, and also performed the full calculation of the background energy density of the decaying neutrino using Eq. (3.8). The more precise treatment of the Boltzmann equations and the background energy evolution in our MCMC study improves the coverage of the case when the neutrinos decay early so that their average momenta are close to their masses. As shown in Fig. 5, if neutrinos decay when having T ν m ν /3, the inclusion of higher moment perturbations C ≥2 gives a negligible change to the power spectra as compared to the experimental uncertainties. However, the complete calculation of the neutrino energy does improve the prediction for the power spectrum significantly from the approximate result using Eq. (3.30) when the decays happen semi-relativistically. Nevertheless, we have found that constraints from Planck 2015,  given their limited precision, are unaffected by these considerations. However, we anticipate that these effects will be relevant for future experiments (as well as an essential contribution in the relativistic case, to be considered in the future). In fact, we have shown that the bulk of the improvement in the constraining power compared to Ref. [18] comes from the use of Planck 2018 data. Indeed, we have demonstrated that the improved τ reio measurement from the low-polarization data helps breaking the degeneracy in the CMB power spectrum amplitude and strengthens the bound on the neutrino mass and lifetime. As a result, we have found that neutrinos with m ν > 0.42 eV (2σ) cannot be made consistent with cosmological data if they decay while non-relativistic, a significant improvement from Planck 2015 data for which masses as high as m ν ∼ 0.9 eV were consistent with the non-relativistic decay scenario [18].
We have argued that one notable application of this result is that, if the KATRIN experiment sees an electron neutrino with m ν ≈ 0.2 eV (the advocated sensitivity), our result would constrain Γ ν 10 5.5 km/s/Mpc, i.e. the neutrinos would need to decay between z ≈ 2 × 10 2 − 4 × 10 3 , while they are still relativistic, so that our bounds and the bounds studied in Ref. [13] would not apply. In case of a neutrino mass discovery at KATRIN, a more involved analysis including inversedecays would be necessary to firmly confirm that the decay scenario can reconcile laboratory and cosmological measurements. Additionally, our results show that the tentative exclusion of the inverted mass ordering [44][45][46]52], based solely on the fact that the inverted ordering predicts m ν > 0.1 eV, is highly dependent on the hypothesis that neutrinos are stable on cosmological time-scales. Nonrelativistic decays can still easily reconcile the inverted ordering with cosmological data.
Finally, let us mention that even though current exclusion bounds in Fig. 6 do not set independent constraints on the neutrino mass and lifetime, next generation measurements of the matter power spectrum at different redshifts can help break that degeneracy [19]. It will be interesting to revisit the forecast on the sensitivity of future cosmological data to the sum of neutrino masses and their lifetime in light of our improved formalism.