# Improved cosmological constraints on the neutrino mass and lifetime

Guillermo Franco Abellán, Zackaria Chacko, Abhish Dev, Peizhi Du, Vivian

Poulin, Yuhsin Tsai

► To cite this version:

Guillermo Franco Abellán, Zackaria Chacko, Abhish Dev, Peizhi Du, Vivian Poulin, et al. Improved cosmological constraints on the neutrino mass and lifetime. Journal of High Energy Physics, 2022, 08, pp.076. 10.1007/JHEP08(2022)076. hal-03520954

HAL Id: hal-03520954

https://hal.science/hal-03520954

Submitted on 17 Nov 2023

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
Improved cosmological constraints on the neutrino mass and lifetime

Guillermo Franco Abellán,ª Zackaria Chacko,ᵇ Abhish Dev,ᶜ Peizhi Du,ᵈ Vivian Poulin,ᵃ and Yuhsin Tsaiᵉ

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

ᵇMaryland Center for Fundamental Physics, Department of Physics, University of Maryland, College Park, MD 20742-4111 USA

ᶜTheoretical Physics Department, Fermilab, P.O. Box 500, Batavia, IL 60510, USA

ᵈC.N. Yang Institute for Theoretical Physics, Stony Brook University, Stony Brook, NY, 11794, USA

ᵉDepartment of Physics, University of Notre Dame, IN 46556, USA

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

1 Introduction 1

2 Parameter space of decaying neutrinos 3

3 Boltzmann equations for massive neutrinos decaying into radiation 5

3.1 Derivation of the equations 6

3.2 Signatures of the non-relativistic neutrino decay on the CMB spectra 10

3.3 Consistency of the implementation of Boltzmann equations 12

4 Updated Monte Carlo analysis of the decaying neutrino scenario 15

4.1 Details of the analysis 15

4.2 Main results: updated limit on the neutrino mass and lifetime 16

4.3 Comparison with former results and the impact of Planck 2018 data 17

5 Conclusions 18

A Excluding the relativistic decay regime from the MCMC analysis 20

Bibliography 20

## 1 Introduction

Even though neutrinos were first detected more than six decades ago, they remain among the most mysterious particles in nature, with many of their fundamental properties still to be determined. In particular, although oscillation experiments have provided convincing evidence that neutrinos have non-vanishing masses, these measurements are only sensitive to the mass-squared splittings and consequently the spectrum of neutrino masses remains unknown. The lifetimes of the neutrinos are also poorly constrained, especially in comparison to the other particles in the Standard Model (SM). The determination of the masses and the lifetimes of these mysterious particles remain some of the most important open problems in fundamental physics.

The fact that cosmic neutrinos are among the most abundant particles in the universe, contributing significantly to the total energy density at early times, provides an opportunity to measure their properties. In particular, the evolution of the cosmological density fluctuations depends on $\sum m_\nu$, the sum of neutrino masses. This translates into characteristic effects on the cosmic microwave background (CMB) and large-scale structure (LSS) [1, 2] (for reviews see [3–6]), that are large enough to allow the sum of neutrino masses to be determined in the near future. This determination is based
on the observation that massive neutrinos contribute differently to cosmological observables than either massless neutrinos or cold dark matter (CDM). At early times, while still relativistic, massive neutrinos contribute to the energy density in radiation, just as in the case of massless neutrinos. However, after neutrinos become non-relativistic, their energy density redshifts as matter and therefore contributes more to the expansion rate than massless neutrinos, which would continue to redshift as radiation. As a result, over a given redshift span, the higher expansion rate reduces the time available for the growth of matter density perturbations. However, since massive neutrinos retain pressure until late times, their contribution to the density perturbations on scales below their free streaming lengths is too small to compensate for the shorter structure formation time. Therefore, if neutrinos become non-relativistic after recombination, the net effect of non-vanishing neutrino masses is a suppression of the matter power spectrum and the CMB lensing potential. Based on this, current observations are able to place a bound on the sum of neutrino masses, $\sum m_\nu \lesssim 0.12$ eV [7]. It is important to note that this result assumes that neutrinos are stable on timescales of order the age of the universe. In scenarios in which the neutrinos decay [8, 9], or annihilate away into lighter species [10, 11] on timescales shorter than the age of the universe, this bound is no longer valid and must be reconsidered.

Cosmological observations can also be used to place limits on the neutrino lifetime. In the case of neutrinos that decay to final states containing photons, the bounds on spectral distortions in the cosmic microwave background (CMB) can be translated into limits on the neutrino lifetime, $\tau_\nu \gtrsim 10^{19}$ s for the larger mass splitting and $\tau_\nu \gtrsim 4 \times 10^{21}$ s for the smaller one [12]. In the case of decays to invisible final states, the limits are much weaker. For neutrinos that decay while still relativistic, the decay and inverse decay processes can prevent neutrinos from free streaming. Measurements of the CMB power spectra set a lower bound on the neutrino lifetime, $\tau_\nu \ge 4 \times 10^6$ s ($m_\nu/0.05\text{eV})^5$, in the case of decay into dark radiation [13] (for earlier work see [14–17]). In the case of non-relativistic neutrino decays into dark radiation, the energy density of the decay products redshifts faster than that of stable massive neutrinos. Unstable neutrinos therefore have less of an effect on structure formation than stable neutrinos of the same mass. Consequently, cosmological observables depend both on the masses of the neutrinos and their lifetimes, and heavier values of $\sum m_\nu$ may still be allowed by the data provided the neutrino lifetime is short enough. In Ref. [18], Planck 2015 and LSS data were used to place constraints on the neutrino mass as a function of the lifetime, and found that values of $\sum m_\nu$ as large as 0.9 eV were still allowed by the data. Future LSS measurements at higher redshifts may be able to break the degeneracy between the neutrino mass and lifetime and measure these parameters independently [19]. It is worth noting that there are also bounds on the neutrino lifetime from Supernova 1987A [20], solar neutrinos [21–24], astrophysical neutrinos measured at IceCube [25–30], atmospheric neutrinos and long baseline experiments [31–34]. However, these constraints are in general much weaker than the limits from cosmology.

In this paper we revisit the scenario in which neutrinos decay into dark radiation after becoming non-relativistic and obtain updated limits based on the newer data from Planck 2018. In order to take advantage of the greater precision of the new data, the analysis we perform is also more accurate. We find that, under the assumption that neutrinos decay after becoming non-relativistic, the neutrino mass bound from Planck 2018 data (in combination with BOSS baryon acoustic oscillation (BAO))
data and Pantheon SN1a data) is relaxed to $\sum m_\nu \lesssim 0.42$ eV (95% C.L.).<fnref n="1" /> While this represents a remarkable relaxation of the constraints as compared to the case of stable neutrinos, we note that it is much stronger than the limit derived from Planck 2015 data for the same decaying neutrino scenario, $\sum m_\nu \lesssim 0.9$ eV at (95% C.L.). We show that the improvement of the bound arises primarily from the more precise low-$\ell$ polarization data from Planck 2018, which allows an improved determination of the optical depth to reionization $\tau_{reio}$, thereby breaking the correlation with $\sum m_\nu$ that appears (for relatively high neutrino masses) through the impact of neutrinos on the overall height of the acoustic peaks (i.e. the “early integrated Sachs-Wolfe effect”) [4].

Besides using up-to-date cosmological data, we also improve the analysis from Ref. [18] by incorporating higher order corrections due to neutrino decays into the Boltzmann equations that describe the evolution of Universe’s energy and metric fluctuations. Recently, Ref. [13] provided a complete set of Boltzmann equations for the neutrino decay, but did not conduct Markov Chain Monte Carlo (MCMC) runs necessary to calculate updated neutrino bounds. In this work, we derive Boltzmann equations exactly valid in the absence of ‘inverse-decays’ and quantum statistics. For the numerical implementation, we follow a consistent $T_{dec}/m_\nu$ expansion, where $T_{dec}$ is the temperature at the time of the decay, so that the analysis is under control when neutrinos decay after become non-relativistic.

This paper is organized as follows. In section 2, we present a summary of constraints on the parameter space of decaying neutrinos. In section 3, we derive the set of Boltzmann equations to describe neutrino decay that are valid in the non-relativistic regime and compare our improved analysis to past work. In section 4, we present a MCMC analysis of the decaying neutrino scenario against up-to-date cosmological data. Finally, we conclude in section 5.

## 2 Parameter space of decaying neutrinos

In this section we outline the constraints on the mass and lifetime of neutrinos decaying into dark radiation. As explained in the introduction, current cosmological observables only place limits on a combination of the sum of neutrino masses and their lifetime. Therefore, in this study we will map out the constraints in the two-dimensional parameter space spanned by the sum of neutrino masses ($\sum m_\nu$) and the neutrino decay width ($\Gamma_\nu$), as shown in Fig. 1. In our analysis we assume that all three neutrinos are degenerate in mass. This is a good approximation because the current bounds on $\sum m_\nu$ are larger than the observed mass splittings (see Fig. 1). 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. While this is a simple parameterization of neutrino decays, our bounds can easily be applied to specific models, as done in great details in Ref. [36].

