# Cosmological Limits on the Neutrino Mass and Lifetime

Zackaria Chacko, Abhish Dev, Peizhi Du, Vivian Poulin, Yuhsin Tsai

► To cite this version:

Zackaria Chacko, Abhish Dev, Peizhi Du, Vivian Poulin, Yuhsin Tsai. Cosmological Limits on the Neutrino Mass and Lifetime. Journal of High Energy Physics, 2020, 04, pp.020. 10.1007/JHEP04(2020)020. hal-02317338

HAL Id: hal-02317338

https://hal.science/hal-02317338

Submitted on 26 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Cosmological limits on the neutrino mass and lifetime

Zackaria Chacko,ª Abhish Dev,ª Peizhi Du,ª Vivian Poulinᵇ and Yuhsin Tsaiª

**ABSTRACT:** 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.

**KEYWORDS:** Cosmology of Theories beyond the SM, Neutrino Physics

ARXIV ePRINT: 1909.05275

<footnote>ªMaryland Center for Fundamental Physics, Department of Physics,
University of Maryland, College Park, MD 20742-4111, U.S.A.</footnote>

<footnote>ᵇLaboratoire Univers & Particules de Montpellier (LUPM),
CNRS & Université de Montpellier (UMR-5299),
Place Eugène Bataillon, F-34095 Montpellier Cedex 05, France

E-mail: zchacko@umd.edu, abhishdev92@gmail.com,
peizhi.du@stonybrook.edu, vivian.poulin@umontpellier.fr,
yhtsai@umd.edu</footnote>
# Contents

1 Introduction 1

2 Parameter space of the unstable neutrino 5

3 Evolution of perturbations in the decay of non-relativistic particles into radiation 8

3.1 Background: zeroth order 9

3.2 Perturbations: first order 11

4 Cosmological signals of neutrino decay 12

4.1 Numerical results 13

4.2 Analytic understanding 15

5 Current limits on the neutrino mass and lifetime from Monte Carlo analysis 20

5.1 The data and analysis pipeline 21

5.2 Current limits on the neutrino mass and lifetime 22

6 Conclusions 23

A A model of massive neutrino decay into dark radiation 23

## 1 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, $\sum m_\nu < 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–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–16], (useful discussions may also be found in [17, 18]). In these scenarios, the lifetime of the massive neutrino is given by $\tau_{\nu} \sim 10^{50}\text{s} (0.05 \text{ eV}/m_{\nu})^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–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–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–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, $\tau_{\nu} \gtrsim 10^{19}\text{s}$ for the larger mass splitting and $\tau_{\nu} \gtrsim 4 \times 10^{21}\text{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–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–44]. Current measurements of the CMB power spectra require neutrinos to
free stream from redshifts $z \approx 8000$ until recombination, $z \approx 1100$ [45–48].<fnref n="1" /> This can be used to set a lower bound on the neutrino lifetime $\tau_\nu \ge 4 \times 10^8$ s ($m_\nu/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, $\tau_{\nu_e}/m_{\nu_e} \ge 5.7 \times 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 $\nu_2$, $\tau_\nu/m_\nu \gtrsim 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–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, 5] 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 non-relativistic earlier and also contribute a greater fraction of the total energy density after becoming non-relativistic, a larger neutrino mass results in 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 non-relativistic 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

<footnote n="1">¹Also see the more recent discussion in [49] for the effects of interacting neutrinos on the CMB.</footnote>
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 $\sum m_\nu$ 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<fnref n="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].<fnref n="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 $\sum m_\nu$ 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 $\sum m_\nu$ than the decay of heavier neutrinos to lighter ones. In particular, at present the cosmological limits on $\sum m_\nu$ 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 section 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 section 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 section 4.2 we derive analytical expressions for these effects. In section 5, we perform a Monte Carlo scan of the parameter space and derive constraints on the mass and lifetime of the neutrino

<footnote n="2">²http://www.class-code.net</footnote>

<footnote n="3">³While this analysis was being finalized, the Planck 2018 data became public [59]. We leave the analysis using Planck 2018 data to future work.</footnote>
<figure><img src="image_2.png" /></figure>

**Figure 1.** The plot shows the current constraints in the $\sum m_\nu - \Gamma_\nu$ 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 dash-dotted 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.

from current data. Our conclusions are in section 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.

## 2 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 ($\sum m_\nu$) and the neutrino decay width ($\Gamma_\nu$), as displayed in figure 1. In our analysis we make the simplifying assumption that all three neutrinos are degenerate in mass. As we shall see, the bounds on $\sum m_\nu$ 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 $\Gamma_{\nu}$. 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 $\sum m_{\nu} \ge \sqrt{\Delta m_{31}^2 + \Delta m_{21}^2} = 0.06$ eV in the case of normal ordering and $\sum m_{\nu} \ge 2 \times \sqrt{\Delta m_{31}^2} = 0.1$ eV in the case of inverted ordering [6]. Therefore, we present the parameter space starting from $\sum m_{\nu} = 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_{\text{eff}}$, during the epoch of acoustic oscillations to be $2.99 \pm 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, $\sum m_{\nu} \lesssim 3T_{\text{rec}} \approx 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_{\text{eff}}$ to be free streaming at redshifts below $z \approx 8000$ until recombination, $z \approx 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 $\tau_{\nu} \ge 4 \times 10^8$ s ($m_{\nu}/0.05$ eV)$^3$ [48]. This bound excludes the grey region at the top of figure 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, $\sum m_{\nu} \le 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 $\sum m_{\nu}$ and the neutrino lifetime. Based on the Monte Carlo study presented in section 5, CMB and LSS data (Planck+BAO+Pantheon+LSS) exclude the blue region labelled as “CMB+LSS (this work)” in figure 1. We have scanned the region between $0 \le \log_{10} \frac{\Gamma_{\nu}}{\text{km/s/Mpc}} \le 5.5$. In figure 1, we simply extrapolate the bound at $\log_{10} \frac{\Gamma_{\nu}}{\text{km/s/Mpc}} = 0$ to $\Gamma_{\nu} = 0$, because the constraint on $\sum m_{\nu}$ is independent of $\Gamma_{\nu}$ when $\Gamma_{\nu} \ll 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 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$, 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 $\sum m_\nu \gtrsim 0.9$ eV, where the CMB limits on $N_{\text{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 ($\Gamma = H(z_{\text{nr}})$) separates the region where neutrinos decay when non-relativistic from the region where they decay while still relativistic. Here $z_{\text{nr}}$, the approximate redshift at which neutrinos become non-relativistic, is defined implicitly from the relation $3T_\nu(z_{\text{nr}}) = m_\nu$. This definition is based on the fact that for relativistic neutrinos at temperature $T_\nu$, the average energy per neutrino is approximately $3T_\nu$. The Hubble scale at $z_{\text{nr}}$ is given by,

$$ H(z_{\text{nr}}) = H_0 \sqrt{\Omega_m} \left( \frac{\sum m_\nu}{9 T_\nu^0} \right)^{3/2} \qquad (2.1) \\ \approx 7.5 \times 10^5 \text{km/s/Mpc} \left( \frac{H_0}{68 \text{km/s/Mpc}} \right) \left( \frac{\Omega_m}{0.3} \right)^{1/2} \left( \frac{\sum m_\nu}{1 \text{eV}} \right)^{3/2} \left( \frac{1.5 \times 10^{-4} \text{eV}}{T_\nu^0} \right)^{3/2}. $$

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 figure 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_{\nu_e}$ as low as 0.2 eV, corresponding to $\sum m_\nu$ of order 0.6 eV. A signal in these experiments would conflict with the current cosmological bound, $\sum m_\nu < 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 $\sum m_\nu$ 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 figure 1, we display a brown vertical line $\sum m_\nu \approx 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 $\sum m_\nu \gtrsim 0.6$ eV (brown vertical
line) [61, 63]. An updated version of KLZ, the KLZ-800, is currently probing $\sum m_\nu$ 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 $\sum m_\nu < 0.12$ eV, since the decaying neutrino paradigm would still admit a normal hierarchy.

### 3 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–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 $\vec{x}$, the comoving momentum $\vec{q} \equiv q\hat{n}$, and the comoving time $\tau$. The evolution of this distribution is determined by the Boltzmann equation,

$$ \frac{df}{d\tau} = \frac{\partial f}{\partial \tau} + \frac{dx^i}{d\tau} \frac{\partial f}{\partial x^i} + \frac{dq}{d\tau} \frac{\partial f}{\partial q} + \frac{d\hat{n}}{d\tau} \cdot \frac{\partial f}{\partial \hat{n}} = C[f], \quad (3.1) $$

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,D_i} \ll 1$. The collision terms for the mother and daughter particles are then given by,

$$ C_M = -\frac{a^2}{2\epsilon_M} \int \prod_i \frac{d^3\vec{q}_i}{2\epsilon_{D_i}} |\mathcal{M}|^2 (2\pi)^4 \delta^{(4)}(\vec{q}_M - \Sigma_i \vec{q}_{D_i}) f_M(q_M), \quad (3.2) $$

