# Linear cosmological constraints on 2-body decaying dark matter scenarios
and the $S_{8}$ tension

Guillermo F. Abellán  Electronic address: guillermo.franco-
<EMAIL_ADDRESS>Riccardo Murgia  Electronic address:
<EMAIL_ADDRESS>Vivian Poulin  Electronic address:
<EMAIL_ADDRESS>Laboratoire Univers & Particules de Montpellier
(LUPM), Université de Montpellier (UMR-5299)
Place Eugène Bataillon, F-34095 Montpellier Cedex 05, France

###### Abstract

The ‘$S_{8}$ tension’ is a longstanding discrepancy between the cosmic
microwave background (CMB) and weak gravitational lensing determination of the
amplitude of matter fluctuations, parametrized as
$S_{8}\equiv\sigma_{8}(\Omega_{m}/0.3)^{0.5}$, where $\sigma_{8}$ is the root
mean square of matter fluctuations on a 8 $h^{-1}$Mpc scale, and $\Omega_{m}$
is the total matter abundance. It was recently shown that dark matter (DM)
decaying into a massless (dark radiation) and a massive (warm DM) species,
with a lifetime $\Gamma^{-1}\simeq 55~{}(\varepsilon/0.007)^{1.4}$ Gyrs –
where $\varepsilon$ represent the mass-energy fraction transferred to the
massless component – can ease the tension. Thanks to a fast and accurate fluid
approximation scheme for the warm species, we perform a comprehensive study of
this 2-body decaying DM scenario, discussing in detail its dynamics and its
impact on the CMB and linear matter power spectra. We then investigate the
implications for the ‘$S_{8}$ tension’ against a number of changes in the
analysis: different $S_{8}$ priors, marginalization over the lensing
information in  _Planck_ data, trading _Planck_ high$-\ell$ polarization data
for those from the SPTpol and ACTPol surveys, and the inclusion of the recent
results from the Xenon1T collaboration. We conclude that the preference for
decaying DM, apparent only when the $S_{8}$ value determined from weak lensing
data is added to the analysis, does not sensibly degrade the fit to any of the
cosmological data-sets considered, and that the model could potentially
explain the anomalous electron recoil excess reported by the Xenon1T
collaboration. Furthermore, we explictly show that while current CMB data
alone are not sensitive enough to distinguish between standard CDM and
decaying DM, next-generation CMB observations (CMB-S4) can unambiguously
detect its signature.

††preprint: APS/123-QED

##  I Introduction

In the last couple of decades, the so-called standard $\Lambda$ Cold Dark
Matter ($\Lambda$CDM) model of cosmology has been firmly established as the
most successful framework to interpret numerous independent experimental
observations up to a very high degree of accuracy. Still, the nature of its
dominant components - Cold Dark Matter (CDM) and Dark Energy (DE) – is yet to
be unveiled. The $\Lambda$CDM model provides indeed a remarkable fit to a wide
variety of early universe data, such as Cosmic Microwave Background (CMB) and
Big Bang Nucleosynthesis (BBN), as well as late universe observations such as
Baryon Acoustic Oscillation (BAO), and luminosity distance to SuperNovae of
type Ia (SNIa). However, as the accuracy of the measurements has increased
over the past few years, a number of intriguing discrepancies between the
values of some cosmological parameter as inferred within $\Lambda$CDM and
their direct measurements at low redshift, has emerged. At the heart of this
study is the mild ($\sim 2-3\sigma$) yet longstanding tension between the CMB
[1] and weak gravitational lensing [2, 3, 4] determination of the amplitude of
the matter fluctuations, tipically parametrized – in the context of weak
lensing surveys – as $S_{8}\equiv\sigma_{8}(\Omega_{m}/0.3)^{0.5}$, where
$\sigma_{8}$ is the root mean square of matter fluctuations on a 8 $h^{-1}$Mpc
scale, and $\Omega_{m}$ is the total matter abundance. Another example is the
so-called Hubble tension [5, 6, 7, 8, 9] , i.e., the significant discrepancy
($\sim 4~{}-~{}5~{}\sigma$ C.L.) between the value of the current expansion
rate of the universe $H_{0}$, directly measured using SNIa data as a cosmic
distance ladder [10], and that inferred from CMB data [1], assuming
$\Lambda$CDM. Despite scrupulous efforts to minimize the systematic errors at
play in the cepheid calibration of the SNIa, the statistical significance of
the disagreement is steadily rising [11, 12, 13, 14, 9]. However, it is still
debated whether the Hubble constant resulting from a calibration of the SNIa
on the ‘tip of the red giant branch’ shows a similar degree of tension with
respect to _Planck_ $\Lambda$CDM model [8, 14, 15, 16]. Alternative methods
have been proposed, but are currently not at the accuracy level required to
unambiguously weigh in the Hubble tension [5].

In the absence of convincing solutions within the standard cosmological model,
and driven by the fact that the nature of the dark sector is unknown,
throughout the years many alternative scenarios have been proposed to explain
these discrepancies. It is not yet clear if introducing new physics in the
pre-recombination era to decrease the sound horizon at recombination could
fully resolve the Hubble tension without spoiling other bounds or exacerbating
the $S_{8}$ tension [17, 18, 19, 20, 21]. On the other hand, resolving the
$S_{8}$ tension requires to decrease the amplitude of matter fluctuations on
scales $k\sim 0.1-1~{}h$/Mpc, which can be easily achieved in a variety of
models departing from $\Lambda$CDM only at late times, often related to new DM
properties [22, 23, 24, 25], including the possibility that DM decays [26, 27,
28, 29, 30].

In this paper, we reassess the phenomenology of a 2-body Decaying Cold Dark
Matter (DCDM) scenario, where the decay products are one massive Warm DM (WDM)
particle and one (massless) DR component, interacting only through gravitation
with the standard model particles. We will refer to the full model as
$\Lambda$DDM. From the point of view of particle physics model building, the
stability over cosmological timescales is one of the most peculiar property of
the dark matter particle, reviewed e.g. in Ref. [31]. Often, an additional
symmetry (typically a discrete $Z_{2}$ symmetry) has to be assumed to make the
DM candidate stable. Nevertheless, DM decays at late-times are known
signatures of many models in the literature such as (for instance) models with
R-parity violation [32, 33], super Weakly Interacting Massive particles (super
WIMPs) [34, 35, 36, 37], sterile neutrinos [38, 39] or models with an
additional U(1) gauge symmetry [40, 41, 42, 43]. Recently, the authors of Ref.
[44] engineered a model in the context of supergravity that explicitely
reproduces the kind of late 2-body decays considered here, and also provides a
natural explanation for the small mass splitting that seems to be favoured by
cosmic data.

Decays to electromagnetically charged particles are tightly constrained by
_Planck_ data [45, 46], $\gamma$-ray [47, 48] and cosmic-ray searches [49,
50], typically requiring $\Gamma^{-1}\\!>\\!{\cal O}(10^{26}$s), with some
level of dependence on the decay channel. Still, a purely gravitational
constraints, although weaker, is very interesting in the spirit of being
‘model-independent’, while applying to models with decay to a dark sector, or
to (non-interacting) neutrinos. The canonical example is perhaps that of the
keV majoron [51, 52] decaying into relativistic neutrinos. Models of CDM
decays with massive daughters have also been invoked as a potential solution
to the observational discrepancies with CDM on small (sub-galactic) scales
after structure formation (e.g. [53, 54, 55, 56, 57, 58, 59, 60, 42]). Even
more recently, decaying dark matter models were proposed [61, 62, 43, 63] as a
way to explain the excess of events in the electronic recoils reported by the
Xenon1T collaboration [64]. In the literature, most studies restricted
themselves to massless daughter particles [65, 26, 66, 46, 67], with the
benefit of greatly simplifying the cosmological analysis, but limiting the
true ‘model-independence’ of the bound, and therefore its robustness.
Nevertheless, some studies have attempted at including the effect of massive
daughters in a cosmological context but either neglected cosmological
perturbations of the daughter particles [28, 30, 29] or were limited by
computational power to perform a complete analysis against a host of
cosmological data [68, 69].

It has originally been suggested that DM decaying into massless daughters
could help with cosmological tensions [26, 66], but careful analysis of this
scenario in light of _Planck_ 2015 data has excluded this possibility [70, 27]
(although see Ref. [71] for a different conclusion if the decay rate is not
constant). Attempting to go beyond these studies, the authors of Ref. [28]
suggested that considering a non-zero mass for (at least one of) the decay
product would affect the phenomenology and allow for a resolution of the
Hubble tension. However, a recent study based on a combination of both BAO and
uncalibrated SNIa data-sets has been carried out in Ref. [29]. As opposed to
Ref. [28], they conclude that a $\Lambda$DDM scenario does not predict higher
$H_{0}$ values. This is in good agreement with model-independent analyses
existing in the literature in which it has been established that a combination
of BAO and uncalibrated SNIa data strongly constrain any late-time
modification as a resolution to the Hubble tension (see, e.g, Refs. [72, 17]).
A similar conclusion is also reached when CMB data are considered [30]. Yet,
these recent analyses were limited to the study of the effects of $\Lambda$DDM
on the background evolution of the universe.

In this paper, we perform a thorough analysis of the $\Lambda$DDM model in
light of up-to-date low- and high-redshift data-sets, including the effects of
linear perturbations. To that end, we introduce a new approximation scheme
that allows to accurately and quickly compute the dynamics of the WDM linear
perturbations by treating the WDM species as a viscous fluid. In a companion
paper [73], thanks to this new fluid approximation, we have shown that the
$\Lambda$DDM, while unable to ease the Hubble tension, can fully explain the
low-$S_{8}$ measurement from recent weak lensing surveys. The aim of this
paper is to: i) introduce the cosmological formalism of the $\Lambda$DDM model
and a new approximation scheme that we developed to accurately describe linear
perturbations of the warm daughter; ii) discuss the background and
perturbation dynamics of the $\Lambda$DDM model and its impact on the CMB and
linear matter power spectrum; iii) compare the constraints obtained with the
inclusion of perturbations to those obtained when neglecting them (as was done
in the past literature); iv) test the robustness of the $\Lambda$DDM
resolution to the $S_{8}-$tension to a number of changes in the analysis
(different $S_{8}$ priors, different CMB datasets, marginalization over the
lensing information in  _Planck_ , including constraints from the Xenon1T
experiment on the model).

This work is structured as follows: in Section II we introduce the formalism
and the novel approximation scheme; in Section III we illustrate the
$\Lambda$DDM impact on the relevant cosmological observables; in Section IV we
discuss the results of our data analyses; in Section V we show that a next
generation CMB experiment (CMB-S4) can detect DDM; finally, in Section VI we
draw our conclusions.

##  II  Formalism of 2-body $\Lambda$DDM

Hereafter, we adopt the Boltzmann formalism by Ref. [69], where the time-
evolution of the Phase-Space Distribution (PSD) for both the mother and the
daughter particles were derived. However, in Section II.1 we explicitly show
that, at the background level, such formalism is equivalent to the one by
[74].

While the (cold) parent particle can be safely described as a perfect fluid,
computing the density perturbation evolution for the daughter particles
requires a more sophisticated treatment. The central role in the game is
played by $\ell_{\rm max}$, i.e. the highest multiple to consider when drawing
up the hierarchy of equations describing the PSD of the daughter particles. In
the massless case, the degrees of freedom associated to momentum can be
removed after the PSD multipole decomposition [27]. Due to its non-trivial
energy-momentum relation, this approach is not possible for the warm daughter
particle. One has to study the full PSD evolution, which would be
computationally prohibitive when performing MCMC analyses. For this reason, in
Section II.2 we provide a detailed description of a novel approximation
scheme, devised in Ref. [73], based on describing the WDM component as a
viscous fluid on sub-Hubble scales. This allows us to integrate out the
momentum degrees of freedom and the hierarchy of equations to be truncated at
$\ell_{\rm max}=1$. We will show that the new, computationally faster scheme
is accurate enough to be used for cosmological analyses, allowing to establish
accurate and robust CMB limits on this class of models.

Our framework is characterized by two additional free parameters with respect
to $\Lambda$CDM: the DCDM lifetime, $\Gamma^{-1}$, and the fraction of DCDM
rest mass energy converted into DR, defined as follows [74]:

$\varepsilon=\frac{1}{2}\left(1-\frac{m^{2}_{\rm{wdm}}}{m^{2}_{\rm{dcdm}}}\right),$
(1)

where $0\leq\varepsilon\leq 1/2$. The lower limit corresponds to the standard
CDM case, so that $\Omega_{\rm{cdm}}=\Omega_{\rm{dcdm}}+\Omega_{\rm{wdm}}$,
whereas $\varepsilon=1/2$ corresponds to DM decaying solely into DR. In
general, small $\varepsilon$ values (i.e. heavy massive daughters) and small
$\Gamma$ values (i.e. lifetimes much longer than the age of the universe)
induce little departures from $\Lambda$CDM.

Let us decompose the PSD function of the $j$-th dark component into a
background contribution $\bar{f}_{j}$ plus a linear perturbation $\Delta
f_{j}$ as

$f_{j}(k,q,\mu,\tau)=\bar{f}_{j}(q,\tau)+\Delta f_{j}(k,q,\mu,\tau),$ (2)

where $\tau$ is the conformal time111We use dots to indicate derivatives with
respect to conformal time and $j=$ $\\{{\rm{dcdm}},{\rm{dr}},{\rm{wdm}}\\}$.
The mean energy density and pressure $\bar{\rho}_{j}$ are obtained by
integrating the background PSD, i.e.

$\displaystyle\bar{\rho_{j}}$
$\displaystyle=\frac{1}{a^{4}}\int_{0}^{\infty}dq\ 4\pi
q^{2}\mathcal{E}_{j}\bar{f}_{j},$ (3) $\displaystyle\bar{P}_{j}$
$\displaystyle=\frac{1}{3a^{4}}\int_{0}^{\infty}dq\ 4\pi
q^{2}\frac{q^{2}}{\mathcal{E}_{j}}\bar{f}_{j},$ (4)

where $\mathcal{E}_{j}=\sqrt{m_{j}^{2}a^{2}+q^{2}}$ is the comoving energy of
$j$-th species. The linear perturbation term $\Delta f_{j}$ is generally
expanded over Legendre polynomials:

$\Delta f_{j}(k,q,\mu,\tau)=\sum_{\ell=0}^{\infty}(-i)^{\ell}(2\ell+1)\Delta
f_{j,\ell}(k,q,\tau)P_{\ell}(\mu).$ (5)

The perturbed energy density, pressure, energy flux and shear stress are thus
given by:

$\displaystyle\delta\rho_{j}=\bar{\rho_{j}}\delta_{j}$
$\displaystyle=\frac{4\pi}{a^{4}}\int_{0}^{\infty}dqq^{2}\mathcal{E}_{j}\Delta
f_{j,0},$ (6) $\displaystyle\delta P_{j}=\bar{\rho}_{j}\Pi_{j}$
$\displaystyle=\frac{4\pi}{3a^{4}}\int_{0}^{\infty}dqq^{2}\frac{q^{2}}{\mathcal{E}_{j}}\Delta
f_{j,0},$ (7) $\displaystyle(\bar{\rho_{j}}+\bar{P}_{j})\theta_{j}$
$\displaystyle=\frac{4\pi k}{a^{4}}\int_{0}^{\infty}dqq^{2}q\Delta f_{j,1},$
(8) $\displaystyle(\bar{\rho_{j}}+\bar{P}_{j})\sigma_{j}$
$\displaystyle=\frac{8\pi}{3a^{4}}\int_{0}^{\infty}dqq^{2}\frac{q^{2}}{\mathcal{E}_{j}}\Delta
f_{j,2}.$ (9)

We choose to work in the synchronous gauge, co-moving with the parent
particle, that means that each $\Delta f_{{\rm{dcdm}},\ell}$ vanishes except
for the term $\Delta f_{{\rm{dcdm}},0}$. In addition, since the parent
particle is non-relativistic, the only relevant dynamical DCDM variables are
$\bar{\rho}_{\rm{dcdm}}$ and $\delta_{\rm{dcdm}}$. We follow the standard
convention of referring to the scalar metric perturbations in this gauge as
$h$ and $\eta$. We implement the equations describing the $\Lambda$DDM model
in a modified version of the Boltzmann code CLASS222http://class-code.net/[75,
76]. Our code is available at https://github.com/PoulinV/class_majoron.

Throughout the rest of this paper and unless stated otherwise, we compare
$\Lambda$DDM models at fixed $\omega_{\rm{dcdm}}^{\rm
ini}\equiv\Omega_{\rm{dcdm}}^{\rm{ini}}h^{2}$ with $\Lambda$CDM models having
the same $\omega_{\rm cdm}$. All other cosmological parameters are fixed to
$\\{H_{0}=67.7\ \rm{km}/\rm{s}/\rm{Mpc}$, $\omega_{\rm
cdm}=\omega_{\rm{dcdm}}^{\rm ini}=0.1194$, $\omega_{\rm b}=0.0224$,
$n_{s}=0.9673$, $\text{ln}(10^{10}A_{s})=3.052$, $\tau_{\rm reio}=0.0582\\}$,
which constitutes our baseline $\Lambda$CDM model. These values correspond to
the best-fit from the combined analysis (including $S_{8}$ data from weak
lensing) of Ref. [73], which are however very similar to the $\Lambda$CDM
best-fit parameters from _Planck_ 2018 [1].

###  II.1 Background dynamics

The background evolution of the mother and daughter particles is described at
the phase-space level by their Boltzmann equations [69]:

$\displaystyle\dot{\bar{f}}_{\rm{dcdm}}$
$\displaystyle=-a\Gamma\bar{f}_{\rm{dcdm}},$ (10)
$\displaystyle\dot{\bar{f}}_{\rm{dr}}=\dot{\bar{f}}_{\rm{wdm}}$
$\displaystyle=\frac{a\Gamma\bar{N}_{\rm{dcdm}}(\tau)}{4\pi
q^{2}}\delta(q-ap_{\rm max}).$ (11)

$\bar{N}_{\rm{dcdm}}=(\Omega_{\rm{dcdm}}^{\rm{ini}}\rho_{c,0}/m_{\rm{dcdm}})e^{-\Gamma
t}$ is the DCDM (comoving) number
density333$\Omega^{0}_{\rm{dcdm}}\equiv\Omega_{\rm{dcdm}}^{\rm{ini}}e^{-\Gamma
t_{0}}$; $\Omega_{\rm{dcdm}}^{\rm{ini}}$ is the initial DCDM abundance,
$t_{0}$ the age of the universe.;
$p_{\rm{max}}=|\vec{p}_{\rm{wdm}}|=|\vec{p}_{\rm{dr}}|=m_{\rm{wdm}}\varepsilon/\sqrt{1-2\varepsilon}$
is the initial momentum of the daughter particles.