The CMB can be used to constrain the masses and decay widths of neutrinos that decay prior to recombination. When neutrinos decay while still relativistic, decay and inverse decay can prevent neutrinos from free-streaming. If this happens before recombination, it can alter the well-known ‘neutrino drag’ effect that manifests as a phase-shift at high-$\ell$’s in the CMB power spectrum [37–40].

<footnote n="1">¹It is a factor of two weaker than the constraints advocated in Ref. [35], which used a model-independent approach to constrain the neutrino mass as a function of redshift, but neglected the effect of the daughter particles.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The plot shows the current constraints on decaying neutrinos 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 represents $\Gamma_\nu = H(a_{nr})$. Our study focuses on the region below this line, meaning decay happens after neutrinos have become non-relativistic. The grey region shows current constraints on neutrino mass and lifetime coming from the requirement that neutrinos are free streaming close to recombination [13]. The light grey region indicates that this bound may not be applicable when neutrino mass is larger than the temperature of recombination: $m_\nu > T_* \sim 0.2$ eV [13]. Our analysis excludes the red (blue) region labelled “Planck 2015” (“Planck 2018”) based on the data (Planck+BAO+Pantheon). The vertical brown line shows the projected KATRIN sensitivity.</figcaption></figure>

Therefore, CMB data can place a constraint on the decay width of neutrinos. The resulting bound depends on neutrino masses, and was recently updated in Ref. [13], $\tau_\nu \ge 4 \times 10^6 \text{s} (m_\nu/0.05 \text{ eV})^5$. This bound excludes the grey region at the top of Fig. 1.

In addition, based on the analysis in this paper, part of the ‘late-decay’ parameter space can also be excluded based on the gravitational impacts of massive neutrinos on the CMB and LSS. Through the Monte Carlo study presented in section 4, the blue (red) shaded region in Fig. 1 is excluded by the data combination Planck 2018(2015)+BAO+Pantheon.<fnref n="2" /> The orange dashed line in the figure ($\Gamma_\nu = H(a_{nr})$) separates the region where neutrinos decay when non-relativistic from the region where they decay while still relativistic. Here $a_{nr}$ corresponds to the approximate scale factor at

<footnote n="2">²Note that in our analysis we scanned the region between $0 \le \log_{10} \frac{\Gamma_\nu}{\text{km/s/Mpc}} \le 6$. In Fig. 1, we have extrapolated 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$.</footnote>
the time that neutrinos transition to non-relativistic, and is defined as $3T_ν(a_{nr}) = m_ν$. This simple definition is based on the fact that for relativistic neutrinos at temperature $T_ν$, the average energy per neutrino is approximately $3T_ν$. The Hubble scale at $a_{nr}$ is given by,

