Linear cosmological constraints on 2-body decaying dark matter scenarios and robustness of the resolution to the S_8 tension
LLinear cosmological constraints on 2-body decaying dark matter scenarios androbustness of the resolution to the S tension Guillermo F. Abell´an, ∗ Riccardo Murgia, † and Vivian Poulin ‡ Laboratoire Univers & Particules de Montpellier (LUPM), Universit´e de Montpellier (UMR-5299)Place Eug`ene Bataillon, F-34095 Montpellier Cedex 05, France (Dated: February 26, 2021)The ‘ S tension’ is a longstanding discrepancy between the cosmological and local determinationof the amplitude of matter fluctuations, parameterized as S ≡ σ (Ω m / . . , where σ is the rootmean square of matter fluctuations on a 8 h − Mpc scale, and Ω m is the total matter abundance. Itwas recently shown that dark matter (DM) decaying into a massless (dark radiation) and a massive(warm DM) species, with a lifetime Γ − (cid:39)
55 ( ε/ . . Gyrs – where ε represent the mass-energyfraction transferred to the massless component – can resolve the tension. Thanks to a new, fastand accurate approximation scheme for the warm species, we perform a comprehensive study of this2-body decaying DM scenario, discussing in details its dynamics and its impact on the CMB andlinear matter power spectra. We then confront the robustness of the resolution to the ‘ S tension’against a number of changes in the analysis: different S priors, marginalization over the lensinginformation in Planck data, trading
Planck high − (cid:96) polarization data for those from the SPTpolcollaboration, and the inclusion of the recent results from the Xenon1T collaboration. We concludethat the preference for decaying DM, while entirely driven by the local S measurements, does notsensibly degrade the fit to any of the cosmological data-sets considered, and that the model couldexplain the anomalous electron recoil excess reported by the Xenon1T collaboration. I. INTRODUCTION
In the last couple of decades, the so-called standardΛ Cold Dark Matter (ΛCDM) model of cosmology hasbeen firmly established as the most successful frameworkto interpret numerous independent experimental obser-vations up to a very high degree of accuracy. Still, thenature of its dominant components - Cold Dark Matter(CDM) and Dark Energy (DE) – is yet to be unveiled.The ΛCDM model provides indeed a remarkable fit to awide variety of early universe data, such as Cosmic Mi-crowave Background (CMB) and Big Bang Nucleosyn-thesis (BBN), as well as late universe observations suchas Large Scale Structure (LSS), Baryon Acoustic Oscil-lation (BAO), and luminosity distance to SuperNovae oftype Ia (SNIa). However, as the accuracy of the measure-ments has increased over the past few years, a number ofintriguing discrepancies between the values of some cos-mological parameter as inferred within ΛCDM and theirdirect measurements at low redshift, has emerged. At theheart of this study is the mild ( ∼ − σ ) yet longstand-ing tension between the cosmological [1] and local [2–4]determination of the amplitude of the matter fluctua-tions, typically parameterized as S ≡ σ (Ω m / . . ,where σ is the root mean square of matter fluctuationson a 8 h − Mpc scale, and Ω m is the total matter abun-dance. Another example is the so-called Hubble tension[5–9] , i.e., the significant discrepancy ( ∼ − σ C.L.)between the value of the current expansion rate of theuniverse H , directly measured using SNIa data as a cos- ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] mic distance ladder [10], and that inferred from CMBdata [1], assuming ΛCDM. Despite scrupulous efforts tominimize the systematic errors at play in the cepheidcalibration of the SNIa, the statistical significance of thedisagreement is steadily rising [9, 11–14]. However, it isstill debated whether the Hubble constant resulting froma calibration of the SNIa on the ‘tip of the red giantbranch’ shows a similar degree of tension with respectto
Planck
ΛCDM model [8, 14–16]. Alternative methodshave been proposed, but are currently not at the accu-racy level required to unambiguously weigh in the Hubbletension [5].In the absence of convincing solutions within the stan-dard cosmological model, and driven by the fact that thenature of the dark sector is unknown, throughout theyears many alternative scenarios have been proposed toexplain these discrepancies. It is now understood thatresolving the Hubble tension most likely requires newphysics in the pre-recombination era, in order to decreasethe sound horizon at recombination before recombination[17]. On the other hand, resolving the S tension requiresto decrease the amplitude of matter fluctuations on scales k ∼ . − h /Mpc, which can be easily achieved in a vari-ety of models often related to new DM properties [18–21],including the possibility that DM decays [22–26].In this paper, we reassess the phenomenology of a2-body Decaying Cold Dark Matter (DCDM) scenario,where the decay products are one massive Warm DM(WDM) particle and one (massless) DR component, in-teracting only through gravitation with the standardmodel particles. We will refer to the full model as ΛDDM.From the point of view of particle physics model build-ing, the stability over cosmological timescales is one ofthe most peculiar property of the dark matter particle,reviewed e.g. in [27]. Often, an additional symmetry a r X i v : . [ a s t r o - ph . C O ] F e b (typically a discrete Z symmetry) has to be assumed tomake the DM candidate stable. Nevertheless, DM de-cays at late-times are known signatures of many modelsin the literature such as (for instance) models with R-parity violation [28, 29], super Weakly Interacting Mas-sive particles (super WIMPs) [30–33], sterile neutrinos[34, 35] or models with an additional U(1) gauge sym-metry [36–39]. Decays to electromagnetically chargedparticles are tightly constrained by Planck data [40, 41], γ -ray [42, 43] and cosmic-ray searches [44, 45], typicallyrequiring Γ − > O (10 s), with some level of dependenceon the decay channel. Still, a purely gravitational con-straints, although weaker, is very interesting in the spiritof being ‘model-independent’, while applying to modelswith decay to a dark sector, or to (non-interacting) neu-trinos. The canonical example is perhaps that of thekeV majoron [46, 47] decaying into relativistic neutrinos.Models of CDM decays with massive daughters have alsobeen invoked as a potential solution to the observationaldiscrepancies with CDM on small (sub-galactic) scalesafter structure formation (e.g. [38, 48–55]). Even morerecently, decaying dark matter models were proposed[39, 56, 57] as a way to explain the excess of events inthe electronic recoils reported by the Xenon1T collabora-tion [58]. In the literature, most studies restricted them-selves to massless daugther particles [22, 41, 59–61], withthe benefit of greatly simplifying the cosmological anal-ysis, but limiting the true ‘model-independence’ of thebound, and therefore its robustness. Nevertheless, somestudies have attempted at including the effect of massivedaughters in a cosmological context but either neglectedcosmological perturbations of the daughter particles [24–26] or were limited by computational power to performa complete analysis against a host of cosmological data[62, 63].It has originally been suggested that DM decaying intomassless daughters could help with cosmological tensions[22, 60], but careful analysis of this scenario in light of Planck H values. This is in good agree-ment with model-independent analyses existing in theliterature in which it has been established that a combi-nation of BAO and uncalibrated SNIa data strongly con-strain any late-time modification as a resolution to theHubble tension (see, e.g, Ref. [17, 66]). A similar conclu-sion is also reached when CMB data are considered [26].Yet, these recent analyses were limited to the study ofthe effects of ΛDDM on the background evolution of theuniverse. In this paper, we perform a thorough analysis ofthe ΛDDM model in light of up-to-date low- and high-redshift data-sets, including the effects of linear pertur-bations. To that end, we introduce a new approximationscheme that allows to accurately and quickly compute thedynamics of the WDM linear perturbations by treatingthe WDM species as a viscous fluid. In a companion pa-per [67], thanks to this new fluid approximation, we haveshown that the ΛDDM, while unable to ease the Hubbletension, can fully explain the low- S measurement fromrecent weak lensing surveys. The aim of this paper is to:i) introduce the cosmological formalism of the ΛDDMmodel and a new approximation scheme that we devel-oped to accurately describe linear perturbations of thewarm daughter; ii) discuss the background and perturba-tion dynamics of the ΛDDM model and its impact on theCMB and linear matter power spectrum; iii) compare theconstraints obtained with the inclusion of perturbationsto those obtained when neglecting them (as was done inthe past literature); iv) test the robustness of the ΛDDMresolution to the S − tension to a number of changes inthe analysis (different S priors, different CMB datasets,marginalization over the lensing information in Planck ,including constraints from the Xenon1T experiment onthe model).This work is structured as follows: in Section II weintroduce the formalism and the novel approximationscheme; in Section III we illustrate the ΛDDM impacton the relevant cosmological observables; in Section IVwe discuss the results of our data analyses; finally, inSection V we draw our conclusions.
II. FORMALISM OF 2-BODY Λ DDM
Hereafter, we adopt the Boltzmann formalism byRef. [63], where the time-evolution of the Phase-SpaceDistribution (PSD) for both the mother and the daugh-ter particles were derived. However, in Section II A weexplicitly show that, at the background level, such for-malism is equivalent to the one by [68].While the (cold) parent particle can be safely describedas a perfect fluid, computing the density perturbationevolution for the daughter particles requires a more so-phisticated treatment. The central role in the game isplayed by (cid:96) max , i.e. the highest multiple to considerwhen drawing up the hierarchy of equations describingthe PSD of the daughter particles. In the massless case,the degrees of freedom associated to momentum can beremoved after the PSD multipole decomposition [23].Due to its non-trivial energy-momentum relation, thisapproach is not possible for the warm daughter particle.One has to study the full PSD evolution, which wouldbe computationally prohibitive when performing MCMCanalyses. For this reason, in Section II B we introducea novel approximation scheme based on describing theWDM component as a viscous fluid on sub-Hubble scales.This allows us to integrate out the momentum degrees offreedom and the hierarchy of equations to be truncatedat (cid:96) max = 1. We will show that the new, computationallyfaster scheme is accurate enough to be used for cosmolog-ical analyses, allowing to establish accurate and robustCMB limits on this class of models.Our framework is characterized by two additional freeparameters with respect to ΛCDM: the DCDM lifetime,Γ − , and the fraction of DCDM rest mass energy con-verted into DR, defined as follows [68]: ε = 12 (cid:18) − m m (cid:19) , (1)where 0 ≤ ε ≤ /
2. The lower limit corresponds tothe standard CDM case, so that Ω cdm = Ω dcdm + Ω wdm ,whereas ε = 1 / ε values (i.e. heavy massive daugh-ters) and small Γ values (i.e. lifetimes much longer thanthe age of the universe) induce little departures fromΛCDM.Let us decompose the PSD function of the j -th darkcomponent into a background contribution ¯ f j plus a lin-ear perturbation ∆ f j as f j ( k, q, µ, τ ) = ¯ f j ( q, τ ) + ∆ f j ( k, q, µ, τ ) , (2)where τ is the conformal time and j = { dcdm , dr , wdm } .The mean energy density and pressure ¯ ρ j are obtainedby integrating the background PSD, i.e.¯ ρ j = 1 a (cid:90) ∞ dq πq E j ¯ f j , (3)¯ P j = 13 a (cid:90) ∞ dq πq q E j ¯ f j , (4)where E j = (cid:113) m j a + q is the comoving energy of j -thspecies.The linear perturbation term ∆ f j is generally ex-panded over Legendre polynomials:∆ f j ( k, q, µ, τ ) = ∞ (cid:88) (cid:96) =0 ( − i ) (cid:96) (2 (cid:96) +1)∆ f j,(cid:96) ( k, q, τ ) P (cid:96) ( µ ) . (5)The perturbed energy density, pressure, energy flux andshear stress are thus given by: δρ j = ¯ ρ j δ j = 4 πa (cid:90) ∞ dqq E j ∆ f j, , (6) δP j = ¯ ρ j Π j = 4 π a (cid:90) ∞ dqq q E j ∆ f j, , (7)( ¯ ρ j + ¯ P j ) θ j = 4 πka (cid:90) ∞ dqq q ∆ f j, , (8)( ¯ ρ j + ¯ P j ) σ j = 8 π a (cid:90) ∞ dqq q E j ∆ f j, . (9) We use dots to indicate derivatives with respect to conformaltime
We choose to work in the synchronous gauge, co-movingwith the parent particle, that means that each ∆ f dcdm ,(cid:96) vanishes except for the term ∆ f dcdm , . In addition, sincethe parent particle is non-relativistic, the only relevantdynamical DCDM variables are ¯ ρ dcdm and δ dcdm . Wefollow the standard convention of referring to the scalarmetric perturbations in this gauge as h and η . We im-plement the equations describing the ΛDDM model in amodified version of the Boltzmann code CLASS [69, 70].Our code is available at https://github.com/PoulinV/class_majoron .Throughout the rest of this paper and unless statedotherwise, we compare ΛDDM models at fixed ω inidcdm ≡ Ω inidcdm h with ΛCDM models having the same ω cdm .All other cosmological parameters are fixed to { H =67 . / s / Mpc, ω cdm = ω inidcdm = 0 . ω b = 0 . n s = 0 . A s ) = 3 . τ reio = 0 . } , whichconstitutes our baseline ΛCDM model. These values cor-respond to the best-fit from the combined analysis (in-cluding S data from weak lensing) of Ref. [67], whichare however very similar to the ΛCDM best-fit parame-ters from Planck
A. Background dynamics
The background evolution of the mother and daughterparticles is described at the phase-space level by theirBoltzmann equations [63]:˙¯ f dcdm = − a Γ ¯ f dcdm , (10)˙¯ f dr = ˙¯ f wdm = a Γ ¯ N dcdm ( τ )4 πq δ ( q − ap max ) . (11)¯ N dcdm = (Ω inidcdm ρ c, /m dcdm ) e − Γ t is the DCDM (co-moving) number density ; p max = | (cid:126)p wdm | = | (cid:126)p dr | = m wdm ε/ √ − ε is the initial momentum of the daughterparticles.By combining Eqs. (3) and (4) with Eqs. (10) and(11) we obtain:˙¯ ρ dcdm = − H ¯ ρ dcdm − a Γ¯ ρ dcdm , (12)˙¯ ρ dr = − H ¯ ρ dr + εa Γ¯ ρ dcdm , (13)˙¯ ρ wdm = − w ) H ¯ ρ wdm + (1 − ε ) a Γ¯ ρ dcdm . (14) H ≡ ˙ a/a is the conformal Hubble parameter, and w ( τ ) ≡ ¯ P wdm / ¯ ρ wdm is the WDM Equation of State (EoS) .Eq. (14) will be useful to analytically derive the fluidequations that we present in Section II C, but for a nu-merical resolution it is much simpler to use an integral http://class-code.net/ Ω ≡ Ω inidcdm e − Γ t ; Ω inidcdm is the initial DCDM abundance, t the age of the universe. Notice that this expression does not coincide with the EoS usedin [24] formula for ¯ ρ wdm , as it was done in Ref. [68]. Thisformula can be obtained by integrating Eq. (11) firstlywith respect to τ , and then with respect to q . The firstintegration requires using the relation δ ( q − ap max ) = δ ( τ − τ q ) /q H , where τ q represents the conformal timewhen daughter particles with co-moving momentum q are born, q = a ( τ q ) p max [63]. For the second integral,changing the integration variable from q to a q = a ( τ q )leads to¯ ρ wdm ( a ) = Ca (cid:90) a da q e − Γ t q H q (cid:113) ε a q + (1 − ε ) a , (15)where C ≡ ρ c, Ω inidcdm Γ, H q ≡ H ( a q ), t q ≡ t ( a q ). Notethe equivalence between Eq. (15) and the analogous ex-pression derived in Ref. [68] with a different formalism.Concerning the massless DR species, we simply take thelimit ε → / H de-pends on ¯ ρ wdm and ¯ ρ dr through the Friedmann equation.For a flat universe, this equation reads H ( a ) = 8 πGa (cid:88) i ¯ ρ i ( a ) , (16)where (cid:88) i ¯ ρ i ( a ) = ¯ ρ dcdm ( a ) + ¯ ρ dr ( a ) + ¯ ρ wdm ( a )+ ¯ ρ γ ( a ) + ¯ ρ ν ( a ) + ¯ ρ b ( a ) + ¯ ρ Λ . (17)Here ¯ ρ γ , ¯ ρ ν , ¯ ρ b and ¯ ρ Λ denote the mean densitiesof photons, neutrinos, baryons and dark energy, respec-tively.We show the evolution of the Hubble parameter forΛCDM and several ΛDDM models in lower panel of Fig.1. In the upper panel of Fig. 1, we also show the evo-lution of the fractional densities Ω dcdm ( z ), Ω wdm ( z ) andΩ dr ( z ) for a particular ΛDDM model (Γ − = 30 Gyrsand ε = 0 .
1) as well as Ω cdm ( z ) with the same initialamount of dark matter, namely Ω = Ω inidcdm . Thetwo new parameters affect H ( a ) as follows: at fixed ε ,a shorter lifetime Γ − implies a lower Hubble parame-ter. This is clearly visible at z (cid:38)
1. The behavior be-low z ∼ H : to fulfill the budgetequation, a lower DM abundance requires a larger Ω Λ ,meaning that the period of accelerated expansion be-gins earlier with respect to ΛCDM. At fixed lifetime Γ − ,a larger ε also induces a lower H ( a ), because more en-ergy is converted into radiation, which dilutes faster. Wethus anticipate a negative correlation between Γ − and ε at the background level. The degeneracy can be cap-tured by the matter contribution from the WDM species,Ω mwdm ≡ Ω (1 − w ), which is well constrained by data.This quantity is roughly approximated byΩ mwdm (cid:39) Ω inidcdm (1 − e − Γ t ) √ − ε, (18) FIG. 1.
Upper − Redshift evolution of the abundances ofthe DCDM, WDM and DR species, assuming Γ − = 30 Gyrsand ε = 0 .
1. We also show the abundance for a standardCDM species with Ω = Ω inidcdm . Lower − Redshift evolu-tion of the Hubble parameter for the ΛCDM and four differ-ent ΛDDM models. The Hubble parameter today is fixed to H = 67 . / s / Mpc.
For small Γ and small ε , Ω mwdm ∝ Γ(1 − ε ) and we expectdata to constrain a parameter combination of the form ε ∝ Γ − . B. Linear perturbation equations
The continuity equation for the DCDM density per-turbation can be obtained either by using the covariantconservation of the stress-energy tensor T µν , or by inte-grating the corresponding first order Boltzmann equation(see Eq. (2.20) in Ref. [63]). As in the case of standardCDM, the resulting expression is˙ δ dcdm = − ˙ h . (19)To obtain the dynamical equations for the daughter par-ticles, we combine Eq. (5) with the corresponding firstorder Boltzmann equation (see Eq. (2.31) in Ref. [63]),so that we have the following hierarchy of equations: ∂ (∆ f j, ) ∂τ = − qk E j ∆ f j, + q ∂ ¯ f j ∂q ˙ h f j δ dcdm , (20) ∂ (∆ f j, ) ∂τ = qk E j [∆ f j, − f j, ] , (21) ∂ (∆ f j, ) ∂τ = qk E j [2∆ f j, − f j, ] − q ∂ ¯ f j ∂q ( ˙ h + 6 ˙ η )15 , (22) ∂ (∆ f j,(cid:96) ) ∂τ = qk (2 l + 1) E j [ l ∆ f j,(cid:96) − − ( (cid:96) + 1)∆ f j,(cid:96) +1 ]( (cid:96) ≥ , (23)where j = { dr, wdm } . The system of Eqs. (20)-(23) isin the same form as the one for massive neutrinos [71],except for the last term in Eq. (20) and the fact that thepartial derivative ∂ ¯ f j /∂q in Eqs. (20) and (22) is nowtime-dependent.Given that the DR species satisfies the condition q/ E dr = 1, the hierarchy of equations can be simplifiedby taking the moments F dr ,(cid:96) ≡ ρ c, (cid:90) ∞ dq πq q ∆ f dr ,(cid:96) (24)of Eqs. (20)-(23), and integrating over all the momentumdegrees of freedom, so that˙ F dr , = − kF dr , − r dr ˙ h + ˙ r dr δ dcdm , (25)˙ F dr . = k F dr , − k F dr . , (26)˙ F dr , = 2 k F dr , − k F dr , + 415 r dr ( ˙ h + 6 ˙ η ) , (27)˙ F dr ,(cid:96) = k (2 (cid:96) + 1) [ (cid:96)F dr ,(cid:96) − − ( (cid:96) + 1) F dr ,(cid:96) +1 ] ( (cid:96) ≥ . (28)We have adopted the convention r dr ≡ a ¯ ρ dr /ρ c, , asin Ref. [23], which in the 2-body decay scenario understudy leads to: ˙ r dr = εa Γ(¯ ρ dcdm / ¯ ρ dr ) r dr . (29)The first three multipoles are given by: F dr , = r dr δ dr , F dr , = 4 r dr k θ dr , F dr , = 2 σ dr r dr . (30)We choose the maximum multipole (cid:96) max to truncate thehierarchies of equations according to the scheme pro-posed in Ref. [71] for both massless and massive neu-trinos, i.e.˙ F dr ,(cid:96) max = kF dr ,(cid:96) max − − (cid:96) max + 1 τ F dr ,(cid:96) max , (31)∆ ˙ f wdm ,(cid:96) max = qk ∆ f wdm ,(cid:96) max − E wdm − (cid:96) max + 1 τ ∆ f wdm ,(cid:96) max . (32) C. A novel WDM fluid approximation scheme
In order to compute the WDM dynamics one cannotintegrate over the momentum degrees of freedom, as wedid in the DR case. Indeed, when taking the momentsof the hierarchy of Eqs. (20)-(23), higher velocity-weightintegrals appearing at (cid:96) = 2 cannot be computed fromthe system of equations itself.Therefore, one has to follow the evolution of thefull time-dependent PSD to obtain the elements of theperturbed stress-energy tensor δ wdm , θ wdm and σ wdm through Eqs. (6)-(9). A typical set-up for CMB analy-ses requires roughly 500 wavenumbers, 50 multipoles and10 momentum bins, i.e. O (10 ) linear differential equa-tions to be computed. On a single processor, this leadsto runs with a CPU time of 1 − O (10 ), with a CPU time per single run ∼ −
40 s. Thenovel approximation scheme is based on the treatment ofmassive neutrinos as a viscous fluid by Ref. [70], and itis only valid at scales deeply inside the Hubble radius,where high- and low- (cid:96) modes are effectively decoupled.In App. A we explicitly demonstrate the accuracy ofthis approximation.Similarly to the DR case, the fluid equations can bederived by multiplying both sides of Eqs. (20) and (21)by 4 πq E wdm a − and 4 πq ka − , respectively, and inte-grating over q . Then, by using Eqs. (6)-(9) and (14), onecan write down the continuity equation,˙ δ wdm = − H ( c − w ) δ wdm − (1 + w ) (cid:32) θ wdm + ˙ h (cid:33) + (1 − ε ) a Γ ¯ ρ dcdm ¯ ρ wdm ( δ dcdm − δ wdm ) , (33)and the Euler equation,˙ θ wdm = −H (1 − c g ) θ wdm + c w k δ wdm − k σ wdm − (1 − ε ) a Γ 1 + c g w ¯ ρ dcdm ¯ ρ wdm θ wdm . (34)Notice that we have introduced the WDM sound speedin the synchronous gauge, c ≡ δP wdm /δρ wdm , and theWDM adiabatic sound speed, c g ≡ ˙¯ P wdm / ˙¯ ρ wdm .The latter can be written as c g = w ( ˙¯ P wdm / ¯ P wdm )( ˙¯ ρ wdm / ¯ ρ wdm ) − , (35)and it can be computed as follows c g = w (cid:18) − p wdm ¯ P wdm − ¯ ρ dcdm ¯ ρ wdm a Γ3 w H ε − ε (cid:19) × (cid:20) w ) − ¯ ρ dcdm ¯ ρ wdm a Γ H (1 − ε ) (cid:21) − . (36)Here p wdm denotes the so-called pseudo-pressure, ahigher momenta integral of ¯ f wdm which is reduced to thestandard pressure in the relativistic limit [70].Obtaining an analytical expression for c is lessstraightforward, since we do not have a dynamical equa-tion for the pressure perturbation δP wdm . In Ref. [70] it isassumed that c is scale-independent and approximatelyequal to c g . For the WDM species, we have found thatthis assumption leads to accurate results for the CMBpower spectrum, but not for the matter power spectrum.In fact, calculations using the full Boltzmann hierarchyof Eqs. (20)-(23) reveal that c exhibits a particular k -dependence that cannot be captured by a backgroundquantity such as c g . In particular, c gets slightly en-hanced on scales k > k fs , where k fs is the free-streamingscale of the WDM species, defined as k fs ( τ ) ≡ (cid:114) H ( τ ) c g ( τ ) . (37)It is possible to gain a semi-analytic understanding of thisbehaviour by building a formal equation for the evolutionof c , as detailed in App. B. To account for such anenhancement, we adopt the following prescription c ( k, τ ) = c g ( τ ) [1 + (1 − ε ) T ( k/k fs )] , (38)where the function T ( x ) = 0 . √ x has been fitted to thesound speed obtained using the full Boltzmann hierar-chy, for the parameter values ε = 0 . , . , . , .
001 andΓ /H = 0 . , ,
10. The factor (1 − ε ) is inserted to makethe k -dependent correction vanishingly small close to therelativistic limit, where c (cid:39) c g (cid:39) / σ wdm onecould follow a similar approach to that of Ref. [70], wherethe authors obtained a dynamical equation for the neu-trino shear by means of an improved truncation schemeat (cid:96) max = 2. We tested the implementation of a general-ization of that equation suitable to the decaying case, butwe found it to be only relevant close to the relativisticcase ε (cid:39) /
2, when it reduces to the DR shear equa-tion from Ref. [22]. In this regime, the dynamics of thedaughter particles do not significantly impact the CMBand matter power spectra, we thus decided to not in-clude any dynamical equation for the shear of the WDMspecies, when switching to the fluid approximation. Inpractice, we simply set σ wdm to a constant value, ob-tained via integration of the second PSD multipole in theBoltzmann hierarchy (see Eq. (9)). We explicitly checkedthat this approach yields better results rather than sim-ply setting σ wdm = 0 in the fluid equations, or than usingthe DR shear equation of Ref. [22] when ε = 1 / CLASS [69] is coded with respect to the CDM. Nev-ertheless, for adiabatic initial conditions, one can choose θ ini , dcdm = θ ini , cdm = 0, such that the gauge co-movingwith CDM is also co-moving with DCDM at all times [59].Hence, hereinafter we neglect this irrelevant complica-tion. D. Dynamics of perturbations
Before discussing the signatures of varying the param-eters Γ and ε on the relevant cosmological observables,it is worth having a look at the linear perturbations ofthe mother and daughter particles. In Fig. 2 we showthe evolution of the linear density perturbations for theDCDM, WDM and DR species, corresponding to a modethat enters the horizon very early (i.e., k = 1 Mpc − ).To clarify the impact of the two extra free parameters,the perturbations are shown for DCDM lifetimes similarand smaller than the age of the universe (Γ − = H − ∼ . − = (10 H ) − ∼ . ε = 10 − and ε = 10 − ).One can see that initially, the perturbations of thedaughter species always track those of the mother, be-cause the coupling term dominates the dynamics (i.e.,the ratios ρ dcdm /ρ wdm and ρ dcdm /ρ dr are large). When amode crosses the free-streaming scale, the pressure sup-port of the daughter particles becomes important and theperturbations develop oscillatory features. For the DRspecies, the free-streaming scale simply corresponds tothe horizon ( k fs ∼ H ,) while for the WDM it correspondsto a larger value, given by k fs ∼ H /c g (see Eq. (37)).This time- and scale-dependent power suppression, to-gether with the different background dynamics, lead tokey signatures on the CMB and matter power spectra, aswe discuss later. By comparing upper and lower panelsof Fig. 2, it is clear that the value of the WDM free-streaming scale is essentially determined by the value of ε , as expected. On the other hand, by comparing left andright panels, one notices that the intensity of the oscilla-tions due to the pressure support can be compensated bythe coupling to the DCDM if the lifetime is long enough,as the daughter particles keep being produced.Interestingly, the decoupling time of the daughter per-turbations from the mother perturbations is always setby the free-streaming crossing time, and not by the char-acteristic decay time. In order to illustrate that, in Fig. 3we show perturbations corresponding to a smaller wave-number, k = 10 − Mpc − , that enters the horizon muchlater. For this mode, the WDM species does not havetime to cross the free-streaming length (the crossing willoccur in the future), so that δ dcdm and δ wdm remainequal, even if the lifetime is smaller than the age of theuniverse. This can be understood from the fluid Eqs. FIG. 2. Time evolution of the linear density perturbationsof the DCDM, WDM and DR species, corresponding to awavenumber k = 1 Mpc − . Each panel displays the pertur-bations for a different combination of the parameters Γ and ε . The black dashed and solid lines indicate the times of hori-zon crossing (at H ( τ ) = k ) and WDM free-streaming scalecrossing (at k fs ( τ ) = k ), respectively. The purple dotted lineindicates the characteristic decay time, given by t ( τ ) = Γ − .FIG. 3. Same as in Fig. 2, but corresponding to a wavenum-ber k = 10 − Mpc − . In this case, the WDM perturbationsnever cross the free-streaming scale. (33)-(34): when the decay term that includes Γ is rel-evant, δ wdm is driven by δ dcdm . Therefore, the WDMdensity perturbation δ wdm will continue to track the be-haviour of δ dcdm , as long as the free-streaming scale isnot crossed, i.e., as long as the pressure term, containing c s k δ wdm , is small compared to ¨ δ wdm ∼ δ wdm H . III. OBSERVABLE IMPACTS OF DARKMATTER DECAYS WITH WARM DAUGHTERSA. The linear matter power spectrum
Firstly, we focus on describing the effects of the 2-body decay on the linear matter power spectrum, sincethis will allow to better understand some of the effectson the CMB spectra. In this Section, we use the sameΛCDM parameters as in Section II A, except for fix-ing 100 θ s = 1 . H , to better connectwith CMB observations, that accurately pin down θ s . InFig. 4, we compare the residual differences in the linearpower spectra (at redshifts z = 0 and z = 2) with respectto our baseline ΛCDM. The left panel shows several life-times and a fixed DR energy fraction ε = 0 .
1, while theright panel shows a fixed lifetime Γ − = 30 Gyrs andseveral values of the DR energy fraction ε .One important feature of the C+WDM scenarios, suchas the one considered in this work, is that they are ex-pected to produce a suppression in the linear matterpower spectrum at scales smaller than a ‘cutoff scale’,with a non-trivial shape [72–75]. The cut-off scale isdetermined by the free-streaming scale of the WDMspecies, k fs , given by Eq. (37). On scales k > k fs , pres-sure becomes important and WDM particles cannot stayconfined in gravitational potential wells, which inhibitsstructure formation.Fig. 4 clearly illustrates that, while the parameter ε fixes the value of the cut-off k fs5 , the lifetime Γ − es-sentially determines the depth of the suppression at verysmall scales. This is to be expected, since the amountof power suppression grows with the WDM abundance,which increase for smaller lifetimes. In a similar way, theeffects of decay become less important when consideringthe matter spectrum at a higher redshift z = 2, since theabundance of WDM was smaller in the past.In general, for late-time decay scenarios (well after re-combination) as the one studied in this work, it is pos-sible to distinguish three different regimes depending onthe value of ε , as it is shown in the right panel of Fig. 4: • Non-relativistic decay : if ε (cid:46) .
001 (blackcurve), the WDM leaves the expansion rate un-affected, since its contribution to the matter den- By looking at Eqs. (36) and (37), we see that for small valuesof ε , the cut-off approximately satisfies the scalings k fs ∝ c − g ∝ w − / ∝ ε − . FIG. 4.
Left − Residuals of the linear matter power spectrum at z = 0 (upper) and z = 2 (lower) for several values of thelifetime Γ − = 10 , , ,
300 Gyrs and a fixed DR energy fraction ε = 0 .
1. Residuals are taken with respect our baselineΛCDM model.
Right − Same as in the left, but for several DR energy fractions ε = 0 . , . , . , .
001 and a fixed lifetimeΓ − = 30 Gyrs. sity, ¯ ρ wdm (1 − ω ) (cid:39) ¯ ρ wdm , compensates the reduc-tion in the DCDM density, ¯ ρ dcdm . In addition, theWDM free-streaming length is very small, induc-ing a power suppression at k (cid:38) h Mpc − . Suchscales are beyond the range of scales probed bythe observables considered in this work, so in thisregime the WDM is almost degenerate with stan-dard CDM. • Relativistic decay : if ε (cid:39) . . Moreover, thefree-streaming length k fs gets as large as the hori-zon, so that the WDM does not cluster at all. Thereduction in the Hubble friction is balanced by areduction in the clustering density of the daugh-ter particles, δρ wdm (cid:39) δρ dr (cid:39)
0, inducing a verylittle overall suppression of the growth of fluctua-tions. However, there is an impact coming fromanother background effect: the significant decreasein the co-moving matter density, a ρ m , leading inturn to an increase of the angular diameter dis-tance, thereby a reduction of θ s . This effect canbe compensated by increasing H , which shifts thelocation of the peak, k eq / ( a H ), towards smaller k in the matter power spectrum – since we are keep-ing the matter-radiation equality era fixed. Thenet effect on the residuals is twofold, a large-scaleenhancement and small-scale suppression of power. • Warm decay : For intermediate values of ε , Note that since now we are fixing 100 θ s instead of H , the Hubblerate H ( z ) can increase with respect to ΛCDM at z (cid:46)
1, once darkenergy starts to dominate. However, this effect is small for longlifetimes, and at early times the Hubble rate is still smaller thanin ΛCDM. namely 0 . (cid:46) ε (cid:46) . k fs thatdetermine the cut-off scale in the matter powerspectrum are not as small as in the case of non-relativistic decay, leading to δρ wdm (cid:39) k > k fs .Hence, on scales k > k fs the Hubble friction gets en-hanced with respect to the clustering density δρ tot ,slowing down the growth of DCDM perturbations.The net impact on the matter power spectrum isthus a suppression on intermediate scales, some-what similar to that induced by massive active neu-trinos [76–78]. B. The CMB angular power spectrum
We now discuss the impact of the 2-body decay sce-nario on the CMB anisotropy temperature and polar-ization angular power spectra, as well as on the lensingpotential power spectrum reconstructed from the CMB4-point correlation function. In the left panel of Fig. 5we report the residuals of the (lensed) TT, EE and lens-ing potential power spectra with respect to our baselineΛCDM, for different lifetimes Γ − and a fixed DR energyfraction ε = 0 .
1. The effects, more and more pronouncedas the lifetime decreases, can be understood as follows: • At the background level, the decay decreases thevalue of Ω m with respect to ΛCDM. This is com-pensated by an increase in Ω Λ (earlier beginningof Λ-domination) and thus an enhancement in theLate Integrated Sachs-Wolfe (LISW) effect, leav-ing a signature in the low- (cid:96) TT power spectrum.Furthermore, a modified background history altersquantities integrated along z , such as τ reio , which FIG. 5.
Left − Residuals (with respect our baseline ΛCDM model) of the CMB lensed TT (upper), EE (middle) and lensingpotential (lower) power spectra for several values of the lifetime Γ − = 10 , , ,
300 Gyrs and a fixed DR energy fraction ε = 0 . Right − Same as in the left, but for several values of the DR energy fraction ε = 0 . , . , . , .
001 and a fixed lifetimeΓ − = 30 Gyrs. impacts the multipoles (cid:96) ∼
10 in the EE powerspectrum. • At the perturbation level, the late-time reductionof a ¯ ρ dcdm implies a reduction of the quantity a δρ dcdm , which acts as a source of gravity throughthe Poisson equation. This induces a damping inthe metric fluctuations, and hence yields a fur-ther enhancement of the LISW effect. Further-more, the suppression in the matter power spec-trum and in Ω m lowers the amplitude of the lens-ing potential power spectrum, consequently reduc-ing the smoothing of the peaks in the high- (cid:96) part ofboth the TT and EE spectra, as one can see fromthe ‘wiggles’ in the corresponding plots.In the right panel of Fig. 5 we show the CMB residualsfor a fixed Γ − and various values of ε . The effects canbe readily understood: • At the background level, smaller values of ε weakenthe effects previously discussed, because the decayproduct dilute in a way similar to dark matter.Namely, the decrease in Ω m is less prominent dueto the significant WDM contribution (i.e., the in-crease in Ω Λ is shallower), and the impact on τ reio is smaller. Therefore, the signatures in the low- (cid:96) part of the TT and EE spectra become less visible. • At the level of perturbations, ε leads to some in-teresting signatures on the LISW effect and on thelensing potential. Since the LISW effect is only rel-evant for small multipoles (cid:96) (i.e., very large scales),one just needs to look at wavenumbers such that k < k fs . On these scales, the growth suppressiondoes not play any role, and the decrease in a δρ dcdm due to the decay gets partially compensated by theincrease in a δρ wdm , which is more significant forsmaller values of ε . Thus, the damping in the met-ric fluctuations is less relevant for smaller ε , reduc-ing the LISW enhancement. • Regarding the effects on the lensing potential, onecan see that the suppression in the correspond-ing power spectrum monotonically decreases forsmaller ε . Naively, one might expect the opposite,since we have argued that the matter power sup-pression increases for small ε . This can be under-stood by looking at the CMB lensing kernel W ( z ),given by the following expression [79] W ( z ) = 3Ω m H H ( z ) (1 + z ) χ ( z ) χ ∗ − χ ( z ) χ ∗ , (39)where χ ∗ is the co-moving distance to the last-scattering surface. Firstly, the CMB lensing kernelpeaks at z ∼
2, where the suppression is less impor-tant (see bottom panels of Fig. 4). Secondly, it getshighly suppressed for higher values of ε , and thiseffect can dominate over the effect on the matterpower spectrum. This suppression happens mainlydue to the smaller Ω m , as we have verified by com-puting W ( z ) for several ΛDDM models, with andwithout including the factor Ω m . Note that coinci-dentally, the effects on the lensing power spectrumat high (cid:96) are very similar for ε = 0 . ε = 0 . m is highly re-duced and the small scales power spectrum is al-most unaffected, while the opposite occurs in thelatter case.Let us finally remark that, even if the effects of vary-ing either ε or Γ on the observables are different, onecan easily exploit the degeneracy mentioned at the back-ground level to get different couples of values (large Γand small ε or vice-versa) with a similar cosmologicalsignature, especially on the CMB. We do indeed expectthe MCMC analysis to show a negative correlation in thereconstructed 2D posteriors for ε and Γ.0 IV. DATA ANALYSISA. Data and method
We now confront the ΛDDM cosmology to a host ofrecent cosmological observations. Our goal is to set con-straints on the lifetime of DM and the mass-ratio of themother and daugther particles. We also wish to check towhat extent the ΛDDM model can play a role in resolvingcosmological tensions as claimed in the recent literature[24, 67]. To that end, we perform comprehensive MCMCanalyses with the
MontePython-v3 [80, 81] code inter-faced with our modified version of CLASS , consideringvarious combinations of the following data-sets: • The BAO measurements from 6dFGS at z =0 .
106 [82], SDSS DR7 at z = 0 .
15 [83], BOSSDR12 at z = 0 . , .
51 and 0 .
61 [84], and thejoint constraints from eBOSS DR14 Ly- α auto-correlation at z = 2 .
34 [85] and cross-correlationat z = 2 .
35 [86]. • The measurements of the growth function f σ ( z )(FS) from the CMASS and LOWZ galaxy samplesof BOSS DR12 at z = 0 .
38, 0 .
51, and 0 .
61 [84]. • The Pantheon SNIa catalogue, spanning redshifts0 . < z < . • The low- (cid:96)
CMB TT, EE, and the high- (cid:96)
TT, TE,EE data + the gravitational lensing potential re-construction from Planck • The high- (cid:96)
CMB EE and TE (50 ≤ (cid:96) ≤ ≤ (cid:96) ≤ • The KIDS-1000+BOSS+2dFLens [91], DES-Y1 [2]and KIDS-1000+Viking+DES-Y1 [4] weak lens-ing data, compressed as a a split-normal likeli-hood, i.e., S = 0 . +0 . − . , S = 0 . +0 . − . , S = 0 . +0 . − . , respectively. • The local measurement of the Hubble constantfrom SH0ES , modelled with a Gaussian likelihoodcentered on H = 74 . ± .
42 km/s/Mpc [10].We start by performing two distinct sets of studies inorder to illustrate the importance of taking CMB datainto account when studying the ΛDDM scenarios, even https://github.com/brinckmann/montepython_public In our main analysis, we use the ‘lite’ version of the clik like-lihood. We verify that this leads to negligible differences withrespect to the full likelihood in App. C. A new version of the SH0ES measurement [9] was published dur-ing completion of this work. We do not expect it to have anyimpact on our conclusions. in the long-lived regime (i.e. when the DCDM decaysafter CMB decoupling):(i) a background-only analysis against BAO+FS andPantheon SNIa data; (ii) full analyses including linearperturbations, where we combine the data-set used in(i) with CMB TT, TE, EE + lensing data, with andwithout including the aforementioned informative priorson S and H . In the case of (i), the parameter space isfully characterized by the following free parameters: (cid:8) Ω inidcdm , H , Γ , ε (cid:9) , whereas in (ii) the whole parameter space is describedby: (cid:8) Ω b h , ln (cid:0) A s (cid:1) , n s , τ reio , Ω inidcdm , H , Γ , ε (cid:9) . For both (i) and (ii) we adopt logarithmic priors on ε and Γ , namely, − ≤ Log ε ≤ Log (0 . − ≤ Log (Γ / Gyrs − ) ≤ ≤ Ω inidcdm ≤ . Secondly, we explore the possibility of resolving the in-famous Hubble and S tensions, and the ‘ A lens ’ anomalythat exists within Planck data. We then test the robust-ness of our results to various changes in the pipeline, andin particular to trading the high- (cid:96)
Planck
CMB data forthose from the SPT collaboration, which are known to beless in tension with local S measurements. Finally, webriefly discuss the viability of the 2-body decay scenarioas solution for the Xenon1T anomaly [58].We follow Planck conventions for all other priors andfor the treatment of neutrinos, i.e. we set two masslessand one massive species with m ν = 0 .
06 eV [1]. Weassume our MCMC chains to be converged when theGelman-Rubin criterion R − < .
05 [92]. To extractthe best-fit parameters, we make use of the
Minuit al-gorithm [93] through the iMinuit python package . InApp. D we report all individual χ ’s per each of the anal-yses performed. B. General constraints: background vs. linearperturbations
In Fig. 6 we show the 2 σ posterior distributionsof the parameters characterizing the ΛDDM model, As discussed in Section I, we calibrate BAO data by imposing aGaussian prior on the sound horizon at recombination r s ( z rec ) =144 . ± . For comparisons with previous works, an useful conversion is thefollowing Log (Γ / Gyrs − ) (cid:39) Log (Γ / km s − Mpc − ) − . https://iminuit.readthedocs.io/ − − − Log (Γ / Gyrs − ) . . . . Ω m . . . S H − − − L og ( ε ) − − − Log ( ε )
66 68 70 72 H .
75 0 . S .
25 0 . Ω m BAO+SNIa+Planck
FIG. 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 measurement from KiDS-1000+BOSS+2dFLens. Log ( ε ) and Log (Γ / Gyrs − ). Red contours refer tothe background-only analysis (i), while blue contours re-fer to the full analysis (ii), at the perturbation level, whenCMB data are also added. From BAO+SNIa and Planck data, the ΛDDM model is not detected. We confirm ourexpection from Sec. II A and III that there is a negativecorrelation between Log ( ε ) and Log (Γ / Gyrs − ), ap-parent within all data sets: as ε decreases, the 95% C.L.upper limit on Γ relaxes. In practice, we find that for de-cay rates Γ ∼ − − − Gyrs − , our 2 σ exclusion curveis roughly described by ε (cid:39) . × − (Γ / Gyrs − ) − . .For large Γ and small ε , the factor ‘1 − e − Γ t ’ in Eq. 18reaches 1 faster than ‘ √ − ε ’ (assuming t (cid:39) . (Γ / Gyrs − ) (cid:46) − ( ε ) (cid:46) − . H slightly higher than the oneinferred assumed ΛCDM, but the full analysis pulls H back to its standard value. This suggests that this kindof models is not suitable for relieving the Hubble ten-sion, contrarily to earlier claims [24], and in agreementwith Refs. [25, 26]. However, in contrast to Ref. [26],we observe a significant decrease in the S contours for ε ∼ .
01 and Γ − ∼ Gyrs. We attribute this dis-agreement to the fact that Ref. [26] does not include atreatment of WDM perturbations, which are responsiblefor the suppression in the matter power spectrum. Hence,the 2-body decay presented here could potentially recon-cile the inferred value of S with its direct measurementsfrom LSS observations, as pointed out in Ref. [67]. Wepresent an explicit comparison of our constraints withthose from Ref. [26] in App. E.Finally, in order to compare our constraints on Γ withprevious literature, we have carried out a MCMC anal-2 Parameter BAO+SNIa +
Planck ω b – 2 . . +0 . − . Ω inidcdm . . +0 . − . . . +0 . − . H / [km / s / Mpc] 68 . . +0 . − . . . +0 . − . ln(10 A s ) – 3 . . +0 . − . n s – 0 . . +0 . − . τ reio – 0 . . +0 . − . log (Γ / [Gyr − ]) unconstrained(0.09) unconstrained( − . ( ε ) unconstrained(-2.89) − . − . +0 . − . Ω m . . +0 . − . . . +0 . − . S – 0 . . +0 . − . χ TABLE I. The mean (best-fit) ± σ errors of the cosmologicalparameters from our ΛDDM analyses against BAO + SNIaand BAO + SNIa + Planck . For each data-set, we also reportthe best-fit χ . ysis including BAO + SNIa + Planck data, but fixing ε = 0 . σ upper limit onthe DCDM decay rate of log (Γ / [Gyr − ]) (cid:46) − .
67, cor-responding to Γ − (cid:38)
468 Gyrs. Our constraints on DMdecays to DR are two to three times tighter than thosefound in previous works [23, 80] due to newer data setsand are the strongest to date. Part of the difference alsocomes from using a logarithmic prior on Γ as opposed tolinear.
C. Implications for cosmological tensions andXenon1T
1. The H tension In order to test the implications of the 2-body decayfor cosmological tensions, we conduct a run that includethe local measurement of H from SH0ES [10], CMB,BAO and SNIa data. For the sake of brevity we donot report the results of the runs here. We find thatthe shape of the posterior probabilities is almost un-changed, except for a tiny shift in H to a higher value, H = 68 . ± . / s / Mpc. We thus confirm the in-ability of this model to resolve the Hubble tension. Thiswas expected since it had already been shown througha model independent reconstruction of the late-time dy-namics of the dark sector that any late-time solution thatdoes not modify the sound horizon at recombination isexpected to fail when combining BAO with SNIa data(e.g. [17, 66]). This is not equivalent to directly reading the constraints on Γ at ε = 0 . ε contours, since the 95% C.L. derived froma χ distribution with different degrees of freedom correspond todifferent ∆ χ . − − − Log ( Γ /Gyrs − ) Ω m S − − − L o g ( e ) − − − Log ( e ) S Ω m KiDS-1000+BOSS+2dFLensDES-Y1KiDS-1000+Viking+DES-Y1
FIG. 7. 2D posterior distribution of a subset of parametersfrom analyses performed by imposing three different priorson S , based on the three sets of weak lensing measure-ments considered, whose mean values are shown as a solid,dashed and dotted line, for KiDS+BOSS+2dFLens, DES andKiDS+Viking+DES respectively. Parameter w/ DES w/ KiDS+Viking+DES ω b . . +0 . − . . . +0 . − . Ω inidcdm . . +0 . − . . . +0 . − . H / [km / s / Mpc] 67 . . +0 . − . . . +0 . − . ln(10 A s ) 3 . . +0 . − . . . +0 . − . n s . . +0 . − . . . +0 . − . τ reio . . +0 . − . . . +0 . − . log (Γ / [Gyr − ]) − . − . +0 . − . − . − . +0 . − log ( ε ) − . − . +0 . − . − . − . +0 . − . Ω m . . +0 . − . . . +0 . − . S . . +0 . − . . . +0 . − . χ TABLE II. The mean (best-fit) ± σ errors of the cosmologicalparameters from our BAO + SNIa + Planck + S , being thelatter from DES-Y1 and KiDS-1000+Viking+DES-Y1 analy-ses. For each data-set we also report the best-fit χ .
2. The S tension: robustness of the resolution and the roleof priors In Ref. [67], it was shown by including a S -prior asmeasured by KiDS+BOSS+2dFLens, that the ΛDDMmodel can fully resolve the tension if DM decays witha lifetime of Γ − (cid:39)
55 Gyrs and converts a fraction ε (cid:39) . σ , without changing significantly thevalue of Ω m . This results in a smaller S without run-ning into disagreement with measurements from Ω m fromBAO and SNIa data. Moreover, the fact that the decaysoccurs at late-times also leaves the CMB power spectralargely unchanged.In this former study, a split-normal likelihood on S as determined by KiDS-1000+BOSS+2dFLens [91], i.e. S = 0 . +0 . − . was used. A more accurate approachwould have been to include the full galaxy shear andclustering power spectra. Making use of the full likeli-hood would however require the ability to compute thematter power spectrum on non-linear scales in ΛDDMuniverse, a task that is beyond the scope of this paper.Let us note that it has been established in various casesthat the reconstructed S value only mildly vary fromone model to another. In particular, the KiDS collabo-ration has established that the reconstructed value of S is insensitive to the neutrino mass [3] – a model that hasphysical effects very similar to the ΛDDM model. Thisprovides confidence in making use of a prior on S derivedin the ΛCDM context. Nevertheless, to highlight the im-pact of a different S measurement, we replace the S prior from KiDS-1000+BOSS+2dFLens by the one de-termined in the combined analysis KiDS+Viking+DES-Y1 [4], S = 0 . +0 . − . and with the DES-Y1 dataonly [2], S = 0 . +0 . − . . The results of these analysisis reported in Tab. II, and shown in Fig. 7, where wedirectly compare the results of the three different analy-ses . From the reconstructed parameterslog (Γ / [Gyr − ]) = − . +0 . − . KiDS1000log ( ε ) = − . +0 . − . +BOSS+ 2dFLens S = 0 . +0 . − . , log (Γ / [Gyr − ]) = − . +0 . − KiDS+ Vikinglog ( ε ) = − . +0 . − . +DES − Y1 S = 0 . +0 . − . , log (Γ / [Gyr − ]) = − . +0 . − . DES − onlylog ( ε ) = − . +0 . − . S = 0 . +0 . − . , one can see that the level of detection is higher for theKiDS+Viking+DESY1 case while it is lower in the DES-only case. In the former case, we nevertheless note a Note that this combined analysis includes a photo-metric redshiftcorrection applied to DES result, slightly lowering the S valuecompared to what is advocated by the DES collaboration [2]. The mean ± σ errors, and the best-fit of the cosmological param-eters in the KiDS-1000+BOSS+2dFLens analysis are reported inthe first column of Tab. III Parameter w/o A lens w/ A lens ω b . . ± .
013 2 . . +0 . − . Ω inidcdm . . +0 . − . . . +0 . − . H / [km / s / Mpc] 67 . . +0 . − . . . +0 . − . ln(10 A s ) 3 . . +0 . − . . . +0 . − . n s . . ± . . . +0 . − . τ reio . . +0 . − . . . +0 . − . log (Γ / [Gyr − ]) − . − . +0 . − . − . − . +0 . − . log ( ε ) − . − . +0 . − . unconstrained (-2.47)Ω m . . +0 . − . . . +0 . − . S . . +0 . − . . . +0 . − . A TTTEEElens – 1 . . +0 . − . A φφ lens – 1 . . +0 . − . χ TABLE III. The mean (best-fit) ± σ errors of the cosmo-logical parameters from our BAO + SNIa + Planck + S (from KiDS+BOSS+2dFLens) analysis, with and withoutmarginalizing over the amplitude of the lensing potential A lens . For each analysis, we also report the best-fit χ . slight degradation of χ ∼ +1 compared to the ‘base-line’ analysis of Ref. [67], while in the latter case the χ does not sensibly change. This explicitly demonstratesthat the statistical significance of the ΛDDM detectionis strongly driven by the level of tension of the S valueused in the analysis. If the S tension increases in the fu-ture, the preference for ΛDDM over ΛCDM would likelyincrease. However, if the S tension disappears, cosmo-logical data would likely tightly constrain the ΛDDM sce-nario.
3. The A lens anomaly Another well studied ‘curiosity’ in the recent litera-ture consists in the anomalous amount of lensing esti-mated from the smoothing of the acoustic peaks at high- (cid:96) ’s within
Planck data, as quantified by the ‘ A lens ’ pa-rameter [94–96]. However, this anomalous ‘lensing’ is notsupported by the lensing power spectrum reconstruction,such that it is now commonly admitted that this tension(oscillating between the 2 − σ statistical level) cannotoriginate from a true lensing effect. On the other hand,it has been understood that this anomaly can be eas-ily resolved if the universe is closed [97–99], in certainmodified gravity theories [100] and in certain alterna-tive inflation scenarios [101]. Here instead, we wish tocheck whether this anomaly could impact constraints onthe ΛDDM model, and conversely if the ΛDDM modelcould help explaining the existence of such anomalies.In fact, it has already been noted that this anomalycould be related to the S tension: indeed, once includ-ing A lens as an extra free-parameter in the analysis, ithas been shown that the reconstructed cosmology hasa smaller A s and ω cdm (as well as a higher H ), show-ing no S tension, but a remnant ∼ . σ Hubble tension[102, 103]. To do so, we follow the approach of the SPT-pol collaboration and implement two new parameters in
CLASS that allows to (roughly) marginalize over the lens-4 − − − Log ( Γ /Gyrs − ) A φφ l e n s A TTTEEE l e n s S − − − L o g ( e ) − − − Log ( e ) S A TTTEEE lens A φφ lens w/ A lens w/o A lens FIG. 8. 2D posterior distribution of a subset of pa-rameters from our BAO + SNIa +
Planck + S (fromKiDS+BOSS+2dFLens) analysis, with and without includ-ing the extra A TTTEEElens and A φφ lens to marginalize over CMBlensing information. ing information in Planck . The parameter A TTTEEElens re-scales the amplitude of the lensing power spectrum en-tering in the high- (cid:96) part of the CMB TT,TE,EE spec-tra, while the parameter A φφ lens re-scales the amplitude ofthe lensing power spectrum reconstruction. We presentthe result of a MCMC analysis including Planck high- (cid:96)
TT,TE,EE+lensing+BAO+SNIa+ S data in Tab. IIIand in Fig. 8, with and without the two extra lensingparameters. One can see that the S parameter recon-structed once marginalizing over the ‘ A lens ’ anomaly islower by ∼ . σ than in the baseline analysis. How-ever, the preference for ΛDDM decreases, with log ( ε )now unconstrained. Note also that the A TTTEEElens is stillmore than 2 σ away than the fiducial value A TTTEEElens = 1.We can therefore conclude that the ΛDDM model cannotexplain this anomaly and that the preference for ΛDDMwould likely disappear if the S tension turns out to beexplained by a systematic in Planck data leading to theanomalous value of the A lens parameters.
4. Confronting
ΛDDM to SPTpol data
It is also interesting to test the robustness of theΛDDM detection to a change of CMB data sets, espe-cially given the impact of marginalizing over the ‘ A lens ’anomaly as discussed in previous section. We thus con-front the ΛDDM scenario under study against a set ofCMB data constituted by low- (cid:96) temperature and polar- − − − Log ( Γ /Gyrs − ) Ω m S − − − L o g ( e ) − − − Log ( e ) S Ω m Λ DDM w/o S Λ CDM w/o S Λ DDM w/ S FIG. 9. 2D posterior distribution of a subset of parameters inthe joint BAO + SNIa +
Planck + SPTpol analysis, with andwithout imposing a prior on S from KiDS+BOSS+2dFLens,compared to the ΛCDM scenario. Parameter w/o S w/ S ω b . . +0 . − . . . +0 . − . Ω inidcdm . . +0 . − . . . +0 . − . H / [km / s / Mpc] 68 . . +0 . − . . . +0 . − . ln(10 A s ) 3 . . +0 . − . . . +0 . − . n s . . +0 . − . . . +0 . − . τ reio . . +0 . − . . . +0 . − . log (Γ / [Gyr − ]) − . − . +0 . − . − . − . +0 . − . log ( ε ) unconstrained(-2.73) unconstrained(-2.57)Ω m . . +0 . − . . . +0 . − . S . . +0 . − . . . +0 . − . χ TABLE IV. The mean (best-fit) ± σ errors of the cosmologi-cal parameters from our BAO + SNIa + Planck
TT + SPTpolanalysis, with and without imposing a split-normal likelihoodon S (from KiDS+BOSS+2dFLens). For each data-set, wealso report the best-fit χ . ization as well as high- (cid:96) temperature data from Planck , incombination with high- (cid:96) polarization data from SPTpol(see section IV A for further details and references). Ithas been shown indeed that such a joint analysis predictsan amount of CMB lensing consistent with the ΛCDMexpectation [90], i.e. no ‘ A lens ’ anomaly, and no S ten-sion. This is manifest in Fig. 9, where we compare pre-dictions from the ΛCDM and the ΛDDM models, thelatter both with and without including information on S from KIDS1000+BOSS+2dFLens. As one can easilysee, both cosmological models predict a S value in ex-cellent agreement with the KIDS1000+BOSS+2dFLensmeasurement, displayed as a green horizontal band. Our5 − − − Log ( Γ /Gyrs − ) Ω m S − − − L o g ( e ) − − − Log ( e ) S Ω m BAO + SNIa + PlanckTT + SPTpol + S BAO + SNIa + Planck + S FIG. 10. 2D posterior distribution of a subset of parame-ters in our from our BAO + SNIa +
Planck + S (fromKiDS+BOSS+2dFLens) analysis, compared to the joint BAO+ SNIa + Planck + SPTpol + S analysis. results regarding the ΛDDM model are also reported inTab. IV.In Fig. 10, instead, we report a comparison betweenour baseline ΛDDM analysis and the Planck +SPT one.First and foremost, SPTpol appears in very good agree-ment with the ΛDDM model resolution of the S tensionrequired by Planck . However, the ΛDDM parametersare largely unconstrained in that case and no deviationsfrom ΛCDM are visible, which further establishes thatif the S tension turns out to be explained by a system-atic in Planck high- (cid:96) polarization data, the preference forΛDDM is likely to vanish.
5. Connections with the recent Xenon1T anomaly
Finally, following our work [67], let us further studythe implications of the ΛDDM model for the excess ofevents in the electronic recoils recently reported by theXenon1T Collaboration [58]. It has been pointed outthat this excess could be explained by the elastic in-teractions of electrons with a fast DM component ofmass m (cid:38) . . (cid:46) v/c (cid:46) . σ range for the speed of the daughter particle extendsup to v/c (cid:39) ε (cid:39) .
05. In order to test this hypothe-sis further, we perform another MCMC analysis includ-ing
Planck high- (cid:96)
TT,TE,EE+lensing+BAO+SNIa+ S ,with the DR energy fraction now fixed to ε = 0 .
05. This − − Log (Γ / Gyrs − ) . . . . Ω m . . . . S .
75 0 . S .
30 0 . Ω m BAO + SNIa + Planck + S FIG. 11. 2D posterior distribution of a subset of pa-rameters in our BAO + SNIa +
Planck + S (fromKiDS+BOSS+2dFLens) analysis when fixing the DR energyfraction to a value motivated by the recent Xenon1T anomaly( ε = 0 . ε = 0 . ω b . . +0 . − . Ω inidcdm . . +0 . − . H / [km / s / Mpc] 67 . . +0 . − . ln(10 A s ) 3 . . +0 . − . n s . . +0 . − . τ reio . . +0 . − . log (Γ / [Gyr − ]) − . − . +0 . − . Ω m . . +0 . − . S . . +0 . − . χ ± σ errors of thecosmological parameters from our Planck high- (cid:96)
TT,TE,EE+lensing+BAO+SNIa+ S analysis, when fix-ing the DR energy fraction to ε = 0 .
05. We also report thebest-fit χ . serves as a proxy for taking into account Xenon1T mea-surement (alternatively, one could enforce ε > . (Γ / Gyrs − ) (cid:39) − . S (cid:39) . χ (cid:39) +1 . S and Xenon1T anomalies. We leave theconstruction of a realistic model and study of the sub-sequent cosmological implications beyond the effect thedecay to another study (see [39, 57] for examples).6 V. CONCLUSIONS
In this paper, we have performed a comprehensivecosmological study of the CDM 2-body decay scenariodubbed ‘ΛDDM’, whereby decays are characterized bothby the decay rate Γ and energy fraction converted toradiation ε , including for the first time a fully consis-tent treatment of the linear perturbations of the WDMdaughter component. To that end, we have devised anew approximation scheme that allows to accurately andquickly compute the dynamics of the WDM linear pertur-bations by treating the WDM species as a viscous fluid.Our approximation scheme is accurate at the O (0 . O (1%) level in thelinear matter power spectrum (see App. A).We have then discussed in details the dynamics of lin-ear density perturbations of the mother and daughterparticles, as well as the physical effects of the ΛDDMmodel on the CMB and matter power spectra. We haveshown that accurate CMB (lensing in particular) andmatter power spectrum measurements have the poten-tial to detect both the decay rate Γ and energy fractionconverted to radiation ε .In a second part, we have performed a set of MCMCanalyses of the ΛDDM model against a suite of up-to-date low- and high-redshift data-sets. We have com-pared the constraints obtained from BAO and SNIa data,thereby solely based on background effects, to those ob-tained from the full Planck data-set, that requires in-stead an accurate description of the WDM linear per-turbations. We find that
Planck
CMB data constrainthe ΛDDM model ∼ S tension, as de-tailed in a previous paper [67].We have then tested the robustness of the ΛDDMresolution to the S tension to a number of change inthe analysis. We show that the mild preference for theΛDDM model over ΛCDM is tied to the S value cho-sen in the analysis. Concretely, the S value from theKiDS+Viking+DES analysis, which has a higher level oftension with the Planck
ΛCDM prediction than the base-line KiDS+BOSS+2dFLens value, leads to a strongerpreference in favor of the ΛDDM model. However, theDES-only result, which is in reasonable agreement with
Planck , leads to a weaker detection of the ΛDDM model.Similarly, once marginalizing over the lensing informationin
Planck through the A lens parameter (we used two ex-tra parameters in practice describing the normalization ofthe lensing power spectrum and the normalization of thelensing smoothing effect in the high- (cid:96) TT,TE,EE powerspectra), or when trading the
Planck high- (cid:96)
TT,TE,EEpower spectra for the SPTpol ones, the preference forthe ΛDDM model decreases. This is because in thesetwo cases the inferred ΛCDM model has a smaller S value, showing less of a tension with the weak lensingsurveys. This indicates that if the S tension ends up coming from an unknown systematic within weak lens-ing surveys or within Planck data, the preference for theΛDDM model would likely disappear.We have finally tested the intriguing possibility thatthe recent Xenon1T excess is due to the ΛDDM model.To that end, we have performed a additional MCMCanalysis fixing ε = 0 .
05 as required by Xenon1T. Wefind that it is easy to resolve the S tension in thatcase, pointing to a DCDM lifetime of log (Γ / [Gyr − ]) = − . +0 . − . . Interestingly, this comes at the cost of a veryminor degradation in Planck fit (∆ χ (cid:39) +1 . Planck , BAO and SNIa data are in good agreementwith this model.It will be very interesting to go beyond the linear as-pects discussed in this work and study the non-linearevolution of density perturbations, in order to be able tomake use of the full power of the KiDS, BOSS and DESlikelihoods. This could for instance be done with N-bodysimulations, as in Refs. [104–106], or via the EffectiveTheory of LSS [107, 108]. This will be even more crucialwith upcoming surveys such as Euclid [109], LSST [110],and DESI [111], which will measure the matter powerspectrum and the growth factor with great accuracy upto z ∼
2. It might also be possible to test the ΛDDMmodel with current Lyman- α forest flux power spectrumdata [20, 72–74, 105, 112]. We plan to study non-linearaspects and the discovery potential of these surveys in anupcoming work. ACKNOWLEDGMENTS
The authors are thankful to Julien Lavalle and Ro-drigo Calder´on for many useful comments and discus-sions. The authors acknowledge the use of computationalresources from the Ulysses SISSA/ICTP super-computerin Trieste, the CNRS/IN2P3 Computing Centre (CC-IN2P3) in Lyon, the IN2P3/CNRS and the Dark En-ergy computing Center funded by the OCEVU Labex(ANR-11-LABX-0060) and the Excellence Initiative ofAix-Marseille University (A*MIDEX) of the “Investisse-ments d’Avenir” programme.
Appendix A: Numerical implementation andaccuracy of the WDM fluid approximation
In order to solve the cosmic evolution equations de-rived in section II, we modified the publicly available nu-merical Boltzmann solver
CLASS [69, 70]. We now brieflyillustrate our implementation of the ΛDDM model.First let us notice that, when solving the backgroundequations for all cosmological species, the DE abundanceis iteratively derived through the budget equation, Ω Λ =1 − (cid:80) i Ω i , where the sum includes the current abundanceof all other components, which are not known a priori.We thus applied a shooting method for the aforemen-tioned parameter, i.e. we guess an initial Ω Λ , we solve7 FIG. 12. Residuals of the lensed CMB TT power spectrum(upper pannel) and EE power spectrum (lower pannel) in theWDM fluid approximation, with respect to the full hierarchycalculation, for a grid of values covering most of the parameterspace: Γ /H = 0 . , ,
10 and ε = 0 . , . , . , . Planck σ errors.FIG. 13. Residuals of the linear matter power spectrum (at z = 0), with respect to the baseline ΛCDM model, for thesame grid of parameter values considered in Fig. 12, bothfrom the full hierarchy calculation (solid lines) and the WDMfluid approximation(dashed lines). the system of background equations to obtain (cid:80) i Ω i , andre-compute Ω Λ . The procedure is iterated until conver-gence is achieved. The WDM density is computed bysolving Eq. (15) in 2800 momentum bins, approximatelyas many as the time-steps used to describe its backgroundevolution.At the linear perturbation level, we truncate the hier-archies of PSD multipole equations for both the daughterparticles at a (cid:96) max = 17. We set the initial conditions forthe WDM species following the same procedure of [63].On conformal times τ < τ q , we set all ∆ f wdm ,(cid:96) = 0,since no daughter particle with comoving momentum q could have been produced. On the crossing time τ = τ q ,one should be more careful, as the terms with ¯ f wdm inEqs. (20) and (22) contain a Dirac delta and, when in-tegrated, a Heaviside function. Thus, the correspondinginitial conditions for ∆ f wdm , ( τ q ) and ∆ f wdm , ( τ q ) arenot-vanishing. We set them according to the analyticalformulas (A.5) and (A.7) from [63].Finally, on times τ > τ q , we treat the WDM com-ponent as a massive neutrino species, and we solve thecorresponding hierarchy of equations in 300 momentum-bins. Such value ensures the needed accuracy, due to thefact that WDM hierarchy is only computed on super-Hubble and Hubble-crossing scales, where the dynamicsis simple and less computationally expensive. On sub-Hubble scales, when kτ is larger than a threshold value( kτ ) fluid , we switch-on the fluid approximation describedin section II C. The WDM dynamics is now described byEqs. (33)-(36) and (38). We have chosen ( kτ ) fluid = 25to provide the speed yet accurate enough for the purposesof the current analyses.In Figs. 12 and 13 we explicitly compare the novel ap-proximation scheme with the results of the “exact” com-putation for the WDM species. For the latter, we solvethe full Boltzmann hierarchy using 10 momentum-binsand (cid:96) max = 17. In Fig. 12 we show the residuals of thelensed CMB TT and EE power spectra in the WDM fluidapproximation, with respect to the full computation, fora grid of parameter values given by Γ /H = 0 . , ,
10 and ε = 0 . , . , . , . Planck uncertainties,which are indicated by the gray shaded regions.The predictions for the matter power spectrum are lessaccurate than for the anisotropy spectra, because it ismore sensible to the dynamics of the daughter particles.However, we are mainly interested in getting accuratepredictions for the departures from ΛCDM, rather thanthe exact shape of the matter spectrum itself. For thisreason, we have computed the residuals of the linear mat-ter power spectrum (at z = 0) with respect to our base-line ΛCDM model, for both the fluid and the full hier-archy calculations. The results are shown in Fig. 13,where we have spanned the same parameter values as inFig. 12. We can see that, for all the ΛDDM models, theshape of the power suppression (that is, the depth andthe cut-off scale) is excellently well-captured by our fluidprescription. In addition, we verified that the residuals inthe structure growth parameter S ≡ σ (Ω m / . . arealways smaller than the ∼ .
8% relative error of the S measurement from [91]. We thus conclude that the newWDM viscous fluid approximation is accurate enough forour analyses.8 Appendix B: Semi-analytic understanding of theWDM sound speed
Here we obtain a formal equation that dictates theevolution of the WDM sound speed in the synchronousgauge, c . The first natural step is to write a dynam-ical equation for the normalized pressure perturbation,Π wdm = δP wdm / ¯ ρ wdm . This can be achieved by multi-plying Eq. (20) by 4 πq q E wdm a − , integrating over q andthen using Eq. (14). By doing so, two higher velocity-weight integrals appear, namely δ P wdm ≡ π a (cid:90) ∞ dq q E ∆ f wdm , , (B1)and(¯ ρ wdm + ¯ P wdm )Θ wdm ≡ πka (cid:90) ∞ dq q E ∆ f wdm , . (B2)The variable Θ was already discussed in Ref. [70] inthe context of massive neutrinos. In the relativisticlimit δ P wdm and Θ wdm become equal to the standardvariables δP wdm and θ wdm (Eqs. (7)-(8)), while in thenon-relativistic limit they are suppressed by a factor( q/ E wdm ) with respect to them. This means that onecan write δ P wdm = δP wdm ω p and Θ wdm = θ wdm ω θ ,where ω p and ω θ are arbitrary functions, going from 1 / wdm reads˙Π wdm = − H Π wdm (cid:18) − ω p − ω (cid:19) − ω θ (1 + ω ) θ wdm − ˙ hω (cid:20) − p wdm p wdm (cid:21) + a Γ ¯ ρ dcdm ¯ ρ wdm (cid:20) ε (1 − ε ) δ dcdm − (1 − ε )Π wdm (cid:21) (B3)One can convert this into an equation for the soundspeed by noting that Π wdm = c δ wdm , and using thecontinuity equation (33). The final result reads: ∂c ∂τ = − H c (cid:18) − ω p − c (cid:19) − (1 + ω ) θ wdm δ wdm ( ω θ − c ) − ˙ h δ wdm (cid:20) ω (cid:18) − p wdm ¯ P wdm (cid:19) − c (1 + ω ) (cid:21) + a Γ ¯ ρ dcdm ¯ ρ wdm δ dcdm δ wdm (cid:20) ε − ε ) − (1 − ε ) c (cid:21) . (B4)We remark that the previous equation is highly non-linear in the perturbed quantities, so it can easily giverise to numerical instabilities. In addition, there is noclosed expression for computing ω p and ω θ . If these func-tions were scale independent, one possible approximationwould be to trade them for some background functions,such as w or c g . However, calculations using the fullhierarchy show that ω p and ω θ exhibit a k -dependence FIG. 14.
Left − The sound speed of the WDM species inthe synchronous gauge, as a function of conformal time andwavenumber.
Right − The ratio between the perturbed densi-ties of the DCDM and the WDM components, as a function ofconformal time and wavenumber. The black solid and dashedlines indicate the horizon and free-streaming crossing scales,respectively. We have set ε = 0 .
007 and Γ − = 55 Gyrs. similar to that of c s . For these reasons, we do not imple-ment Eq. (B4) in our code.Nonetheless, by making some simplifying assumptions,Eq. (B4) allows to qualitatively understand why there isa particular k -dependence of c in the decaying scenario,that is not present in the case of massive neutrinos. Letus consider the non-relativistic limit of Eq. (B4), sincedata favors in general very small DR energy fractions, ε (cid:28)
1. This also implies that w p , c , w (cid:28)
1. Let usfurther assume that ω θ and c behave similarly, so thatthe difference ω θ − c can be neglected. Finally, let us alsorestrict to sub-Hubble scales, for which the term ˙ h/δ wdm is very small and can be also neglected. In this case,Eq. (B4) simplifies to ∂c ∂τ = − H c − a Γ ¯ ρ dcdm ¯ ρ wdm δ dcdm δ wdm c . (B5)In absence of the decay term, we see that the sound speeddilutes as c ∝ a − , which is a well-known result formassive neutrinos. This dilution can be compensatedby the presence of the decay term, leading to a c ∼ cte, as long as the ratio δ dcdm /δ wdm doesn’t change . Inpractice, for scales and times such that k < k fs ( τ ), wehave δ dcdm /δ wdm = 1. In this regime, the sound speed c is well approximated by the adiabatic sound speed c g .However, when k > k fs ( τ ), δ wdm oscillates and starts tobecome suppressed with respect to δ dcdm , which leadsto oscillatory features and a small enhancement in theevolution of c . This is visible in Fig. 14, where we haveplotted c and δ dcdm /δ wdm in the k − τ plane using thefull Boltzmann hierarchy, for the best-fit parameters fromthe combined analysis of Ref. [67] (i.e., setting ε = 0 . − − − Log ( Γ /Gyrs − ) Ω m S − − − L o g ( e ) − − − Log ( e ) S Ω m Planck ”full”Planck ”lite”
FIG. 15. 2D posterior distribution for a subset of parametersobtained by using the
Planck ‘full’ likelihood, compared tothose from the
Planck ‘lite’ analysis. and Γ − = 55 Gyrs).One can see that this k -dependent effect appears onlybecause of the coupling term in Eq. (B5), which is notpresent for massive neutrinos. This also justifies why thesound speed c is well approximated by a backgroundfunction such as c g in the case of massive neutrinos. Inthe case of the WDM daughter species, the approxima-tion c (cid:39) c g will only work when k < k fs ( τ ). This moti-vates the use of the fitting formula introduced in Eq. (38),that accounts for the small enhancement at scales smallerthan the free-streaming scale. While this simple fittingformula is not able to capture the oscillatory features de-scribed previously, it leads to results that are accurateenough for all the observables analysed in this work. Appendix C: Consistency between the use of the
Planck
TTTEEE ‘lite’ and ‘full’ likelihoods
Given that for all our analyses we made use of the ‘lite’version of the
Planck likelihood, in Fig. 15 and Tab. VIwe show that our conclusions are not affected by such achoice, by comparing the results from our main run withthe predictions from a ‘full’
Planck analysis.
Appendix D: Best fit χ per experiment In Tab. VII we report all χ ’s obtained with the Mi-nuit algorithm [93] through the iMinuit python package
Parameter
P lanck ‘full’
P lanck ‘lite’100 ω b . +0 . − . . ± . inidcdm . +0 . − . . +0 . − . H / [km / s / Mpc] 67 . +0 . − . . +0 . − . ln(10 A s ) 3 . +0 . − . . +0 . − . n s . +0 . − . . ± . τ reio . +0 . − . . +0 . − . log (Γ / [Gyr − ]) − . +1 − . − . +0 . − . log ( ε ) − . +0 . − . − . +0 . − . Ω m . +0 . − . . +0 . − . S . +0 . − . . +0 . − . TABLE VI. The mean ± σ errors of the cosmological param-eters from our BAO + SNIa + Planck + S analysis, obtainedby using both the Planck ‘full’ and ‘lite’ likelihoods. for the various model and data-set combinations consid-ered in this work.
Appendix E: Comparison with the Planckconstraints from Ref. [26]
Here we carry out an explicit comparison of our con-straints with those of Ref. [26], which performed an anal-ysis of the ΛDDM model against
Planck data, neglectingthe perturbations of the warm daughter particles. Asshown in Fig. 16, we find that the constraints on theΛDDM models are up to (roughly) one order of magni-tude stronger when our improved treatment is consid-ered.
FIG. 16. Comparison between the 2 σ exclusion bounds(solid lines) from the Planck analysis of Ref. [26] and our
Planck +BAO+SNIa analysis. In each chase, the dashed lineindicates a fit that roughly describes the 2 σ limit in the rangeΓ ∼ − − − Gyrs − . BAO+SNIa +
Planck + Planck w/ S +SPTpolw/o S w/ A lens ε = 0 .
05 w/o S w/ S Pantheon SNIa 1026.9 1028.2 1027.5 1026.8 1029.2 1028.0 1027.1 1027.0 1026.9BAO+FS BOSS DR12 – 6.63 7.06 6.93 7.11 6.14 6.59 5.94 6.17BAO BOSS DR12 3.52 – – – – – – – –BAO BOSS low − z Planck high − (cid:96) TT,TE,EE ‘lite’ – 584.8 585.9 585.3 586.9 577.7 587.5 – –
Planck high − (cid:96) TT ‘lite’ – – – – – – – 207.8 207.8
Planck low − (cid:96) EE – 396.9 396.9 397.2 396.3 395.7 396.3 395.8 396.1
Planck low − (cid:96) TT – 23.1 23.1 23.2 23.0 22.1 23.0 22.3 22.1
Planck lensing – 8.78 9.12 8.88 9.47 8.53 9.83 – –SPTpol high − (cid:96) TE,EE – – – – – – – 145.9 145.6SPTpol lensing – – – – – – – 5.43 5.93KiDS+BOSS+2dFLens – – 0.0003 – – 0.0097 0.98 – 0.0015DES – – – 0.19 – – – – –KiDS+Viking+DES – – – – 0.20 – – – –total 1036.6 2053.4 2055.0 2054.8 2055.9 2043.2 2057.6 1816.3 1816.8TABLE VII. Best-fit χ per experiment (and total) for all the ΛDDM analyses performed in this work.[1] N. Aghanim et al. (Planck), (2018), 1807.06209.[2] T. Abbott et al. (DES), Phys. Rev. D , 043526 (2018),1708.01530.[3] H. Hildebrandt et al. , Astron. Astrophys. , A69(2020), 1812.06076.[4] S. Joudaki et al. , Astron. Astrophys. , L1 (2020),1906.09262.[5] L. Verde, T. Treu, and A. G. Riess, in Nature Astronomy2019 (2019) 1907.10625.[6] K. Aylor, M. Joy, L. Knox, M. Millea, S. Raghunathan,and W. K. Wu, Astrophys. J. , 4 (2019), 1811.00537.[7] K. C. Wong et al. et al. , 85 (2019), 1903.07603.[11] T. Shanks, L. Hogarth, and N. Metcalfe, Mon. Not. Roy.Astron. Soc. , L64 (2019), 1810.02595.[12] A. G. Riess, S. Casertano, D. Kenworthy, D. Scolnic,and L. Macri, (2018), 1810.03526.[13] T. M. Davis, S. R. Hinton, C. Howlett, and J. Cal-cino, Mon. Not. Roy. Astron. Soc. , 2948 (2019),1907.12639.[14] W. Yuan, A. G. Riess, L. M. Macri, S. Casertano, andD. Scolnic, Astrophys. J. , 61 (2019), 1908.00993.[15] W. Cerny, W. L. Freedman, B. F. Madore, F. Ashmead,T. Hoyt, E. Oakes, N. Q. H. Tran, and B. Moss, (2020),2012.09701.[16] J. Soltis, S. Casertano, and A. G. Riess, (2020),2012.09196.[17] L. Knox and M. Millea, (2019), 1908.03663.[18] S. Kumar and R. C. Nunes, Phys. Rev.
D94 , 123511(2016), 1608.02454. [19] R. Murgia, S. Gariazzo, and N. Fornengo, JCAP ,014.[20] M. Archidiacono, D. C. Hooper, R. Murgia, S. Bohr,J. Lesgourgues, and M. Viel, JCAP , 055, 1907.01496.[21] E. Di Valentino et al. , (2020), 2008.11285.[22] K. Enqvist, S. Nadathur, T. Sekiguchi, and T. Taka-hashi, JCAP , 067, 1505.05511.[23] V. Poulin, P. D. Serpico, and J. Lesgourgues, JCAP (08), 036, 1606.02073.[24] K. Vattis, S. M. Koushiappas, and A. Loeb, Phys. Rev. D99 , 121302 (2019), 1903.06220.[25] B. S. Haridasu and M. Viel, (2020), 2004.07709.[26] S. J. Clark, K. Vattis, and S. M. Koushiappas, (2020),2006.03678.[27] T. Hambye, PoS
IDM2010 , 098 (2011), 1012.4587.[28] V. Berezinsky, A. Masiero, and J. Valle, Physics LettersB , 382 (1991).[29] H. B. Kim and J. E. Kim, Physics Letters B , 18(2002).[30] L. Covi, J. E. Kim, and L. Roszkowski, Phys. Rev. Lett. , 4180 (1999), hep-ph/9905212.[31] J. L. Feng, A. Rajaraman, and F. Takayama, Phys. Rev.Lett. , 011302 (2003), hep-ph/0302215.[32] J. L. Feng, A. Rajaraman, and F. Takayama, Phys. Rev.D , 063504 (2003), hep-ph/0306024 [hep-ph].[33] R. Allahverdi, B. Dutta, F. S. Queiroz, L. E. Stri-gari, and M.-Y. Wang, Phys. Rev. D , 055033 (2015),1412.4391 [hep-ph].[34] K. N. Abazajian et al. , Light Sterile Neutrinos: A WhitePaper, (2012), arXiv:1204.5379 [hep-ph].[35] M. Drewes et al. , JCAP , 025, 1602.04816.[36] C.-R. Chen, F. Takahashi, and T. T. Yanagida, Phys.Lett. B , 71 (2009), 0809.0792.[37] G. Choi, M. Suzuki, and T. T. Yanagida, Phys. Rev. D , 075031 (2020), 2002.00036.[38] G. Choi, M. Suzuki, and T. T. Yanagida, Phys. Rev. D , 035022 (2020), 2004.07863. [39] G. Choi, M. Suzuki, and T. T. Yanagida, Phys. Lett. B , 135976 (2020), 2006.12348.[40] T. R. Slatyer and C.-L. Wu, Phys. Rev. D , 023010(2017), 1610.06933.[41] V. Poulin, J. Lesgourgues, and P. D. Serpico, JCAP ,043, 1610.10051.[42] M. Cirelli, E. Moulin, P. Panci, P. D. Serpico, andA. Viana, Phys. Rev. D , 083506 (2012), 1205.5283.[43] R. Essig, E. Kuflik, S. D. McDermott, T. Volansky, andK. M. Zurek, JHEP , 193, 1309.4091.[44] H.-B. Jin, Y.-L. Wu, and Y.-F. Zhou, JCAP , 026,1304.1997.[45] G. Giesen, M. Boudaud, Y. G´enolini, V. Poulin,M. Cirelli, P. Salati, and P. D. Serpico, JCAP , 023,1504.04276.[46] V. Berezinsky and J. W. F. Valle, Phys. Lett. B ,360 (1993), hep-ph/9309214.[47] M. Lattanzi, AIP Conf. Proc. , 163 (2008),arXiv:0802.3155 [astro-ph].[48] W. B. Lin, D. H. Huang, X. Zhang, and R. Bran-denberger, Phys. Rev. Lett. , 954 (2001), astro-ph/0009003 [astro-ph].[49] F. J. S´anchez-Salcedo, ApJL , L107 (2003), astro-ph/0305496 [astro-ph].[50] J. A. R. Cembranos, J. L. Feng, A. Rajaraman, andF. Takayama, Phys. Rev. Lett. , 181301 (2005), hep-ph/0507150.[51] M. Kaplinghat, Phys. Rev. D , 063510 (2005), astro-ph/0507300 [astro-ph].[52] L. E. Strigari, M. Kaplinghat, and J. S. Bullock, Phys.Rev. Lett. , 061303 (2007), astro-ph/0606281 [astro-ph].[53] F. Borzumati, T. Bringmann, and P. Ullio, Phys. Rev.D , 063514 (2008), hep-ph/0701007.[54] A. H. G. Peter, C. E. Moody, and M. Kamionkowski,Phys. Rev. D , 103501 (2010), 1003.0419 [astro-ph.CO].[55] A. H. G. Peter and A. J. Benson, Phys. Rev. D ,123521 (2010), 1009.1912.[56] K. Kannike, M. Raidal, H. Veerm¨ae, A. Strumia, andD. Teresi, (2020), 2006.10735.[57] S. Xu and S. Zheng, (2020), 2012.10827.[58] The Xenon-1t Collaboration, E. Aprile, J. Aalbers, andF. a. Agostini, (2020), 2006.09721.[59] B. Audren, J. Lesgourgues, G. Mangano, P. D. Serpico,and T. Tram, JCAP , 028, 1407.2418.[60] Z. Berezhiani, A. D. Dolgov, and I. I. Tkachev, Phys.Rev. D92 , 061303 (2015), 1505.03644.[61] A. Nygaard, T. Tram, and S. Hannestad, Updatedconstraints on decaying cold dark matter, (2020),arXiv:2011.01632 [astro-ph.CO].[62] S. Aoyama, K. Ichiki, D. Nitta, and N. Sugiyama, JCAP , 025, 1106.1984.[63] S. Aoyama, T. Sekiguchi, K. Ichiki, and N. Sugiyama,JCAP , 021, 1402.2972.[64] A. Chudaykin, D. Gorbunov, and I. Tkachev, Phys. Rev. D94 , 023528 (2016), 1602.08121.[65] T. Bringmann, F. Kahlhoefer, K. Schmidt-Hoberg, andP. Walia, Phys. Rev.
D98 , 023543 (2018), 1803.03644.[66] V. Poulin, K. K. Boddy, S. Bird, and M. Kamionkowski,Phys. Rev.
D97 , 123504 (2018), 1803.02474.[67] G. F. Abellan, R. Murgia, V. Poulin, and J. Lavalle,(2020), 2008.09615. [68] G. Blackadder and S. M. Koushiappas, Phys. Rev.
D90 ,103527 (2014), 1410.0683.[69] D. Blas, J. Lesgourgues, and T. Tram, JCAP , 034,1104.2933.[70] J. Lesgourgues and T. Tram, JCAP , 032, 1104.2935.[71] C.-P. Ma and E. Bertschinger, Astrophys. J. , 7(1995), astro-ph/9506072.[72] R. Murgia, A. Merle, M. Viel, M. Totzauer, andA. Schneider, JCAP , 046, 1704.07838.[73] R. Murgia, V. Irsic, and M. Viel, Phys. Rev. D ,083540 (2018), 1806.08371.[74] C. Miller, A. L. Erickcek, and R. Murgia, Phys. Rev. D , 123520 (2019), 1908.10369.[75] S. Bohr, J. Zavala, F.-Y. Cyr-Racine, M. Vogelsberger,T. Bringmann, and C. Pfrommer, (2020), 2006.01842.[76] J. Lesgourgues and S. Pastor, Massive neutrinosand cosmology, Phys. Rept. , 307 (2006), astro-ph/0603494.[77] E. Giusarma, M. Gerbino, O. Mena, S. Vagnozzi,S. Ho, and K. Freese, Phys. Rev. D , 083522 (2016),1605.04320.[78] G. Parimbelli, M. Viel, and E. Sefusatti, JCAP , 010,1809.06634.[79] A. Manzotti, Phys. Rev. D , 10.1103/Phys-RevD.97.043527 (2018), 1710.11038.[80] B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet,JCAP , 001, 1210.7183.[81] T. Brinckmann and J. Lesgourgues, Phys. Dark Univ. , 100260 (2019), 1804.07261.[82] F. Beutler, C. Blake, M. Colless, D. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, andF. Watson, Mon. Not. Roy. Astron. Soc. , 3017(2011), 1106.3366.[83] A. J. Ross, L. Samushia, C. Howlett, W. J. Percival,A. Burden, and M. Manera, Mon. Not. Roy. Astron.Soc. , 835 (2015), 1409.3242.[84] S. Alam et al. (BOSS), Mon. Not. Roy. Astron. Soc. , 2617 (2017), 1607.03155.[85] V. de Sainte Agathe et al. , Astron. Astrophys. , A85(2019), 1904.03400.[86] M. Blomqvist et al. , Astron. Astrophys. , A86(2019), 1904.03430.[87] D. M. Scolnic et al. , Astrophys. J. , 101 (2018),1710.00845.[88] J. Henning et al. (SPT), Astrophys. J. , 97 (2018),1707.09353.[89] F. Bianchini et al. (SPT), Astrophys. J. , 119 (2020),1910.07157.[90] A. Chudaykin, D. Gorbunov, and N. Nedelko, JCAP ,013, 2004.13046.[91] C. Heymans et al. (2020) 2007.15632 [astro-ph.CO].[92] A. Gelman and D. B. Rubin, Statist. Sci. , 457 (1992).[93] F. James and M. Roos, Comput. Phys. Commun. ,343 (1975).[94] N. Aghanim et al. (Planck), Astron. Astrophys. ,A95 (2017), 1608.02487.[95] N. Aghanim et al. (Planck), Astron. Astrophys. , A6(2020), 1807.06209.[96] G. Efstathiou and S. Gratton, (2019), 1910.00483.[97] W. Handley, Phys. Rev. D , L041301 (2021),1908.09139.[98] E. Di Valentino, A. Melchiorri, and J. Silk, Nature As-tron. , 196 (2019), 1911.02087. [99] G. Efstathiou and S. Gratton, Mon. Not. Roy. Astron.Soc. , L91 (2020), 2002.06892.[100] H. Moshafi, S. Baghram, and N. Khosravi, (2020),2012.14377.[101] G. Dom`enech, X. Chen, M. Kamionkowski, andA. Loeb, JCAP , 005, 2005.08998.[102] P. Motloch and W. Hu, Phys. Rev. D , 103536 (2018),1803.11526.[103] P. Motloch and W. Hu, Phys. Rev. D , 083515(2020), 1912.06601.[104] M.-Y. Wang and A. R. Zentner, Phys. Rev. D ,043514 (2012), 1201.2426.[105] M.-Y. Wang, R. A. C. Croft, A. H. G. Peter, A. R.Zentner, and C. W. Purcell, Phys. Rev. D , 123515 (2013), 1309.7354.[106] M.-Y. Wang, A. H. G. Peter, L. E. Strigari, A. R. Zent-ner, B. Arant, S. Garrison-Kimmel, and M. Rocha, Mon.Not. Roy. Astron. Soc. , 614 (2014), 1406.0527.[107] G. D’Amico, L. Senatore, and P. Zhang, JCAP , 006,2003.07956.[108] A. Chudaykin, M. M. Ivanov, and M. Simonovi´c,(2020), arXiv:2004.10607 [astro-ph.CO].[109] L. Amendola et al. , Living Rev. Rel. , 2 (2018),1606.00180.[110] D. Alonso et al. (LSST Dark Energy Science), (2018),1809.01669.[111] A. Aghamousa et al. (DESI), (2016), 1611.00036.[112] W. Enzi et al.et al.