By combining Eqs. (3) and (4) with Eqs. (10) and (11) we obtain:

$\displaystyle\dot{\bar{\rho}}_{\rm{dcdm}}$
$\displaystyle=-3\mathcal{H}\bar{\rho}_{\rm{dcdm}}-a\Gamma\bar{\rho}_{\rm{dcdm}},$
(12) $\displaystyle\dot{\bar{\rho}}_{\rm{dr}}$
$\displaystyle=-4\mathcal{H}\bar{\rho}_{\rm{dr}}+\varepsilon
a\Gamma\bar{\rho}_{\rm{dcdm}},$ (13)
$\displaystyle\dot{\bar{\rho}}_{\rm{wdm}}$
$\displaystyle=-3(1+w)\mathcal{H}\bar{\rho}_{\rm{wdm}}+(1-\varepsilon)a\Gamma\bar{\rho}_{\rm{dcdm}}.$
(14)

$\mathcal{H}\equiv\dot{a}/a$ is the conformal Hubble parameter, and
$w(\tau)\equiv\bar{P}_{\rm{wdm}}/\bar{\rho}_{\rm{wdm}}$ is the WDM Equation of
State (EoS)444Notice that this expression does not coincide with the EoS used
in [28].

Eq. (14) will be useful to analytically derive the fluid equations that we
present in Section II.3, but for a numerical resolution it is much simpler to
use an integral formula for $\bar{\rho}_{\rm{wdm}}$, as it was done in Ref.
[74]. This formula can be obtained by integrating Eq. (11) firstly with
respect to $\tau$, and then with respect to $q$. The first integration
requires using the relation $\delta(q-ap_{\rm
max})=\delta(\tau-\tau_{q})/q\mathcal{H}$, where $\tau_{q}$ represents the
conformal time when daughter particles with co-moving momentum $q$ are born,
$q=a(\tau_{q})p_{\rm{max}}$ [69]. For the second integral, changing the
integration variable from $q$ to $a_{q}=a(\tau_{q})$ leads to

$\displaystyle\bar{\rho}_{\rm{wdm}}(a)$
$\displaystyle=\frac{C}{a^{4}}\int_{0}^{a}da_{q}\frac{e^{-\Gamma
t_{q}}}{\mathcal{H}_{q}}\sqrt{\varepsilon^{2}a_{q}^{2}+(1-2\varepsilon)a^{2}},$
(15)

where $C\equiv\rho_{c,0}\Omega_{\rm{dcdm}}^{\rm{ini}}\Gamma$,
$\mathcal{H}_{q}\equiv\mathcal{H}(a_{q})$, $t_{q}\equiv t(a_{q})$. Note the
equivalence between Eq. (15) and the analogous expression derived in Ref. [74]
with a different formalism. Concerning the massless DR species, we simply take
the limit $\varepsilon\rightarrow 1/2$ of Eq. (15).

In both cases, the background evolution Eq. (15) needs to be solved
iteratively, as the Hubble parameter $\mathcal{H}$ depends on
$\bar{\rho}_{\rm{wdm}}$ and $\bar{\rho}_{\rm{dr}}$ through the Friedmann
equation. For a flat universe, this equation reads

$\mathcal{H}^{2}(a)=\frac{8\pi Ga^{2}}{3}\sum_{i}\bar{\rho}_{i}(a),$ (16)

where

$\displaystyle\sum_{i}\bar{\rho}_{i}(a)$
$\displaystyle=\bar{\rho}_{\rm{dcdm}}(a)+\bar{\rho}_{\rm{dr}}(a)+\bar{\rho}_{\rm{wdm}}(a)$
$\displaystyle+\bar{\rho}_{\gamma}(a)+\bar{\rho}_{\nu}(a)+\bar{\rho}_{b}(a)+\bar{\rho}_{\Lambda}.$
(17)

Here $\bar{\rho}_{\gamma}$, $\bar{\rho}_{\nu}$, $\bar{\rho}_{b}$ and
$\bar{\rho}_{\Lambda}$ denote the mean densities of photons, neutrinos,
baryons and dark energy, respectively.

Figure 1: Upper $-$ Redshift evolution of the abundances of the DCDM, WDM and
DR species, assuming $\Gamma^{-1}=30\ \text{Gyrs}$ and $\varepsilon=0.1$. We
also show the abundance for a standard CDM species with $\Omega^{0}_{\rm
cdm}=\Omega_{\rm{dcdm}}^{\rm{ini}}$. Lower $-$ Redshift evolution of the
Hubble parameter for the $\Lambda$CDM and four different $\Lambda$DDM models.
The Hubble parameter today is fixed to $H_{0}=67.7\ \rm{km}/\rm{s}/\rm{Mpc}$.

We show the evolution of the Hubble parameter for $\Lambda$CDM and several
$\Lambda$DDM models in lower panel of Fig. 1. In the upper panel of Fig. 1, we
also show the evolution of the fractional densities $\Omega_{{\rm{dcdm}}}(z)$,
$\Omega_{{\rm{wdm}}}(z)$ and $\Omega_{{\rm{dr}}}(z)$ for a particular
$\Lambda$DDM model ($\Gamma^{-1}=30\ \text{Gyrs}$ and $\varepsilon=0.1$) as
well as $\Omega_{\rm cdm}(z)$ with the same initial amount of dark matter,
namely $\Omega^{0}_{\rm cdm}=\Omega_{\rm{dcdm}}^{\rm{ini}}$. The two new
parameters affect $\mathcal{H}(a)$ as follows: at fixed $\varepsilon$, a
shorter lifetime $\Gamma^{-1}$ implies a lower Hubble parameter. This is
clearly visible at $z\gtrsim 1$. The behavior below $z\sim 1$ is due to our
assumption of a flat universe and the requirement of fixing $H_{0}$: to
fulfill the budget equation, a lower DM abundance requires a larger
$\Omega_{\Lambda}$, meaning that the period of accelerated expansion begins
earlier with respect to $\Lambda$CDM. At fixed lifetime $\Gamma^{-1}$, a
larger $\varepsilon$ also induces a lower $\mathcal{H}(a)$, because more
energy is converted into radiation, which dilutes faster. We thus anticipate a
negative correlation between $\Gamma^{-1}$ and $\varepsilon$ at the background
level. The degeneracy can be captured by the matter contribution from the WDM
species, $\Omega_{\rm{wdm}}^{\rm m}\equiv\Omega_{\rm{wdm}}^{0}(1-3w)$, which
is well constrained by data. This quantity is roughly approximated by

$\Omega_{\rm{wdm}}^{\rm m}\simeq\Omega_{\rm{dcdm}}^{\rm{ini}}(1-e^{-\Gamma
t_{0}})\sqrt{1-2\varepsilon},$ (18)

For small $\Gamma$ and small $\varepsilon$, $\Omega_{\rm{wdm}}^{\rm
m}\propto\Gamma(1-\varepsilon)$ and we expect data to constrain a parameter
combination of the form $\varepsilon\propto\Gamma^{-1}$.

###  II.2 Linear perturbation equations

The continuity equation for the DCDM density perturbation can be obtained
either by using the covariant conservation of the stress-energy tensor
$T_{\mu\nu}$, or by integrating the corresponding first order Boltzmann
equation (see Eq. (2.20) in Ref. [69]). As in the case of standard CDM, the
resulting expression is

$\dot{\delta}_{\rm{dcdm}}=-\frac{\dot{h}}{2}.$ (19)

To obtain the dynamical equations for the daughter particles, we combine Eq.
(5) with the corresponding first order Boltzmann equation (see Eq. (2.31) in
Ref. [69]), so that we have the following hierarchy of equations:

$\displaystyle\frac{\partial\left(\Delta f_{j,0}\right)}{\partial\tau}$
$\displaystyle=-\frac{qk}{\mathcal{E}_{j}}\Delta
f_{j,1}+q\frac{\partial\bar{f}_{j}}{\partial
q}\frac{\dot{h}}{6}+\dot{\bar{f}}_{j}\delta_{\rm{dcdm}},$ (20)
$\displaystyle\frac{\partial\left(\Delta f_{j,1}\right)}{\partial\tau}$
$\displaystyle=\frac{qk}{3\mathcal{E}_{j}}\left[\Delta f_{j,0}-2\Delta
f_{j,2}\right],$ (21) $\displaystyle\frac{\partial\left(\Delta
f_{j,2}\right)}{\partial\tau}$
$\displaystyle=\frac{qk}{5\mathcal{E}_{j}}\left[2\Delta f_{j,1}-3\Delta
f_{j,3}\right]-q\frac{\partial\bar{f}_{j}}{\partial
q}\frac{(\dot{h}+6\dot{\eta})}{15},\ $ (22)
$\displaystyle\frac{\partial\left(\Delta f_{j,\ell}\right)}{\partial\tau}$
$\displaystyle=\frac{qk}{(2l+1)\mathcal{E}_{j}}\left[l\Delta
f_{j,\ell-1}-(\ell+1)\Delta f_{j,\ell+1}\right]$ $\displaystyle\hskip
142.26378pt(\ell\geq 3),$ (23)

where $j=\\{{\rm{dr}}$, ${\rm{wdm}}\\}$. The system of Eqs. (20)-(23) is in
the same form as the one for massive neutrinos [77], except for the last term
in Eq. (20) and the fact that the partial derivative
$\partial\bar{f}_{j}/\partial q$ in Eqs. (20) and (22) is now time-dependent.

Given that the DR species satisfies the condition $q/\mathcal{E}_{\rm{dr}}=1$,
the hierarchy of equations can be simplified by taking the moments

$F_{{\rm{dr}},\ell}\equiv\frac{1}{\rho_{c,0}}\int_{0}^{\infty}dq\ 4\pi
q^{2}q\Delta f_{{\rm{dr}},\ell}$ (24)

of Eqs. (20)-(23), and integrating over all the momentum degrees of freedom,
so that

$\displaystyle\dot{F}_{{\rm{dr}},0}$
$\displaystyle=-kF_{{\rm{dr}},1}-\frac{2}{3}r_{\rm dr}\dot{h}+\dot{r}_{\rm
dr}\delta_{\rm{dcdm}},$ (25) $\displaystyle\dot{F}_{{\rm{dr}}.1}$
$\displaystyle=\frac{k}{3}F_{{\rm{dr}},0}-\frac{2k}{3}F_{{\rm{dr}}.2},$ (26)
$\displaystyle\dot{F}_{{\rm{dr}},2}$
$\displaystyle=\frac{2k}{5}F_{{\rm{dr}},1}-\frac{3k}{5}F_{{\rm{dr}},3}+\frac{4}{15}r_{\rm
dr}(\dot{h}+6\dot{\eta}),$ (27) $\displaystyle\dot{F}_{{\rm{dr}},\ell}$
$\displaystyle=\frac{k}{(2\ell+1)}\left[\ell
F_{{\rm{dr}},\ell-1}-(\ell+1)F_{{\rm{dr}},\ell+1}\right]\ \ \ \ \ \ (\ell\geq
3).$ (28)

We have adopted the convention $r_{\rm{dr}}\equiv
a^{4}\bar{\rho}_{\rm{dr}}/\rho_{c,0}$, as in Ref. [27], which in the 2-body
decay scenario under study leads to:

$\dot{r}_{\rm{dr}}=\varepsilon
a\Gamma(\bar{\rho}_{\rm{dcdm}}/\bar{\rho}_{\rm{dr}})r_{\rm{dr}}.$ (29)

The first three multipoles are given by:

$F_{{\rm{dr}},0}=r_{\rm{dr}}\delta_{\rm{dr}},\hskip
5.69054ptF_{{\rm{dr}},1}=\frac{4r_{\rm{dr}}}{3k}\theta_{\rm{dr}},\hskip
5.69054ptF_{{\rm{dr}},2}=2\sigma_{\rm{dr}}r_{\rm{dr}}.$ (30)

We choose the maximum multipole $\ell_{\rm max}$ to truncate the hierarchies
of equations according to the scheme proposed in Ref. [77] for both massless
and massive neutrinos, i.e.

$\displaystyle\dot{F}_{{\rm{dr}},\ell_{\rm max}}$
$\displaystyle=kF_{{\rm{dr}},\ell_{\rm max}-1}-\frac{\ell_{\rm
max}+1}{\tau}F_{{\rm{dr}},\ell_{\rm max}},$ (31)
$\displaystyle\Delta\dot{f}_{{\rm{wdm}},\ell_{\rm max}}$
$\displaystyle=\frac{qk\Delta f_{{\rm{wdm}},\ell_{\rm
max}-1}}{\mathcal{E}_{\rm{wdm}}}-\frac{\ell_{\rm max}+1}{\tau}\Delta
f_{{\rm{wdm}},\ell_{\rm max}}.$ (32)

###  II.3 A novel WDM fluid approximation scheme

In order to compute the WDM dynamics one cannot integrate over the momentum
degrees of freedom, as we did in the DR case. Indeed, when taking the moments
of the hierarchy of Eqs. (20)-(23), higher velocity-weight integrals appearing
at $\ell=2$ cannot be computed from the system of equations itself.

Therefore, one has to follow the evolution of the full time-dependent PSD to
obtain the elements of the perturbed stress-energy tensor $\delta_{\rm{wdm}}$,
$\theta_{\rm{wdm}}$ and $\sigma_{\rm{wdm}}$ through Eqs. (6)-(9). A typical
set-up for CMB analyses requires roughly $500$ wavenumbers, $50$ multipoles
and $10^{4}$ momentum bins, i.e. $\mathcal{O}(10^{8})$ linear differential
equations to be computed. On a single processor, this leads to runs with a CPU
time of $1-2$ days per each parameter choice, making a systematic scan of the
parameter space computationally prohibitive.

To overcome the problem, we make use of a new fluid approximation for the WDM
species, introduced in Ref. [73], where the momentum dependence is removed,
and one only needs to track the evolution of the first two multipoles. The
number of linear differential equations to be solved is now reduced to
$\mathcal{O}(10^{3})$, with a CPU time per single run $\sim 30-40$ s. The
novel approximation scheme is based on the treatment of massive neutrinos as a
viscous fluid by Ref. [76], and it is only valid at scales deeply inside the
Hubble radius, where high- and low-$\ell$ modes are effectively decoupled. In
App. A we explicitly demonstrate the accuracy of this approximation.

Similarly to the DR case, the fluid equations can be derived by multiplying
both sides of Eqs. (20) and (21) by $4\pi q^{2}\mathcal{E}_{\rm{wdm}}a^{-4}$
and $4\pi q^{3}ka^{-4}$, respectively, and integrating over $q$. Then, by
using Eqs. (6)-(9) and (14), one can write down the continuity equation,

$\displaystyle\dot{\delta}_{\rm{wdm}}$ $\displaystyle=-3\mathcal{H}(c^{2}_{\rm
s}-w)\delta_{\rm{wdm}}-(1+w)\left(\theta_{\rm{wdm}}+\frac{\dot{h}}{2}\right)$
$\displaystyle+(1-\varepsilon)a\Gamma\frac{\bar{\rho}_{\rm{dcdm}}}{\bar{\rho}_{\rm{wdm}}}(\delta_{\rm{dcdm}}-\delta_{\rm{wdm}}),$
(33)

and the Euler equation,

$\displaystyle\dot{\theta}_{\rm{wdm}}$
$\displaystyle=-\mathcal{H}(1-3c_{g}^{2})\theta_{\rm{wdm}}+\frac{c^{2}_{\rm
s}}{1+w}k^{2}\delta_{\rm{wdm}}-k^{2}\sigma_{\rm{wdm}}$
$\displaystyle-(1-\varepsilon)a\Gamma\frac{1+c_{g}^{2}}{1+w}\frac{\bar{\rho}_{\rm{dcdm}}}{\bar{\rho}_{\rm{wdm}}}\theta_{\rm{wdm}}.$
(34)

Notice that we have introduced the WDM sound speed in the synchronous gauge,
$c^{2}_{\rm s}\equiv\delta P_{\rm{wdm}}/\delta\rho_{\rm{wdm}}$, and the WDM
adiabatic sound speed,
$c_{g}^{2}\equiv\dot{\bar{P}}_{\rm{wdm}}/\dot{\bar{\rho}}_{\rm{wdm}}$.

The latter can be written as

$c_{g}^{2}=w(\dot{\bar{P}}_{\rm{wdm}}/\bar{P}_{\rm{wdm}})(\dot{\bar{\rho}}_{\rm{wdm}}/\bar{\rho}_{\rm{wdm}})^{-1},$
(35)

and it can be computed as follows

$\displaystyle c_{g}^{2}$
$\displaystyle=w\left(5-\frac{\mathfrak{p}_{\rm{wdm}}}{\bar{P}_{\rm{wdm}}}-\frac{\bar{\rho}_{\rm{dcdm}}}{\bar{\rho}_{\rm{wdm}}}\frac{a\Gamma}{3w\mathcal{H}}\frac{\varepsilon^{2}}{1-\varepsilon}\right)$
$\displaystyle\times\left[3(1+w)-\frac{\bar{\rho}_{\rm{dcdm}}}{\bar{\rho}_{\rm{wdm}}}\frac{a\Gamma}{\mathcal{H}}(1-\varepsilon)\right]^{-1}.$
(36)

Here $\mathfrak{p}_{\rm{wdm}}$ denotes the so-called pseudo-pressure, a higher
momenta integral of $\bar{f}_{\rm{wdm}}$ which is reduced to the standard
pressure in the relativistic limit [76].

Obtaining an analytical expression for $c^{2}_{\rm s}$ is less
straightforward, since we do not have a dynamical equation for the pressure
perturbation $\delta P_{\rm wdm}$. In Ref. [76] it is assumed that $c^{2}_{\rm
s}$ is scale-independent and approximately equal to $c_{g}^{2}$. For the WDM
species, we have found that this assumption leads to accurate results for the
CMB power spectrum, but not for the matter power spectrum. In fact,
calculations using the full Boltzmann hierarchy of Eqs. (20)-(23) reveal that
$c_{\rm s}^{2}$ exhibits a particular $k$-dependence that cannot be captured
by a background quantity such as $c_{g}^{2}$. In particular, $c_{\rm s}^{2}$
gets slightly enhanced on scales $k>k_{\rm fs}$, where $k_{\rm fs}$ is the
free-streaming scale of the WDM species, defined as