$$ H(a_{nr}) = H_0 \sqrt{\Omega_m} \left(\frac{\sum m_\nu}{9T_{\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}, $$

where $T_{ν0}$ is the present neutrino temperature. Since our study focuses on the decay of neutrinos after they become non-relativistic, we only present constraints below the orange dashed line. Our analysis shows that $\sum m_ν$ as large as 0.42 eV is still allowed by the data.

Our results have important implications for current and future laboratory experiments designed to detect neutrino masses. Next generation tritium decay experiments such as KATRIN [41] are expected to be sensitive to values of $m_{νe}$ as low as 0.2 eV, corresponding to $\sum m_ν$ of order 0.6 eV. Naively, a signal in these experiments would conflict with the current cosmological bound for stable neutrinos, $\sum m_ν < 0.12$ eV. However, since the unstable neutrino paradigm greatly expands the range of neutrino masses allowed by current cosmological data, it is interesting to explore whether this scenario can accommodate a potential signal at KATRIN. In Fig. 1, we display a brown vertical line $\sum m_ν = 0.6$ eV that corresponds to the expected KATRIN sensitivity. We see that this value of $\sum m_ν$ is too large to be accommodated in the non-relativistic decay regime, where our analysis is valid. However, our result, in combination with those from the ‘relativistic decay’ scenario studied in Ref. [13], leaves open the interesting possibility that neutrinos decaying with a decay width between $\log_{10} \frac{\Gamma_ν}{\text{km/s/Mpc}} \sim 5.5 - 9$ could reconcile cosmological observations with a potential detection at KATRIN, thereby opening a large discovery potential for laboratory experiments. To confirm this conjecture, more work needs to be done to cover the ‘intermediate’ decay regime (i.e. where neutrinos are neither fully relativistic nor fully non-relativistic). We leave this for future work.

In recent years, a number of studies have attempted to constrain the neutrino mass ordering, showing that under the assumption of stable neutrinos, the inverted ordering is now disfavored by constraints from joint analysis of cosmological and oscillation data [42–47] (see also Refs. [48–51] for a different take) as well as from Ly-α observations [52]. However, these arguments are centered on the fact that these analysis lead to a constraint on $\sum m_ν$ at odds with the lower bound on the sum of neutrino masses in the case of inverted ordering, $\sum m_ν \gtrsim 0.1$ eV. Our result suggests that these constraints are strongly dependent on the assumption of neutrino stability over cosmological timescales, and therefore that the inverted ordering is not robustly excluded. It would be very interesting to extend our analysis to the inclusion of Ly-α data to confirm this conclusion.

### 3 Boltzmann equations for massive neutrinos decaying into radiation

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

## 3.1 Derivation of the equations

We denote the phase space distribution of each species as $f(q, \hat{n}, \vec{x}, \tau)$, which is a function of the comoving momentum $q\hat{n}$, coordinates $\vec{x}$ and conformal time $\tau$. The general time evolution of $f$ is controlled by the Boltzmann equations,

$$ \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], \qquad (3.1) $$

where $C[f]$ is the collision term that includes all the processes involving the species.

This phase space distribution has the leading order contribution $\bar{f}(q, \tau)$ that only depends on $q$ and $\tau$, while perturbations are encoded in $\Delta f(q, \hat{n}, \vec{x}, \tau)$,

$$ f(q, \hat{n}, \vec{x}, \tau) \equiv \bar{f}(q, \tau) + \Delta f(q, \hat{n}, \vec{x}, \tau). \qquad (3.2) $$

Treating $\Delta f$ fluctuations about the homogeneous background as higher order perturbations, the zeroth order Boltzmann equations for $\bar{f}$ take the form

$$ \frac{\partial \bar{f}}{\partial \tau} = C[\bar{f}]. \qquad (3.3) $$

In this work, our focus is on the case in which neutrinos decay after turning non-relativistic. In this scenario, we can neglect the effects of inverse decay processes and quantum statistics. The collision term for the neutrino and its daughters are respectively given by [18]

$$ C_{\nu} = - \frac{a^2}{2\epsilon_{\nu}} \int \prod_{i} \frac{d^3 q_i}{2\epsilon_i} |\mathcal{M}|^2 (2\pi)^4 \delta^{(4)}(q - \Sigma_i q_i) f_{\nu}(q), \qquad (3.4) $$

$$ C_{Dj} = + \frac{a^2}{2\epsilon_j} \int \frac{d^3 q}{2\epsilon_\nu} \prod_{i \neq j} \frac{d q_i^3}{2\epsilon_i} |\mathcal{M}|^2 (2\pi)^4 \delta^{(4)}(q - \Sigma_i q_i) f_\nu(q). \qquad (3.5) $$

Here $d^3 q \equiv d^3 q/(2\pi)^3$, $\epsilon \equiv \sqrt{q^2 + a^2 m^2}$ represents the comoving energy and $a$ is the scale factor. The label $i(j)$ denotes the $i$th(jth) daughter. In the case of two body decays to massless daughters, the amplitude squared $|\mathcal{M}|^2$ is simply related to the rest-frame decay width of the neutrino as $|\mathcal{M}|^2 = 16\pi\Gamma_\nu m_\nu$. From the collision terms above, the background evolution for decaying neutrinos is given by

$$ \frac{\partial \bar{f}_{\nu}}{\partial \tau} = -a \frac{\Gamma_{\nu}}{\gamma} \bar{f}_{\nu}, \qquad (3.6) $$

where $\Gamma_\nu$ is the neutrino decay width and $\gamma$ is the Lorentz boost factor,

$$ \gamma = \frac{\sqrt{q^2 + a^2 m_{\nu}^2}}{(a m_{\nu})}. \qquad (3.7) $$
The formal solution to $\bar{f}_{\nu}(q, \tau)$ from the differential equation Eq. (3.6) is

$$ \bar{f}_{\nu}(q, \tau) = \bar{f}_{\text{ini}}(q) e^{-\Gamma_{\nu} \int_{\tau_{\text{ini}}}^{\tau} \frac{a}{\gamma(a)} d\tau'}, \quad (3.8) $$

where $\tau_{\text{ini}}$ denotes the initial conformal time and $\bar{f}_{\text{ini}}(q)$ represents the initial momentum distribution, which we take to be of the Fermi-Dirac form, $\bar{f}_{\text{ini}} = 1/(e^{q/T_{\nu 0}} + 1)$.

The Boltzmann equations for the individual daughter particles do not have a simple form, especially when the daughters consist of more than two species. However, since the daughter particles are taken to be massless in this study the total background density of daughter radiation can be defined as

$$ \bar{\rho}_D \equiv 4\pi a^{-4} \sum_i \int dq \, q^3 \bar{f}_{Di}(q), \quad (3.9) $$

where $\bar{f}_{Di}$ is the background phase space distribution of the $i$th daughter particle. With the definition in Eq. (3.9), regardless of the number of daughter particles and their spins, the Boltzmann equation for the total background daughter density $\bar{\rho}_D$ has the simple form

$$ \frac{\partial \bar{\rho}_D}{\partial \tau} + 4aH\bar{\rho}_D = a\Gamma_{\nu}m_{\nu}\bar{n}_{\nu}, \quad (3.10) $$

where $\bar{n}_{\nu} \equiv 4\pi a^{-3} \int dq \, q^2 \bar{f}_{\nu}(q)$.

We now turn to the Boltzmann equations describing the perturbations of the phase space distribution of decaying neutrinos and their decay products. We work in the synchronous gauge for which the metric perturbations can be parametrized as [53]

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

In Fourier space, $H_{ij}$ is given by

$$ 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.12) $$

where $\vec{k}$ is conjugate to $\vec{x}$ and $h$ and $\eta$ are the two independent scalar metric perturbations. To obtain the Boltzmann hierarchy, we expand the angular dependence of the perturbations as a series in Legendre polynomials,

$$ \Delta f(q, \hat{n}, \vec{k}, \tau) = \sum_{\ell=0}^{\infty} (-i)^{\ell} (2\ell + 1) \Delta f_{\ell}(q, k, \tau) P_{\ell}(\hat{k} \cdot \hat{n}), \quad (3.13) $$

where $P_{\ell}$ represents the $\ell$th Legendre polynomial. The Boltzmann hierarchy for the perturbations of the decaying massive neutrinos $\Delta f_{\nu(\ell)}$ read [18]

$$ \Delta \dot{f}_{\nu(0)} = -\frac{qk}{\varepsilon_{\nu}} \Delta f_{\nu(1)} + q \frac{\partial \bar{f}_{\nu}}{\partial q} \frac{\dot{h}}{6} - \frac{a^2 \Gamma_{\nu} m_{\nu}}{\varepsilon_{\nu}} \Delta f_{\nu(0)}, \quad (3.14) $$

$$ \Delta \dot{f}_{\nu(1)} = \frac{qk}{3\varepsilon_{\nu}} [\Delta f_{\nu(0)} - 2\Delta f_{\nu(2)}] - \frac{a^2 \Gamma_{\nu} m_{\nu}}{\varepsilon_{\nu}} \Delta f_{\nu(1)}, \quad (3.15) $$

$$ \Delta \dot{f}_{\nu(2)} = \frac{qk}{5\varepsilon_{\nu}} [2\Delta f_{\nu(1)} - 3\Delta f_{\nu(3)}] - q \frac{\partial \bar{f}_{\nu}}{\partial q} \frac{(\dot{h} + 6\dot{\eta})}{15} - \frac{a^2 \Gamma_{\nu} m_{\nu}}{\varepsilon_{\nu}} \Delta f_{\nu(2)}, \quad (3.16) $$

$$ \Delta \dot{f}_{\nu(\ell>2)} = \frac{qk}{(2\ell+1)\varepsilon_{\nu}} [\ell \Delta f_{\nu(\ell-1)} - (\ell+1) \Delta f_{\nu(\ell+1)}] - \frac{a^2 \Gamma_{\nu} m_{\nu}}{\varepsilon_{\nu}} \Delta f_{\nu(\ell)}. \quad (3.17) $$
Here $\epsilon_ν = \sqrt{q^2 + a^2 m_ν^2}$ indicates the comoving energy of the neutrinos.

To study the perturbations of the daughter radiation, we focus on the case of two-body decay. In this case, the phase space distributions of the two massless particles are basically identical and they can be considered effectively as one species with a single $f_D$. We can therefore define multipoles $F_{D(ℓ)}$ as in Ref. [54],

$$F_{D(\ell)} \equiv \frac{4\pi}{\rho_c} \int dq q^3 \Delta f_{D(\ell)}, \quad (3.18)$$

where $\rho_c$ is the critical density of Universe today. The Boltzmann hierarchy of the $F_{D(ℓ)}$ can be written as,

$$ \begin{aligned} \dot{F}_{D(0)} &= -k F_{D(1)} - \frac{2}{3} \dot{h} r_D + C_0, \\ \dot{F}_{D(1)} &= \frac{k}{3} F_{D(0)} - \frac{2k}{3} F_{D(2)} + C_1, \\ \dot{F}_{D(2)} &= \frac{2k}{5} F_{D(1)} - \frac{3k}{5} F_{D(3)} + \frac{4(\dot{h} + 6\dot{\eta})}{15} r_D + C_2, \\ \dot{F}_{D(\ell>2)} &= \frac{k}{(2\ell+1)} [\ell F_{D(\ell-1)} - (\ell+1)F_{D(\ell+1)}] + C_\ell, \end{aligned} \quad (3.19) $$

where $r_D \equiv a^4 \bar{\rho}_D / \rho_c$. The terms $C_ℓ$ appearing in Eq. (3.19) arise from the integrated daughter collision term in Eq. (3.5) expanded in terms of Legendre polynomials. The expression for $C_ℓ$ is given by,

$$ \begin{aligned} C_{\ell} &= 2i^{\ell} \int \frac{d\Omega_k}{4\pi} P_{\ell}(\hat{q}_1 \cdot \hat{k}) \left( \frac{4\pi}{\rho_c} \int dq_1 q_1^3 C_{D1}[q_1, \hat{q}_1 \cdot \hat{k}] \right), \\ &= i^{\ell} \left( \frac{32\pi m_\nu \Gamma_\nu a^2}{\rho_c} \right) \int d\Omega_k P_{\ell}(\hat{q}_1 \cdot \hat{k}) \int \frac{dq_1}{2\epsilon_1} q_1^3 \int \frac{d^3 q_2}{2\epsilon_2} \frac{d^3 q}{2\epsilon_\nu} \Delta f_\nu(q, \hat{q} \cdot \hat{k}) (2\pi)^4 \delta^{(4)}(q - q_1 - q_2). \end{aligned} \quad (3.20) $$

The overall factor of two in the equation above arises because we are adding the collision integrals of the two massless daughters, which are of the same form. In this expression $d\Omega_k$ represents the differential solid angle along the direction $\hat{k}$, while $q_{1,2}$ are the momenta of daughter particles. The $d^3q_2$ integral can be easily evaluated using the delta function corresponding to momentum conservation. In order to perform the integral over $d\Omega_k$, we notice that the direction of $\hat{k}$ enters only via $P_\ell(\hat{q}_1 \cdot \hat{k})$ and $\Delta f_\nu(q, \hat{q} \cdot \hat{k})$. Now, using the Legendre expansion of $\Delta f_\nu(q, \hat{q} \cdot \hat{k})$ in Eq. (3.13) and employing the identity

$$ \int d\Omega_k P_\ell(\hat{k} \cdot \hat{q}) P_{\ell'}(\hat{k} \cdot \hat{q}_1) = \left( \frac{4\pi}{2\ell+1} \right) P_\ell(\hat{q} \cdot \hat{q}_1) \delta_{\ell\ell'}, \quad (3.21) $$

we can evaluate the $d\Omega_k$ integral to obtain

$$ C_{\ell} = \left( \frac{128\pi^2 m_{\nu} \Gamma_{\nu} a^2}{\rho_c} \right) \int \frac{d^3 q dq_1}{8\epsilon_{\nu} \epsilon_1 \epsilon_2} q_1^3 P_{\ell}(\hat{q}_1 \cdot \hat{q}) \Delta f_{\nu\ell}(q) (2\pi)^4 \delta(\epsilon_{\nu} - \epsilon_1 - \epsilon_2). \quad (3.22) $$
Now, notice that the direction of the neutrino momentum only enters the integrand via the angle between the neutrino momentum $q$ and the daughter momentum $q_1$, defined as $\cos \theta_1 \equiv \hat{q} \cdot \hat{q}_1$. The energy conserving delta function can be expressed in terms of this angle as

$$ \delta(\epsilon_{\nu} - \epsilon_{1} - \epsilon_{2}) = \frac{\epsilon_{2}}{qq_{1}}\delta(\cos \theta_{1} - \cos \theta_{1}^{*}), \qquad (3.23) $$

where

$$ \cos \theta_{1}^{*} = \frac{2\epsilon_{\nu}q_{1} - a^{2}m_{\nu}^{2}}{2qq_{1}}. \qquad (3.24) $$

The energy conservation restricts the daughter momentum to a range of values ($q_1^+, q_1^-$). The edges of this range occur when the extreme values, $\cos \theta_1^* = \pm 1$, are reached. For these values,

$$ q_{1}^{\pm} = \frac{m_{\nu}^{2}a^{2}}{2(\epsilon_{\nu} \pm q)}. \qquad (3.25) $$

After integrating over the delta function corresponding to energy conservation, this reduces to the simpler form,

$$ C_{\ell} = \left( \frac{8\pi m_{\nu} \Gamma_{\nu} a^2}{\rho_c} \right) \int \frac{dq}{\epsilon_{\nu}} q \Delta f_{\nu(\ell)} \int_{q_1^+}^{q_1^-} dq_1 q_1 P_{\ell} \left( \frac{2\epsilon_{\nu}q_1 - a^2m_{\nu}^2}{2qq_1} \right). \qquad (3.26) $$

Eq. (3.26) may also be obtained by taking the appropriate limit of the more general expression in Ref. [13]. The same Boltzmann hierarchy has been derived in the context of warm matter decaying into dark radiation [55].

Performing the integral over $q_1$, we can obtain the following expressions for the first few $C_\ell$'s,

$$ \begin{aligned} C_0 &= \frac{4\pi a^2 \Gamma_\nu m_\nu}{\rho_c} \int dq q^2 \Delta f_{\nu(0)}, \\ C_1 &= \frac{4\pi a^2 \Gamma_\nu m_\nu}{\rho_c} \int dq \frac{q^3}{\epsilon_\nu} \Delta f_{\nu(1)}, \\ C_2 &= \frac{4\pi a^2 \Gamma_\nu m_\nu}{\rho_c} \int dq q^2 g_2(q, \epsilon_\nu) \Delta f_{\nu(2)}, \\ C_3 &= \frac{4\pi a^2 \Gamma_\nu m_\nu}{\rho_c} \int dq q^2 g_3(q, \epsilon_\nu) \Delta f_{\nu(3)}. \end{aligned} \qquad (3.27) $$

Here the functions $g_2(q, \epsilon_\nu)$ and $g_3(q, \epsilon_\nu)$ are given by,

$$ \begin{aligned} g_2(q, \epsilon_\nu) &\equiv \frac{5}{2} - \frac{3}{2} \frac{\epsilon_\nu^2}{q^2} + \frac{3}{4} \frac{(\epsilon_\nu^2 - q^2)^2}{\epsilon_\nu q^3} \ln \left( \frac{\epsilon_\nu + q}{\epsilon_\nu - q} \right), \\ g_3(q, \epsilon_\nu) &\equiv \frac{25}{2} \frac{\epsilon_\nu}{q} - \frac{4q}{\epsilon_\nu} - \frac{15}{2} \frac{\epsilon_\nu^3}{q^3} + \frac{15}{4} \frac{(\epsilon_\nu^2 - q^2)^2}{q^4} \ln \left( \frac{\epsilon_\nu + q}{\epsilon_\nu - q} \right). \end{aligned} \qquad (3.28) $$

Given the complicated integrals in Eq. (3.26), it is technically challenging to keep track of all the collision terms in the Boltzmann hierarchy. Instead, we choose to keep just the first few $C_\ell$'s for $\ell \le \ell_{max}$. The idea behind this approach is that $C_\ell$ is of $O((T_{dec}/m_\nu)^\ell)$ around the time of decay.
Therefore, for non-relativistic decay ($T_{dec}/m_ν \ll 1$), it is self-consistent to set $C_{ℓ>ℓ_{max}} = 0$ because those terms only have negligible effect on physical observables. To understand the scaling of $C_ℓ$, we first note that the integral over $q$ in Eq. (3.27) receives most of its support from the region around $q \sim T_{ν0}$ because $\Delta f_{ν(ℓ)}$ inherits features of the Fermi-Dirac distribution from $\bar{f}_{ini} = 1/(e^{q/T_{ν0}} + 1)$. Deep in the non-relativistic region, $q \ll \epsilon_ν$ and $T_{ν0} \ll a m_ν$. In this regime, we can employ a Taylor expansion for the functions $g_2$ and $g_3$ in powers of $q/\epsilon_ν$ to obtain,

$$g_2(q, \epsilon_ν) \approx \frac{4 q^2}{5 \epsilon_ν^2}, \quad g_3(q, \epsilon_ν) \approx \frac{4 q^3}{7 \epsilon_ν^3} \quad \text{for } (q \ll \epsilon_ν). \qquad (3.29)$$

Inserting Eq. (3.29) above into Eq. (3.27), it is straightforward to see that $C_ℓ \propto (T_{ν0}/a m_ν)^ℓ$. Moreover, if we assume decay happens deep in the non-relativistic region, we will get $C_ℓ \propto (T_{dec}/m_ν)^ℓ$ when decay happens, where $T_{dec} = T_{ν0}/a_{dec}$. Therefore, $C_ℓ$ is suppressed by powers of $T_{dec}/m_ν \ll 1$ for higher $ℓ$. To further justify this argument, we show in section 3.3 that setting $ℓ_{max} = 2$ or $ℓ_{max} = 3$ makes negligible difference to cosmological observables (see Fig. 4). Therefore, we only keep $C_{ℓ≤3}$ and set $C_{ℓ>3} = 0$ in our numerical study for simplicity.

Physically, the expansion in the small parameter $T_{dec}/m_ν$ corresponds to perturbing about the ultra-nonrelativistic limit in which the momentum of the mother particle has completely redshifted away, so that it has come to rest in the cosmic frame. Energy and momentum conservation is respected order by order in this expansion. The earlier work [18] approximated the Boltzmann hierarchy for daughter radiation (Eq. 3.19) by just keeping $C_0$ and setting all the $C_{ℓ≥1} = 0$. It is clear from the above discussion that this is a consistent approximation to zeroth order in an expansion in the small parameter $T_{dec}/m_ν$. The authors in Ref. [13] argued that the Boltzmann hierarchy for daughter radiation in Ref. [18] does not reproduce the standard decaying CDM scenario and does not respect momentum conservation. Both criticisms can be addressed by considering the term $C_1$. Since $C_1$ begins at $O(T_{dec}/m_ν)$, we see that the Boltzmann hierarchy in Ref. [18] does in fact reproduce the decaying CDM scenario and respects momentum conservation up to $O(T_{dec}/m_ν)$ corrections, consistent with the approximation. In this limit, the momenta of the daughter particles arise entirely from the rest mass of the mother. In practice, since the contributions of neutrinos to the density perturbations are small, we will see that the higher order terms do not significantly affect the constraints derived in Ref. [18] with Planck 2015 data.

## 3.2 Signatures of the non-relativistic neutrino decay on the CMB spectra

To make this work fully self-contained, we briefly summarize the impact of the non-relativistic invisible neutrino decays on the CMB spectra, following the discussion in Ref. [18]. In Fig. 2, we display the residuals in the CMB (lensed) TT, EE and lensing power spectra, for the sum of neutrino masses $\sum m_ν = 0.6$ eV and several decay widths $\text{Log}_{10}(\Gamma_ν/\text{km/s/Mpc}) = 0, 2, 4, 6$. In all cases, the $\Lambda$CDM parameters are set to their best-fit values from Planck 2018, that is, $\{100\theta_s = 1.04089, \omega_{cdm} = 0.1198, \omega_b = 0.02233, n_s = 0.9652, \ln(10^{10}A_s) = 3.043, \tau_{reio} = 0.0540\}$. Our reference $\Lambda$CDM model makes use of the same parameters and assumes standard massless neutrinos.

For the value of the mass considered ($\sum m_ν = 0.6$ eV) and at fixed angular size of the sound horizon $\theta_s$, neutrino masses primarily impact the lensing spectrum. Indeed, as they reduce power
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Residuals in the CMB lensed TT (upper), EE (middle) and lensing (lower) spectrum for a fixed value of the neutrino mass and several decay widths. The residuals are taken with respect to the ΛCDM best-fit parameters from Planck 2018. The ΛCDM parameters are kept fixed in all cases.</figcaption></figure>

below the free-streaming scale, they produce a significant matter power suppression at small scales, which leads to a ~20% reduction in the $C_\ell^{\phi\phi}$ at large $\ell$ (blue curve in Fig. 2). Consequently, this power suppression decreases the smoothing in the high-$\ell$ part of the TT and EE spectra, which can be seen as ‘wiggles’ in the corresponding plots.

In addition, stable neutrinos dilute like non-relativistic matter at late times ($\bar{\rho}_\nu \sim a^{-3}$), which increases the value of $\Omega_m$. As we impose the closure relation $\Omega_m + \Omega_\Lambda = 1$ at late-times, this is compensated for by a decrease in $\Omega_\Lambda$ (later beginning of Λ-domination), and thus a reduction in the Late Integrated Sachs-Wolfe effect (LISW), leaving a signature in the low-$\ell$ TT spectrum. Furthermore, the modified expansion history $H(z)$ changes quantities integrated along $z$, such as $\tau_{reio}$, which affects the multipoles at $\ell \sim 10$ in the EE spectrum.

When a non-negligible $\Gamma_\nu$ is considered (orange, green and red curves in Fig. 2), one can see that the aforementioned effects typically become less prominent for earlier decays. This is particularly true for the high-$\ell$ part of the lensing spectrum (and consequently the smoothing at high-$\ell$ in TT and EE spectra).
EE) since decay of neutrinos reduce their impact on structure formation. The reduction of the effect in the low-$\ell$ part of the TT and EE spectra is not entirely monotonic, as intermediate values of $\Gamma_\nu$ can induce additional time variation in the gravitational potentials (thereby affecting the LISW effect), as well as time variations in $H(z)$ (thereby affecting $\tau_{reio}$). As a result, the $\Lambda$CDM limit is reached not only for small values of $\sum m_\nu$, but also for high values of $\Gamma_\nu$. This will be reflected in the MCMC analysis in section 4, which shows a large positive correlation between both parameters. It is precisely this degeneracy which relaxes the neutrino mass bounds.