$$ C_{Dj} = +\frac{a^2}{2\epsilon_{Dj}} \int \frac{d^3\vec{q}_M}{2\epsilon_M} \prod_{i \neq j} \frac{d^3\vec{q}_i}{2\epsilon_{D_i}} |\mathcal{M}|^2 (2\pi)^4 \delta^{(4)}(\vec{q}_M - \Sigma_i \vec{q}_{D_i}) f_M(q_M). \quad (3.3) $$

where $\epsilon_S \equiv (q_S^2 + m_S^2 a^2)^{1/2}$ represents the comoving energy of the species $S(\equiv M, D_i)$ and $d^3\vec{q} \equiv d^3\vec{q}/(2\pi)^3$. From the definition of the decay width, the collision term for the mother particle can be simplified to

$$ C_M = -\frac{a\Gamma}{\gamma} f_M, \quad (3.4) $$
where $\Gamma$ denotes the decay width in the rest frame of the decaying particle, and the relativistic boost factor $\gamma \equiv \sqrt{q_M^2 + M^2a^2/(Ma)}$ 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, \hat{n}, \vec{x}, \tau) = f_S^0(q_S, \tau) + \Delta f_S(q_S, \hat{n}, \vec{x}, \tau), \quad S = M, D_i. \qquad (3.5)$$

## 3.1 Background: zeroth order

Treating $\Delta f_M$ and fluctuations about the gravitational background as higher order perturbations, the zeroth order Boltzmann equations for $f_M^0$ arising from eq. (3.1) take the form,

$$\frac{\partial f_M^0}{\partial \tau} = -a \frac{\Gamma}{\gamma} f_M^0. \qquad (3.6)$$

The formal solution to $f_M^0(q, \tau)$ from the differential equations in eq. (3.6) is given by,

$$f_M^0(q, \tau) = f_i(q) e^{-\Gamma \int_{\tau_i}^{\tau} \frac{a}{\gamma(a)} d\tau'}, \qquad (3.7)$$

where $\tau_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,

$$\Gamma \int_{\tau_i}^{\tau} \frac{a}{\gamma(a)} \frac{d\tau'}{\gamma(a)} \Big|_{\tau_i}^{t} - \Gamma \int_{\tau_i}^{t} dt' \frac{d}{dt'} \left( \frac{1}{\gamma(a)} \right), \qquad (3.8)$$

where we have used $ad\tau = dt$. It is computationally demanding to solve the integral for general $a(\tau)$. However, the behavior of the exponential factor is rather simple: the exponential is close to 1 when $\tau$ is smaller than the mother lifetime $\sim \gamma/\Gamma a$, and $f_M$ no longer contributes when $\tau$ is much larger than the mother lifetime. The only time that the exponential factor exhibits a non-trivial $a$-dependence is when $\tau \sim \gamma/\Gamma a$. Since our focus is on decays in the non-relativistic regime, so that $\gamma(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 $\gamma(a)$, can be neglected in favor of the first term. This allows us to approximate the exponent as

$$\Gamma \int_{\tau_i}^{\tau} \frac{a}{\gamma(a)} \frac{d\tau'}{\gamma(a)} \approx \frac{\Gamma t}{\gamma(a)}. \qquad (3.9)$$

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

$$f_M^0(q, \tau) \approx f_i(q) e^{-\frac{\Gamma}{\gamma(a)}t}. \qquad (3.10)$$

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) = \delta(q_M)N_{Mi}/(4\pi q_M^2)$, where $N_{Mi}$ represents the initial
comoving number density of mother particles. Since this distribution is localized entirely at $q_M = 0$, the boost factor $\gamma(a) = 1$. Then eq. (3.10) reduces to the known result for decaying cold dark matter [70–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_{\nu 0}} + 1)$, leading to

$$f_M^0 = \frac{1}{e^{q/T_{\nu 0}} + 1} \exp \left( -\frac{\Gamma}{\gamma} t \right). \qquad (3.11)$$

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 $d^3\vec{q}_{Di}\epsilon_{Di}$ and adding them up. The resulting total integrated collision term for the daughter species is given by,

$$\begin{align}
\sum_j \int d^3 \vec{q}_{Dj} \epsilon_{Dj} C_{Dj}^0 &= a^2 \int \frac{d^3 \vec{q}_M}{2\epsilon_M} \prod_i \frac{d^3 \vec{q}_{Di}}{2\epsilon_{Di}} \left( \sum_j \epsilon_{Dj} \right) |\mathcal{M}|^2 (2\pi)^4 \delta^{(4)}(\vec{q}_M - \sum_i \vec{q}_{Di}) f_M^0(q_M), \nonumber \\
&= a^2 \Gamma M \int d^3 \vec{q}_M f_M^0. \tag{3.12}
\end{align}$$

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 $\epsilon_M/\gamma = Ma$ 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 $N_M$

$$E_D \equiv \sum_i \int dq_{Di} q_{Di}^3 f_{Di}, \qquad N_M \equiv \int dq_M q_M^2 f_M. \qquad (3.13)$$

Since the daughter particles constitute massless radiation, we can rewrite the expression for the evolution of the daughter distribution in eq. (3.1) in terms of the background daughter energy density $\bar{\rho}_D \equiv 4\pi a^{-4} E_D^0$ and the background mother number density $\bar{n}_M \equiv 4\pi a^{-3} N_M^0$, where $E_D^0$ and $N_M^0$ are defined as in eq. (3.13) after expanding out $f_M$ and $f_{Di}$ as in eq. (3.5),

$$\frac{\partial \bar{\rho}_D}{\partial \tau} + 4aH\bar{\rho}_D = a\Gamma M \bar{n}_M. \qquad (3.14)$$

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 $\epsilon_M/\gamma = Ma$ that was used in obtaining eq. (3.12).
## 3.2 Perturbations: first order

In the synchronous gauge, the metric perturbations can be parametrized as,

$$ds^2 = a(\tau)^2 [-d\tau^2 + (\delta_{ij} + H_{ij}) dx^i dx^j], \quad (3.15)$$

where $d\tau = dt/a(\tau)$ 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,

$$H_{ij}(\vec{k}, \tau) = \hat{k}_i \hat{k}_j h(\vec{k}, \tau) + \left( \hat{k}_i \hat{k}_j - \frac{1}{3} \delta_{ij} \right) 6\eta(\vec{k}, \tau), \quad (3.16)$$

where $\vec{k}$ is conjugate to $\vec{x}$ and $\hat{k}$ is the unit vector. In Fourier space the first order terms in eq. (3.1) for the mother particle can be collected as,

$$\Delta f'_{M} + i \frac{qk}{\epsilon_M} P_1(\mu) \Delta f_M + q \frac{\partial f_M^0}{\partial q} \left[ -\frac{h'}{6} - \frac{P_2(\mu)}{3}(h' + 6\eta') \right] = -a^2 \frac{\Gamma M}{\epsilon_M} \Delta f_M, \quad (3.17)$$

where $\mu \equiv \hat{k} \cdot \hat{n}$ and $P_l(\mu)$ are the Legendre polynomials.

As usual, we can expand the angular dependence of the perturbations as a series in Legendre polynomials,

$$X(\dots, \vec{k}, \hat{n}) = \sum_{l=0}^{\infty} (-i)^l (2l+1) X_l(\dots, k) P_l(\hat{k} \cdot \hat{n}). \quad (3.18)$$

Here $X$ represents any of the perturbations $\Delta f_{M,D_j}$, $\Delta E_D$ or $\Delta 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,

$$\begin{align*}
\Delta f'_{M(0)} &= -\frac{qk}{\epsilon_M} \Delta f_{M(1)} + \frac{h'}{6} q \frac{\partial f_M^0}{\partial q} - \frac{a^2 \Gamma M}{\epsilon_M} \Delta f_{M(0)}, \\
\Delta f'_{M(1)} &= \frac{qk}{3\epsilon_M} (\Delta f_{M(0)} - 2\Delta f_{M(2)}) - \frac{a^2 \Gamma M}{\epsilon_M} \Delta f_{M(1)}, \\
\Delta f'_{M(2)} &= \frac{qk}{5\epsilon_M} (2\Delta f_{M(1)} - 3\Delta f_{M(3)}) - \left(\frac{1}{15}h' + \frac{2}{5}\eta'\right) q \frac{\partial f_M^0}{\partial q} - \frac{a^2 \Gamma M}{\epsilon_M} \Delta f_{M(2)}, \\
\Delta f'_{M(l)} &= \frac{qk}{(2l+1)\epsilon_M} [l\Delta f_{M(l-1)} - (l+1)\Delta f_{M(l+1)}] - \frac{a^2 \Gamma M}{\epsilon_M} \Delta f_{M(l)}, \quad l \ge 3. \quad (3.19)
\end{align*}$$

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 $\int d^3\vec{q}_{D_j} q_{D_j} P_l(\mu_{D_j})$ on both sides of eq. (3.1) for each daughter particle and add them
up. The collision term becomes