$k_{\rm
fs}(\tau)\equiv\sqrt{\frac{3}{2}}\frac{\mathcal{H}(\tau)}{c_{g}(\tau)}.$ (37)

It is possible to gain a semi-analytic understanding of this behaviour by
building a formal equation for the evolution of $c^{2}_{\rm s}$, as detailed
in App. B. To account for such an enhancement, we adopt the following
prescription

$\displaystyle c_{\rm
s}^{2}(k,\tau)=c_{g}^{2}(\tau)\left[1+(1-2\varepsilon)T(k/k_{\rm fs})\right],$
(38)

where the function $T(x)=0.2\sqrt{x}$ has been fitted to the sound speed
obtained using the full Boltzmann hierarchy, for the parameter values
$\varepsilon=0.5,0.1,0.01,0.001$ and $\Gamma/H_{0}=0.1,1,10$. The factor
$(1-2\varepsilon)$ is inserted to make the $k$-dependent correction
vanishingly small close to the relativistic limit, where $c^{2}_{\rm s}\simeq
c_{g}^{2}\simeq 1/3$.

In order to trace the evolution of the shear $\sigma_{\rm{wdm}}$ one could
follow a similar approach to that of Ref. [76], where the authors obtained a
dynamical equation for the neutrino shear by means of an improved truncation
scheme at $\ell_{\rm max}=2$. We tested the implementation of a generalization
of that equation suitable to the decaying case, but we found it to be only
relevant close to the relativistic case $\varepsilon\simeq 1/2$, when it
reduces to the DR shear equation from Ref. [26]. In this regime, the dynamics
of the daughter particles do not significantly impact the CMB and matter power
spectra, we thus decided to not include any dynamical equation for the shear
of the WDM species, when switching to the fluid approximation. In practice, we
simply set $\sigma_{\rm{wdm}}$ to a constant value, obtained via integration
of the second PSD multipole in the Boltzmann hierarchy (see Eq. (9)). We
explicitly checked that this approach yields better results rather than simply
setting $\sigma_{\rm{wdm}}=0$ in the fluid equations, or than using the DR
shear equation of Ref. [26] when $\varepsilon=1/2$

Finally, let us recall that, strictly speaking, these equations are valid in
the synchronous gauge co-moving with the DCDM. In practice, however, the
synchronous gauge in CLASS [75] is coded with respect to the CDM.
Nevertheless, for adiabatic initial conditions, one can choose $\theta_{\rm
ini,dcdm}=\theta_{\rm ini,cdm}=0$, such that the gauge co-moving with CDM is
also co-moving with DCDM at all times [65]. Hence, hereinafter we neglect this
irrelevant complication.

###  II.4 Dynamics of perturbations

Before discussing the signatures of varying the parameters $\Gamma$ and
$\varepsilon$ on the relevant cosmological observables, it is worth having a
look at the linear perturbations of the mother and daughter particles. In Fig.
2 we show the evolution of the linear density perturbations for the DCDM, WDM
and DR species, corresponding to a mode that enters the horizon very early
(i.e., $k=1\ \rm{Mpc}^{-1}$). To clarify the impact of the two extra free
parameters, the perturbations are shown for DCDM lifetimes similar and smaller
than the age of the universe ($\Gamma^{-1}=H_{0}^{-1}\sim 14.5\ \text{Gyrs}$
and $\Gamma^{-1}=(10H_{0})^{-1}\sim 1.5\ \text{Gyrs}$) , and for massive
daughters behaving either as warm or cold particles ($\varepsilon=10^{-2}$ and
$\varepsilon=10^{-3}$).

One can see that initially, the perturbations of the daughter species always
track those of the mother, because the coupling term dominates the dynamics
(i.e., the ratios $\rho_{\rm{dcdm}}/\rho_{\rm{wdm}}$ and
$\rho_{\rm{dcdm}}/\rho_{\rm{dr}}$ are large). When a mode crosses the free-
streaming scale, the pressure support of the daughter particles becomes
important and the perturbations develop oscillatory features. For the DR
species, the free-streaming scale simply corresponds to the horizon ($k_{\rm
fs}\sim\mathcal{H}$,) while for the WDM it corresponds to a larger value,
given by $k_{\rm fs}\sim\mathcal{H}/c_{g}$ (see Eq. (37)). This time- and
scale-dependent power suppression, together with the different background
dynamics, lead to key signatures on the CMB and matter power spectra, as we
discuss later. By comparing upper and lower panels of Fig. 2, it is clear that
the value of the WDM free-streaming scale is essentially determined by the
value of $\varepsilon$, as expected. On the other hand, by comparing left and
right panels, one notices that the intensity of the oscillations due to the
pressure support can be compensated by the coupling to the DCDM if the
lifetime is long enough, as the daughter particles keep being produced.

Figure 2: Time evolution of the linear density perturbations of the DCDM, WDM
and DR species, corresponding to a wavenumber $k=1\ \text{Mpc}^{-1}$. Each
panel displays the perturbations for a different combination of the parameters
$\Gamma$ and $\varepsilon$. The black dashed and solid lines indicate the
times of horizon crossing (at $\mathcal{H}(\tau)=k$) and WDM free-streaming
scale crossing (at $k_{\rm fs}(\tau)=k$), respectively. The purple dotted line
indicates the characteristic decay time, given by $t(\tau)=\Gamma^{-1}$.
Figure 3: Same as in Fig. 2, but corresponding to a wavenumber $k=10^{-2}\
\text{Mpc}^{-1}$. In this case, the WDM perturbations never cross the free-
streaming scale.

Interestingly, the decoupling time of the daughter perturbations from the
mother perturbations is always set by the free-streaming crossing time, and
not by the characteristic decay time. In order to illustrate that, in Fig. 3
we show perturbations corresponding to a smaller wave-number, $k=10^{-2}\
\text{Mpc}^{-1}$, that enters the horizon much later. For this mode, the WDM
species does not have time to cross the free-streaming length (the crossing
will occur in the future), so that $\delta_{\rm{dcdm}}$ and
$\delta_{\rm{wdm}}$ remain equal, even if the lifetime is smaller than the age
of the universe. This can be understood from the fluid Eqs. (33)-(34): when
the decay term that includes $\Gamma$ is relevant, $\delta_{{\rm{wdm}}}$ is
driven by $\delta_{\rm dcdm}$. Therefore, the WDM density perturbation
$\delta_{\rm{wdm}}$ will continue to track the behaviour of
$\delta_{\rm{dcdm}}$, as long as the free-streaming scale is not crossed,
i.e., as long as the pressure term, containing
$c_{s}^{2}k^{2}\delta_{\rm{wdm}}$, is small compared to
$\ddot{\delta}_{\rm{wdm}}\sim\delta_{\rm{wdm}}\mathcal{H}^{2}$.

##  III  Observable impact of Dark Matter decays with warm daughters

###  III.1 The linear matter power spectrum

Figure 4:  Left $-$ Residuals of the linear matter power spectrum at $z=0$
(upper) and $z=2$ (lower) for several values of the lifetime
$\Gamma^{-1}=10,30,100,300\ \text{Gyrs}$ and a fixed DR energy fraction
$\varepsilon=0.1$. Residuals are taken with respect our baseline $\Lambda$CDM
model. Right $-$ Same as in the left, but for several DR energy fractions
$\varepsilon=0.5,0.1,0.01,0.001$ and a fixed lifetime $\Gamma^{-1}=30$ Gyrs.

Firstly, we focus on describing the effects of the 2-body decay on the linear
matter power spectrum, since this will allow to better understand some of the
effects on the CMB spectra. In this Section, we use the same $\Lambda$CDM
parameters as in Section II.1, except for fixing $100\theta_{s}=1.04217$
instead of $H_{0}$, to better connect with CMB observations, that accurately
pin down $\theta_{s}$. In Fig. 4, we compare the residual differences in the
linear power spectra (at redshifts $z=0$ and $z=2$) with respect to our
baseline $\Lambda$CDM. The left panel shows several lifetimes and a fixed DR
energy fraction $\varepsilon=0.1$, while the right panel shows a fixed
lifetime $\Gamma^{-1}=30$ Gyrs and several values of the DR energy fraction
$\varepsilon$.

One important feature of the C+WDM scenarios, such as the one considered in
this work, is that they are expected to produce a suppression in the linear
matter power spectrum at scales smaller than a ‘cutoff scale’, with a non-
trivial shape [78, 79, 80, 81]. The cut-off scale is determined by the free-
streaming scale of the WDM species, $k_{\rm fs}$, given by Eq. (37). On scales
$k>k_{\rm fs}$, pressure becomes important and WDM particles cannot stay
confined in gravitational potential wells, which inhibits structure formation.

Fig. 4 clearly illustrates that, while the parameter $\varepsilon$ fixes the
value of the cut-off $k_{\rm fs}$555By looking at Eqs. (36) and (37), we see
that for small values of $\varepsilon$, the cut-off approximately satisfies
the scalings $k_{\rm fs}\propto c_{g}^{-1}\propto
w^{-1/2}\propto\varepsilon^{-1}$., the lifetime $\Gamma^{-1}$ essentially
determines the depth of the suppression at very small scales. This is to be
expected, since the amount of power suppression grows with the WDM abundance,
which increase for smaller lifetimes. In a similar way, the effects of decay
become less important when considering the matter spectrum at a higher
redshift $z=2$, since the abundance of WDM was smaller in the past.

In general, for late-time decay scenarios (well after recombination) as the
one studied in this work, it is possible to distinguish three different
regimes depending on the value of $\varepsilon$, as it is shown in the right
panel of Fig. 4:

  * •

Non-relativistic decay: if $\varepsilon\lesssim 0.001$ (black curve), the WDM
leaves the expansion rate unaffected, since its contribution to the matter
density, $\bar{\rho}_{\rm{wdm}}(1-3\omega)\simeq\bar{\rho}_{\rm{wdm}}$,
compensates the reduction in the DCDM density, $\bar{\rho}_{\rm dcdm}$. In
addition, the WDM free-streaming length is very small, inducing a power
suppression at $k\gtrsim 1\ h\mathrm{Mpc}^{-1}$. Such scales are beyond the
range of scales probed by the observables considered in this work, so in this
regime the WDM is almost degenerate with standard CDM.

  * •

Relativistic decay: if $\varepsilon\simeq 0.5$ (red curve), the WDM component
acts as DR, which can appreciably reduce the expansion rate666Note that since
now we are fixing $100\theta_{s}$ instead of $H_{0}$, the Hubble rate $H(z)$
can increase with respect to $\Lambda$CDM at $z\lesssim 1$, once dark energy
starts to dominate. However, this effect is small for long lifetimes, and at
early times the Hubble rate is still smaller than in $\Lambda$CDM. . Moreover,
the free-streaming length $k_{\rm fs}$ gets as large as the horizon, so that
the WDM does not cluster at all. The reduction in the Hubble friction is
balanced by a reduction in the clustering density of the daughter particles,
$\delta\rho_{\rm{wdm}}\simeq\delta\rho_{\rm{dr}}\simeq 0$, inducing a very
little overall suppression of the growth of fluctuations. However, there is an
impact coming from another background effect: the significant decrease in the
co-moving matter density, $a^{3}\rho_{\rm m}$, leading in turn to an increase
of the angular diameter distance, thereby a reduction of $\theta_{s}$. This
effect can be compensated by increasing $H_{0}$, which shifts the location of
the peak, $k_{\rm eq}/(a_{0}H_{0})$, towards smaller $k$ in the matter power
spectrum – since we are keeping the matter-radiation equality era fixed. The
net effect on the residuals is twofold, a large-scale enhancement and small-
scale suppression of power.

  * •

Warm decay: For intermediate values of $\varepsilon$, namely
$0.001\lesssim\varepsilon\lesssim 0.5$ (green and blue curves), the WDM
component partially contributes to the matter energy density, leaving to an
expansion rate almost unchanged. However, the values of $k_{\rm fs}$ that
determine the cut-off scale in the matter power spectrum are not as small as
in the case of non-relativistic decay, leading to
$\delta\rho_{{\rm{wdm}}}\simeq 0$ for $k>k_{\rm fs}$. Hence, on scales
$k>k_{\rm fs}$ the Hubble friction gets enhanced with respect to the
clustering density $\delta\rho_{\rm tot}$, slowing down the growth of DCDM
perturbations. The net impact on the matter power spectrum is thus a
suppression on intermediate scales, somewhat similar to that induced by
massive active neutrinos [82, 83, 84].

###  III.2 The CMB angular power spectrum

Figure 5:  Left $-$ Residuals (with respect our baseline $\Lambda$CDM model)
of the CMB lensed TT (upper), EE (middle) and lensing potential (lower) power
spectra for several values of the lifetime $\Gamma^{-1}=10,30,100,300\
\text{Gyrs}$ and a fixed DR energy fraction $\varepsilon=0.1$. Right $-$ Same
as in the left, but for several values of the DR energy fraction
$\varepsilon=0.5,0.1,0.01,0.001$ and a fixed lifetime $\Gamma^{-1}=30\
\text{Gyrs}$.

We now discuss the impact of the 2-body decay scenario on the CMB anisotropy
temperature and polarization angular power spectra, as well as on the lensing
potential power spectrum reconstructed from the CMB 4-point correlation
function. In the left panel of Fig. 5 we report the residuals of the (lensed)
TT, EE and lensing potential power spectra with respect to our baseline
$\Lambda$CDM, for different lifetimes $\Gamma^{-1}$ and a fixed DR energy
fraction $\varepsilon=0.1$. The effects, more and more pronounced as the
lifetime decreases, can be understood as follows:

  * •

At the background level, the decay decreases the value of $\Omega_{\rm m}$
with respect to $\Lambda$CDM. This is compensated by an increase in
$\Omega_{\Lambda}$ (earlier beginning of $\Lambda$-domination) and thus an
enhancement in the Late Integrated Sachs-Wolfe (LISW) effect, leaving a
signature in the low-$\ell$ TT power spectrum. Furthermore, a modified
background history alters quantities integrated along $z$, such as $\tau_{\rm
reio}$, which impacts the multipoles $\ell\sim 10$ in the EE power spectrum.

  * •

At the perturbation level, the late-time reduction of
$a^{3}\bar{\rho}_{\rm{dcdm}}$ implies a reduction of the quantity
$a^{2}\delta\rho_{\rm{dcdm}}$, which acts as a source of gravity through the
Poisson equation. This induces a damping in the metric fluctuations, and hence
yields a further enhancement of the LISW effect. Furthermore, the suppression
in the matter power spectrum and in $\Omega_{\rm m}$ lowers the amplitude of
the lensing potential power spectrum, consequently reducing the smoothing of
the peaks in the high-$\ell$ part of both the TT and EE spectra, as one can
see from the ‘wiggles’ in the corresponding plots.

In the right panel of Fig. 5 we show the CMB residuals for a fixed
$\Gamma^{-1}$ and various values of $\varepsilon$. The effects can be readily
understood:

  * •

At the background level, smaller values of $\varepsilon$ weaken the effects
previously discussed, because the decay product dilute in a way similar to
dark matter. Namely, the decrease in $\Omega_{\rm m}$ is less prominent due to
the significant WDM contribution (i.e., the increase in $\Omega_{\Lambda}$ is
shallower), and the impact on $\tau_{\rm reio}$ is smaller. Therefore, the
signatures in the low-$\ell$ part of the TT and EE spectra become less
visible.

  * •

At the level of perturbations, $\varepsilon$ leads to some interesting
signatures on the LISW effect and on the lensing potential. Since the LISW
effect is only relevant for small multipoles $\ell$ (i.e., very large scales),
one just needs to look at wavenumbers such that $k<k_{\rm fs}$. On these
scales, the growth suppression does not play any role, and the decrease in
$a^{2}\delta\rho_{\rm{dcdm}}$ due to the decay gets partially compensated by
the increase in $a^{2}\delta\rho_{\rm{wdm}}$, which is more significant for
smaller values of $\varepsilon$. Thus, the damping in the metric fluctuations
is less relevant for smaller $\varepsilon$, reducing the LISW enhancement.

  * •

Regarding the effects on the lensing potential, one can see that the
suppression in the corresponding power spectrum monotonically decreases for
smaller $\varepsilon$. Naively, one might expect the opposite, since we have
argued that the matter power suppression increases for small $\varepsilon$.
This can be understood by looking at the CMB lensing kernel $W(z)$, given by
the following expression [85]

$\centering W(z)=\frac{3\Omega_{\rm
m}}{2}\frac{H_{0}^{2}}{H(z)}(1+z)\chi(z)\frac{\chi_{\ast}-\chi(z)}{\chi_{\ast}},\@add@centering$
(39)

where $\chi_{\ast}$ is the co-moving distance to the last-scattering surface.
Firstly, the CMB lensing kernel peaks at $z\sim 2$, where the suppression is
less important (see bottom panels of Fig. 4). Secondly, it gets highly
suppressed for higher values of $\varepsilon$, and this effect can dominate
over the effect on the matter power spectrum. This suppression happens mainly
due to the smaller $\Omega_{\rm m}$, as we have verified by computing $W(z)$
for several $\Lambda$DDM models, with and without including the factor
$\Omega_{\rm m}$. Note that coincidentally, the effects on the lensing power
spectrum at high $\ell$ are very similar for $\varepsilon=0.5$ and
$\varepsilon=0.1$: this is because in the former case $\Omega_{\rm m}$ is
highly reduced and the small scales power spectrum is almost unaffected, while
the opposite occurs in the latter case.

Let us finally remark that, even if the effects of varying either
$\varepsilon$ or $\Gamma$ on the observables are different, one can easily
exploit the degeneracy mentioned at the background level to get different
couples of values (large $\Gamma$ and small $\varepsilon$ or vice-versa) with
a similar cosmological signature, especially on the CMB. We do indeed expect
the MCMC analysis to show a negative correlation in the reconstructed 2D
posteriors for $\varepsilon$ and $\Gamma$.

##  IV Data analysis

###  IV.1 Data and method

We now confront the $\Lambda$DDM cosmology to a host of recent cosmological
observations. Our goal is to set constraints on the lifetime of DM and the
mass-ratio of the mother and daugther particles. We also wish to check to what
extent the $\Lambda$DDM model can play a role in resolving cosmological
tensions as claimed in the recent literature [28, 73]. To that end, we perform
comprehensive MCMC analyses with the
MontePython-v3777https://github.com/brinckmann/montepython_public [86, 87]
code interfaced with our modified version of CLASS, considering various
combinations of the following data-sets:

  * •