### 3.3 Consistency of the implementation of Boltzmann equations

We begin by comparing the approximation used in Ref. [18] for the background energy density of decaying massive neutrinos to the more accurate results obtained by evaluating the integral in Eq. (3.8) numerically. In Ref. [18], the phase space distribution of neutrinos in Eq. (3.8) is approximated through the following analytic formula,

$$ \bar{f}_{\nu}(q, \tau) = \bar{f}_{\text{ini}}(q)e^{-\Gamma_{\nu}t/\gamma}. \qquad (3.30) $$

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Redshift evolution of the quantity $(\bar{\rho}_{\nu} + \bar{\rho}_{D})/\bar{\rho}_{ur}$ (where $\bar{\rho}_{ur}$ denotes the energy density of stable massless neutrinos), which should be equal to 1 in the limit of relativistic decays. We consider a very small value of the neutrino mass sum, $\sum m_\nu = 0.06$ eV, and several values for the decay width, $\text{Log}_{10}(\Gamma_\nu/\text{km/s/Mpc})$. “approx. PSD” refers to the approximated phase space distribution in Eq. (3.30) while “Full PSD” refers to the exact solutions of Eq. (3.8).</figcaption></figure>

As argued in Ref. [18], this approximation is valid under the assumption that the decay happens deep in the non-relativistic regime. To see the difference between the approximation and the full
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Fractional change in the CMB TT (upper), EE (middle) and lensing (lower) spectrum, when imposing different prescriptions for the background energy density distribution and Boltzmann hierarchies. “approx. PSD” refers to the approximate phase space distribution in Eq. (3.30) while “Full PSD” refers to the exact solution of Eq. (3.8). “$C_\ell$” in the plot means we only keep those collision terms in Eq. (3.19). The chosen values of the neutrino mass ($\sum m_\nu = 0.6$ eV) and decay width ($\text{Log}_{10}(\Gamma_\nu/[\text{km/s/Mpc}]) = 5.5$) correspond to the case when neutrinos decay close to non-relativistic transition ($T_{\text{dec}}/m_\nu \approx 0.3$). The gray shaded region indicates Planck 2018 1-$\sigma$ uncertainties, while the pink boxes indicate the (binned) cosmic variance.</figcaption></figure>