$$ \sum_{l'} (-i)^{l'} (2l' + 1) a^2 \int \frac{d^3 \vec{q}_M}{2\epsilon_M} \prod_i \frac{d^3 \vec{q}_{Di}}{2\epsilon_{Di}} \left( \sum_j \epsilon_{Dj} \right) (2\pi)^4 \\ \times \delta^{(4)}(\vec{q}_M - \sum_i \vec{q}_{Di}) \Delta f_{M(l')} P_l(\mu_{D_j}) P_{l'}(\mu_M). $$

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/(Ma)$ 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 ($\hat{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/(Ma)$ and assume that $\mu_M$ and $\mu_{D_j}$ are uncorrelated. In this case, the angular integrals over the Legendre polynomials can be performed independently, so that

$$ \sum_{l'} (-i)^{l'} (2l' + 1) \int_{-1}^{1} d\mu_{D_j} P_l(\mu_{D_j}) \int_{-1}^{1} d\mu_M P_{l'}(\mu_M) = 0 \quad (\text{if } l \text{ or } l' > 0). \quad (3.20) $$

This implies that in the daughter equations, only the zeroth moment of the source term from the mother particle decay ($\Delta f_{M(0)}$) survives in the limit of non-relativistic decay. The source term shows up in the equation for $\Delta 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_M^0$ replaced by the perturbation $\Delta f_{M(0)}$. Therefore, the Boltzmann hierarchy for the daughter energy perturbations, $\Delta E_{D(l)}$, in terms of the $\Delta N_{M(l)}$ and the metric perturbations $h$ and $\eta$ is given by,

$$ \begin{aligned} \Delta E'_{D(0)} &= -k \Delta E_{D(1)} - \frac{2}{3} h' E_D^0 + a^2 M \Gamma \Delta N_{M(0)}, \\ \Delta E'_{D(1)} &= \frac{k}{3} \Delta E_{D(0)} - \frac{2k}{3} \Delta E_{D(2)}, \\ \Delta E'_{D(2)} &= \frac{2k}{5} \Delta E_{D(1)} - \frac{3k}{5} \Delta E_{D(3)} + \frac{4}{15} E_D^0 (h' + 6\eta'), \\ \Delta E'_{D(l)} &= \frac{k}{2l+1} [l\Delta E_{D(l-1)} - (l+1)\Delta E_{D(l+1)}], \quad l \ge 3. \end{aligned} \quad (3.21) $$

Similar equations can also be found in [65, 66, 68, 71]. Again, we neglect the source terms with $\Delta N_{M(l>0)}$ due to the additional $q_M/(Ma)$ 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].

## 4 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 section 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_{\ell}^{\phi\phi}$ 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 section 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.

## 4.1 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, $\sum m_{\nu}$, and the logarithm of the decay width, $\log_{10}\Gamma_{\nu}$. In our analysis, we fix the following cosmological parameters to their central values from the Planck 2015 TT, TE, EE+low-P data: {$\omega_b = 0.022032$, $\omega_{cdm} = 0.12038$, $\ln(10^{10}A_s) = 3.052$, $n_s = 0.96229$, $\tau_{reio} = 0.0648$}. The impact of neutrino masses on the matter power spectrum looks different depending on whether $\theta_s$ or $H_0$ is kept fixed [75]. This is because, to keep $\theta_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 $\theta_s$ gives a better reflection of the constraining effects of a combined analysis of CMB+LSS data, since CMB data pins $\theta_s$ down very precisely. In the following, we will show results with either $H_0 = 67.56$ km/s/Mpc or $100 \times \theta_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

$$P_{cb}(k) = \left\langle \frac{\delta\rho_{cb}}{\bar{\rho}_{cb}} \frac{\delta\rho_{cb}}{\bar{\rho}_{cb}} \right\rangle, \qquad (4.1)$$

where $\bar{\rho}_{cb}$ ($\delta\rho_{cb}$) is the average (perturbation) of the sum of CDM and baryon energy densities.<fnref n="4" /> In figure 2, we display the residuals of $P_{cb}$ (left) and the CMB lensing potential (right) with respect to the case of massless neutrinos for $\sum m_{\nu}$ fixed at 0.25 eV, keeping the value of $H_0$ fixed. We compare three different values of $\Gamma_{\nu}$ and the limiting case of stable neutrinos. The curves run from top to bottom in order of decreasing $\Gamma_{\nu}$. The analytic results are shown as dashed lines in the plot, and are seen to agree reasonably well at large $k$ or $\ell$ 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 figure 2 that in this region there are no qualitative features in $P_{cb}|_{z=0}$ or $C_{\ell}^{\phi\phi}$ that would allow unstable neutrinos to be distinguished from

<footnote n="4">⁴Note that this is different from the matter power spectrum conventionally defined as $P_m = \langle[(\delta\rho_{cb} + \delta\rho_{\nu})/(\bar{\rho}_{cb} + \bar{\rho}_{\nu})]^2\rangle$, which differs from $P_{cb}$ by an extra factor $[\bar{\rho}_{cb}/(\bar{\rho}_{cb} + \bar{\rho}_{\nu})]^2$.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Plots of the fractional difference in the CDM+Baryon power spectrum $P_{cb}$ (left) and CMB-lensing potential $C_\ell^{\phi\phi}$ (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, $\Gamma_\nu = 10^{4.0}, 3.5, 3.0$ (km/s/Mpc) (top to bottom), and also the stable neutrino scenario, holding $\sum m_\nu = 0.25$ eV and $H_0$ fixed. The dashed lines represent the corresponding analytic estimates from section 4.2.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** The fractional differences in the CMB-lensing potential $C_\ell^{\phi\phi}$ (left), CDM+Baryon power spectrum $P_{cb}$ (right) for an unstable (purple) and a stable (blue) neutrino scenario with respect to the case of massless neutrinos (black) at fixed $H_0$. The grey regions show the $1\sigma$ uncertainties from Planck and SDSS DR7 respectively.</figcaption></figure>

stable ones. Although $P_{cb}|_{z=0}$ and $C_\ell^{\phi\phi}$ are more suppressed in the stable case, as expected, this effect can be mimicked if the stable 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_\ell^{\phi\phi}$.

In figure 3 we show an explicit example of the degeneracy between mass and lifetime in the values of $P_{cb}$ and $C_\ell^{\phi\phi}$ at fixed $H_0$. We consider a model with stable neutrinos of mass $\sum m_\nu = 0.2$ eV, and a different model with unstable neutrinos of mass $\sum m_\nu = 0.36$ eV and width $\Gamma_\nu = 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 section 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 \approx 3$ [76]. We will explore the possibility of breaking the degeneracy by using next generation measurements at different redshifts in future work.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** The fractional differences in the CMB-lensing potential $C_{\ell}^{\phi\phi}$ (top left), CDM+Baryon power spectrum $P_{cb}$ (top right), $C_{\ell}^{TT}$ (bottom left), and $C_{\ell}^{EE}$ (bottom right) for an unstable (purple) and a stable (blue) neutrino scenario with respect to the case of massless neutrinos (black) at fixed $\theta_s$. The grey regions show the $1\sigma$ uncertainties from Planck and SDSS DR7 respectively.</figcaption></figure>

Finally, we show in figures 4 the effects of neutrino masses and decay at fixed $\theta_s$ on $P_{cb}$, $C_{\ell}^{\phi\phi}$ and $C_{\ell}^{TT,EE}$. This fixes the peak locations in the CMB power spectra and only generates negligible deviations away from the massless neutrino case in $C_{\ell}^{TT,EE}$ [75]. The same choices of parameters, however, do generate sizeable deviations in $C_{\ell}^{\phi\phi}$ and $P_{cb}$ away from the massless neutrino case that are close to the current sensitivities. This demonstrates that as expected, for sub-eV $\sum m_\nu$, 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 $\theta_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.

## 4.2 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_{\ell}^{\phi\phi}$) 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 $\delta_i = \delta\rho_i/\bar{\rho}_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

$$k^2\phi \approx -4\pi Ga^2(\delta_{cb}\bar{\rho}_{cb} + \delta_{\nu}\bar{\rho}_{\nu}). \quad (4.2)$$

Here $\phi$ is the metric perturbation in the conformal Newtonian gauge [73].<fnref n="5" /> We 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 $\delta_{\nu} \ll \delta_{cb}$ for perturbation modes that enter the horizon before $z_{nr}$, we can write,

$$k^2\phi \approx -\frac{6}{\tau} \left( 1 - \frac{\bar{\rho}_{\nu}(\tau)}{\bar{\rho}_{\text{tot}}(\tau)} \right) \delta_{cb}, \quad (4.3)$$

where $\tau$ is the comoving time and $\bar{\rho}_{\text{tot}} \equiv \bar{\rho}_{cb} + \bar{\rho}_{\nu}$. Inserting this expression into the Boltzmann equation for CDM perturbations yields,

$$\ddot{\delta}_{cb} + \frac{2}{\tau}\dot{\delta}_{cb} - \frac{6}{\tau^2}(1 - f_{\nu}(\tau))\delta_{cb} = 0, \qquad f_{\nu}(\tau) = \frac{\bar{\rho}_{\nu}(\tau)}{\bar{\rho}_{\text{tot}}(\tau)}. \quad (4.4)$$

where the dots represent derivatives with respect to $\tau$. 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_{\nu}$ ($\ll 1$). We look for a solution of the form,

$$\delta_{cb} = \delta_{cb,i}h(\tau) \left(\frac{\tau}{\tau_i}\right)^2 \exp\left[-\frac{6}{5}\int_{\tau_i}^{\tau} \frac{d\hat{\tau}}{\hat{\tau}} f_{\nu}(\hat{\tau})\right] \quad (4.5)$$

where now the function $h(\tau)$ is to be determined. Inserting this expression into eq. (4.4) and dropping the term proportional to $f_{\nu}^2$, we obtain the following differential equation for $h(\tau)$.

$$\tau \ddot{h} + 6\dot{h} - \frac{6}{5}h \dot{f} = 0. \quad (4.6)$$

Thus far we have not made any assumption about the redshift dependence of $f_{\nu}$. For massless or ultrarelativistic neutrinos in the matter dominated era, we have

$$f_{\nu}(\tau) = f_{\nu,i} \left(\frac{\tau_i}{\tau}\right)^2, \quad (4.7)$$

In this case we can solve for the function $h(\tau)$ as,

$$h^{\prime\prime\nu}(\tau) = \exp\left[k \int_{\tau_i}^{\tau} d\hat{\tau} \dot{f}_{\nu}^{\prime\prime\nu}(\hat{\tau})\right] = \exp\left[\frac{2}{5} \left(f_{\nu}^{\prime\prime\nu}(\tau_i) - f_{\nu}^{\prime\prime\nu}(\tau)\right)\right]. \quad (4.8)$$

<footnote n="5">⁵We use the metric $ds^2 = a^2(\tau)[-(1+2\psi)d\tau^2 + (1-2\phi)\delta_{ij}dx^i dx^j]$ and approximate $\psi = -\phi$, ignoring the small correction arising from the presence of free streaming radiation.</footnote>
This leads to the following approximate solution for perturbations in the case of massless or ultrarelativistic neutrinos,

$$ \delta_{cb}^{\eta\eta\nu}(\tau) = \delta_{cb,i} \left(\frac{\tau}{\tau_i}\right)^2 \exp\left[-\frac{6}{5} \int_{\tau_i}^{\tau} \frac{d\hat{\tau}}{\hat{\tau}} f_{\nu}^{\eta\eta\nu}(\hat{\tau})\right] h^{\eta\eta\nu}(\tau). \quad (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

$$ h^{m\nu}(\tau) = \exp\left\{\frac{2}{5}\left[f_{\nu}^{\eta\eta\nu}(\tau_i) - f_{\nu}^{\eta\eta\nu}(\min(\tau, \tau_{nr}))\right]\right\}. \quad (4.10) $$

The result is almost identical to eq. (4.8) since in both cases the exponent is dominated by $f_ν^{\eta\eta\nu}(\tau_i)$ (and $f_ν^{\eta\eta\nu}$ after $τ > τ_{nr}$ is much smaller than the expansion parameter $f_ν^{m\nu}$)

$$ h^{m\nu}(\tau) \approx h^{\eta\eta\nu}(\tau). \quad (4.11) $$

This means that the solution in the case of massive neutrinos can be approximated as,

$$ \delta_{cb}^{m\nu}(\tau) = \delta_{cb,i} \left(\frac{\tau}{\tau_i}\right)^2 \exp\left[-\frac{6}{5} \int_{\tau_i}^{\tau} \frac{d\hat{\tau}}{\hat{\tau}} f_{\nu}^{m\nu}(\hat{\tau})\right] h^{\eta\eta\nu}(\tau). \quad (4.12) $$

Then the ratio of the perturbations in the two cases is given by,

$$ \frac{\delta_{cb}^{m\nu}(\tau)}{\delta_{cb}^{\eta\eta\nu}(\tau)} = \exp\left[-\frac{6}{5} \int_{\tau_i}^{\tau} \frac{d\hat{\tau}}{\hat{\tau}} \left(f_{\nu}^{m\nu}(\hat{\tau}) - f_{\nu}^{\eta\eta\nu}(\hat{\tau})\right)\right]. \quad (4.13) $$

This ratio can be expressed in terms of the scale factor as,

$$ \frac{\delta_{cb}^{m\nu}(a)}{\delta_{cb}^{\eta\eta\nu}(a)} \approx \frac{\delta_{cb}^{m\nu}(a_i)}{\delta_{cb}^{\eta\eta\nu}(a_i)} \exp\left[-\frac{3}{5} \int_{a_i}^{a} \frac{da}{a} \frac{\hat{\rho}_{\nu}(a)}{\bar{\rho}_{\text{tot}}(a)}\right], \quad (4.14) $$

where $\hat{\rho}_{\nu}(a) \equiv \bar{\rho}_{\nu,m_{\nu}}(a) - \bar{\rho}_{\nu,\eta\eta\nu}(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$, $\hat{\rho}_{\nu}(a)/\bar{\rho}_{\text{tot}}(a) = \bar{\rho}_{\nu,m_{\nu}}/\bar{\rho}_{\text{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,

$$ \frac{\delta_{cb}^{m\nu}(a)}{\delta_{cb}^{\eta\eta\nu}(a)} \propto \left(\frac{a}{a_i}\right)^{-\frac{3}{5} \frac{\hat{\rho}_{\nu,m\nu}}{\bar{\rho}_{\text{tot}}}}. \quad (4.15) $$

We can improve on this estimate by incorporating a more precise expression for the neutrino energy in eq. (4.14),

$$ \hat{\rho}_{\nu}(a) = 4\pi a^{-4} \int_{0}^{\infty} dq \, q^2 \left( \sqrt{q^2 + m_{\nu}^2 a^2} - q \right) f(q). \quad (4.16) $$

Here $q = a p_{\nu}$ denotes the neutrino's conformal momentum, and $f(q) = [e^{q/T_{\nu 0}} + 1]^{-1}$ represents the momentum distribution of neutrinos. $\hat{\rho}_{\nu}(a)$ exhibits non-trivial redshift
dependence since the neutrino energy goes from being radiation-like to being matter-like. In figure 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 \times 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 $\delta_{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,

$$ \frac{P_{cb,m_\nu}}{P_{cb,\bar{\eta}\eta_\nu}} \approx \left( \frac{\delta_{cb}^{m_\nu}(a_f)}{\delta_{cb}^{\bar{\eta}\eta_\nu}(a_f)} \right)^2 . \qquad (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.<fnref n="6" />

We now turn our attention to the effects of massive neutrinos on CMB lensing. The difference in the density perturbation $\delta\rho_{cb}$ between the massive and massless neutrino scenarios results in a change in the gravity perturbation $\phi$. The photons are therefore deflected differently in the CMB lensing process. The correlation function of the lensing potential, $C_{\ell}^{\phi\phi} \sim \langle\phi\phi\rangle$, parameterizes the size of angular deflection of CMB photons. The ratio of $C_{\ell}^{\phi\phi}$ in the massive neutrino case to that in the massless case can be approximated using Limber's formula [77, 78]

$$ \frac{C_{\ell, m_\nu}^{\phi\phi}}{C_{\ell, \bar{\eta}\bar{\eta}_\nu}^{\phi\phi}} \approx \frac{\chi_*^{m_\nu} \int_0^1 dx \, \phi_{m_\nu}^2 \left( \frac{l}{x^{m_\nu} \chi_*} \right) (1-x)^2}{\chi_*^{\bar{\eta}\bar{\eta}_\nu} \int_0^1 dx \, \phi_{\bar{\eta}\bar{\eta}_\nu}^2 \left( \frac{l}{x^{\bar{\eta}\bar{\eta}_\nu} \chi_*} \right) (1-x)^2}, \quad x \equiv \frac{\tau_f - \tau}{\tau_f - \tau_*}, \quad \chi_* \equiv \tau_f - \tau_* . \qquad (4.18) $$

Here $\tau_* \approx 2.8 \times 10^2$ Mpc is the conformal time at last scattering, while $\tau_f \approx 1.4 \times 10^4$ Mpc is the conformal time today. The value of $\tau_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 $\chi_*$ between the two scenarios can be neglected. Then, the difference between $C_{\ell}^{\phi\phi}$ in the two cases primarily arises from differences in the evolution of $\phi$.

According to the Einstein eq. (4.2), the ratio of $\phi$ between the two scenarios for large $l$ modes at a given value of the scale factor is,

$$ \frac{\phi_{m_\nu}(a)}{\phi_{\bar{\eta}\bar{\eta}_\nu}(a)} \approx \frac{\delta_{cb}^{m_\nu}(a)}{\delta_{cb}^{\bar{\eta}\bar{\eta}_\nu}(a)}. \qquad (4.19) $$

<footnote n="6">⁶We can also use $a_f = a g(a)$ with the growth function $g(a)$ for a reasonable approximation [5].</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Evolution of the ratio of the CDM+baryon density perturbation with respect to the case of a massless neutrino, $\delta_{cb}^{m_\nu}/\delta_{cb}^{\bar{\nu}_\nu}$. The results are shown for the case of a single massive neutrino with $m_\nu = 60$ meV. All the solid curves are obtained from numerical calculations using the modified CLASS code discussed in section 3. The black curve is for the stable neutrino scenario, and the blue (orange) curve is for the neutrino with decay rate $\Gamma_\nu = 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_\nu$ and $a_f = a_{\text{dec}}$. As we see, eq. (4.14) provides a much better approximation to the full numerical result.</figcaption></figure>

Since $C_\ell^{\phi\phi}$ receives its dominant contribution close to $z \approx 3$ [76], we can estimate the ratio of the $C_\ell^{\phi\phi}$ as,

$$ \frac{C_{\ell, m_\nu}^{\phi\phi}}{C_{\ell, \bar{\nu}_\nu}^{\phi\phi}} \approx \left( \frac{\delta_{cb}^{m_\nu}}{\delta_{cb}^{\bar{\nu}_\nu}} \right)^2 \bigg|_{z=3} . \qquad (4.20) $$

Based on a very similar analysis, we can predict the suppression of $P_{cb}(k)$ and $C_\ell^{\phi\phi}$ for large $k$ and $\ell$ 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_{\text{dec}}$, which is implicitly determined by the equation,

$$ \Gamma_{\nu} = H(a_{\text{dec}}). \qquad (4.21) $$

The difference in energy density $\hat{\rho}_{\nu}$ 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,

$$ \hat{\rho}_{\nu}(a) = 4\pi a^{-4} \int_{0}^{\infty} dq \, q^2 \left( \sqrt{q^2 + m_{\nu}^2 a^2} - q \right) f(q), \quad a < a_{\text{dec}}. \qquad (4.22) $$

In the instantaneous decay approximation, the energy density in non-relativistic neutrinos is immediately transferred into radiation energy at $a_{\text{dec}}$. It subsequently redshifts with an extra $(a_{\text{dec}}/a)$ factor as compared to a non-relativistic neutrino, so that

$$ \hat{\rho}_{\nu}(a) = 4\pi a^{-4} \int_{0}^{\infty} dq \, q^2 \left[ \sqrt{q^2 + m_{\nu}^2 a^2} \left( \frac{a_{\text{dec}}}{a} \right) - q \right] f(q), \quad a \ge a_{\text{dec}}. \qquad (4.23) $$

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_{\ell}^{\phi\phi}$ in the limit of large $k$ and $\ell$ can be obtained from eqs. (4.17) and (4.20)

In figure 5, we show the ratio of $\delta_{cb}$ from the numerical calculation described in section 3 for both the decaying (blue and orange) and stable (black) neutrinos. The plots are for a single massive neutrino with $m_{\nu} = 60$ meV (upper) and 80 meV (lower), and a decay rate $\Gamma_{\nu} = 10^4 (10^3)$ km/s/Mpc for the blue (orange) curves. In this scenario, more than 80% of the neutrinos have momenta $p_{\nu} < m_{\nu}$ after $a > 0.012$ ($a > 0.0096$) for $m_{\nu} = 60$ (80) meV neutrino. It is at this point, when most of the neutrinos have become non-relativistic, that the major suppression of $\delta_{cb}$ begins. During this period the $\delta_{cb}$ ratio drops with the power described in eq. (4.15) (grey line). The blue (orange) dotted lines give the value of the $\delta_{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 $\delta_{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 $\delta_{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 $\delta_{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.

# 5 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.
## 5.1 The data and analysis pipeline

Our analysis makes use of various combinations of the following datasets.

* CMB: we include Planck 2015 CMB high-$\ell$ TT, TE, and EE and low-$\ell$ 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

$$f\sigma_8 \equiv \frac{[\sigma_8^{(vd)}(z)]^2}{\sigma_8^{(dd)}(z)}, \qquad (5.1)$$

where $\sigma_8^{(vd)}$ measures the smoothed density-velocity correlation, analogous to $\sigma_8 \equiv \sigma_8^{(dd)}$ 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-V3$^8$ [87, 88], we run Monte Carlo Markov chain analyses using the Metropolis-Hastings algorithm assuming flat priors on all parameters. Our ΛCDM parameters are,

$$\{\omega_{cdm}, \omega_b, \theta_s, \ln(10^{10} A_s), n_s, \tau_{reio}\},$$

to which we add the sum of neutrino masses $\sum m_\nu$ and the logarithm of the neutrino lifetime $\text{Log}_{10}\Gamma_\nu$. 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.

<footnote>$^7$More recent measurements are not yet available in MONTEPYTHON-V3. These could naturally make the bounds presented here slightly stronger.</footnote>

<footnote>$^8$https://github.com/brinckmann/montepython_public</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Posterior distributions of $\sum m_\nu$ and $\log_{10}\Gamma_\nu$ for each dataset. Small decay rate $\log_{10}\Gamma_\nu/(\text{km/s/Mpc}) \in [0, 3]$ are shown in the left panel, while large decay rate $\log_{10}\Gamma_\nu/(\text{km/s/Mpc}) \in [3, 5.5]$ are shown in the right panel.</figcaption></figure>

## 5.2 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 $\sum m_\nu < 0.28 \text{ 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 $\sum m_\nu < 0.25 \text{ eV (95\% C.L.)}$. This constraint when LSS data is included is also in good agreement with what is reported in ref. [90].

In figure 6 we show the 1D and 2D marginalized posterior distribution of $\sum m_\nu$ and $\log_{10}\Gamma_\nu$ for both datasets, cutting the parameter space between small decay rate $\log_{10}\Gamma_\nu/(\text{km/s/Mpc}) \in [0, 3]$ (left panel) and large decay rate $\log_{10}\Gamma_\nu/(\text{km/s/Mpc}) \in [3, 5.5]$ (right panel) to accelerate convergence. Strikingly, once the neutrino lifetime is let free to vary, the constraint on $\sum m_\nu$ is driven by our prior on $\log_{10}\Gamma_\nu$. 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}\Gamma_\nu/(\text{km/s/Mpc}) < 2.5$, but relaxes to $\sum m_\nu < 0.9 \text{ 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 $\sum m_\nu$ 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 refs. [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.
# 6 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 well-motivated 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 $\sum m_\nu$ 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].

# Acknowledgments

We would like to thank warmly Thejs Brinckmann for his help with MONTEPYTHON-V3, and Marilena Loverde and Gustavo Marques-Tavares for useful discussions. ZC, AD, PD and YT are supported in part by the National Science Foundation under Grant Number PHY-1620074. ZC would like to thank the Fermilab Theory Group for hospitality during the completion of this work. ZC's stay at Fermilab was supported by the Fermilab Intensity Frontier Fellowship and by the Visiting Scholars Award #17-S-02 from the Universities Research Association. VP thanks the Johns Hopkins University and particularly Marc Kamionkowski for his hospitality during completion of this work. This work has been made possible thanks to the facilities offered by the Université Savoie Mont-Blanc MUST computing center.

# A A model of massive neutrino decay into dark radiation

In this appendix we present a simple, realistic model in which massive neutrinos decay into invisible dark radiation on timescales of order the age of the universe. To illustrate
the main features of the model, we first consider a simplified version with just a single flavor of SM neutrino, denoted by $\nu$, and two singlet right-handed neutrinos, labelled as $n$ and $n'$. The model also contains two complex scalars, labelled as $\Phi$ and $\Phi'$. We introduce $U(1)_n \times U(1)_{n'}$ global symmetries that act on the right-handed neutrinos. While $n$ and $\Phi$ carry equal and opposite charges under $U(1)_n$, $n'$ and $\Phi'$ are neutral under this symmetry. Similarly, $n'$ and $\Phi'$ carry equal and opposite charges under $U(1)_{n'}$, while $n$ and $\Phi$ are neutral. Then the part of the Lagrangian responsible for generating the neutrino masses takes the form,

$$ -\mathcal{L} \supset \frac{y}{\Lambda} \bar{L} \tilde{H} n \Phi + \frac{y'}{\Lambda} \bar{L} \tilde{H} n' \Phi' + \text{H.c.} \qquad (A.1) $$

Here $L$ represents the SM lepton doublet and $\tilde{H} = i\sigma_2 H^*$, where $H$ denotes the SM Higgs doublet. $\Lambda$ is a UV mass scale while $y$ and $y'$ are coupling constants. Although this Lagrangian is nonrenormalizable, it can be interpreted as the low energy description of a renormalizable theory after particles with masses of $\Lambda$ have been integrated out. For example, consider the renormalizable Lagrangian,

$$ -\mathcal{L} = \tilde{y} \bar{L} \tilde{H} N + M_N N N^C + \tilde{\lambda} n N^c \Phi + \tilde{y}' \bar{L} \tilde{H} N' + M'_N N' N^{c'} + \tilde{\lambda}' n' N^{c'} \Phi' + \text{H.c.} \qquad (A.2) $$

Terms of the form shown in eq. (A.1) are obtained after the heavy fermions $N, N^c, N'$ and $N^{c'}$ have been integrated out.

Once the scalars $\Phi, \Phi'$ and the SM Higgs each acquire a vacuum expectation value (VEV), we obtain Dirac masses for the SM neutrino,

$$ \begin{aligned} -\mathcal{L} & \supset \frac{y f v}{2 \Lambda} \bar{\nu} n + \frac{y' f' v}{2 \Lambda} \bar{\nu} n' + \text{H.c.} \\ & = m \bar{\nu} n_h + \text{H.c.} \end{aligned} \qquad (A.3) $$

Here $\frac{f}{\sqrt{2}}, \frac{f'}{\sqrt{2}}$ and $\frac{v}{\sqrt{2}}$ denote the VEVs of $\Phi, \Phi'$ and $H$ respectively. The SM neutrino acquires a mass $m = \sqrt{(yf)^2 + (y'f')^2}v/(2\Lambda)$. Its Dirac partner $n_h$ is one linear combination of $n$ and $n'$,

$$ \begin{pmatrix} n_h \\ n_l \end{pmatrix} = \begin{pmatrix} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{pmatrix} \begin{pmatrix} n \\ n' \end{pmatrix} ; \quad \cos \theta = \frac{yf}{\sqrt{(yf)^2 + (y'f')^2}}. \qquad (A.4) $$

It is clear from eq. (A.3) that the spectrum contains one massive Dirac neutrino and one massless singlet neutrino $n_l$.

Below the spontaneous symmetry breaking scales $f$ and $f'$, the Goldstone bosons can be parametrized as

$$ \Phi = \frac{f}{\sqrt{2}} e^{i\phi/f}, \quad \bar{\Phi}' = \frac{f'}{\sqrt{2}} e^{i\phi'/f'}, \qquad (A.5) $$

where $\phi$ and $\phi'$ denote the Goldstone bosons from $U(1)_n$ and $U(1)_{n'}$ respectively. The couplings of the Goldstone bosons in the low energy effective theory are dictated by the non-linearly realized global symmetries. To leading order in $1/f$ and $1/f'$, they are given by,

$$ -\mathcal{L} \supset i \frac{y f v \phi}{2 \Lambda} \bar{\nu} n + i \frac{y f' v \phi'}{2 \Lambda} \bar{\nu} n' + \text{H.c.} \qquad (A.6) $$
In the mass basis these interactions take the form,

$$ -\mathcal{L} \supset i m \bar{\nu} \left[ \left( \frac{\phi}{f} \cos^2 \theta + \frac{\phi'}{f'} \sin^2 \theta \right) n_h + \left( \frac{\phi'}{f'} - \frac{\phi}{f} \right) \sin \theta \cos \theta n_l \right] + \text{H.c.} \quad (A.7) $$

We see from this that the massive neutrino can decay into $n_l$ and either $\phi$ or $\phi'$. Its partial widths into these decay modes are given by,

$$ \Gamma(\nu \rightarrow n_l \phi) = \frac{m^3}{32\pi f^2}, \quad \Gamma(\nu \rightarrow n_l \phi') = \frac{m^3}{32\pi f'^2}, \quad (A.8) $$

where $\bar{f} \equiv f/(\cos\theta \sin\theta)$ and $\bar{f}' \equiv f'/(\cos\theta \sin\theta)$.

Now we move on to discuss the realistic case in which there are three flavors of SM neutrinos $\nu_\alpha$ ($\alpha = e, \mu, \tau$). We also introduce three flavors of the sterile neutrinos $n_\alpha$ and $n'_\alpha$, as well as a new scalar field $\Sigma_{\alpha\beta}$. The global symmetry in the neutrino sector is now extended to $SU(3)_L \times SU(3)_R \times U(1)_n \times U(1)_{n'}$. The charge assignments under $U(1)_n \times U(1)_{n'}$ are the same as before, but with all 3 flavors of $n_\alpha$ and $n'_\alpha$ now being charged under $U(1)_n$ and $U(1)_{n'}$ respectively. Under $SU(3)_L \times SU(3)_R$, the various fields transform as

$$ L \rightarrow U_L L \quad n \rightarrow U_R n \quad n' \rightarrow U_R n' \quad \Sigma \rightarrow U_L \Sigma U_R^\dagger, \quad (A.9) $$

where $U_L$ and $U_R$ are the rotation matrices of $SU(3)_L$ and $SU(3)_R$ respectively. The neutrino masses now arise from terms in the Lagrangian of the form,

$$ -\mathcal{L} \supset \frac{y}{\Lambda^2} \bar{L}_\alpha \tilde{H} \Sigma_{\alpha\beta} n_\beta \Phi + \frac{y'}{\Lambda^2} \bar{L}_\alpha \tilde{H} \Sigma_{\alpha\beta} n'_\beta \Phi' + \text{H.c.} \quad (A.10) $$

Once the $\Sigma$ field acquires a VEV, we can diagonalize its VEV $\langle\Sigma\rangle$ to obtain,

$$ -\mathcal{L} \supset \sum_i \left( \frac{y}{\Lambda^2} \bar{L}_i \tilde{H} \langle\Sigma\rangle_i n_i \Phi + \frac{y'}{\Lambda^2} \bar{L}_i \tilde{H} \langle\Sigma\rangle_i n'_i \Phi' \right) + \text{H.c.} \quad (A.11) $$

where the index $i$ runs over $i = 1, 2, 3$ and $\langle\Sigma\rangle_i$ denotes the $i$th eigenvalue of $\langle\Sigma\rangle$. The Lagrangian in eq. (A.11) can be viewed as three copies of eq. (A.1). After the scalars $\Phi, \Phi'$ and $H$ acquire VEVs, all three generations of $(n_i, n'_i)$ can be simultaneously transformed to the mass basis $(n_{hi}, n_{li})$ using the same orthogonal matrix,

$$ \begin{pmatrix} n_{hi} \\ n_{li} \end{pmatrix} = \begin{pmatrix} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{pmatrix} \begin{pmatrix} n_i \\ n'_i \end{pmatrix}, \quad (A.12) $$

where $\cos \theta$ is exactly the same as in eq. (A.4). Now the neutrino masses are given by,

$$ m_i = \sqrt{(yf)^2 + (y'f')^2} \frac{\langle\Sigma\rangle_i v}{2\Lambda^2}. \quad (A.13) $$

Assuming that the Goldstone bosons from $\Sigma$ 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 $\phi$ or $\phi'$. Following the discussion above, the total neutrino decay width is given by

$$ \Gamma_{\nu_i} = \Gamma(\nu_i \to n_{li}\phi) + \Gamma(\nu_i \to n_{li}\phi') = \frac{m_i^3}{32\pi f^2} + \frac{m_i^3}{32\pi f'^2}, \quad (A.14) $$

where $\bar{f}$ and $\bar{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, $\Gamma_{\nu_i}/\Gamma_{\nu_j} = m_i^3/m_j^3$. In the case of quasi-degenerate neutrinos, $m_1 \approx m_2 \approx m_3$, it is clear that all neutrinos have almost the same total width. Assuming $\bar{f} = \bar{f}'$, we find that the total width is of order $H_0$ for $\bar{f} \sim 10^5$ GeV and neutrino masses of order 0.1 eV,

$$ \frac{\Gamma_{\nu_i}}{H_0} \approx 1.3 \left( \frac{m_i}{0.1 \text{ eV}} \right)^3 \left( \frac{10^5 \text{ GeV}}{\bar{f}} \right)^2. \quad (A.15) $$

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–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' \gtrsim 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' \gtrsim 100$ keV [93–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' \sim 10^5$ GeV are easily consistent with all current bounds.

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.
References

[1] PLANCK collaboration, Planck 2018 results. VI. Cosmological parameters, arXiv:1807.06209 [INSPIRE].

[2] J.R. Bond, G. Efstathiou and J. Silk, Massive neutrinos and the large scale structure of the universe, Phys. Rev. Lett. 45 (1980) 1980 [INSPIRE].

[3] W. Hu, D.J. Eisenstein and M. Tegmark, Weighing neutrinos with galaxy surveys, Phys. Rev. Lett. 80 (1998) 5255 [astro-ph/9712057] [INSPIRE].

[4] Y.Y.Y. Wong, Neutrino mass in cosmology: status and prospects, Ann. Rev. Nucl. Part. Sci. 61 (2011) 69 [arXiv:1111.1436] [INSPIRE].

[5] J. Lesgourgues, G. Mangano, G. Miele and S. Pastor, Neutrino cosmology, Cambridge University Press, Cambridge U.K. (2018).

[6] PARTICLE DATA GROUP collaboration, Review of particle physics, Phys. Rev. D 98 (2018) 030001.

[7] M. Lattanzi and M. Gerbino, Status of neutrino properties and future prospects - Cosmological and astrophysical constraints, Front. in Phys. 5 (2018) 70 [arXiv:1712.07109] [INSPIRE].

[8] P.D. Serpico, Cosmological neutrino mass detection: the best probe of neutrino lifetime, Phys. Rev. Lett. 98 (2007) 171301 [astro-ph/0701699] [INSPIRE].

[9] P.D. Serpico, Neutrinos and cosmology: a lifetime relationship, J. Phys. Conf. Ser. 173 (2009) 012018 [INSPIRE].

[10] J.F. Beacom, N.F. Bell and S. Dodelson, Neutrinoless universe, Phys. Rev. Lett. 93 (2004) 121302 [astro-ph/0404585] [INSPIRE].

[11] Y. Farzan and S. Hannestad, Neutrinos secretly converting to lighter particles to please both KATRIN and the cosmos, JCAP 02 (2016) 058 [arXiv:1510.02201] [INSPIRE].

[12] S.T. Petcov, The Processes $\mu \to e\Gamma$, $\mu \to ee\bar{e}$, neutrino' $\to$ neutrino gamma in the Weinberg-Salam model with neutrino mixing, Sov. J. Nucl. Phys. 25 (1977) 340 [Erratum ibid. 25 (1977) 698] [INSPIRE].

[13] J.T. Goldman and G.J. Stephenson, Jr., Limits on the mass of the muon-neutrino in the absence of muon lepton number conservation, Phys. Rev. D 16 (1977) 2256 [INSPIRE].

[14] W.J. Marciano and A.I. Sanda, Exotic decays of the muon and heavy leptons in gauge theories, Phys. Lett. B 67 (1977) 303.

[15] B.W. Lee and R.E. Shrock, Natural suppression of symmetry violation in gauge theories: muon-lepton and electron lepton number nonconservation, Phys. Rev. D 16 (1977) 1444 [INSPIRE].

[16] P.B. Pal and L. Wolfenstein, Radiative decays of massive neutrinos, Phys. Rev. D 25 (1982) 766 [INSPIRE].

[17] R.N. Mohapatra and P.B. Pal, Massive neutrinos in physics and astrophysics, Lecture Notes in Physics volume 41, World Scientific, Singapore (1991).

[18] F. Boehm and P. Vogel, Physics of massive neutrinos, (1992).
[19] G.B. Gelmini and M. Roncadelli, Left-handed neutrino mass scale and spontaneously broken lepton number, Phys. Lett. B 99 (1981) 411.

[20] Y. Chikashige, R. Mohapatra and R. Peccei, Are there real goldstone bosons associated with broken lepton number?, Phys. Lett. B 98 (1981) 265.

[21] H.M. Georgi, S.L. Glashow and S. Nussinov, Unconventional model of neutrino masses, Nucl. Phys. B 193 (1981) 297 [INSPIRE].

[22] J. Valle, Fast neutrino decay in horizontal majoron models, Phys. Lett. B 131 (1983) 87.

[23] G.B. Gelmini and J.W.F. Valle, Fast invisible neutrino decays, Phys. Lett. B 142 (1984) 181.

[24] G. Dvali and L. Funcke, Small neutrino masses from gravitational θ-term, Phys. Rev. D 93 (2016) 113002 [arXiv:1602.03191] [INSPIRE].

[25] L. Funcke, G. Raffelt and E. Vitagliano, Distinguishing Dirac and Majorana neutrinos by their decays via Nambu-Goldstone bosons in the gravitational-anomaly model of neutrino masses, Phys. Rev. D 101 (2020) 015025 [arXiv:1905.01264] [INSPIRE].

[26] J.N. Bahcall, N. Cabibbo and A. Yahil, Are neutrinos stable particles?, Phys. Rev. Lett. 28 (1972) 316 [INSPIRE].

[27] Z.G. Berezhiani, G. Fiorentini, M. Moretti and A. Rossi, Fast neutrino decay and solar neutrino detectors, Z. Phys. C 54 (1992) 581 [INSPIRE].

[28] A. Acker, S. Pakvasa and J.T. Pantaleone, Decaying Dirac neutrinos, Phys. Rev. D 45 (1992) 1 [INSPIRE].

[29] V.D. Barger, J.G. Learned, S. Pakvasa and T.J. Weiler, Neutrino decay as an explanation of atmospheric neutrino observations, Phys. Rev. Lett. 82 (1999) 2640 [astro-ph/9810121] [INSPIRE].

[30] A. Acker and S. Pakvasa, Solar neutrino decay, Phys. Lett. B 320 (1994) 320 [hep-ph/9310207] [INSPIRE].

[31] S. Choubey, S. Goswami and D. Majumdar, Status of the neutrino decay solution to the solar neutrino problem, Phys. Lett. B 484 (2000) 73.

[32] A.S. Joshipura, E. Masso and S. Mohanty, Constraints on decay plus oscillation solutions of the solar neutrino problem, Phys. Rev. D 66 (2002) 113008 [hep-ph/0203181] [INSPIRE].

[33] A.G. Doroshkevich and M.Y. Khlopov, Formation of structure in a universe with unstable neutrinos, Mon. Not. Roy. Astron. Soc. 211 (1984) 277.

[34] A.G. Doroshkevich, A.A. Klyin and M.Y. Khlopov, Cosmological models with unstable neutrinos, Sov. Astron. 32 (1988) 127.

[35] M. Chianese, P. Di Bari, K. Farrag and R. Samanta, Probing relic neutrino radiative decays with 21 cm cosmology, Phys. Lett. B 790 (2019) 64 [arXiv:1805.11717] [INSPIRE].

[36] J.L. Aalberts et al., Precision constraints on radiative neutrino decay with CMB spectral distortion, Phys. Rev. D 98 (2018) 023001 [arXiv:1803.00588] [INSPIRE].

[37] A.G. Beda et al., Gemma experiment: the results of neutrino magnetic moment search, Phys. Part. Nucl. Lett. 10 (2013) 139.

[38] BOREXINO collaboration, Limiting neutrino magnetic moments with Borexino Phase-II solar neutrino data, Phys. Rev. D 96 (2017) 091103 [arXiv:1707.09355] [INSPIRE].
[39] G.G. Raffelt, New bound on neutrino dipole moments from globular cluster stars, Phys. Rev. Lett. 64 (1990) 2856 [INSPIRE].

[40] G.G. Raffelt, Limits on neutrino electromagnetic properties: An update, Phys. Rept. 320 (1999) 319 [INSPIRE].

[41] S. Arceo-Díaz, K.P. Schröder, K. Zuber and D. Jack, Constraint on the magnetic dipole moment of neutrinos by the tip-RGB luminosity in $\omega$-Centauri, Astropart. Phys. 70 (2015) 1 [INSPIRE].

[42] P.J.E. Peebles, The role of neutrinos in the evolution of primeval adiabatic perturbations, Astrophys. J. 180 (1973) 1.

[43] W. Hu and N. Sugiyama, Small scale cosmological perturbations: an analytic approach, Astrophys. J. 471 (1996) 542 [astro-ph/9510117] [INSPIRE].

[44] S. Bashinsky and U. Seljak, Neutrino perturbations in CMB anisotropy and matter clustering, Phys. Rev. D 69 (2004) 083002 [astro-ph/0310198] [INSPIRE].

[45] M. Archidiacono and S. Hannestad, Updated constraints on non-standard neutrino interactions from Planck, JCAP 07 (2014) 046 [arXiv:1311.3873] [INSPIRE].

[46] B. Audren et al., Robustness of cosmic neutrino background detection in the cosmic microwave background, JCAP 03 (2015) 036 [arXiv:1412.5948] [INSPIRE].

[47] B. Follin, L. Knox, M. Millea and Z. Pan, First detection of the acoustic oscillation phase shift expected from the cosmic neutrino background, Phys. Rev. Lett. 115 (2015) 091301 [arXiv:1503.07863] [INSPIRE].

[48] M. Escudero and M. Fairbairn, Cosmological constraints on invisible neutrino decays revisited, Phys. Rev. D 100 (2019) 103531 [arXiv:1907.05425] [INSPIRE].

[49] C.D. Kreisch, F.-Y. Cyr-Racine and O. Doré, The neutrino puzzle: anomalies, interactions and cosmological tensions, arXiv:1902.00534 [INSPIRE].

[50] J.A. Frieman, H.E. Haber and K. Freese, Neutrino mixing, decays and supernova SN1987A, Phys. Lett. B 200 (1988) 115 [INSPIRE].

[51] J.F. Beacom and N.F. Bell, Do solar neutrinos decay?, Phys. Rev. D 65 (2002) 113009 [hep-ph/0204111] [INSPIRE].

[52] A. Bandyopadhyay, S. Choubey and S. Goswami, Neutrino decay confronts the SNO data, Phys. Lett. B 555 (2003) 33 [hep-ph/0204173] [INSPIRE].

[53] M.C. Gonzalez-Garcia and M. Maltoni, Status of oscillation plus decay of atmospheric and long-baseline neutrinos, Phys. Lett. B 663 (2008) 405 [arXiv:0802.3699] [INSPIRE].

[54] R.A. Gomes, A.L.G. Gomes and O.L.G. Peres, Constraints on neutrino decay lifetime using long-baseline charged and neutral current data, Phys. Lett. B 740 (2015) 345 [arXiv:1407.5640] [INSPIRE].

[55] S. Choubey, D. Dutta and D. Pramanik, Invisible neutrino decay in the light of NOvA and T2K data, JHEP 08 (2018) 141 [arXiv:1805.01848] [INSPIRE].

[56] SNO collaboration, Constraints on neutrino lifetime from the sudbury neutrino observatory, Phys. Rev. D 99 (2019) 032013 [arXiv:1812.01088] [INSPIRE].

[57] D. Blas, J. Lesgourgues and T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS) II: approximation schemes, JCAP 07 (2011) 034 [arXiv:1104.2933] [INSPIRE].
[58] PLANCK collaboration, Planck 2015 results. XIII. Cosmological parameters, Astron. Astrophys. 594 (2016) A13 [arXiv:1502.01589] [INSPIRE].

[59] PLANCK collaboration, Planck 2018 results. V. CMB power spectra and likelihoods, arXiv:1907.12875 [INSPIRE].

[60] KATRIN collaboration, KATRIN design report 2004, (2005).

[61] KAMLAND-ZEN collaboration, Search for Majorana neutrinos near the inverted mass hierarchy region with KamLAND-Zen, Phys. Rev. Lett. 117 (2016) 082503 [arXiv:1605.02889] [INSPIRE].

[62] M. Auger et al., The EXO-200 detector, part I: detector design and construction, 2012 JINST 7 P05010 [arXiv:1202.2192] [INSPIRE].

[63] EXO-200 collaboration, Search for Majorana neutrinos with the first two years of EXO-200 data, Nature 510 (2014) 229 [arXiv:1402.6956] [INSPIRE].

[64] T. Brunner and L. Winslow, Searching for $0\nu\beta\beta$ decay in $^{136}$Xe — towards the tonne-scale and beyond, Nucl. Phys. News 27 (2017) 14 [arXiv:1704.01528] [INSPIRE].

[65] B. Audren et al., Strongest model-independent bound on the lifetime of dark matter, JCAP 12 (2014) 028 [arXiv:1407.2418] [INSPIRE].

[66] V. Poulin, P.D. Serpico and J. Lesgourgues, A fresh look at linear cosmological constraints on a decaying dark matter component, JCAP 08 (2016) 036 [arXiv:1606.02073] [INSPIRE].

[67] M. Kawasaki, G. Steigman and H.-S. Kang, Cosmological evolution of an early decaying particle, Nucl. Phys. B 403 (1993) 671 [INSPIRE].

[68] S. Bharadwaj and S.K. Sethi, Decaying neutrinos and large scale structure formation, Astrophys. J. Suppl. 114 (1998) 37 [astro-ph/9707143] [INSPIRE].

[69] M. Kaplinghat, R.E. Lopez, S. Dodelson and R.J. Scherrer, Improved treatment of cosmic microwave background fluctuations induced by a late decaying massive neutrino, Phys. Rev. D 60 (1999) 123508 [astro-ph/9907388] [INSPIRE].

[70] S. Aoyama, K. Ichiki, D. Nitta and N. Sugiyama, Formulation and constraints on decaying dark matter with finite mass daughter particles, JCAP 09 (2011) 025 [arXiv:1106.1984] [INSPIRE].

[71] M.-Y. Wang and A.R. Zentner, Effects of unstable dark matter on large-scale structure and constraints from future surveys, Phys. Rev. D 85 (2012) 043514 [arXiv:1201.2426] [INSPIRE].

[72] S. Aoyama, T. Sekiguchi, K. Ichiki and N. Sugiyama, Evolution of perturbations and cosmological constraints in decaying dark matter models with arbitrary decay mass products, JCAP 07 (2014) 021 [arXiv:1402.2972] [INSPIRE].

[73] C.-P. Ma and E. Bertschinger, Cosmological perturbation theory in the synchronous and conformal Newtonian gauges, Astrophys. J. 455 (1995) 7 [astro-ph/9506072] [INSPIRE].

[74] J. Lesgourgues and T. Tram, Fast and accurate CMB computations in non-flat FLRW universes, JCAP 09 (2014) 032 [arXiv:1312.2697] [INSPIRE].

[75] M. Archidiacono, T. Brinckmann, J. Lesgourgues and V. Poulin, Physical effects involved in the measurements of neutrino masses with future cosmological data, JCAP 02 (2017) 052 [arXiv:1610.09852] [INSPIRE].
[76] F. Bernardeau, *Weak lensing detection in CMB maps*, Astron. Astrophys. **324** (1997) 15 [astro-ph/9611012] [INSPIRE].

[77] D.N. Limber, *The analysis of counts of the extragalactic nebulae in terms of a fluctuating density field*, Astrophys. J. **117** (1953) 134.

[78] Z. Pan, L. Knox and M. White, *Dependence of the cosmic microwave background lensing power spectrum on the matter density*, Mon. Not. Roy. Astron. Soc. **445** (2014) 2941 [arXiv:1406.5459] [INSPIRE].

[79] PLANCK collaboration, *Planck 2015 results. XI. CMB power spectra, likelihoods and robustness of parameters*, Astron. Astrophys. **594** (2016) A11 [arXiv:1507.02704] [INSPIRE].

[80] PLANCK collaboration, *Planck 2015 results. XV. Gravitational lensing*, Astron. Astrophys. **594** (2016) A15 [arXiv:1502.01591] [INSPIRE].

[81] F. Beutler et al., *The 6dF galaxy survey: baryon acoustic oscillations and the local Hubble constant*, Mon. Not. Roy. Astron. Soc. **416** (2011) 3017 [arXiv:1106.3366] [INSPIRE].

[82] A.J. Ross et al., *The clustering of the SDSS DR7 main Galaxy sample — I. A 4 per cent distance measure at z = 0.15*, Mon. Not. Roy. Astron. Soc. **449** (2015) 835 [arXiv:1409.3242] [INSPIRE].

[83] BOSS collaboration, *The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample*, Mon. Not. Roy. Astron. Soc. **470** (2017) 2617 [arXiv:1607.03155] [INSPIRE].

[84] D.M. Scolnic et al., *The complete light-curve sample of spectroscopically confirmed SNe Ia from Pan-STARRS1 and cosmological constraints from the combined pantheon sample*, Astrophys. J. **859** (2018) 101 [arXiv:1710.00845] [INSPIRE].

[85] B.A. Reid et al., *Cosmological constraints from the clustering of the Sloan Digital Sky Survey DR7 luminous red galaxies*, Mon. Not. Roy. Astron. Soc. **404** (2010) 60 [arXiv:0907.1659] [INSPIRE].

[86] F. Köhlinger et al., *KiDS-450: the tomographic weak lensing power spectrum and constraints on cosmological parameters*, Mon. Not. Roy. Astron. Soc. **471** (2017) 4412 [arXiv:1706.02892] [INSPIRE].

[87] B. Audren, J. Lesgourgues, K. Benabed and S. Prunet, *Conservative constraints on early cosmology: an illustration of the Monte Python cosmological parameter inference code*, JCAP **02** (2013) 001 [arXiv:1210.7183] [INSPIRE].

[88] T. Brinckmann and J. Lesgourgues, *MontePython 3: boosted MCMC sampler and other features*, Phys. Dark Univ. **24** (2019) 100260 [arXiv:1804.07261] [INSPIRE].

[89] A. Lewis, *Efficient sampling of fast and slow cosmological parameters*, Phys. Rev. D **87** (2013) 103529 [arXiv:1304.4473] [INSPIRE].

[90] S. Vagnozzi et al., *Unveiling $\nu$ secrets with cosmological data: neutrino masses and mass hierarchy*, Phys. Rev. D **96** (2017) 123503 [arXiv:1701.08172] [INSPIRE].

[91] Z. Chacko, A. Dev, P. Du, V. Poulin and Y. Tsai, *Determining the neutrino lifetime from cosmology*, arXiv:2002.08401 [INSPIRE].

[92] Z. Chacko, L.J. Hall, T. Okui and S.J. Oliver, *CMB signals of neutrino mass generation*, Phys. Rev. D **70** (2004) 085008 [hep-ph/0312267] [INSPIRE].
[93] G.B. Gelmini, S. Nussinov and M. Roncadelli, *Bounds and prospects for the Majoron model of left-handed neutrino masses*, Nucl. Phys. B 209 (1982) 157 [INSPIRE].

[94] E.W. Kolb and M.S. Turner, *Supernova SN 1987a and the secret interactions of neutrinos*, Phys. Rev. D 36 (1987) 2895 [INSPIRE].

[95] K. Choi and A. Santamaria, *Majorons and supernova cooling*, Phys. Rev. D 42 (1990) 293 [INSPIRE].

[96] M. Kachelriess, R. Tomas and J.W.F. Valle, *Supernova bounds on Majoron emitting decays of light neutrinos*, Phys. Rev. D 62 (2000) 023004 [hep-ph/0001039] [INSPIRE].

[97] Y. Farzan, *Bounds on the coupling of the Majoron to light neutrinos from supernova cooling*, Phys. Rev. D 67 (2003) 073015 [hep-ph/0211375] [INSPIRE].

[98] M. Doi, T. Kotani and E. Takasugi, *Double beta decay and Majorana neutrino*, Prog. Theor. Phys. Suppl. 83 (1985) 1 [INSPIRE].

[99] M. Doi, T. Kotani and E. Takasugi, *The neutrinoless double beta decay with Majoron emission*, Phys. Rev. D 37 (1988) 2575 [INSPIRE].

[100] V.D. Barger, W.-Y. Keung and S. Pakvasa, *Majoron emission by neutrinos*, Phys. Rev. D 25 (1982) 907 [INSPIRE].

[101] A.P. Lessa and O.L.G. Peres, *Revising limits on neutrino-Majoron couplings*, Phys. Rev. D 75 (2007) 094001 [hep-ph/0701068] [INSPIRE].

[102] G. Arcadi et al., *Tritium beta decay with additional emission of new light bosons*, JHEP 01 (2019) 206 [arXiv:1811.03530] [INSPIRE].