The BAO measurements from 6dFGS at $z=0.106$ [88], SDSS DR7 at $z=0.15$ [89],
BOSS DR12 at $z=0.38,0.51$ and $0.61$ [90], and the joint constraints from
eBOSS DR14 Ly-$\alpha$ auto-correlation at $z=2.34$ [91] and cross-correlation
at $z=2.35$ [92].

  * •

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$ [90].

  * •

The Pantheon SNIa catalogue, spanning redshifts $0.01<z<2.3$ [93].

  * •

The low-$\ell$ CMB TT, EE, and the high-$\ell$ TT, TE, EE data888In our main
analysis, we use the ‘lite’ version of the clik likelihood. We verify that
this leads to negligible differences with respect to the full likelihood in
App. C. \+ the gravitational lensing potential reconstruction from _Planck_
2018 [1].

  * •

The KIDS-1000+BOSS+2dFLens [94], DES-Y1 [2] and KIDS-1000+Viking+DES-Y1 [4]
weak lensing data, compressed as a a split-normal likelihood, i.e.,
$S_{8}=0.766^{+0.02}_{-0.014}$, $S_{8}=0.773_{-0.02}^{+0.026}$,
$S_{8}=0.755_{-0.021}^{+0.019}$, respectively.

  * •

The local measurement of the Hubble constant from SH0ES999A new version of the
SH0ES measurement [9] was published during completion of this work. We do not
expect it to have any impact on our conclusions., modelled with a Gaussian
likelihood centered on $H_{0}=74.03\pm 1.42$ km/s/Mpc [10].

We start by performing two distinct sets of studies in order to illustrate the
importance of taking CMB data into account when studying the $\Lambda$DDM
scenarios, even in the long-lived regime (i.e. when the DCDM decays after CMB
decoupling):

(i) a background-only analysis against BAO101010As discussed in Section I, we
calibrate BAO data by imposing a Gaussian prior on the sound horizon at
recombination $r_{s}(z_{\rm rec})=144.7\pm 0.5\ \rm{Mpc}$, to not to spoil CMB
data [1]. and Pantheon SNIa data;

(ii) full analyses including linear perturbations, where we combine the data-
set used in (i) with CMB TT, TE, EE + lensing data, with and without including
the aforementioned informative priors on $S_{8}$ and $H_{0}$.

In the case of (i), the parameter space is fully characterized by the
following free parameters:

$\left\\{\Omega_{\rm{dcdm}}^{\rm{ini}},H_{0},\Gamma,\varepsilon\right\\},$

whereas in (ii) the whole parameter space is described by:

$\left\\{\Omega_{b}h^{2},\text{ln}\left(10^{10}A_{s}\right),n_{s},\tau_{\rm
reio},\Omega_{\rm{dcdm}}^{\rm{ini}},H_{0},\Gamma,\varepsilon\right\\}.$

For both (i) and (ii) we adopt logarithmic priors on $\varepsilon$ and
$\Gamma$ 111111For comparisons with previous works, an useful conversion is
the following
$\mathrm{Log}_{10}(\Gamma/\mathrm{Gyrs}^{-1})\simeq\mathrm{Log}_{10}(\Gamma/\mathrm{km}\
\mathrm{s}^{-1}\mathrm{Mpc}^{-1})-2.991$. , namely,

$-4\leq{\mathrm{Log}_{10}}~{}\varepsilon\leq{\mathrm{Log}_{10}}(0.5)$,
---
$-4\leq{\mathrm{Log}_{10}}~{}(\Gamma/{\rm Gyrs}^{-1})\leq 1$,

and a flat prior on the initial DCDM abundance:

$0\leq\Omega_{\rm{dcdm}}^{\rm{ini}}\leq 1.$

Secondly, we explore the possibility of resolving the infamous Hubble and
$S_{8}$ tensions, and the ‘$A_{\rm lens}$’ anomaly that exists within _Planck_
data. We then test the robustness of our results to various changes in the
pipeline, and in particular to trading the high-$\ell$ _Planck_ CMB data for
those from the SPT collaboration [95], which are known to be less in tension
with local $S_{8}$ measurements, as well as ACTPol data [96], which shows a
level of tension with $S_{8}$ measurements similar to Planck. Finally, we
briefly discuss the viability of the 2-body decay scenario as solution for the
Xenon1T anomaly [64].

We adopt flat priors on all other parameters, and we set two massless and one
massive active neutrino species with $m_{\nu}=0.06\ \text{eV}$, following
_Planck_ ’s conventions [1]. We assume our MCMC chains to be converged when
the Gelman-Rubin criterion $R-1<0.02$ [97]. To extract the best-fit
parameters, we make use of the Minuit algorithm [98] through the iMinuit
python package121212https://iminuit.readthedocs.io/. In App. D we report all
individual $\chi^{2}$’s per each of the analyses performed.

###  IV.2  General constraints: background vs. linear perturbations

Figure 6:  1 and 2D marginalized posterior distributions for the cosmological
parameters relevant for our analysis. Hereafter, unless otherwise stated, the
green shaded bands refer to the joint $S_{8}$ measurement from
KiDS-1000+BOSS+2dFLens, while the gray bands stand for the $H_{0}$ measurement
by the SH0ES collaboration. Note that the BAO+SNIa analysis is based on
background evolution only, whereas the BAO+SNIa+CMB analysis also includes
linear perturbations (see the main text for further details).

In Fig. 6 we show the 2$\sigma$ posterior distributions of the parameters
characterizing the $\Lambda$DDM model, ${\mathrm{Log}_{10}}(\varepsilon)$ and
${\mathrm{Log}_{10}}(\Gamma/{\rm Gyrs}^{-1})$. Red contours refer to the
background-only analysis (i), while blue contours refer to the full analysis
(ii), at the perturbation level, when CMB data are also added. From BAO+SNIa
and _Planck_ data, the DDM is not detected. We confirm our expectation from
Sec. II.1 and III that there is a negative correlation between
${\mathrm{Log}_{10}}(\varepsilon)$ and ${\mathrm{Log}_{10}}(\Gamma/{\rm
Gyrs}^{-1})$, apparent within all data sets: as $\varepsilon$ decreases, the
95% C.L. upper limit on $\Gamma$ relaxes. In practice, we find that for decay
rates $\Gamma\sim 10^{-3}-10^{-1}\text{Gyrs}^{-1}$, our $2\sigma$ exclusion
curve is roughly described by $\varepsilon\simeq 1.6\times
10^{-4}(\Gamma/\text{Gyrs}^{-1})^{-1.1}$. For large $\Gamma$ and small
$\varepsilon$, the factor ‘$1-e^{-\Gamma t_{0}}$’ in Eq. 18 reaches 1 faster
than ‘$\sqrt{1-2\varepsilon}$’ (assuming $t_{0}\simeq 13.8$ Gyrs), explaining
why our constraint on the DR energy fraction becomes flat even for very large
decay rates. In particular, for small decay rates
${\mathrm{Log}_{10}}(\Gamma/{\rm Gyrs}^{-1})\lesssim-3$ and very massive
daughters ${\mathrm{Log}_{10}}(\varepsilon)\lesssim-2.7$ the DCDM behaves like
CDM, leading to departures indistinguishable from $\Lambda$CDM. Of utmost
importance, we find that constraints on the $\Lambda$DDM free parameters
become much stronger when CMB data are included, increasing by more than one
order of magnitude over the whole parameter space, contrarily to the naive
expectations that CMB data do not weigh-in on late-time decays.

Interestingly, we see that the background-only analysis predicts a value of
$H_{0}$ slightly higher than the one inferred assumed $\Lambda$CDM, though it
can be noticed that the standard $\Lambda$CDM value for $H_{0}$ is still
perfectly compatible, due to larger uncertainties with respect to the full
analysis, which anyhow pulls $H_{0}$ back to its standard value. This suggests
that this kind of models is not suitable for relieving the Hubble tension,
contrarily to earlier claims [28], and in agreement with Refs. [29, 30].
However, in contrast to Ref. [30], we observe a significant decrease in the
$S_{8}$ contours for $\varepsilon\sim 0.01$ and $\Gamma^{-1}\sim 10^{2}$ Gyrs.
We attribute this disagreement to the fact that Ref. [30] does not include a
treatment of WDM perturbations, which are responsible for the suppression in
the matter power spectrum. Hence, the 2-body decay presented here could
potentially reconcile the inferred value of $S_{8}$ with its direct
measurements from LSS observations, as pointed out in Ref. [73]. We present an
explicit comparison of our constraints with those from Ref. [30] in App. E.

Finally, in order to compare our constraints on $\Gamma$ with previous
literature, we have carried out a MCMC analysis including BAO + SNIa +
_Planck_ data, but fixing131313This is not equivalent to directly reading the
constraints on $\Gamma$ at $\varepsilon=0.5$ from the $\Gamma$ vs.
$\varepsilon$ contours, since the 95% C.L. derived from a $\chi^{2}$
distribution with different degrees of freedom correspond to different
$\Delta\chi^{2}$. $\varepsilon=0.5$ (i.e. in the limit in which the daughter
particle behaves as dark radiation). We find a $2\sigma$ upper limit on the
DCDM decay rate of $\text{log}_{10}(\Gamma/[\text{Gyr}^{-1}])\lesssim-2.67$,
corresponding to $\Gamma^{-1}\gtrsim 468\ \rm{Gyrs}$. Our constraints on DM
decays to DR are three times tighter than those found in works using older
Planck data [65, 27] but also $\sim 40\%$ tighter than Ref. [67] due to the
use of a logarithmic prior on $\Gamma$ as opposed to linear.

Parameter | BAO+SNIa | +_Planck_
---|---|---
$100\ \omega_{b}$ | – | $2.243(2.244)_{-0.013}^{+0.014}$
$\Omega_{\rm dcdm}^{\rm ini}$ | $0.2529(0.2532)_{-0.01}^{+0.0098}$ | $0.2606(0.2619)_{-0.0054}^{+0.0051}$
$H_{0}/[{\rm km/s/Mpc}]$ | $68.48(68.44)_{-0.91}^{+0.88}$ | $67.71(67.71)_{-0.43}^{+0.42}$
$\text{ln}(10^{10}A_{s})$ | – | $3.051(3.052)_{-0.015}^{+0.014}$
$n_{s}$ | – | $0.9674(0.9672)_{-0.0038}^{+0.0038}$
$\tau_{\rm reio}$ | – | $0.0576(0.0582)_{-0.0079}^{+0.0069}$
$\text{log}_{10}(\Gamma/[{\rm Gyr}^{-1}])$ | unconstrained(0.09) | unconstrained($-3.86$)
$\text{log}_{10}(\varepsilon)$ | unconstrained(-2.89) | $-2.69(-2.97)_{-1.3}^{+0.32}$
$\Omega_{\rm m}$ | $0.299(0.2992)_{-0.011}^{+0.013}$ | $0.3102(0.3109)_{-0.0058}^{+0.0056}$
$S_{8}$ | – | $0.821(0.828)_{-0.011}^{+0.017}$
$\chi^{2}_{\rm min}$ | 1036.6 | 2053.4

Table 1: The mean (best-fit) $\pm 1\sigma$ errors of the cosmological
parameters from our $\Lambda$DDM analyses against BAO + SNIa and BAO + SNIa +
_Planck_. For each data-set, we also report the best-fit $\chi^{2}$.

###  IV.3 Implications for cosmological tensions and Xenon1T

####  IV.3.1 The $H_{0}$ tension

In order to test the implications of the 2-body decay for cosmological
tensions, we conduct a run that include the local measurement of $H_{0}$ from
SH0ES [10], CMB, BAO and SNIa data. For the sake of brevity we do not report
the results of the runs here. We find that the shape of the posterior
probabilities is almost unchanged, except for a tiny shift in $H_{0}$ to a
higher value, $H_{0}=68.21\pm 0.4\ \rm{km}/\rm{s}/\rm{Mpc}$. We thus confirm
the inability of this model to resolve the Hubble tension. This was expected
since it had already been shown through a model independent reconstruction of
the late-time dynamics of the dark sector that any late-time solution that
does not modify the sound horizon at recombination is expected to fail when
combining BAO with SNIa data (e.g. [72, 99, 17, 100]). Although this is not of
material importance in the context of the present work, let us note that
recent Refs. [101, 102] emphasized that a more correct way of combining
Pantheon and SH0ES is through a prior on the intrinsic magnitude of SN1a. Yet,
this does not affect our conclusions, as it was explicitly shown in Ref.
[103].

####  IV.3.2 The $S_{8}$ tension and the role of priors

Figure 7: 2D posterior distribution of a subset of parameters from analyses
performed by imposing three different priors on $S_{8}$, based on the three
sets of weak lensing measurements considered, whose mean values are shown as a
solid, dashed and dotted line, for KiDS+BOSS+2dFLens, DES and KiDS+Viking+DES
respectively. Note that all the three analyses also include our baseline data-
set, i.e. Planck+BAO+SNIa.

Parameter | w/ KiDS+BOSS+2dFLens | w/ DES | w/ KiDS+Viking+DES
---|---|---|---
$100\ \omega_{b}$ | $2.246(2.241)\pm 0.013$ | $2.245(2.241)_{-0.013}^{+0.014}$ | $2.246(2.244)_{-0.013}^{+0.014}$
$\Omega_{\rm dcdm}^{\rm ini}$ | $0.2581(0.2606)_{-0.0054}^{+0.005}$ | $0.2585(0.2610)_{-0.0051}^{+0.005}$ | $0.2585(0.2603)_{-0.0056}^{+0.0053}$
$H_{0}/[{\rm km/s/Mpc}]$ | $67.92(67.70)_{-0.42}^{+0.43}$ | $67.88(67.66)_{-0.42}^{+0.41}$ | $67.89(67.73)_{-0.47}^{+0.42}$
$\text{ln}(10^{10}A_{s})$ | $3.048(3.052)_{-0.016}^{+0.014}$ | $3.048(3.054)_{-0.015}^{+0.014}$ | $3.050(3.047)_{-0.015}^{+0.015}$
$n_{s}$ | $0.9682(0.9673)\pm 0.0037$ | $0.9681(0.9671)_{-0.0038}^{+0.0037}$ | $0.9681(0.9671)_{-0.004}^{+0.0036}$
$\tau_{\rm reio}$ | $0.0570(0.0582)_{-0.0077}^{+0.0071}$ | $0.0568(0.0590)_{-0.0078}^{+0.0069}$ | $0.0578(0.0554)_{-0.0077}^{+0.007}$
$\text{log}_{10}(\Gamma/[{\rm Gyr}^{-1}])$ | $-1.89(-1.74)_{-1.5}^{+0.82}$ | $-2.15(-1.97)_{-1.8}^{+0.41}$ | $-1.62(-1.08)_{-1}^{+0.81}$
$\text{log}_{10}(\varepsilon)$ | $-2.28(-2.16)_{-0.78}^{+0.8}$ | $<-2.14(-2.10)$ | $-2.23(-2.52)_{-0.59}^{+0.48}$
$\Omega_{\rm m}$ | $0.3071(0.3099)_{-0.0057}^{+0.0053}$ | $0.3078(0.3107)_{-0.0054}^{+0.0055}$ | $0.3073(0.3085)_{-0.0058}^{+0.0058}$
$S_{8}$ | $0.795(0.767)_{-0.016}^{+0.024}$ | $0.809(0.784)_{-0.01}^{+0.021}$ | $0.778(0.763)_{-0.023}^{+0.025}$
$\chi^{2}_{\rm min}$ | 2055.0 | 2054.8 | 2055.9
$\footnotesize{\chi^{2}_{\rm min}(\Lambda\text{DDM})-\chi^{2}_{\rm min}(\Lambda\text{CDM})}$ | -5.7 | -2.3 | -8.6

Table 2: The mean (best-fit) $\pm 1\sigma$ errors of the cosmological
parameters from our BAO + SNIa + _Planck_ \+ $S_{8}$, being the latter from
KiDS+BOSS+2dFLens, DES-Y1 and KiDS-1000+Viking+DES-Y1 analyses. For each data-
set we also report its best-fit $\chi^{2}$, and the $\Delta\chi^{2}$ with
respect to the analogous $\Lambda$CDM best-fit model.

As illustrated in Sec. IV.2, in a CMB+BAO+SNIa analysis of the $\Lambda$DDM
scenario, not including weak lensing data, the $S_{8}$ tension is only
marginally alleviated and DDM is not detected. On the other hand, in Ref. [73]
it was shown that, by including a $S_{8}$-prior as measured by
KiDS+BOSS+2dFLens, the $\Lambda$DDM model can fully resolve the tension if DM
decays with a lifetime of $\Gamma^{-1}\simeq 55\ \text{Gyrs}$ and converts a
fraction $\varepsilon\simeq 0.7\ \%$ of its rest mass energy into kinetic
energy for the massless component. The reason for this is explained in detail
in the latter reference: it is due to the fact that the presence of WDM at
late-times reduce the amplitude of the matter power spectrum on small scales,
thereby decreasing $\sigma_{8}$, without changing significantly the value of
$\Omega_{m}$. This results in a smaller $S_{8}$ without running into
disagreement with measurements from $\Omega_{m}$ from BAO and SNIa data.
Moreover, the fact that the decays occurs at late-times leaves the CMB power
spectra largely unchanged.

In this former study, a split-normal likelihood on $S_{8}$ as determined by
KiDS-1000+BOSS+2dFLens [94], i.e. $S_{8}=0.766^{+0.02}_{-0.014}$ was used. A
more accurate approach would have been to include the full galaxy shear and
clustering power spectra. Making use of the full likelihood would however
require the ability to compute the matter power spectrum on non-linear scales
in a $\Lambda$DDM universe, a task that is beyond the scope of this paper. Let
us note that it has been established in various cases that the reconstructed
$S_{8}$ value only mildly vary from one model to another. In particular, the
KiDS collaboration has established that the reconstructed value of $S_{8}$ is
insensitive to the neutrino mass [3] – a model that has physical effects very
similar to the $\Lambda$DDM model. This provides confidence in making use of a
prior on $S_{8}$ derived in the $\Lambda$CDM context. Nevertheless, to
highlight the impact of a different $S_{8}$ measurement, we replace the
$S_{8}$ prior from KiDS-1000+BOSS+2dFLens by the one determined in the
combined analysis KiDS+Viking+DES-Y1141414Note that this combined analysis
includes a photo-metric redshift correction applied to DES result, slightly
lowering the $S_{8}$ value compared to what is advocated by the DES
collaboration [2]. [4], $S_{8}=0.755_{-0.021}^{+0.019}$ and with the DES-Y1
data only [2], $S_{8}=0.773_{-0.02}^{+0.026}$. The results of these analyses
are reported in Tab. 2, and shown in Fig. 7, where we directly compare the
results of the three different runs151515The mean $\pm 1\sigma$ errors, and
the best-fit of the cosmological parameters in the KiDS-1000+BOSS+2dFLens
analysis are reported in the first column of Tab. 3. From the reconstructed
parameters