result, we plot the ratio $r = (\bar{\rho}_\nu + \bar{\rho}_\text{D})/\bar{\rho}_\text{ur}$ in Fig. 3, for several values of the decay width $\Gamma_\nu$ and a fixed value of the total neutrino mass $\sum m_\nu = 0.06$ eV. Here $\bar{\rho}_\text{ur}$ denotes the energy density of stable massless neutrinos. If neutrinos decay while relativistic, this ratio always gives $r \simeq 1$. However, if the decay happens when the neutrinos are already non-relativistic ($\bar{\rho}_\nu \sim a^{-3}$), then the ratio evolves from $r \simeq 1$ to $r \sim a$, and will eventually reach a plateau once all the neutrinos have decayed. From Fig. 3, we can see that the approximate formula in Eq. (3.30) gradually improves as we go to smaller decay widths (that is, going deeper into the regime of non-relativistic decays), as expected. The error in the case of neutrinos decaying right around the time of the non-relativistic transition ($\text{Log}_{10}(\Gamma_\nu/[\text{km/s/Mpc}]) \simeq 4$ for $\sum 0.06$ eV) is around 25%. Nevertheless, as we argue below, the
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Same as in Fig. 4, but with a smaller decay width ($\text{Log}_{10}(\Gamma_{\nu}/[\text{km/s/Mpc}]) = 4$), corresponding to a neutrino decay happening deep in the non-relativistic limit ($T_{\text{dec}}/m_{\nu} \simeq 0.03$).</figcaption></figure>

impact on observables is much smaller given that neutrinos only contribute a small fraction of the total energy density for masses considered in this work. Not surprisingly, the approximate formula fails in the relativistic regime, leading to $r < 1$ at late-times. Therefore future work focusing on this regime should make use of the exact formula.

In Figs. 4 and 5, we show the effects of various approximations in dealing with decaying neutrinos (at the background and perturbation level) on the CMB TT, EE and lensing spectra. We compare the impact of using either the approximated or the exact PSD of neutrinos discussed above, as well as the impact of only keeping $C_{\ell \le \ell_{\text{max}}}$ in the Boltzmann hierarchy of daughter particles in Eq. (3.19), where we vary $\ell_{\text{max}}$ from zero to three. We show the residuals of these approximations with respect to the ‘optimal’ case (i.e. including all terms up to $\ell_{\text{max}} = 3$ and the exact background PSD) for a fixed value of the neutrino mass ($\sum m_{\nu} = 0.6$ eV) and two different decay widths ($\text{Log}_{10}(\Gamma_{\nu}/[\text{km/s/Mpc}]) = 5.5$ in Fig. 4 and $\text{Log}_{10}(\Gamma_{\nu}/[\text{km/s/Mpc}]) = 4$ in Fig. 5). Fig. 4 corresponds to decays happening around the time of the non-relativistic transition, $T_{\text{dec}}/m_{\nu} \simeq 0.3$, where
the effects of the approximations are expected to be largest. Fig. 5 on the other hand refers to decays happening deep in the non-relativistic regime, $T_{\text{dec}}/m_{\nu} \simeq 0.03$. We also show the Planck 2018 1-$\sigma$ error bars, as well as the (binned) cosmic variance.

For decays close to the non-relativistic transition $T_{\text{dec}}/m_{\nu} \simeq 0.3$ shown in Fig. 4, we find that the biggest improvement in the CMB TT spectrum occurs when including $C_{\ell \le 1}$ (i.e., the contribution from the decaying neutrino bulk velocity) in the Boltzmann hierarchy of daughter radiation, which impacts the integrated Sachs-Wolfe (ISW) effect at multipoles $\ell \lesssim 100$. On the other hand, the approximate background distribution of neutrinos does not have a significant effect. For the CMB EE spectrum shown in the same figure, which is not sourced by the ISW effect, the impact of the approximate background distribution of neutrinos is comparable to the effect of the approximate perturbed hierarchy. Nevertheless, one can see that for $\ell_{\text{max}} \ge 2$, additional contributions to the daughter hierarchy have negligible impacts, which justifies our choice of cutting the collision term $C_{\ell}$ contribution at $\ell_{\text{max}} = 3$. Finally for the CMB lensing spectrum, the effects due to the approximate treatment of the background PSD dominate over the ones due to including higher order terms in the Boltzmann hierarchy of the dark radiation. This is expected given that the matter power spectrum suppression scales approximately with $\bar{\rho}_{\nu}/\bar{\rho}_{m}$ [2, 4] where $\bar{\rho}_{m}$ is the total matter density, while neutrino perturbations are very small well below the free-streaming scale, so that their detailed dynamics is not as important as on larger scales.

The impact of the various approximations in the case of decays deep in the non-relativistic regime $T_{\text{dec}}/m_{\nu} \simeq 0.03$, displayed in Fig. 5, is much less visible. In that case, one can therefore safely neglect $C_{\ell>0}$ and consider the approximate PSD, as done in Ref. [18].

# 4 Updated Monte Carlo analysis of the decaying neutrino scenario

## 4.1 Details of the analysis

In this section we perform a numerical scan over the parameter space to obtain updated limits on the neutrino mass and lifetime. We perform comprehensive MCMC analyses with the `MontePython-v3`[^3] [56, 57] code interfaced with our modified version of `CLASS`. We fit the decaying neutrino model to a combination of the following data-sets:

* The Planck 2018 high-$\ell$ TT, TE, EE + low-$\ell$ data TT, EE + lensing data [7]. We will also compare these results with the use of Planck 2015 data to disentangle the effects of our improved formalism and that of the new data.

* The BAO measurements from 6dFGS at $z = 0.106$ [58], SDSS DR7 at $z = 0.15$ [59], BOSS DR12 at $z = 0.38, 0.51$ and $0.61$ [60], and the joint constraints from eBOSS DR14 Ly-$\alpha$ autocorrelation at $z = 2.34$ [61] and cross-correlation at $z = 2.35$ [62].

* The measurements of the growth function $f\sigma_8(z)$ (FS) from the CMASS and LOWZ galaxy samples of BOSS DR12 at $z = 0.38, 0.51$, and $0.61$ [60].

<footnote>[^3]: https://github.com/brinckmann/montepython_public</footnote>
• The Pantheon SNIa catalogue, spanning redshifts $0.01 < z < 2.3$ [63].

We adopt wide flat priors on the following six ΛCDM parameters: {$\omega_b, \omega_{cdm}, H_0, n_s, A_s, \tau_{reio}$}. We assume three degenerate neutrinos decaying into massless radiation and consider flat priors on $\sum m_\nu$/eV and Log$_{10}$(Γ$_\nu$/[km/s/Mpc]). To accelerate convergence, we split the parameter space between Log$_{10}$(Γ$_\nu$/[km/s/Mpc]) ∈ [0.1, 2.5] and Log$_{10}$(Γ$_\nu$/[km/s/Mpc]) ∈ [2.5, 6.5]. In both cases we take wide priors on $\sum m_\nu$ ∈ [0.06, 1.5] eV. We assume our MCMC chains to be converged when the Gelman-Rubin criterion $R - 1 < 0.05$ [64]. In our baseline analysis, we do not apply any specific cut to the parameter space, even if neutrinos decay in the relativistic regime (this occurs for low $\sum m_\nu$ and high Γ$_\nu$). In appendix A, we investigate the impact of imposing a prior that excludes the parameter space corresponding to relativistic decay from our analysis and show that the limit at 95% on $\sum m_\nu$ agrees within a few percent.

## 4.2 Main results: updated limit on the neutrino mass and lifetime

The results of our analyses are presented in Figs. 6. For very late decays, Log$_{10}$(Γ$_\nu$/[km/s/Mpc]) ≲ 2.5, no relaxation of the constraints on $\sum m_\nu$/eV is visible, in agreement with what was found in Ref. [18]. The impact of the new Planck data is visible as a significantly improved bound on the sum of neutrino mass, namely we find $\sum m_\nu < 0.127$ eV (95%C.L.), an improvement of about ~35% over 2015 data, in good agreement with Ref. [7]. For Log$_{10}$(Γ$_\nu$/[km/s/Mpc]) ≳ 2.5, one can see that the bound relaxes as expected, although not as much with Planck 2018 data as for Planck 2015 data.

Taking the intersect of the non-relativistic decay line as our $2\sigma$ limit, we find that Planck 2018 allows neutrinos with masses up to $\sum m_\nu = 0.42$ eV. In appendix A, we present an alternative analysis that directly imposes the non-relativistic decay criterion as a prior while performing the scan. Marginalizing over all parameters we find the same result, $\sum m_\nu < 0.42$ eV (95% C.L.). The excellent agreement between these different analyses leads us to conclude with confidence that, within the regime of non-relativistic decay, values of $\sum m_\nu$ as large as 0.4 eV are still allowed by the data. This bound is significantly stronger than the limit from Planck 2015 data, for which $\sum m_\nu \sim 0.9$ eV was still allowed in the non-relativistic decay scenario.

Our result also has implications for laboratory searches. For $\sum m_\nu = 0.6$ eV, the smallest mass scale that the KATRIN experiment is designed to probe, Planck 2018 data requires decay rate Γ$_\nu$ ≳ $10^{5.5}$ km/s/Mpc, a constraint roughly one order of magnitude stronger than from Planck 2015 data. However, this value of the decay rate is now slightly beyond the regime of validity of our work$^4$, indicating that, in the event of a neutrino mass discovery at KATRIN, a more involved analysis including inverse-decays would be necessary to confirm that the decay scenario can reconcile laboratory and cosmological measurements.

<footnote>$^4$For $\sum m_\nu = 0.6$ eV and assuming degenerate neutrino masses, the non-relativistic condition requires $\Gamma_\nu < 10^{5.3}$ km/s/Mpc.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** 2D posterior distribution of the decaying neutrino model reconstructed from the analysis of BAO + FS + Pantheon together with either *Planck* 2015 or *Planck* 2018 data. In the top panel, we show the correlation with other cosmological parameters.</figcaption></figure>

## 4.3 Comparison with former results and the impact of *Planck* 2018 data

Comparing with the constraints presented in Ref. [18] for *Planck* 2015, we find that, while the impact of our improved treatment is clearly visible in the CMB power spectra (and will be relevant for future experiments), it has only a marginal impact on the constraints, and our bounds are in very good agreement with those derived in Ref. [18], which only included the leading order term in the daughter radiation hierarchy<fnref n="5" />. The bulk of the improvement is due to the newest *Planck* 2018 data and can be understood as follows. As shown in Fig. 2, for the masses we consider, the main effect is an almost scale independent suppression of CMB lensing spectrum. This suppression can be compensated for by increasing the primordial amplitude $A_s$ or by adjusting the matter density $\omega_{cdm}$ (see Ref. [65] for

<footnote n="5">⁵Let us note that the implementation of the BAO/$f\sigma_8$ DR12 likelihood used in Ref. [18] within the MontePython code had an issue that led to constraints on $\sum m_\nu$ that were somewhat milder than the true bounds. MontePython has since then been corrected, leading to an improvement on the constraints on the stable/long-lived ($\Gamma_\nu < 10^3$ km/s/Mpc) case by about 20%. However, we have verified that this bug had no impact in the short-lived case ($\Gamma_\nu > 10^{2.5}$ km/s/Mpc).</footnote>
a discussion of the correlation between {$\sum m_\nu, A_s, \tau_{reio}, \omega_{cdm}$}). Due to the well-known degeneracy between $A_s$ and $e^{-2\tau_{reio}}$, Planck 2015 data, which was limited in polarization, were unable to place a tight constraint on $A_s$, and thus the constraining power on the sum of neutrino mass and lifetime was limited. The precise measurements of low-$\ell$ polarization from Planck 2018 leads to constraints on $\tau_{reio}$ that are tighter by a factor of two than those from Planck 2015. As a result, parameters degenerate with $\tau_{reio}$ such as $A_s$ are now much better constrained. Consequently, the constraints on the sum of neutrino mass and lifetime have significantly improved with Planck 2018 data. To confirm this simple argument, we perform another MCMC run with Planck 2015 data and a tight gaussian prior on $\tau_{reio} = 0.0540 \pm 0.0074$, chosen to match the optical depth to reionization reconstructed from Planck 2018. Given that the constraints on $\sum m_\nu$ are independent of $\Gamma_\nu$ below $\Gamma_\nu \lesssim 10^3$, and the scaling above $\Gamma_\nu \lesssim 10^{5.5}$ is monotonic, we focus on the parameter space $\text{Log}_{10}(\Gamma_\nu/[\text{km/s/Mpc}]) \in [3, 5.5]$ to accelerate convergence. Our results are presented in Fig. 7, where one can see that this simple prescription leads to constraints that are very similar to those from the full Planck 2018 data. We attribute the remaining differences to the additional constraining power of Planck 2018 data on the parameters $\omega_{cdm}$ and $\omega_b$, which are mildly correlated with $\sum m_\nu$ (see Fig. 6, top panel). Note that our constraints are a factor of two weaker than those advocated in Ref. [35], which performed a ‘model-independent’ reconstruction of the neutrino mass as a function of redshift, but neglects the decay products. As we show here, including details about the daughter radiation is necessary to accurately compute the effect of neutrino decays even in the non-relativistic regime. Finally, as discussed in Refs. [19, 65], a combination of CMB data with future tomographic measurements of the power spectrum by DESI [66] or Euclid [67], and an improved determination of the optical depth to reionization by 21-cm observations with SKA [68, 69], could greatly increase the sensitivity of cosmological probes to neutrino masses and lifetimes.

# 5 Conclusions

Cosmological observations are known to set the strongest constraints on the sum of neutrino masses. Yet, the existing mass bound from CMB and LSS measurements, which assumes that neutrinos are stable, is significantly weakened if neutrinos decay. In this work, we provide up-to-date limits on the lifetime of massive neutrinos that decay into dark radiation after becoming non-relativistic, from a combination of CMB, BAO, growth factor measurements, and Pantheon SN1a data.

Compared to the earlier analysis [18], we have incorporated higher-order corrections up to $O((T_{\text{dec}}/m_\nu)^3)$ when solving the dark radiation perturbations, and also performed the full calculation of the background energy density of the decaying neutrino using Eq. (3.8). The more precise treatment of the Boltzmann equations and the background energy evolution in our MCMC study improves the coverage of the case when the neutrinos decay early so that their average momenta are close to their masses. As shown in Fig. 5, if neutrinos decay when having $T_\nu \ll m_\nu/3$, the inclusion of higher moment perturbations $C_{\ell \ge 2}$ gives a negligible change to the power spectra as compared to the experimental uncertainties. However, the complete calculation of the neutrino energy does improve the prediction for the power spectrum significantly from the approximate result using Eq. (3.30) when the decays happen semi-relativistically. Nevertheless, we have found that constraints from Planck 2015,
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Posterior distribution of $\sum m_\nu$ and $\text{Log}_{10}(\Gamma_\nu/[\text{km/s/Mpc}])$ with Planck 2018 and Planck 2015 + a primer on $\tau_{\text{reio}}$ from Planck 2018. The agreement of the posteriors shows the dominant constraining power on $\sum m_\nu$ and $\text{Log}_{10}(\Gamma_\nu/[\text{km/s/Mpc}])$ comes from a precise measurement of $\tau_{\text{reio}}$ from Planck 2018.</figcaption></figure>

given their limited precision, are unaffected by these considerations. However, we anticipate that these effects will be relevant for future experiments (as well as an essential contribution in the relativistic case, to be considered in the future).