$\displaystyle\text{log}_{10}(\Gamma/[{\rm Gyr}^{-1}])$ $\displaystyle=$
$\displaystyle-1.89_{-1.5}^{+0.82}~{}~{}~{}~{}~{}~{}{\rm KiDS1000}$
$\displaystyle\text{log}_{10}(\varepsilon)$ $\displaystyle=$
$\displaystyle-2.28_{-0.78}^{+0.8}~{}~{}~{}~{}~{}~{}{\rm+BOSS\\!+2dFLens}$
$\displaystyle S_{8}$ $\displaystyle=$ $\displaystyle
0.795_{-0.016}^{+0.024},$ $\displaystyle\text{log}_{10}(\Gamma/[{\rm
Gyr}^{-1}])$ $\displaystyle=$
$\displaystyle-1.62_{-1}^{+0.81}~{}~{}~{}~{}~{}~{}{\rm KiDS\\!+Viking}$
$\displaystyle\text{log}_{10}(\varepsilon)$ $\displaystyle=$
$\displaystyle-2.23_{-0.59}^{+0.48}~{}~{}~{}~{}~{}~{}{\rm+DES-Y1}$
$\displaystyle S_{8}$ $\displaystyle=$ $\displaystyle
0.778_{-0.023}^{+0.025},$ $\displaystyle\text{log}_{10}(\Gamma/[{\rm
Gyr}^{-1}])$ $\displaystyle=$
$\displaystyle-2.15_{-1.8}^{+0.41}~{}~{}~{}~{}~{}~{}~{}~{}{\rm DES-only}$
$\displaystyle\text{log}_{10}(\varepsilon)$ $\displaystyle=$
$\displaystyle-2.52_{-1.5}^{+0.38}$ $\displaystyle S_{8}$ $\displaystyle=$
$\displaystyle 0.809_{-0.01}^{+0.021},$

one can see that the level of preference is higher for the KiDS+Viking+DESY1
case while it is lower in the DES-only case. In the former case, we
nevertheless note a slight degradation of $\chi^{2}_{\rm min}\sim+1$ compared
to the ‘baseline’ analysis of Ref. [73], while in the latter case the
$\chi^{2}_{\rm min}$ does not sensibly change. This explicitly demonstrates
that the statistical significance of the DDM “detection” is strongly driven by
the level of tension of the $S_{8}$ value used in the analysis (see the values
of $\Delta\chi^{2}$ with respect to $\Lambda$CDM reported in Tab. 2). If the
$S_{8}$ tension increases in the future, the preference for $\Lambda$DDM over
$\Lambda$CDM would likely increase. On the other hand, if the $S_{8}$ tension
disappears, cosmological data would not favour the $\Lambda$DDM scenario
compared to the standard $\Lambda$CDM scenario.

####  IV.3.3 The $A_{\rm lens}$ anomaly

Parameter | w/ $A_{\rm lens}$
---|---
$100\ \omega_{b}$ | $2.262(2.260)_{-0.015}^{+0.016}$
$\Omega_{\rm dcdm}^{\rm ini}$ | $0.2506(0.2526)_{-0.0059}^{+0.0049}$
$H_{0}/[{\rm km/s/Mpc}]$ | $68.56(68.38)_{-0.45}^{+0.5}$
$\text{ln}(10^{10}A_{s})$ | $3.025(3.032)_{-0.018}^{+0.02}$
$n_{s}$ | $0.9725(0.9718)_{-0.0039}^{+0.0042}$
$\tau_{\rm reio}$ | $0.0474(0.0506)_{-0.008}^{+0.0098}$
$\text{log}_{10}(\Gamma/[{\rm Gyr}^{-1}])$ | $-2.10(-1.49)_{-1.9}^{+0.39}$
$\text{log}_{10}(\varepsilon)$ | unconstrained (-2.47)
$\Omega_{\rm m}$ | $0.2991(0.3009)_{-0.0066}^{+0.0053}$
$S_{8}$ | $0.784(0.768)_{-0.014}^{+0.018}$
$A_{\rm lens}^{\rm TTTEEE}$ | $1.208(1.192)_{-0.064}^{+0.066}$
$A_{\rm lens}^{\phi\phi}$ | $1.086(1.072)_{-0.041}^{+0.035}$
$\chi^{2}_{\rm min}$ | 2043.2

Table 3: The mean (best-fit) $\pm 1\sigma$ errors of the cosmological
parameters from our BAO + SNIa + _Planck_ \+ $S_{8}$ (from KiDS+BOSS+2dFLens)
analysis performed by marginalizing over the amplitude of the lensing
potential $A_{\rm lens}$. We also report the best-fit $\chi^{2}$.

Another well studied ‘curiosity’ in the recent literature consists in the
anomalous amount of lensing estimated from the smoothing of the acoustic peaks
at high-$\ell$’s within _Planck_ data, as quantified by the ‘$A_{\rm lens}$’
parameter [104, 105, 1, 106]. However, this anomalous ‘lensing’ is not
supported by the lensing power spectrum reconstruction, such that it is now
commonly admitted that this tension (oscillating between the $2-3\sigma$
statistical level) cannot originate from a true lensing effect. On the other
hand, it has been understood that this anomaly can be easily resolved if the
universe is closed [107, 108, 109], in certain modified gravity theories
[110], or in early-universe scenarios inducing a pattern of primordial
oscillatory features [111, 112, 113, 114]. Here instead, we wish to check
whether this anomaly could impact constraints on the $\Lambda$DDM model, and
conversely if the $\Lambda$DDM model could help explaining the existence of
such anomalies. In fact, it has already been noted that this anomaly could be
related to the $S_{8}$ tension: indeed, once including $A_{\rm lens}$ as an
extra free-parameter in the analysis, it has been shown that the reconstructed
cosmology has a smaller $A_{s}$ and $\omega_{\rm cdm}$ (as well as a higher
$H_{0}$), showing no $S_{8}$ tension, but a remnant $\sim 3.5\sigma$ Hubble
tension [115, 116, 117]. To do so, we follow the approach of the SPTpol
collaboration and implement two new parameters in CLASS that allows to
(roughly) marginalize over the lensing information in _Planck_. The parameter
$A_{\rm lens}^{\rm TTTEEE}$ re-scales the amplitude of the lensing power
spectrum entering in the high-$\ell$ part of the CMB TT,TE,EE spectra, while
the parameter $A_{\rm lens}^{\phi\phi}$ re-scales the amplitude of the lensing
power spectrum reconstruction. We present the result of a MCMC analysis
including _Planck_ high-$\ell$ TT,TE,EE+lensing+BAO+SNIa+$S_{8}$ data – with
the two extra lensing parameters – in Tab. 3 and in Fig. 8. One can see that
the $S_{8}$ parameter reconstructed once marginalizing over the ‘$A_{\rm
lens}$’ anomaly is lower by $\sim 0.5\sigma$ than in the baseline analysis.
However, the preference for $\Lambda$DDM decreases, with
$\text{log}_{10}(\varepsilon)$ now unconstrained. Note also that the $A_{\rm
lens}^{\rm TTTEEE}$ is still more than 2$\sigma$ away than the fiducial value
$A_{\rm lens}^{\rm TTTEEE}=1$. We can therefore conclude that the $\Lambda$DDM
model cannot explain this anomaly and that the preference for $\Lambda$DDM
would likely disappear if the $S_{8}$ tension turns out to be explained by a
systematic in _Planck_ data leading to the anomalous value of the $A_{\rm
lens}$ parameters.

Figure 8: 2D posterior distribution of a subset of parameters from our BAO +
SNIa + _Planck_ \+ $S_{8}$ (from KiDS+BOSS+2dFLens) analysis, with and without
including the extra $A_{\rm lens}^{\rm TTTEEE}$ and $A_{\rm lens}^{\phi\phi}$
to marginalize over CMB lensing information.

####  IV.3.4 Connections with the recent Xenon1T anomaly

Following our work [73], let us further study the implications of the
$\Lambda$DDM model for the excess of events in the electronic recoils recently
reported by the Xenon1T Collaboration [64]. It has been pointed out that this
excess could potentially be explained by the elastic interactions of electrons
with a fast DM component of mass $m\gtrsim 0.1\ \rm{MeV}$ and velocities
$0.05\lesssim v/c\lesssim 0.13$ [61, 43, 62] (see also [63] for an alternative
decaying scenario). Interestingly, the WDM daughter species in the
$\Lambda$DDM scenario could in principle play the role of such a fast
component, since our results indicate that the 1-$\sigma$ range for the speed
of the daughter particle extends up to $v/c\simeq\varepsilon\simeq 0.05$. In
order to test this hypothesis further, we perform another MCMC analysis
including _Planck_ high-$\ell$ TT,TE,EE+lensing+BAO+SNIa+$S_{8}$, with the DR
energy fraction now fixed to $\varepsilon=0.05$. This serves as a proxy for
taking into account Xenon1T measurement (alternatively, one could enforce
$\varepsilon>0.05$). The results are summarized in Fig. 9 and Tab. 4. We find
best-fit values $\mathrm{Log}_{10}(\Gamma/\mathrm{Gyrs}^{-1})\simeq-2.4$ and
$S_{8}\simeq 0.784$, at the cost of a mild degradation in the fit to _Planck_
data ($\Delta\chi^{2}\simeq+1.7$), indicating that the 2-body decaying
scenario has indeed the potential to provide a common resolution to the
$S_{8}$ and Xenon1T anomalies. We leave the construction of a realistic model
and study of the subsequent cosmological implications beyond the effect the
decay to another study (see [43, 62] for examples).

Figure 9: 2D posterior distribution of a subset of parameters in our BAO + SNIa + _Planck_ \+ $S_{8}$ (from KiDS+BOSS+2dFLens) analysis when fixing the DR energy fraction to a value motivated by the recent Xenon1T anomaly ($\varepsilon=0.05$). Parameter |  $\Lambda$DDM ($\varepsilon=0.05$)
---|---
$100\ \omega_{b}$ | $2.244(2.245)_{-0.015}^{+0.014}$
$\Omega_{\rm{dcdm}}^{\rm{ini}}$ | $0.2597(0.2595)_{-0.0067}^{+0.0054}$
$H_{0}/[{\rm km/s/Mpc}]$ | $67.76(67.78)_{-0.45}^{+0.54}$
$\text{ln}(10^{10}A_{s})$ | $3.050(3.047)_{-0.016}^{+0.015}$
$n_{s}$ | $0.9674(0.9674)_{-0.0039}^{+0.0042}$
$\tau_{\rm reio}$ | $0.0573(0.0559)_{-0.0079}^{+0.0073}$
$\text{log}_{10}(\Gamma/[{\rm Gyr}^{-1}])$ | $-2.72(-2.44)_{-0.21}^{+0.61}$
$\Omega_{\rm m}$ | $0.3093(0.3090)_{-0.007}^{+0.0057}$
$S_{8}$ | $0.794(0.786)_{-0.015}^{+0.021}$
$\chi^{2}_{\rm min}$ | 2057.6
Table 4: The mean (best-fit) $\pm 1\sigma$ errors of the cosmological
parameters from our _Planck_ high-$\ell$ TT,TE,EE+lensing+BAO+SNIa+$S_{8}$
analysis, when fixing the DR energy fraction to $\varepsilon=0.05$. We also
report the best-fit $\chi^{2}$.

##  V Detecting DDM in the CMB: impact of current and future data

In this section, we confront the DDM model to additional CMB data from current
ground based surveys and perform forecast for future surveys. In addition to
Planck, we consider first, the high-$\ell$ CMB EE and TE ($50\leq\ell\leq
8000$) [95] measurements and the reconstructed gravitational lensing potential
($100\leq\ell\leq 8000$) [118] from the 500deg SPTpol survey [119]. We then
include the high-$\ell$ CMB TT, EE and TE ($350\leq\ell\leq 4125$) data from
the DR4 of the ACTPol survey [96, 120]. Finally, we demonstrate that an
experiment like CMB-S4 can unambiguously detect the DDM model.

####  V.0.1 Confronting $\Lambda$DDM to SPTpol data

Figure 10: 2D posterior distribution of a subset of parameters in the joint
BAO + SNIa + _Planck_ \+ SPTpol analysis, with and without imposing a prior on
S8 from KiDS+BOSS+2dFLens, compared to the $\Lambda$CDM scenario. Figure 11:
2D posterior distribution of a subset of parameters in our from our BAO + SNIa
+ _Planck_ \+ $S_{8}$ (from KiDS+BOSS+2dFLens) analysis, compared to the joint
BAO + SNIa + _Planck_ \+ SPTpol + $S_{8}$ analysis.

Data | BAO + SNIa + _Planck_ TT + SPTpol
---|---
Parameter | w/o $S_{8}$ | w/ $S_{8}$
$100\ \omega_{b}$ | $2.239(2.2378)_{-0.018}^{+0.019}$ | $2.241(2.247)_{-0.019}^{+0.016}$
$\Omega_{\rm dcdm}^{\rm ini}$ | $0.2544(0.2557)_{-0.0061}^{+0.0057}$ | $0.2514(0.2532)_{-0.0055}^{+0.0054}$
$H_{0}/[{\rm km/s/Mpc}]$ | $68.15(68.03)_{-0.48}^{+0.51}$ | $68.39(68.25)_{-0.46}^{+0.47}$
$\text{ln}(10^{10}A_{s})$ | $3.031(3.026)_{-0.017}^{+0.016}$ | $3.026(3.018)_{-0.014}^{+0.017}$
$n_{s}$ | $0.9701(0.9695)_{-0.0042}^{+0.0041}$ | $0.9712(0.9707)_{-0.004}^{+0.0037}$
$\tau_{\rm reio}$ | $0.0509(0.0481)_{-0.0081}^{+0.0074}$ | $0.0494(0.0457)_{-0.0073}^{+0.0082}$
$\text{log}_{10}(\Gamma/[{\rm Gyr}^{-1}])$ | $-2.38(-1.73)_{-1.6}^{+0.38}$ | $-2.25(-1.35)_{-1.7}^{+0.42}$
$\text{log}_{10}(\varepsilon)$ | unconstrained(-2.73) | unconstrained(-2.57)
$\Omega_{\rm m}$ | $0.3033(0.3046)_{-0.0067}^{+0.0064}$ | $0.2999(0.3013)_{-0.0062}^{+0.0059}$
$S_{8}$ | $0.799(0.798)_{-0.015}^{+0.022}$ | $0.787(0.767)_{-0.013}^{+0.016}$
$\chi^{2}_{\rm min}$ | 1816.3 | 1816.8

Table 5: The mean (best-fit) $\pm 1\sigma$ errors of the cosmological
parameters from our BAO + SNIa + _Planck_ TT + SPTpol analysis, with and
without imposing a split-normal likelihood on $S_{8}$ (from
KiDS+BOSS+2dFLens). For each data-set, we also report the best-fit $\chi^{2}$.

It is interesting to test the robustness of the DDM “detection” to a change of
CMB data sets, especially given the impact of marginalizing over the ‘$A_{\rm
lens}$’ anomaly as discussed in previous section. We thus confront the
$\Lambda$DDM scenario under study against a set of CMB data constituted by
low-$\ell$ temperature and polarization as well as high-$\ell$ temperature
data from _Planck_ , in combination with high-$\ell$ polarization data from
SPTpol (see section IV.1 for further details and references). It has been
shown indeed that such a joint analysis predicts an amount of CMB lensing
consistent with the $\Lambda$CDM expectation [119], i.e. no ‘$A_{\rm lens}$’
anomaly, and no $S_{8}$ tension. This is manifest in Fig. 10, where we compare
predictions from the $\Lambda$CDM and the $\Lambda$DDM models, the latter both
with and without including information on $S_{8}$ from KIDS1000+BOSS+2dFLens.
As one can easily see, both cosmological models predict a $S_{8}$ value in
excellent agreement with the KIDS1000+BOSS+2dFLens measurement, displayed as a
green horizontal band. Our results regarding the $\Lambda$DDM model are also
reported in Tab. 5.

In Fig. 11, instead, we report a comparison between our baseline $\Lambda$DDM
analysis and the _Planck_ +SPT one. First and foremost, SPTpol appears in very
good agreement with the $\Lambda$DDM model resolution of the $S_{8}$ tension
required by Planck. However, the $\Lambda$DDM parameters are largely
unconstrained in that case and no deviations from $\Lambda$CDM are visible,
which further establishes that if the $S_{8}$ tension turns out to be
explained by a systematic in _Planck_ high-$\ell$ polarization data, the
preference for $\Lambda$DDM is likely to vanish.

####  V.0.2 Confronting $\Lambda$DDM to ACTPol data

Second, we confront the $\Lambda$DDM model to the combination of Planck and
ACTPol data, to test whether more accurate measurements at high-$\ell$ can
further constrain the model. Within $\Lambda$CDM, it has been found that
ACTPol data (when combined with WMAP) also favor relatively high $S_{8}$, in
$2.1\sigma$ tension with KiDS1000+BOSS+2dFLens [96]. To limit double counting
of information, we follow the procedure of the ACT collaboration and truncate
multipoles $\ell<1800$ in the ACT TT data. The results of this analysis are
presented in Fig. 12. and Table 6. Interestingly, the $\Lambda$DDM parameter
$\varepsilon$ is more precisely measured with the inclusion of the ACTPol
data, while the mean value is barely affected. Compared with $\Lambda$CDM, the
$\Delta\chi^{2}$ in favor of $\Lambda$DDM is now $-6.7$, and the level of
tension between Planck+ACT+BAO+SN1a and $S_{8}$ from KIDS1000+BOSS+2dFLens is
$1.3\sigma$. We conclude that ACT data are in very good agreement with the
$\Lambda$DDM model, slightly increasing its preference over $\Lambda$CDM.

Figure 12: 2D posterior distribution of a subset of parameters in the case BAO
+ SNIa + _Planck_ \+ $S_{8}$ (from KiDS+BOSS+2dFLens), with and without the
inclusion of ACT data.