In fact, we have shown that the bulk of the improvement in the constraining power compared to Ref. [18] comes from the use of Planck 2018 data. Indeed, we have demonstrated that the improved $\tau_{\text{reio}}$ measurement from the low-$\ell$ polarization data helps breaking the degeneracy in the CMB power spectrum amplitude and strengthens the bound on the neutrino mass and lifetime. As a result, we have found that neutrinos with $\sum m_\nu > 0.42 \text{ eV} (2\sigma)$ cannot be made consistent with cosmological data if they decay while non-relativistic, a significant improvement from Planck 2015 data for which masses as high as $\sum m_\nu \sim 0.9 \text{ eV}$ were consistent with the non-relativistic decay scenario [18].

We have argued that one notable application of this result is that, if the KATRIN experiment sees an electron neutrino with $m_\nu \approx 0.2 \text{ eV}$ (the advocated sensitivity), our result would constrain $\Gamma_\nu \gtrsim 10^{5.5} \text{ km/s/Mpc}$, i.e. the neutrinos would need to decay between $z \approx 2 \times 10^2 - 4 \times 10^3$, while they are still relativistic, so that our bounds and the bounds studied in Ref. [13] would not apply. In case of a neutrino mass discovery at KATRIN, a more involved analysis including inverse-decays would be necessary to firmly confirm that the decay scenario can reconcile laboratory and cosmological measurements. Additionally, our results show that the tentative exclusion of the inverted mass ordering [44–46, 52], based solely on the fact that the inverted ordering predicts $\sum m_\nu > 0.1 \text{ eV}$, is highly dependent on the hypothesis that neutrinos are stable on cosmological time-scales. Non-relativistic decays can still easily reconcile the inverted ordering with cosmological data.

Finally, let us mention that even though current exclusion bounds in Fig. 6 do not set independent
constraints on the neutrino mass and lifetime, next generation measurements of the matter power
spectrum at different redshifts can help break that degeneracy [19]. It will be interesting to revisit the
forecast on the sensitivity of future cosmological data to the sum of neutrino masses and their lifetime
in light of our improved formalism.

Acknowledgements

PD is supported in part by NSF grant PHY-1915093. ZC is supported in part by the National Science Foundation under Grant Number PHY-1914731. ZC is also supported in part by the US-Israeli BSF grant 2018236. YT is supported by the NSF grant PHY-2014165 and PHY-2112540. Fermilab is operated by Fermi Research Alliance, LLC under contract number DE-AC02-07CH11359 with the United States Department of Energy. This work has been partly supported by the CNRS-IN2P3 grant Dark21. The authors acknowledge the use of computational resources from the Dark Energy computing Center funded by the Excellence Initiative of Aix-Marseille University - A*MIDEX, a French "Investissements d'Avenir" programme (AMX-19-IET-008 - IPhU). PD and YT thank the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611, where part of this work was performed. PD was also partially supported by a grant from the Simons Foundation during the stay at the Aspen Center for Physics.

A Excluding the relativistic decay regime from the MCMC analysis

In our baseline analysis, we have extrapolated our scans to the (mildly-)relativistic decay regime,
despite the fact that the equations do not include inverse decays. We have then interpreted the bound
on the sum of neutrino masses when considering non-relativistic decays as the intersect between the
non-relativistic decay condition $\Gamma_{\nu} > H(T_{\nu} = m_{\nu}/3)$ and the $2\sigma$ limit derived from our analysis.

In this appendix, we investigate how excluding the relativistic decay regime of parameter space
from the scan can affect the bounds on $\sum m_{\nu}/\text{eV}$ and $\text{Log}_{10}(\Gamma_{\nu}/[\text{km/s/Mpc}])$. As we are interested
in (semi-)relativistic decays, we focus on the parameter space $\text{Log}_{10}(\Gamma_{\nu}/[\text{km/s/Mpc}]) \in [3, 6.5]$.
Our results are presented in Fig. 8. In the 2D plane $\{\text{Log}_{10}\Gamma_{\nu}, \sum m_{\nu}\}$ and below the non-relativistic
line $\Gamma_{\nu} = H(T_{\nu} = m_{\nu}/3)$, we find that imposing the condition directly within the MCMC prior
relaxes the bound by $\sim 10-20\%$. Nevertheless, after marginalizing over $\text{Log}_{10}(\Gamma_{\nu})$, we find that the
'naive' bound coming from the intersect between the non-relativistic line ($\Gamma_{\nu} > H(T_{\nu} = m_{\nu}/3)$) and
the $2\sigma$ limit without priors is in excellent agreement with that coming from imposing this condition
as a prior in the analysis, both yielding $\sum m_{\nu} < 0.42 \text{ eV}$.

References

[1] J. R. Bond, G. Efstathiou, and J. Silk, “Massive Neutrinos and the Large Scale Structure of the Universe,” *Phys. Rev. Lett.*, vol. 45, pp. 1980–1984, 1980. [,61(1980)].

[2] W. Hu, D. J. Eisenstein, and M. Tegmark, “Weighing neutrinos with galaxy surveys,” *Phys. Rev. Lett.*, vol. 80, pp. 5255–5258, 1998.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Posterior distribution of $\sum m_\nu$ and $\text{Log}_{10}(\Gamma_\nu/[\text{km/s/Mpc}])$ when confronted to Planck 2018 + BAO + FS + Pantheon for two different choices of priors on $\Gamma_\nu$ (see legend).</figcaption></figure>

[3] Y. Y. Y. Wong, “Neutrino mass in cosmology: status and prospects,” *Ann. Rev. Nucl. Part. Sci.*, vol. 61, pp. 69–98, 2011.

[4] J. Lesgourgues, G. Mangano, G. Miele, and S. Pastor, *Neutrino Cosmology*. Cambridge University Press, 2 2013.

[5] M. Tanabashi et al., “Review of Particle Physics,” *Phys. Rev. D*, vol. 98, no. 3, p. 030001. 1898 p. 2018.

[6] M. Lattanzi and M. Gerbino, “Status of neutrino properties and future prospects - Cosmological and astrophysical constraints,” *Front.in Phys.*, vol. 5, p. 70, 2018.

[7] N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” *Astron. Astrophys.*, vol. 641, p. A6, 2020.

[8] P. D. Serpico, “Cosmological neutrino mass detection: The best probe of neutrino lifetime,” *Phys. Rev. Lett.*, vol. 98, p. 171301, 2007.

[9] P. D. Serpico, “Neutrinos and cosmology: a lifetime relationship,” *J. Phys. Conf. Ser.*, vol. 173, p. 012018, 2009.

[10] J. F. Beacom, N. F. Bell, and S. Dodelson, “Neutrinoless universe,” *Phys. Rev. Lett.*, vol. 93, p. 121302, 2004.

[11] Y. Farzan and S. Hannestad, “Neutrinos secretly converting to lighter particles to please both KATRIN and the cosmos,” *JCAP*, vol. 1602, no. 02, p. 058, 2016.

[12] J. L. Aalberts et al., “Precision constraints on radiative neutrino decay with CMB spectral distortion,” *Phys. Rev.*, vol. D98, p. 023001, 2018.

[13] G. Barenboim, J. Z. Chen, S. Hannestad, I. M. Oldengott, T. Tram, and Y. Y. Y. Wong, “Invisible neutrino decay in precision cosmology,” *JCAP*, vol. 03, p. 087, 2021.
[14] S. Hannestad and G. Raffelt, “Constraining invisible neutrino decays with the cosmic microwave background,” *Phys. Rev. D*, vol. 72, p. 103514, 2005.

[15] A. Basboll, O. E. Bjaelde, S. Hannestad, and G. G. Raffelt, “Are cosmological neutrinos free-streaming?,” *Phys. Rev. D*, vol. 79, p. 043512, 2009.

[16] M. Archidiacono and S. Hannestad, “Updated constraints on non-standard neutrino interactions from Planck,” *JCAP*, vol. 07, p. 046, 2014.

[17] M. Escudero and M. Fairbairn, “Cosmological Constraints on Invisible Neutrino Decays Revisited,” *Phys. Rev. D*, vol. 100, no. 10, p. 103531, 2019.

[18] Z. Chacko, A. Dev, P. Du, V. Poulin, and Y. Tsai, “Cosmological Limits on the Neutrino Mass and Lifetime,” *JHEP*, vol. 04, p. 020, 2020.

[19] Z. Chacko, A. Dev, P. Du, V. Poulin, and Y. Tsai, “Determining the Neutrino Lifetime from Cosmology,” *Phys. Rev. D*, vol. 103, no. 4, p. 043519, 2021.

[20] J. A. Frieman, H. E. Haber, and K. Freese, “Neutrino Mixing, Decays and Supernova Sn1987a,” *Phys. Lett.*, vol. B200, pp. 115–121, 1988.

[21] A. S. Joshipura, E. Masso, and S. Mohanty, “Constraints on decay plus oscillation solutions of the solar neutrino problem,” *Phys. Rev.*, vol. D66, p. 113008, 2002.

[22] J. F. Beacom and N. F. Bell, “Do solar neutrinos decay?,” *Phys. Rev.*, vol. D65, p. 113009, 2002.