Data | BAO + SNIa + _Planck_ \+ ACTPol
---|---
Parameter | w/o $S_{8}$ | w/ $S_{8}$
$100\ \omega_{b}$ | $2.245(2.242)\pm 0.013$ | $2.243(2.243)\pm 0.012$
$\Omega_{\rm dcdm}^{\rm ini}$ | $0.26(0.261)\pm 0.005$ | $0.2583(0.2591)_{-0.0051}^{+0.0048}$
$H_{0}/[{\rm km/s/Mpc}]$ | $67.77(67.73)_{-0.42}^{+0.33}$ | $67.92(67.84)\pm 0.40$
$\text{ln}(10^{10}A_{s})$ | $3.058(3.054)\pm 0.015$ | $3.056(3.061)\pm 0.015$
$n_{s}$ | $0.9703(0.9697)_{-0.004}^{+0.0034}$ | $0.9708(0.9706)\pm 0.0036$
$\tau_{\rm reio}$ | $0.0571(0.0549)_{-0.0084}^{+0.0069}$ | $0.0567(0.0592)_{-0.0075}^{+0.007}$
$\text{log}_{10}(\Gamma/[{\rm Gyr}^{-1}])$ | $0.92(-2.92)_{-1.90}^{+0.55}$ | $-1.56(-1.19)_{-1.50}^{+1.10}$
$\text{log}_{10}(\varepsilon)$ | unconstrained (-3.88) | $-2.34(-2.54)\pm 0.65$
$\Omega_{\rm m}$ | $0.3094(0.3108)_{-0.0054}^{+0.0051}$ | $0.3069(0.3068)_{-0.0056}^{+0.0050}$
$S_{8}$ | $0.822(0.829)_{-0.011}^{+0.018}$ | $0.792(0.772)_{-0.019}^{+0.025}$
$\chi^{2}_{\rm min}$ | 2294.54 | 2296.20
$\Delta\chi^{2}(\Lambda\text{CDM})$ | -0.12 | -6.7

Table 6: The mean (best-fit) $\pm 1\sigma$ errors of the cosmological
parameters from our BAO + SNIa + _Planck_ \+ ACTPol analysis, with and without
imposing a split-normal likelihood on $S_{8}$ (from KiDS+BOSS+2dFLens). For
each data-set, we also report the best-fit $\chi^{2}$ and the
$\Delta\chi^{2}(\Lambda\text{CDM})\equiv\chi^{2}_{\rm
min}(\Lambda\text{DDM})-\chi^{2}_{\rm min}(\Lambda\text{CDM})$.

####  V.0.3 Towards detecting the $\Lambda$DDM model with CMB-S4

As we have extensively discussed, current CMB data are not sensitive enough to
detect DDM, so that the preference for non-standard values for $\varepsilon$
and $\Gamma$ is fully driven by the inclusion of the $S_{8}$ measurement from
weak lensing data in the analysis. To further stress this aspect we have
performed additional analyses fitting a set of mock CMB data generated
starting from our reference best-fit $\Lambda$DDM model (i.e., that from the
$Planck$+BAO+SNIa+KiDS+BOSS+2dFLens run, reported in the first column of Tab.
2). The resulting contour plots are shown in Fig. 13, where we compare the
constraints that $Planck$ would obtain if the “true” cosmological model
actually contained DDM, with those that a future generation CMB survey
(CMB-S4) would get. Concretely, this task was pursued by using the ‘Planck-
fake-realistic’ and ‘CMB-S4’ likelihoods available in MontePython-v3. The
former allowed us to generate synthetic Planck data, whereas the latter
includes multipoles $\ell$ from 30 to 3000, assuming a sky coverage of 40$\%$,
uncorrelated Gaussian error on each $a_{\ell m}$’s, uncorrelated temperature
and polarization noise, and perfect foreground cleaning up to $\ell_{\rm
max}$161616To overcome the lack of low-$\ell$ data in the CMB-S4 analysis, we
have imposed a Gaussian prior on the optical depth to reionization, centred on
its best-fit value from our reference $Planck$+BAO+SNIa+KiDS+BOSS+2dFLens
analysis, namely $\tau_{\rm reio}=0.0582\pm 0.008$.. All details about the
likelihood can be found in Tab. 1 of Ref. [121].

As expected, from Fig. 13 it is manifest that $Planck$ alone could not detect
DDM even if its signature was truly imprinted in CMB data. Note indeed that
the red contours barely features an overlap at 2$\sigma$ in $S_{8}$, and only
upper limits on ${\rm Log}_{10}\Gamma$, ${\rm Log}_{10}\varepsilon$. This, as
we explained earlier, is a consequence of the degeneracy that exists within
$\Lambda$CDM and leads to a bias in the Bayesian analysis. Therefore, the
information that matters in quantifying the success of the resolution is
rather contained in the $\chi^{2}$ values: just like in the analysis of real
data, we find that, when testing the $\Lambda$CDM model against the mock data
that contains the $\Lambda$DDM signal, the $\chi^{2}$ in the $\Lambda$CDM
model is identical to that of the $\Lambda$DDM model. This clearly shows that
Planck cannot disentangle between $\Lambda$CDM and $\Lambda$DDM, while the
$S_{8}$ measurements favors $\Lambda$DDM (in terms of $\chi^{2}$). This is not
the case for CMB-S4, which contours are over-plotted in blue: if the real
Universe contains DDM, CMB-S4 would unequivocally detect its signature,
finding $\varepsilon\neq 0$ at $\gtrsim$ 2 $\sigma$ level. The
$\Delta\chi^{2}$ in favour of the $\Lambda$DDM model from (mock) CMB-S4 data
alone is in fact $\simeq+8$.

Figure 13: 2D posterior distribution of a subset of parameters reconstructed
from a fit to simulated $Planck$ and CMB-S4 data. The fiducial model has ${\rm
log}_{10}(\epsilon)=-2.16$ and ${\rm log}_{10}(\Gamma/[{\rm
Gyr^{-1}}])=-1.74$, as denoted by the gray dotted lines.

##  VI Conclusions

In this paper, we have performed a comprehensive cosmological study of the CDM
2-body decay scenario dubbed ‘$\Lambda$DDM’, whereby decays are characterized
both by the decay rate $\Gamma$ and energy fraction converted to radiation
$\varepsilon$, including a fully consistent treatment of the linear
perturbations of the WDM daughter component.

To that end, we have made use of a new approximation scheme, introduced in
Ref. [73], that allows to accurately and quickly compute the dynamics of the
WDM linear perturbations by treating the WDM species as a viscous fluid. Close
to the best-fit values, our approximation scheme is accurate at the ${\cal
O}(0.1\%)$ level in the CMB power spectra and ${\cal O}(1\%)$ level in the
linear matter power spectrum (see App. A).

We have then discussed in detail the dynamics of linear density perturbations
of the mother and daughter particles, as well as the physical effects of the
$\Lambda$DDM model on the CMB and matter power spectra. We have shown that
accurate CMB (lensing in particular) and matter power spectrum measurements
have the potential to detect both the decay rate $\Gamma$ and energy fraction
converted to radiation $\varepsilon$.

In a second part, we have performed a set of MCMC analyses of the $\Lambda$DDM
model against a suite of up-to-date low- and high-redshift data-sets. We have
compared the constraints obtained from BAO and SNIa data, thereby solely based
on background effects, to those obtained from the full Planck data-set, that
requires instead an accurate description of the WDM linear perturbations. We
find that _Planck_ CMB data constrain the $\Lambda$DDM model $\sim 1$ order of
magnitude better than current BAO+SNIa data. However, we also show that
despite these constraints, the $\Lambda$DDM model provides a promising
possibility to resolve the $S_{8}$ tension, as detailed in a previous paper
[73].

We have then tested the robustness of the $\Lambda$DDM resolution to the
$S_{8}$ tension to a number of change in the analysis. We show that the mild
preference for the $\Lambda$DDM model over $\Lambda$CDM is tied to the $S_{8}$
value chosen in the analysis. Concretely, the $S_{8}$ value from the
KiDS+Viking+DES analysis, which has a higher level of tension with the
_Planck_ $\Lambda$CDM prediction than the baseline KiDS+BOSS+2dFLens value,
leads to a stronger preference in favor of the $\Lambda$DDM model. However,
the DES-only result, which is in reasonable agreement with _Planck_ , leads to
a weaker preference of the $\Lambda$DDM model. Similarly, once marginalizing
over the lensing information in _Planck_ through the $A_{\rm lens}$ parameter
(we used two extra parameters in practice describing the normalization of the
lensing power spectrum and the normalization of the lensing smoothing effect
in the high-$\ell$ TT,TE,EE power spectra), or when trading the _Planck_
high-$\ell$ TE,EE power spectra for the SPTpol ones, the preference for the
$\Lambda$DDM model decreases. This is because in these two cases the inferred
$\Lambda$CDM model has a smaller $S_{8}$ value, showing less of a tension with
the weak lensing surveys. This indicates that if the $S_{8}$ tension ends up
coming from an unknown systematic within weak lensing surveys or within
_Planck_ data, the preference for the $\Lambda$DDM model would likely
disappear. On the other hand, when combining _Planck_ with ACTPol the mild
preference for $\Lambda$DDM increases, and the remaining ‘tension’ with
$S_{8}$ is now only $\sim 1.3\sigma$.

We have also tested the intriguing possibility that the recent Xenon1T excess
is due to the $\Lambda$DDM model. To that end, we have performed a additional
MCMC analysis fixing $\varepsilon=0.05$ as required by Xenon1T. We find that
it is easy to resolve the $S_{8}$ tension in that case, pointing to a DCDM
lifetime of $\text{log}_{10}(\Gamma/[{\rm Gyr}^{-1}])=-2.72_{-0.21}^{+0.61}$.
Interestingly, this comes at the cost of a very minor degradation in Planck
fit ($\Delta\chi^{2}\simeq+1.7$), indicating that Planck, BAO and SNIa data
are in good agreement with this model.

Finally, by performing an analysis where we artificially introduce a DDM
signal in a set of mock CMB data, we explicitly demonstrate that while current
CMB data alone are not sensitive enough to distinguish between standard CDM
and DDM, next-generation CMB experiments (CMB-S4) can unambiguously detect its
signature.

It will be very interesting to go beyond the linear aspects discussed in this
work and study the non-linear evolution of density perturbations, in order to
be able to make use of the full power of the KiDS, BOSS and DES likelihoods.
This could for instance be done with N-body simulations, as in Refs. [122,
123, 124], or via the Effective Theory of LSS [125, 126]. This will be even
more crucial with upcoming surveys such as Euclid [127], LSST [128], and DESI
[129], which will measure the matter power spectrum and the growth factor with
great accuracy up to $z\sim 2$. It might also be possible to test the
$\Lambda$DDM model with current Lyman-$\alpha$ forest flux power spectrum data
[123, 78, 79, 24, 80, 130]. We plan to study non-linear aspects and the
discovery potential of these surveys in an upcoming work.

###### Acknowledgements.

The authors are thankful to Julien Lavalle and Rodrigo Calderón for many
useful comments and discussions. The authors acknowledge the use of
computational resources from the CNRS/IN2P3 Computing Centre (CC-IN2P3) in
Lyon, the IN2P3/CNRS and the Dark Energy computing Center funded by the OCEVU
Labex (ANR-11-LABX-0060) and the Excellence Initiative of Aix-Marseille
University (A*MIDEX) of the “Investissements d’Avenir” programme. This project
has received support from the European Union’s Horizon 2020 research and
innovation program under the Marie Skodowska-Curie grant agreement No
860881-HIDDeN.

##  Appendix A Numerical implementation and accuracy of the WDM fluid
approximation

Figure 14:  Residuals of the lensed CMB TT power spectrum (upper) and EE power
spectrum (lower) in the WDM fluid approximation, with respect to the full
hierarchy calculation, for a grid of values covering most of the parameter
space: $\Gamma/H_{0}=0.1,1,10$ and $\varepsilon=0.5,0.1,0.01,0.001$. The gray
shaded regions indicate _Planck_ 1$\sigma$ errors, while the pink shaded areas
indicate cosmic variance up to $\ell=3000$. Figure 15:  Residuals of the
linear matter power spectrum (at $z=0$), with respect to the baseline
$\Lambda$CDM model, for the same grid of parameter values considered in Fig.
14, both from the full hierarchy calculation (solid lines) and the WDM fluid
approximation(dashed lines).

In order to solve the cosmic evolution equations derived in section II, we
modified the publicly available numerical Boltzmann solver CLASS [75, 76]. We
now briefly illustrate our implementation of the $\Lambda$DDM model.

First let us notice that, when solving the background equations for all
cosmological species, the DE abundance is iteratively derived through the
budget equation, $\Omega_{\Lambda}=1-\sum_{i}\Omega_{i}$, where the sum
includes the current abundance of all other components, which are not known a
priori.

We thus applied a shooting method for the aforementioned parameter, i.e. we
guess an initial $\Omega_{\Lambda}$, we solve the system of background
equations to obtain $\sum_{i}\Omega_{i}$, and re-compute $\Omega_{\Lambda}$.
The procedure is iterated until convergence is achieved. The WDM density is
computed by solving Eq. (15) in 2800 momentum bins, approximately as many as
the time-steps used to describe its background evolution.

At the linear perturbation level, we truncate the hierarchy of the PSD
multipole equations for both the daughter particles at a $\ell_{\rm max}=17$.
We set the initial conditions for the WDM species following the same procedure
of [69].

On conformal times $\tau<\tau_{q}$, we set all $\Delta f_{{\rm{wdm}},\ell}=0$,
since no daughter particle with comoving momentum $q$ could have been
produced. On the crossing time $\tau=\tau_{q}$, one should be more careful, as
the terms with $\bar{f}_{\rm{wdm}}$ in Eqs. (20) and (22) contain a Dirac
delta and, when integrated, a Heaviside function. Thus, the corresponding
initial conditions for $\Delta f_{{\rm{wdm}},0}(\tau_{q})$ and $\Delta
f_{{\rm{wdm}},2}(\tau_{q})$ are not-vanishing. We set them according to the
analytical formulas (A.5) and (A.7) from [69].

Finally, on times $\tau>\tau_{q}$, we treat the WDM component as a massive
neutrino species, and we solve the corresponding hierarchy of equations in 300
momentum-bins. This number of bins is chosen simply because it gives a good
compromise between speed and accuracy: it is large enough to accurately
describe the super-Hubble and Hubble-crossing scales, where the dynamics is
relatively simple, and small enough to not become too computationally
expensive 171717 Note that the number of bins used at both the background and
perturbation level is much larger than the one typically used in standard
CLASS analysis for massive active neutrinos, given that the time-dependence of
the background PSD of the WDM requires a finer momentum resolution. Regarding
the momentum spacing, we have considered a logarithmic Simpson quadrature
instead of the Gauss-Laguerre quadrature typically used in standard CLASS
analysis. . On sub-Hubble scales, when $k\tau$ is larger than a threshold
value $(k\tau)_{\rm fluid}$, we switch-on the fluid approximation described in
section II.3. The WDM dynamics is now described by Eqs. (33)-(36) and (38). We
have chosen $(k\tau)_{\rm fluid}=25$ to provide the speed yet accurate enough
for the purposes of the current analyses.

In Figs. 14 and 15 we explicitly compare the novel approximation scheme with
the results of the “exact” computation for the WDM species. For the latter, we
solve the full Boltzmann hierarchy using $10^{4}$ momentum-bins and $\ell_{\rm
max}=17$. In Fig. 14 we show the residuals of the lensed CMB TT and EE power
spectra in the WDM fluid approximation, with respect to the full computation,
for a grid of parameter values given by $\Gamma/H_{0}=0.1,1,10$ and
$\varepsilon=0.5,0.1,0.01,0.001$. These values span most of the parameter
space in the $\Lambda$DDM framework, and for none of them the residuals exceed
the _Planck_ $1\sigma$ uncertainties, which are indicated by the gray shaded
regions, nor the error bars for a cosmic-variance-limited experiment (close to
CMB-S4 errors), indicated by the pink shaded regions.

The predictions for the linear matter power spectrum $P(k)$ are less accurate
than for the anisotropy spectra, because the former is more sensitive to the
dynamics of the daughter particles. Close to the best-fit parameter values,
and in general close to $\Lambda$CDM, the residual errors between the full and
the fluid calculations in $P(k)$ are $\mathcal{O}(1\%)$, but they can become
higher far away from the best-fit. In particular, we have verified that inside
the parameter region delimited by $\rm{Log}_{10}\varepsilon\in[-2.3,-0.7]$ and
$\mathbf{\rm{Log}_{10}(\Gamma/\rm{Gyrs}^{-1})\in[-1.3,1]}$, the residual
errors are typically larger than $10\%$, so the fluid approximation should be
used with caution in this region. However, this portion of the parameter space
is deeply inside the $2\sigma$ exclusion region, as one can check by looking
at Fig. 6. In addition, given that current data are mostly sensitive to
integrals over $P(k)$ (e.g. $S_{8}$, CMB lensing), we are mainly interested in
getting accurate predictions for the departures from $\Lambda$CDM (rather than
the exact shape of the matter spectrum itself). To illustrate that, we have
computed the residuals of the linear matter power spectrum (at $z=0$) with
respect to our baseline $\Lambda$CDM model, for both the fluid and the full
hierarchy calculations. The results are shown in Fig. 15, where we have
spanned the same parameter values as in Fig. 14. We can see that, for all the
$\Lambda$DDM models, the shape of the power suppression (that is, the depth
and the cut-off scale) is excellently well-captured by our fluid prescription.
Furthermore, we verified that the residuals in the structure growth parameter
$S_{8}\equiv\sigma_{8}(\Omega_{m}/0.3)^{0.5}$ are always smaller than the
$\sim 1.8\%$ relative error of the $S_{8}$ measurement from [94]. We thus
conclude that the new WDM viscous fluid approximation is accurate enough for
our analyses.

##  Appendix B Semi-analytic understanding of the WDM sound speed

Here we obtain a formal equation that dictates the evolution of the WDM sound
speed in the synchronous gauge, $c^{2}_{\rm s}$. The first natural step is to
write a dynamical equation for the normalized pressure perturbation,
$\Pi_{\rm{wdm}}=\delta P_{\rm{wdm}}/\bar{\rho}_{\rm{wdm}}$. This can be
achieved by multiplying Eq. (20) by $4\pi
q^{2}\frac{q^{2}}{3\mathcal{E}_{\rm{wdm}}}a^{-4}$, integrating over $q$ and
then using Eq. (14). By doing so, two higher velocity-weight integrals appear,
namely

$\delta\mathcal{P}_{\rm{wdm}}\equiv\frac{4\pi}{3a^{4}}\int_{0}^{\infty}dq\frac{q^{6}}{\mathcal{E}_{\rm{wdm}}^{3}}\Delta
f_{{\rm{wdm}},0},$ (40)

and

$(\bar{\rho}_{\rm{wdm}}+\bar{P}_{\rm{wdm}})\Theta_{\rm{wdm}}\equiv\frac{4\pi
k}{a^{4}}\int_{0}^{\infty}dq\frac{q^{5}}{\mathcal{E}_{\rm{wdm}}^{2}}\Delta
f_{{\rm{wdm}},1}.$ (41)

The variable $\Theta$ was already discussed in Ref. [76] in the context of
massive neutrinos. In the relativistic limit $\delta\mathcal{P}_{\rm{wdm}}$
and $\Theta_{\rm{wdm}}$ become equal to the standard variables $\delta
P_{\rm{wdm}}$ and $\theta_{\rm{wdm}}$ (Eqs. (7)-(8)), while in the non-
relativistic limit they are suppressed by a factor
$(q/\mathcal{E}_{\rm{wdm}})^{2}$ with respect to them. This means that one can
write $\delta\mathcal{P}_{\rm{wdm}}=\delta P_{\rm{wdm}}3\omega_{p}$ and
$\Theta_{\rm{wdm}}=\theta_{\rm{wdm}}3\omega_{\theta}$, where $\omega_{p}$ and
$\omega_{\theta}$ are arbitrary functions, going from $1/3$ in the
relativistic limit, to $0$ in the non-relativistic case. In terms of these
functions, the equation for $\Pi_{\rm{wdm}}$ reads

$\displaystyle\dot{\Pi}_{\rm{wdm}}$
$\displaystyle=-3\mathcal{H}\Pi_{\rm{wdm}}\left(\frac{2}{3}-\omega_{p}-\omega\right)$
$\displaystyle-\omega_{\theta}(1+\omega)\theta_{\rm{wdm}}-\frac{\dot{h}\omega}{6}\left[5-\frac{\mathfrak{p}_{\rm{wdm}}}{p_{\rm{wdm}}}\right]$
$\displaystyle+a\Gamma\frac{\bar{\rho}_{\rm{dcdm}}}{\bar{\rho}_{\rm{wdm}}}\left[\frac{\varepsilon^{2}}{(1-\varepsilon)}\frac{\delta_{\rm{dcdm}}}{3}-(1-\varepsilon)\Pi_{\rm{wdm}}\right]$
(42)

Figure 16:  Left $-$ The sound speed of the WDM species in the synchronous
gauge, as a function of conformal time and wavenumber. Right $-$ The ratio
between the perturbed densities of the DCDM and the WDM components, as a
function of conformal time and wavenumber. The black solid and dashed lines
indicate the horizon and free-streaming crossing scales, respectively. We have
set $\varepsilon=0.007$ and $\Gamma^{-1}=55\ \rm{Gyrs}$.

One can convert this into an equation for the sound speed by noting that
$\Pi_{\rm{wdm}}=c_{\rm s}^{2}\delta_{\rm{wdm}}$, and using the continuity
equation (33). The final result reads:

$\displaystyle\frac{\partial c^{2}_{\rm s}}{\partial\tau}$
$\displaystyle=-3\mathcal{H}c^{2}_{\rm
s}\left(\frac{2}{3}-\omega_{p}-c^{2}_{\rm
s}\right)-(1+\omega)\frac{\theta_{\rm{wdm}}}{\delta_{\rm{wdm}}}(\omega_{\theta}-c^{2}_{\rm
s})$
$\displaystyle-\frac{\dot{h}}{2\delta_{\rm{wdm}}}\left[\frac{\omega}{3}\left(5-\frac{\mathfrak{p}_{\rm{wdm}}}{\bar{P}_{\rm{wdm}}}\right)-c^{2}_{\rm
s}(1+\omega)\right]$
$\displaystyle+a\Gamma\frac{\bar{\rho}_{\rm{dcdm}}}{\bar{\rho}_{\rm{wdm}}}\frac{\delta_{\rm{dcdm}}}{\delta_{\rm{wdm}}}\left[\frac{\varepsilon^{2}}{3(1-\varepsilon)}-(1-\varepsilon)c^{2}_{\rm
s}\right].$ (43)

We remark that the previous equation is highly non-linear in the perturbed
quantities, so it can easily give rise to numerical instabilities. In
addition, there is no closed expression for computing $\omega_{p}$ and
$\omega_{\theta}$. If these functions were scale independent, one possible
approximation would be to trade them for some background functions, such as
$w$ or $c_{g}^{2}$. However, calculations using the full hierarchy show that
$\omega_{p}$ and $\omega_{\theta}$ exhibit a $k$-dependence similar to that of
$c_{s}^{2}$. For these reasons, we do not implement Eq. (43) in our code.

Nonetheless, by making some simplifying assumptions, Eq. (43) allows to
qualitatively understand why there is a particular $k$-dependence of
$c^{2}_{\rm s}$ in the decaying scenario, that is not present in the case of
massive neutrinos. Let us consider the non-relativistic limit of Eq. (43),
since data favors in general very small DR energy fractions, $\varepsilon\ll
1$. This also implies that $w_{p},c^{2}_{\rm s},w\ll 1$. Let us further assume
that $\omega_{\theta}$ and $c^{2}_{\rm s}$ behave similarly, so that the
difference $\omega_{\theta}-c^{2}_{\rm s}$ can be neglected. Finally, let us
also restrict to sub-Hubble scales, for which the term
$\dot{h}/\delta_{\rm{wdm}}$ is very small and can be also neglected. In this
case, Eq. (43) reduces to

$\displaystyle\frac{\partial c^{2}_{\rm s}}{\partial\tau}$
$\displaystyle=-2\mathcal{H}c^{2}_{\rm s}\
-a\Gamma\frac{\bar{\rho}_{\rm{dcdm}}}{\bar{\rho}_{\rm{wdm}}}\frac{\delta_{\rm{dcdm}}}{\delta_{\rm{wdm}}}c^{2}_{\rm
s}.$ (44)

In absence of the decay term, we see that the sound speed dilutes as
$c^{2}_{\rm s}\propto a^{-2}$, which is a well-known result for massive
neutrinos. This dilution can be compensated by the presence of the decay term,
leading to a $c^{2}_{\rm s}\sim\rm{cte}$, as long as the ratio
$\delta_{\rm{dcdm}}/\delta_{\rm{wdm}}$ doesn’t change. In practice, for scales
and times such that $k<k_{\rm fs}(\tau)$, we have
$\delta_{\rm{dcdm}}/\delta_{\rm{wdm}}=1$. In this regime, the sound speed
$c^{2}_{\rm s}$ is well approximated by the adiabatic sound speed $c_{g}^{2}$.
However, when $k>k_{\rm fs}(\tau)$, $\delta_{\rm{wdm}}$ oscillates and starts
to become suppressed with respect to $\delta_{\rm{dcdm}}$, which leads to
oscillatory features and a small enhancement in the evolution of $c^{2}_{\rm
s}$. This is visible in Fig. 16, where we have plotted $c^{2}_{\rm s}$ and
$\delta_{\rm{dcdm}}/\delta_{\rm{wdm}}$ in the $k-\tau$ plane using the full
Boltzmann hierarchy, for the best-fit parameters from the combined analysis of
Ref. [73] (i.e., setting $\varepsilon=0.007$ and $\Gamma^{-1}=55\ \rm{Gyrs}$).

One can see that this $k$-dependent effect appears only because of the
coupling term in Eq. (44), which is not present for massive neutrinos. This
also justifies why the sound speed $c^{2}_{\rm s}$ is well approximated by a
background function such as $c_{g}^{2}$ in the case of massive neutrinos. In
the case of the WDM daughter species, the approximation $c^{2}_{\rm s}\simeq
c_{g}^{2}$ will only work when $k<k_{\rm fs}(\tau)$. This motivates the use of
the fitting formula introduced in Eq. (38), that accounts for the small
enhancement at scales smaller than the free-streaming scale. While this simple
fitting formula is not able to capture the oscillatory features described
previously, it leads to results that are accurate enough for all the
observables analysed in this work.

##  Appendix C Consistency between the use of the _Planck_ TTTEEE ‘lite’ and
‘full’ likelihoods

Figure 17: 2D posterior distribution for a subset of parameters obtained by using the _Planck_ ‘full’ likelihood, compared to those from the _Planck_ ‘lite’ analysis. Parameter |  $Planck$ ‘full’ |  $Planck$ ‘lite’
---|---|---
$100\ \omega_{b}$ | $2.243_{-0.015}^{+0.012}$ | $2.246\pm 0.013$
$\Omega_{\rm dcdm}^{\rm ini}$ | $0.2585_{-0.0048}^{+0.0054}$ | $0.2581_{-0.0054}^{+0.005}$
$H_{0}/[{\rm km/s/Mpc}]$ | $67.88_{-0.44}^{+0.36}$ | $67.92_{-0.42}^{+0.43}$
$\text{ln}(10^{10}A_{s})$ | $3.046_{-0.015}^{+0.015}$ | $3.048_{-0.016}^{+0.014}$
$n_{s}$ | $0.9667_{-0.0039}^{+0.0036}$ | $0.9682\pm 0.0037$
$\tau_{\rm reio}$ | $0.0562_{-0.0072}^{+0.0073}$ | $0.0570_{-0.0077}^{+0.0071}$
$\text{log}_{10}(\Gamma/[{\rm Gyr}^{-1}])$ | $-1.84_{-1.4}^{+1}$ | $-1.89_{-1.5}^{+0.82}$
$\text{log}_{10}(\varepsilon)$ | $-2.26_{-0.82}^{+0.84}$ | $-2.28_{-0.78}^{+0.8}$
$\Omega_{\rm m}$ | $0.3075_{-0.0049}^{+0.0057}$ | $0.3071_{-0.0057}^{+0.0053}$
$S_{8}$ | $0.794_{-0.018}^{+0.024}$ | $0.795_{-0.016}^{+0.024}$
Table 7: The mean $\pm 1\sigma$ errors of the cosmological parameters from our
BAO + SNIa + _Planck_ \+ $S_{8}$ analysis, obtained by using both the _Planck_
‘full’ and ‘lite’ likelihoods.

Given that for all our analyses we made use of the ‘lite’ version of the
_Planck_ likelihood, in Fig. 17 and Tab. 7 we show that our conclusions are
not affected by such a choice, by comparing the results from our main run with
the predictions from a ‘full’ _Planck_ analysis.

##  Appendix D Best-fit $\chi^{2}$ per experiment

In Tab. 8 we report all $\chi^{2}_{\rm min}$’s obtained with the Minuit
algorithm [98] through the iMinuit python package for the various model and
data-set combinations considered in this work.

| BAO+SNIa | +_Planck_ | +_Planck_ w/ $S_{8}$ | +SPTpol | +ACTPol
---|---|---|---|---|---
|  | w/o $S_{8}$ |  | w/ $A_{\rm lens}$ | $\varepsilon=0.05$ | w/o $S_{8}$ | w/ $S_{8}$ | w/o $S_{8}$ | w/ $S_{8}$
Pantheon SNIa | 1026.9 | 1028.2 | 1027.5 | 1026.8 | 1029.2 | 1028.0 | 1027.1 | 1027.0 | 1026.9 | 1027.0 | 1026.84
BAO+FS BOSS DR12 | – | 6.63 | 7.06 | 6.93 | 7.11 | 6.14 | 6.59 | 5.94 | 6.17 | 6.78 | 6.66
BAO BOSS DR12 | 3.52 | – | – | – | – | – | – | – | – | – | –
BAO BOSS low$-z$ | 1.84 | 1.20 | 1.22 | 1.18 | 1.24 | 1.83 | 1.28 | 1.58 | 1.79 | 1.17 | 1.34
BAO eBOSS DR14 | 4.29 | 4.94 | 4.91 | 4.94 | 4.83 | 4.52 | 4.88 | 4.68 | 4.53 | 4.97 | 4.77
_Planck_ high$-\ell$ TT,TE,EE ‘lite’ | – | 584.8 | 585.9 | 585.3 | 586.9 | 577.7 | 587.5 | – | – | 589.216 | 590.545
_Planck_ high$-\ell$ TT ‘lite’ | – | – | – | – | – | – | – | 207.8 | 207.8 | – | –
_Planck_ low$-\ell$ EE | – | 396.9 | 396.9 | 397.2 | 396.3 | 395.7 | 396.3 | 395.8 | 396.1 | 396.2 | 397.15
_Planck_ low$-\ell$ TT | – | 23.1 | 23.1 | 23.2 | 23.0 | 22.1 | 23.0 | 22.3 | 22.1 | 22.6 | 22.46
_Planck_ lensing | – | 8.78 | 9.12 | 8.88 | 9.47 | 8.53 | 9.83 | – | – | 8.8 | 8.94
SPTpol high$-\ell$ TE,EE | – | – | – | – | – | – | – | 145.9 | 145.6 | – | –
SPTpol lensing | – | – | – | – | – | – | – | 5.43 | 5.93 | – | –
ACTPol | – | – | – | – | – | – | – | – | – | 238.235 | 237.359
KiDS+BOSS+2dFLens | – | – | 0.0003 | – | – | 0.0097 | 0.98 | – | 0.0015 | – |
DES | – | – | – | 0.19 | – | – | – | – | – | – | –
KiDS+Viking+DES | – | – | – | – | 0.20 | – | – | – | – | – | –
total $\chi^{2}$ | 1036.6 | 2053.4 | 2055.0 | 2054.8 | 2055.9 | 2043.2 | 2057.6 | 1816.3 | 1816.8 | 2294.8 | 2296.2

Table 8: Best-fit $\chi^{2}$ per experiment (and total) for all the
$\Lambda$DDM analyses performed in this work.

##  Appendix E  Comparison with the Planck constraints from Ref. [30]

Here we carry out an explicit comparison of our constraints with those of Ref.
[30], which performed an analysis of the $\Lambda$DDM model against Planck
data, neglecting the perturbations of the warm daughter particles. As shown in
Fig. 18, we find that the constraints on the $\Lambda$DDM models are up to
(roughly) one order of magnitude stronger when our improved treatment is
considered.

Figure 18: Comparison between the $2\sigma$ exclusion bounds (solid lines)
from the _Planck_ analysis of Ref. [30] and our _Planck_ +BAO+SNIa analysis.
In each case, the dashed line indicates a fit that roughly describes the
$2\sigma$ limit in the range $\Gamma\sim 10^{-3}-10^{-1}\ \text{Gyrs}^{-1}$.

##  Appendix F Results with a linear prior on $\Gamma$ and $\varepsilon$

Figure 19: Reconstructed 2D posteriors of a BAO + SNIa + _Planck_ \+ $S_{8}$
(from KiDS+BOSS+2dFLens) analysis, with linear priors and sampling either with
the original prior range (upper panel) or within a restricted prior range
(lower panel).

In our baseline analysis we have made use of log-prior on $\varepsilon$ and
$\Gamma$, to ease comparison with earlier works [28, 30] who adopted the same
choice. Here we present results using linear priors on the DDM parameters. Let
us however stress that the use of a linear prior is less informative than
adopting a logarithmic one. That is because a linear prior carries a scale
(due to the large error bars used in the proposal distribution of the MCMC
sampler), so that it is hard for the sampler to detect fine structure over 4
orders of magnitude by using a linear scale, in particular at very small
values. In other words, given the difference between the scale of the upper
limit on $\Gamma/\rm{Gyrs}^{-1}$ an $\varepsilon$ ($\sim 10^{-1}$) and that of
the lower limit ($\sim 10^{-3}$)), it is very difficult to accurately
reconstruct the parameter space with a linear prior. Such a difficulty is
illustrated in Fig. 19, where we provide the results of two linear-prior
analyses: the upper panel corresponding to the original prior range, the lower
panel corresponding to a more restricted range: $\varepsilon\in[0.0001,0.015]$
and $\Gamma/\rm{Gyrs}^{-1}\in[0.0001,0.05]$. While in the latter case
$\varepsilon$ and $\Gamma$ are detected at the 2$\sigma$ level – similarly to
the log-prior results – the former case weighs in favor of larger
$\varepsilon$ values, so that one would uncorrectly deduce an upper limit
only.