[23] A. Bandyopadhyay, S. Choubey, and S. Goswami, “Neutrino decay confronts the SNO data,” *Phys. Lett.*, vol. B555, pp. 33–42, 2003.

[24] J. M. Berryman, A. de Gouvea, and D. Hernandez, “Solar Neutrinos and the Decaying Neutrino Hypothesis,” *Phys. Rev. D*, vol. 92, no. 7, p. 073003, 2015.

[25] P. Baerwald, M. Bustamante, and W. Winter, “Neutrino Decays over Cosmological Distances and the Implications for Neutrino Telescopes,” *JCAP*, vol. 10, p. 020, 2012.

[26] G. Pagliaroli, A. Palladino, F. L. Villante, and F. Vissani, “Testing nonradiative neutrino decay scenarios with IceCube data,” *Phys. Rev. D*, vol. 92, no. 11, p. 113008, 2015.

[27] M. Bustamante, J. F. Beacom, and K. Murase, “Testing decay of astrophysical neutrinos with incomplete information,” *Phys. Rev. D*, vol. 95, no. 6, p. 063013, 2017.

[28] P. B. Denton and I. Tamborra, “Invisible Neutrino Decay Could Resolve IceCube’s Track and Cascade Tension,” *Phys. Rev. Lett.*, vol. 121, no. 12, p. 121802, 2018.

[29] A. Abdullahi and P. B. Denton, “Visible Decay of Astrophysical Neutrinos at IceCube,” *Phys. Rev. D*, vol. 102, no. 2, p. 023018, 2020.

[30] M. Bustamante, “New limits on neutrino decay from the Glashow resonance of high-energy cosmic neutrinos,” *4 2020*.

[31] M. C. Gonzalez-Garcia and M. Maltoni, “Status of Oscillation plus Decay of Atmospheric and Long-Baseline Neutrinos,” *Phys. Lett.*, vol. B663, pp. 405–409, 2008.

[32] 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.*, vol. B740, pp. 345–352, 2015.

[33] S. Choubey, D. Dutta, and D. Pramanik, “Invisible neutrino decay in the light of NOvA and T2K data,” *JHEP*, vol. 08, p. 141, 2018.
[34] B. Aharmim et al., "Constraints on Neutrino Lifetime from the Sudbury Neutrino Observatory," Phys. Rev., vol. D99, no. 3, p. 032013, 2019.

[35] C. S. Lorenz, L. Funcke, M. Löffler, and E. Calabrese, "Reconstruction of the neutrino mass as a function of redshift," Phys. Rev. D, vol. 104, no. 12, p. 123518, 2021.

[36] M. Escudero, J. Lopez-Pavon, N. Rius, and S. Sandner, "Relaxing Cosmological Neutrino Mass Bounds with Unstable Neutrinos," JHEP, vol. 12, p. 119, 2020.

[37] S. Bashinsky and U. Seljak, "Neutrino perturbations in CMB anisotropy and matter clustering," Phys. Rev. D, vol. 69, p. 083002, 2004.

[38] B. Audren et al., "Robustness of cosmic neutrino background detection in the cosmic microwave background," JCAP, vol. 1503, p. 036, 2015.

[39] 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., vol. 115, no. 9, p. 091301, 2015.

[40] D. Baumann, D. Green, J. Meyers, and B. Wallisch, "Phases of New Physics in the CMB," JCAP, vol. 01, p. 007, 2016.

[41] J. Angrik et al., "KATRIN design report 2004," 2005.

[42] M. Gerbino, M. Lattanzi, O. Mena, and K. Freese, "A novel approach to quantifying the sensitivity of current and future cosmological datasets to the neutrino mass ordering through Bayesian hierarchical modeling," Phys. Lett. B, vol. 775, pp. 239-250, 2017.

[43] A. Caldwell, M. Ettengruber, A. Merle, O. Schulz, and M. Totzauer, "Global Bayesian analysis of neutrino mass data," Phys. Rev. D, vol. 96, no. 7, p. 073001, 2017.

[44] S. Vagnozzi, E. Giusarma, O. Mena, K. Freese, M. Gerbino, S. Ho, and M. Lattanzi, "Unveiling $\nu$ secrets with cosmological data: neutrino masses and mass hierarchy," Phys. Rev. D, vol. 96, no. 12, p. 123503, 2017.

[45] F. Simpson, R. Jimenez, C. Pena-Garay, and L. Verde, "Strong Bayesian Evidence for the Normal Neutrino Hierarchy," JCAP, vol. 06, p. 029, 2017.

[46] E. Di Valentino, S. Garciazo, and O. Mena, "Most constraining cosmological neutrino mass bounds," Phys. Rev. D, vol. 104, no. 8, p. 083504, 2021.

[47] R. Jimenez, C. Pena-Garay, K. Short, F. Simpson, and L. Verde, "Neutrino Masses and Mass Hierarchy: Evidence for the Normal Hierarchy," 3 2022.

[48] T. Schwetz, K. Freese, M. Gerbino, E. Giusarma, S. Hannestad, M. Lattanzi, O. Mena, and S. Vagnozzi, "Comment on "Strong Evidence for the Normal Neutrino Hierarchy"," 3 2017.

[49] S. Garciazo, M. Archidiacono, P. F. de Salas, O. Mena, C. A. Ternes, and M. Tórtola, "Neutrino masses and their ordering: Global Data, Priors and Models," JCAP, vol. 03, p. 011, 2018.

[50] L. T. Hergt, W. J. Handley, M. P. Hobson, and A. N. Lasenby, "Bayesian evidence for the tensor-to-scalar ratio $r$ and neutrino masses $m_\nu$: Effects of uniform vs logarithmic priors," Phys. Rev. D, vol. 103, p. 123511, 2021.

[51] S. Garciazo et al., "Neutrino mass and mass ordering: No conclusive evidence for normal ordering," 5 2022.
[52] N. Palanque-Delabrouille, C. Yèche, N. Schöneberg, J. Lesgourgues, M. Walther, S. Chabanier, and E. Armengaud, "Hints, neutrino bounds and WDM constraints from SDSS DR14 Lyman-α and Planck full-survey data," JCAP, vol. 04, p. 038, 2020.

[53] C.-P. Ma and E. Bertschinger, “Cosmological perturbation theory in the synchronous versus conformal Newtonian gauge,” *J 1994*.

[54] V. Poulin, P. D. Serpico, and J. Lesgourgues, “A fresh look at linear cosmological constraints on a decaying dark matter component,” *JCAP*, vol. 08, p. 036, 2016.

[55] N. Blinov, C. Keith, and D. Hooper, “Warm Decaying Dark Matter and the Hubble Tension,” *JCAP*, vol. 06, p. 005, 2020.

[56] 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*, vol. 02, p. 001, 2013.

[57] T. Brinckmann and J. Lesgourgues, “MontePython 3: boosted MCMC sampler and other features,” *Phys. Dark Univ.*, vol. 24, p. 100260, 2019.

[58] F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, “The 6dF Galaxy Survey: baryon acoustic oscillations and the local Hubble constant,” *MNRAS*, vol. 416, pp. 3017–3032, Oct. 2011.

[59] A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, “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.*, vol. 449, no. 1, pp. 835–847, 2015.

[60] S. Alam et al. *Mon. Not. Roy. Astron. Soc.*, vol. 470, no. 3, pp. 2617–2652, 2017.

[61] V. de Sainte Agathe et al., “Baryon acoustic oscillations at z = 2.34 from the correlations of Lyα absorption in eBOSS DR14,” *Astron. Astrophys.*, vol. 629, p. A85, 2019.

[62] M. Blomqvist et al., “Baryon acoustic oscillations from the cross-correlation of Lyα absorption and quasars in eBOSS DR14,” *Astron. Astrophys.*, vol. 629, p. A86, 2019.

[63] 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.*, vol. 859, no. 2, p. 101, 2018.

[64] A. Gelman and D. B. Rubin *Statist. Sci.*, vol. 7, pp. 457–472, 1992.

[65] M. Archidiacono, T. Brinckmann, J. Lesgourgues, and V. Poulin, “Physical effects involved in the measurements of neutrino masses with future cosmological data,” *JCAP*, vol. 02, p. 052, 2017.

[66] A. Font-Ribera, P. McDonald, N. Mostek, B. A. Reid, H.-J. Seo, and A. Slosar, “DESI and other dark energy experiments in the era of neutrino mass measurements,” *JCAP*, vol. 05, p. 023, 2014.

[67] L. Amendola et al., “Cosmology and fundamental physics with the Euclid satellite,” *Living Rev. Rel.*, vol. 21, no. 1, p. 2, 2018.

[68] A. Liu, J. R. Pritchard, R. Allison, A. R. Parsons, U. Seljak, and B. D. Sherwin, “Eliminating the optical depth nuisance from the CMB with 21 cm cosmology,” *Phys. Rev. D*, vol. 93, no. 4, p. 043013, 2016.

[69] R. Maartens, F. B. Abdalla, M. Jarvis, and M. G. Santos, “Overview of Cosmology with the SKA,” PoS, vol. AASKA14, p. 016, 2015.