## References

  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), 1807.06209 .
  * Abbott  _et al._ [2018] T. Abbott _et al._ (DES), Phys. Rev. D 98, 043526 (2018), 1708.01530 .
  * Hildebrandt  _et al._ [2020] H. Hildebrandt _et al._ , Astron. Astrophys. 633, A69 (2020), 1812.06076 .
  * Joudaki  _et al._ [2020] S. Joudaki _et al._ , Astron. Astrophys. 638, L1 (2020), 1906.09262 .
  * Verde  _et al._ [2019] L. Verde, T. Treu, and A. G. Riess, in  _Nature Astronomy 2019_ (2019) 1907.10625 .
  * Aylor  _et al._ [2019] K. Aylor, M. Joy, L. Knox, M. Millea, S. Raghunathan, and W. K. Wu, Astrophys. J. 874, 4 (2019), 1811.00537 .
  * Wong  _et al._ [2020] K. C. Wong _et al._ , Mon. Not. Roy. Astron. Soc. 498, 1420 (2020), 1907.04869 .
  * Freedman  _et al._ [2019] W. L. Freedman _et al._ , (2019), 1907.05922 .
  * Riess  _et al._ [2020] A. G. Riess, S. Casertano, W. Yuan, J. B. Bowers, L. Macri, J. C. Zinn, and D. Scolnic, (2020), 2012.08534 .
  * Riess  _et al._ [2019] A. G. Riess, S. Casertano, W. Yuan, L. M. Macri, and D. Scolnic, Astrophys. J. 876, 85 (2019), 1903.07603 .
  * Shanks  _et al._ [2019] T. Shanks, L. Hogarth, and N. Metcalfe, Mon. Not. Roy. Astron. Soc. 484, L64 (2019), 1810.02595 .
  * Riess  _et al._ [2018] A. G. Riess, S. Casertano, D. Kenworthy, D. Scolnic, and L. Macri, (2018), 1810.03526 .
  * Davis  _et al._ [2019] T. M. Davis, S. R. Hinton, C. Howlett, and J. Calcino, Mon. Not. Roy. Astron. Soc. 490, 2948 (2019), 1907.12639 .
  * Yuan  _et al._ [2019] W. Yuan, A. G. Riess, L. M. Macri, S. Casertano, and D. Scolnic, Astrophys. J. 886, 61 (2019), 1908.00993 .
  * Cerny  _et al._ [2020] W. Cerny, W. L. Freedman, B. F. Madore, F. Ashmead, T. Hoyt, E. Oakes, N. Q. H. Tran, and B. Moss, (2020), 2012.09701 .
  * Soltis  _et al._ [2020] J. Soltis, S. Casertano, and A. G. Riess, (2020), 2012.09196 .
  * Knox and Millea [2019] L. Knox and M. Millea, (2019), 1908.03663 .
  * Jedamzik  _et al._ [2020] K. Jedamzik, L. Pogosian, and G.-B. Zhao, (2020), 2010.04158 .
  * Haridasu  _et al._ [2021] B. S. Haridasu, M. Viel, and N. Vittorio, Phys. Rev. D 103, 063539 (2021), 2012.10324 [astro-ph.CO] .
  * Di Valentino  _et al._ [2021] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, (2021), 2103.01183 .
  * Vagnozzi [2021] S. Vagnozzi, (2021), 2105.10425 .
  * Kumar and Nunes [2016] S. Kumar and R. C. Nunes, Phys. Rev. D94, 123511 (2016), 1608.02454 .
  * Murgia  _et al._ [2016] R. Murgia, S. Gariazzo, and N. Fornengo, JCAP 04, 014.
  * Archidiacono  _et al._ [2019] M. Archidiacono, D. C. Hooper, R. Murgia, S. Bohr, J. Lesgourgues, and M. Viel, JCAP 10, 055, 1907.01496 .
  * Di Valentino  _et al._ [2020] E. Di Valentino _et al._ , (2020), 2008.11285 .
  * Enqvist  _et al._ [2015] K. Enqvist, S. Nadathur, T. Sekiguchi, and T. Takahashi, JCAP 09, 067, 1505.05511 .
  * Poulin  _et al._ [2016] V. Poulin, P. D. Serpico, and J. Lesgourgues, JCAP 1608 (08), 036, 1606.02073 .
  * Vattis  _et al._ [2019] K. Vattis, S. M. Koushiappas, and A. Loeb, Phys. Rev. D99, 121302 (2019), 1903.06220 .
  * Haridasu and Viel [2020] B. S. Haridasu and M. Viel, (2020), 2004.07709 .
  * Clark  _et al._ [2020] S. J. Clark, K. Vattis, and S. M. Koushiappas, (2020), 2006.03678 .
  * Hambye [2011] T. Hambye, PoS IDM2010, 098 (2011), 1012.4587 .
  * Berezinsky  _et al._ [1991] V. Berezinsky, A. Masiero, and J. Valle, Physics Letters B 266, 382 (1991).
  * Kim and Kim [2002] H. B. Kim and J. E. Kim, Physics Letters B 527, 18 (2002).
  * Covi  _et al._ [1999] L. Covi, J. E. Kim, and L. Roszkowski, Phys. Rev. Lett. 82, 4180 (1999), hep-ph/9905212 .
  * Feng  _et al._ [2003] J. L. Feng, A. Rajaraman, and F. Takayama, Phys. Rev. Lett. 91, 011302 (2003), hep-ph/0302215 .
  * Feng  _et al._ [2003] J. L. Feng, A. Rajaraman, and F. Takayama, Phys. Rev. D 68, 063504 (2003), hep-ph/0306024 [hep-ph] .
  * Allahverdi  _et al._ [2015] R. Allahverdi, B. Dutta, F. S. Queiroz, L. E. Strigari, and M.-Y. Wang, Phys. Rev. D 91, 055033 (2015), 1412.4391 [hep-ph] .
  * Abazajian  _et al._ [2012] K. N. Abazajian _et al._ , (2012), 1204.5379 .
  * Drewes  _et al._ [2017] M. Drewes _et al._ , JCAP 01, 025, 1602.04816 .
  * Chen  _et al._ [2009] C.-R. Chen, F. Takahashi, and T. T. Yanagida, Phys. Lett. B 671, 71 (2009), 0809.0792 .
  * Choi  _et al._ [2020a] G. Choi, M. Suzuki, and T. T. Yanagida, Phys. Rev. D 101, 075031 (2020a), 2002.00036 .
  * Choi  _et al._ [2020b] G. Choi, M. Suzuki, and T. T. Yanagida, Phys. Rev. D 102, 035022 (2020b), 2004.07863 .
  * Choi  _et al._ [2020c] G. Choi, M. Suzuki, and T. T. Yanagida, Phys. Lett. B 811, 135976 (2020c), 2006.12348 .
  * Choi and Yanagida [2021] G. Choi and T. T. Yanagida, (2021), 2104.02958 .
  * Slatyer and Wu [2017] T. R. Slatyer and C.-L. Wu, Phys. Rev. D 95, 023010 (2017), 1610.06933 .
  * Poulin  _et al._ [2017] V. Poulin, J. Lesgourgues, and P. D. Serpico, JCAP 03, 043, 1610.10051 .
  * Cirelli  _et al._ [2012] M. Cirelli, E. Moulin, P. Panci, P. D. Serpico, and A. Viana, Phys. Rev. D 86, 083506 (2012), 1205.5283 .
  * Essig  _et al._ [2013] R. Essig, E. Kuflik, S. D. McDermott, T. Volansky, and K. M. Zurek, JHEP 11, 193, 1309.4091 .
  * Jin  _et al._ [2013] H.-B. Jin, Y.-L. Wu, and Y.-F. Zhou, JCAP 11, 026, 1304.1997 .
  * Giesen  _et al._ [2015] G. Giesen, M. Boudaud, Y. Génolini, V. Poulin, M. Cirelli, P. Salati, and P. D. Serpico, JCAP 09, 023, 1504.04276 .
  * Berezinsky and Valle [1993] V. Berezinsky and J. W. F. Valle, Phys. Lett. B 318, 360 (1993), hep-ph/9309214 .
  * Lattanzi [2008] M. Lattanzi, AIP Conf. Proc. 966, 163 (2008), arXiv:0802.3155 [astro-ph] .
  * Lin  _et al._ [2001] W. B. Lin, D. H. Huang, X. Zhang, and R. Brandenberger, Phys. Rev. Lett. 86, 954 (2001), astro-ph/0009003 [astro-ph] .
  * Sánchez-Salcedo [2003] F. J. Sánchez-Salcedo, apjl 591, L107 (2003), astro-ph/0305496 .
  * Cembranos  _et al._ [2005] J. A. R. Cembranos, J. L. Feng, A. Rajaraman, and F. Takayama, Phys. Rev. Lett. 95, 181301 (2005), hep-ph/0507150 .
  * Kaplinghat [2005] M. Kaplinghat, Phys. Rev. D 72, 063510 (2005), astro-ph/0507300 [astro-ph] .
  * Strigari  _et al._ [2007] L. E. Strigari, M. Kaplinghat, and J. S. Bullock, Phys. Rev. D 75, 061303 (2007), astro-ph/0606281 [astro-ph] .
  * Borzumati  _et al._ [2008] F. Borzumati, T. Bringmann, and P. Ullio, Phys. Rev. D 77, 063514 (2008), hep-ph/0701007 .
  * Peter  _et al._ [2010] A. H. G. Peter, C. E. Moody, and M. Kamionkowski, Phys. Rev. D 81, 103501 (2010), 1003.0419 [astro-ph.CO] .
  * Peter and Benson [2010] A. H. G. Peter and A. J. Benson, Phys. Rev. D 82, 123521 (2010), 1009.1912 .
  * Kannike  _et al._ [2020] K. Kannike, M. Raidal, H. Veermäe, A. Strumia, and D. Teresi, (2020), 2006.10735 .
  * Xu and Zheng [2020] S. Xu and S. Zheng, (2020), 2012.10827 .
  * Buch  _et al._ [2020] J. Buch, M. A. Buen-Abad, J. Fan, and J. S. C. Leung, Phys. Rev. D 102, 083010 (2020), 2007.13750 .
  * The Xenon-1t Collaboration  _et al._ [2020] The Xenon-1t Collaboration, E. Aprile, J. Aalbers, and F. a. Agostini, (2020), 2006.09721 .
  * Audren  _et al._ [2014] B. Audren, J. Lesgourgues, G. Mangano, P. D. Serpico, and T. Tram, JCAP 12, 028, 1407.2418 .
  * Berezhiani  _et al._ [2015] Z. Berezhiani, A. D. Dolgov, and I. I. Tkachev, Phys. Rev. D92, 061303 (2015), 1505.03644 .
  * Nygaard  _et al._ [2021] A. Nygaard, T. Tram, and S. Hannestad, JCAP 05, 017, 2011.01632 .
  * Aoyama  _et al._ [2011] S. Aoyama, K. Ichiki, D. Nitta, and N. Sugiyama, JCAP 09, 025, 1106.1984 .
  * Aoyama  _et al._ [2014] S. Aoyama, T. Sekiguchi, K. Ichiki, and N. Sugiyama, JCAP 1407, 021, 1402.2972 .
  * Chudaykin  _et al._ [2016] A. Chudaykin, D. Gorbunov, and I. Tkachev, Phys. Rev. D94, 023528 (2016), 1602.08121 .
  * Bringmann  _et al._ [2018] T. Bringmann, F. Kahlhoefer, K. Schmidt-Hoberg, and P. Walia, Phys. Rev. D98, 023543 (2018), 1803.03644 .
  * Poulin  _et al._ [2018] V. Poulin, K. K. Boddy, S. Bird, and M. Kamionkowski, Phys. Rev. D97, 123504 (2018), 1803.02474 .
  * Abellan  _et al._ [2020] G. F. Abellan, R. Murgia, V. Poulin, and J. Lavalle, (2020), 2008.09615 .
  * Blackadder and Koushiappas [2014] G. Blackadder and S. M. Koushiappas, Phys. Rev. D90, 103527 (2014), 1410.0683 .
  * Blas  _et al._ [2011] D. Blas, J. Lesgourgues, and T. Tram, JCAP 1107, 034, 1104.2933 .
  * Lesgourgues and Tram [2011] J. Lesgourgues and T. Tram, JCAP 09, 032, 1104.2935 .
  * Ma and Bertschinger [1995] C.-P. Ma and E. Bertschinger, Astrophys. J. 455, 7 (1995), astro-ph/9506072 .
  * Murgia  _et al._ [2017] R. Murgia, A. Merle, M. Viel, M. Totzauer, and A. Schneider, JCAP 11, 046, 1704.07838 .
  * Murgia  _et al._ [2018] R. Murgia, V. Irsic, and M. Viel, Phys. Rev. D 98, 083540 (2018), 1806.08371 .
  * Miller  _et al._ [2019] C. Miller, A. L. Erickcek, and R. Murgia, Phys. Rev. D 100, 123520 (2019), 1908.10369 .
  * Bohr  _et al._ [2020] S. Bohr, J. Zavala, F.-Y. Cyr-Racine, M. Vogelsberger, T. Bringmann, and C. Pfrommer, (2020), 2006.01842 .
  * Lesgourgues and Pastor [2006] J. Lesgourgues and S. Pastor, Massive neutrinos and cosmology, Phys. Rept. 429, 307 (2006), astro-ph/0603494 .
  * Giusarma  _et al._ [2016] E. Giusarma, M. Gerbino, O. Mena, S. Vagnozzi, S. Ho, and K. Freese, Phys. Rev. D 94, 083522 (2016), 1605.04320 .
  * Parimbelli  _et al._ [2019] G. Parimbelli, M. Viel, and E. Sefusatti, JCAP 01, 010, 1809.06634 .
  * Manzotti [2018] A. Manzotti, Phys. Rev. D 97, 10.1103/PhysRevD.97.043527 (2018), 1710.11038 .
  * Audren  _et al._ [2013] B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, JCAP 1302, 001, 1210.7183 .
  * Brinckmann and Lesgourgues [2019] T. Brinckmann and J. Lesgourgues, Phys. Dark Univ. 24, 100260 (2019), 1804.07261 .
  * Beutler  _et al._ [2011] F. Beutler, C. Blake, M. Colless, D. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, Mon. Not. Roy. Astron. Soc. 416, 3017 (2011), 1106.3366 .
  * Ross  _et al._ [2015] A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, Mon. Not. Roy. Astron. Soc. 449, 835 (2015), 1409.3242 .
  * Alam  _et al._ [2017] S. Alam _et al._ (BOSS), Mon. Not. Roy. Astron. Soc. 470, 2617 (2017), 1607.03155 .
  * de Sainte Agathe  _et al._ [2019] V. de Sainte Agathe _et al._ , Astron. Astrophys. 629, A85 (2019), 1904.03400 .
  * Blomqvist  _et al._ [2019] M. Blomqvist _et al._ , Astron. Astrophys. 629, A86 (2019), 1904.03430 .
  * Scolnic  _et al._ [2018] D. M. Scolnic _et al._ , Astrophys. J. 859, 101 (2018), 1710.00845 .
  * Heymans  _et al._ [2020] C. Heymans _et al._ (2020) 2007.15632 [astro-ph.CO] .
  * Henning  _et al._ [2018] J. Henning _et al._ (SPT), Astrophys. J. 852, 97 (2018), 1707.09353 .
  * Aiola  _et al._ [2020] S. Aiola _et al._ (ACT), (2020), 2007.07288 .
  * Gelman and Rubin [1992] A. Gelman and D. B. Rubin, Statist. Sci. 7, 457 (1992).
  * James and Roos [1975] F. James and M. Roos, Comput. Phys. Commun. 10, 343 (1975).
  * Lemos  _et al._ [2019] P. Lemos, E. Lee, G. Efstathiou, and S. Gratton, Mon. Not. Roy. Astron. Soc. 483, 4803 (2019), 1806.06781 .
  * Benevento  _et al._ [2020] G. Benevento, W. Hu, and M. Raveri, Phys. Rev. D 101, 103517 (2020), 2002.11707 .
  * Camarena and Marra [2021] D. Camarena and V. Marra, (2021), 2101.08641 .
  * Efstathiou [2021] G. Efstathiou, (2021), 2103.08723 .
  * Schöneberg  _et al._ [2021] N. Schöneberg, G. Franco Abellán, A. Pérez Sánchez, S. J. Witte, V. Poulin, and J. Lesgourgues, (2021), 2107.10291 .
  * Calabrese  _et al._ [2008] E. Calabrese, A. Slosar, A. Melchiorri, G. F. Smoot, and O. Zahn, Phys. Rev. D 77, 123531 (2008), 0803.2309 .
  * Aghanim  _et al._ [2017] N. Aghanim _et al._ (Planck), Astron. Astrophys. 607, A95 (2017), 1608.02487 .
  * Efstathiou and Gratton [2019] G. Efstathiou and S. Gratton, (2019), 1910.00483 .
  * Handley [2021] W. Handley, Phys. Rev. D 103, L041301 (2021), 1908.09139 .
  * Di Valentino  _et al._ [2019] E. Di Valentino, A. Melchiorri, and J. Silk, Nature Astron. 4, 196 (2019), 1911.02087 .
  * Efstathiou and Gratton [2020] G. Efstathiou and S. Gratton, Mon. Not. Roy. Astron. Soc. 496, L91 (2020), 2002.06892 .
  * Moshafi  _et al._ [2020] H. Moshafi, S. Baghram, and N. Khosravi, (2020), 2012.14377 .
  * Chen [2010] X. Chen, Primordial Non-Gaussianities from Inflation Models, Adv. Astron. 2010, 638979 (2010), arXiv:1002.1416 [astro-ph.CO] .
  * Chluba  _et al._ [2015] J. Chluba, J. Hamann, and S. P. Patil, Int. J. Mod. Phys. D 24, 1530023 (2015), 1505.01834 .
  * Slosar  _et al._ [2019] A. Slosar, X. Chen, C. Dvorkin, D. Meerburg, B. Wallisch, D. Green, and E. Silverstein, Astro2020: Decadal Survey on Astronomy and Astrophysics, science white papers 51, 98 (2019), 1903.09883 .
  * Domènech  _et al._ [2020] G. Domènech, X. Chen, M. Kamionkowski, and A. Loeb, JCAP 10, 005, 2005.08998 .
  * Motloch and Hu [2018] P. Motloch and W. Hu, Phys. Rev. D 97, 103536 (2018), 1803.11526 .
  * Di Valentino and Bridle [2018] E. Di Valentino and S. Bridle, Symmetry 10, 585 (2018).
  * Motloch and Hu [2020] P. Motloch and W. Hu, Phys. Rev. D 101, 083515 (2020), 1912.06601 .
  * Bianchini  _et al._ [2020] F. Bianchini _et al._ (SPT), Astrophys. J. 888, 119 (2020), 1910.07157 .
  * Chudaykin  _et al._ [2020a] A. Chudaykin, D. Gorbunov, and N. Nedelko, JCAP 08, 013, 2004.13046 .
  * Choi  _et al._ [2020d] S. K. Choi _et al._ (ACT), JCAP 12, 045, 2007.07289 .
  * Brinckmann  _et al._ [2019] T. Brinckmann, D. C. Hooper, M. Archidiacono, J. Lesgourgues, and T. Sprenger, JCAP 01, 059, 1808.05955 .
  * Wang and Zentner [2012] M.-Y. Wang and A. R. Zentner, Phys. Rev. D 85, 043514 (2012), 1201.2426 .
  * Wang  _et al._ [2013] M.-Y. Wang, R. A. C. Croft, A. H. G. Peter, A. R. Zentner, and C. W. Purcell, Phys. Rev. D 88, 123515 (2013), 1309.7354 .
  * Wang  _et al._ [2014] M.-Y. Wang, A. H. G. Peter, L. E. Strigari, A. R. Zentner, B. Arant, S. Garrison-Kimmel, and M. Rocha, Mon. Not. Roy. Astron. Soc. 445, 614 (2014), 1406.0527 .
  * D’Amico  _et al._ [2021] G. D’Amico, L. Senatore, and P. Zhang, JCAP 01, 006, 2003.07956 .
  * Chudaykin  _et al._ [2020b] A. Chudaykin, M. M. Ivanov, and M. Simonović, (2020b), 2004.10607 .
  * Amendola  _et al._ [2018] L. Amendola _et al._ , Living Rev. Rel. 21, 2 (2018), 1606.00180 .
  * Alonso  _et al._ [2018] D. Alonso _et al._ (LSST Dark Energy Science), (2018), 1809.01669 .
  * Aghamousa  _et al._ [2016] A. Aghamousa _et al._ (DESI), (2016), 1611.00036 .
  * Enzi  _et al._ [2021] W. Enzi _et al._ , Mon. Not. Roy. Astron. Soc. 506, 5848 (2021), 2010.13802 .