© The Authors 2023

1. Introduction

There is overwhelming evidence for the existence of dark matter (DM), but we still know very little about its nature and composition. DM most probably consists of one or several new particles, requiring an extension of the standard model (Bertone et al. 2005; Feng 2010). The bulk of these particles has to be rather cold, interacts weakly at most, and is stable over at least one Hubble time. However, small deviations from these assumptions remain possible. Furthermore, a multi-particle DM sector would allow sub-species to evade the requirements mentioned above. They might be hot, interact strongly, or be very unstable, for instance.

In this paper, we focus on the possibility that a fraction or all of the DM fluid decays into radiation via a simple one-body decay channel. The nature of this radiation component is not specified and is not relevant to our analysis. Decay into photons or other standard model particles would lead to constraints from the absence of an observable radiation signal in the sky, however, which would exceed the constraints provided here. We therefore implicitly assume a DM decay into dark radiation.

Recent weak-lensing (WL) surveys such as CFHTLenS¹ (Heymans et al. 2012; Fu et al. 2014), KiDS² (Kuijken et al. 2019; Giblin et al. 2021; Hildebrandt et al. 2021; Asgari et al. 2021, A21), HSC³ (Aihara et al. 2017, 2022; Hamana et al. 2020; Liu et al. 2022), and DES⁴ (The Dark Energy Survey Collaboration 2005; Abbott et al. 2022; Amon et al. 2022) have reported a mild but persistent difference of the clustering amplitude σ₈ of the cosmic microwave background (CMB) as measured by the Planck satellite (Planck Collaboration I 2020; Planck Collaboration VI 2020; Planck Collaboration V 2020). This difference is usually quantified with the combined S₈ parameter, which is defined as $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$, with Ω_m being the total matter budget of the Universe. If a fraction of the DM were allowed to decay, the clustering signal at low redshift would be modified, which might provide a solution to the S₈ difference in principle, as was pointed out by Enqvist et al. (2015, E15), Berezhiani et al. (2015), Chudaykin et al. (2016) and Archidiacono et al. (2019). However, other authors have questioned these conclusions, showing that an agreement of the clustering amplitude between WL and the CMB cannot be easily achieved (Simon et al. 2022, S22; McCarthy & Hill 2022).

Independent of the S₈ difference, several works have focused on providing forecasts and constraints for the one-body decaying dark matter (DDM) model using a variety of data from Milky Way satellite counts (Mau et al. 2022), WL shear observations (E15; Enqvist et al. 2020, E20), and CMB data (S22). Most authors have focused on the assumption that all DM is unstable, while models with decaying sub-species as part of a more complicated DM sector were investigated only little (Poulin et al. 2016; S22).

In the present paper, we study the effect of a one-body DDM fluid on the temperature and polarization spectra from Planck and on the WL band power spectrum from the latest KiDS data release. We use the Boltzmann solver Class (Blas et al. 2011; Lesgourgues & Tram 2011) together with the nonlinear prescription of Hubert et al. (2021) to model the effects of DM decay on the high- and low-redshift Universe. Our goal is on one hand to re-investigate the effect of one-body decay on the S₈ difference, and on the other hand, to provide new constraints on the half-life time of DDM and on the fraction of decaying to total DM.

The paper is structured in the following way: In Sect. 2 we review the theoretical aspects of the one-body DDM model. Section 3 and 4 are dedicated to the presentation of our modelling pipeline, including the specifics of the Bayesian inference or Markov chain Monte Carlo (MCMC) process. In Sect. 5 we present our results, before we conclude in Sect. 6. We benchmark our ΛCDM pipelines in Appendix A and provide more details about our MCMC analyses in Appendix B.

2. Decaying dark matter model

The DDM consists of a minimal extension of the standard ΛCDM model, where DM particles, instead of being stable, decay into massless relativistic particles propagating at the speed of light. A phenomenological description of this model includes two parameters (in addition to those describing the ΛCDM model), namely the decay rate of the DM particles Γ and the fraction f of decaying to total DM budget. As a result, the matter is transformed into radiation affecting the background evolution of the Universe, that is,

$$\begin{array}{cc}& {\rho }_{\mathrm{dcdm}}^{\prime }+3\mathcal{H}{\rho }_{\mathrm{dcdm}}=-a\mathrm{\Gamma }{\rho }_{\mathrm{dcdm}},\hfill \end{array}$$(1)

$$\begin{array}{cc}& {\rho }_{\mathrm{dr}}^{\prime }+4\mathcal{H}{\rho }_{\mathrm{dr}}=a\mathrm{\Gamma }{\rho }_{\mathrm{dcdm}},\hfill \end{array}$$(2)

where derivatives are expressed with respect to conformal time, ℋ is a conformal Hubble parameter, and ρ_dcdm and ρ_dr are background densities of decaying cold DM and dark radiation, respectively (see e.g. Hubert et al. 2021 for more details about the DM decay process). When only a fraction (f) of the total DM is allowed to decay, we define

$$\begin{array}{c}\hfill f={\mathrm{\Omega }}_{\mathrm{ddm},\mathrm{ini}}/{\mathrm{\Omega }}_{\mathrm{dm},\mathrm{ini}},{\mathrm{\Omega }}_{\mathrm{dm},\mathrm{ini}}={\mathrm{\Omega }}_{\mathrm{ddm},\mathrm{ini}}+{\mathrm{\Omega }}_{\mathrm{cdm},\mathrm{ini}},\end{array}$$(3)

where Ω_ddm, ini, Ω_cdm, ini, and Ω_dm, ini are the decaying, stable, and total DM abundances at a time t ≪ τ = Γ⁻¹, that is, before the start of the decay process.

The background evolution of the Universe was modified as described in Eqs. (1) and (2). In particular, the source terms whose amplitudes are set by the decay rate Γ cause a decrease in the DM and an increase in radiation abundance. In Fig. 1 we show the evolution of the DM abundance between redshift 0 and 5 (solid lines). As expected, the DM abundance decreases towards low redshifts, whereas the amplitude of the effect depends on the decay rate (Γ). We also indicate the redshift range of the WL data from KiDS as well as the range of late-time integrated Sachs-Wolfe (ISW) effect as measured by Planck (see e.g. Nishizawa 2014). Both observables overlap with the regime in which the effects of DDM are most prominent, making them promising probes to constrain DM decays.

Fig. 1. Redshift evolution of the DM abundance for different DM decay rates (Ωdm) compared to the corresponding ΛCDM model (Ωdm, ΛCDM). We assume that all DM is unstable (f = 1). The green and pink arrows indicate the sensitivity ranges of the WL data from KiDS and the ISW effect from Planck.

The decay process affects not only the background evolution of the Universe, but also the process of structure formation. Since the scale factor a evolves at a somewhat slower rate (compared to ΛCDM), the Universe is less evolved, and the clustering process is therefore delayed. We are therefore left with suppression of power at small scales at a given redshift (see Fig. 2, described in the next section). This suppression becomes more pronounced with high Γ and with large f. Scenarios with Γ → 0 and f → 0 correspond to the ΛCDM model.

Fig. 2. Suppression of the linear (solid) and nonlinear (dashed) matter power spectrum due to one-body DDM. The left panel shows the effects of varying the decay rate Γ (in Gyr−1) for a models in which all DM is allowed to decay (f = 1.0). In the right panel, the fraction of decaying to total DM is varied, and the decay rate is fixed to 1/Γ = 13.5 Gyr.

3. Modelling pipeline

In this section, we provide details of our modelling pipeline for both the CMB and the WL observables. We specifically focus on nonlinear clustering, including the effects from DM decay, and we discuss our implementations of baryonic feedback and intrinsic alignment.

3.1. Cosmic microwave background modelling

Although originating from the early Universe, the CMB temperature fluctuations provide strong constraints on reduced models, even for half-life times of the order of (or longer than) a Hubble time. The reason for this behaviour is the late-time ISW effect, which causes a modification of the large-scale CMB modes due to the gravitational redshifting of the CMB photons that pass through evolving potential wells. Following Nishizawa (2014), we can write the ISW spectrum as

$$\begin{array}{c}\hfill C_l^{\mathrm{ISW}}=\frac{18}{{\pi }^2}{\mathrm{\Omega }}_{\mathrm{m}}^2{H}_0^4{\displaystyle \int \mathrm{d}kP(k){[\int \mathrm{d}r\mathcal{H}D(f_{\mathrm{lin}}-1)j_l(kr)]}^2,}\end{array}$$(4)

where D is the linear growth factor, and $f_{\mathrm{lin}}=\frac{\mathrm{d}lnD}{\mathrm{d}lna}$ is the velocity growth rate. The cosmological dependence of the ISW effect is governed by the term ℋD(f_lin − 1) and by the linear power spectrum. As discussed in Nishizawa (2014), the late-time ISW kernel starts to become important at z ∼ 3 − 4 steadily increasing towards z → 0 when the Universe becomes dominated by dark energy. In Fig. 1 we indicate the redshift range in which the CMB signal becomes sensitive to the ISW effect with a pink arrow.

To model the CMB including the ISW effect, we relied on the publicly available Boltzmann code Class⁵, which comes with the option to include DM that decays into dark radiation. We modelled high-ℓ TT, TE, and EE power spectra of the Planck 2018 CMB data using the lightweight version of the Planck plik likelihood, called plik_lite (Planck Collaboration XX 2016; Planck Collaboration V 2020), and mimicked SimAll (EE for 2 ≤ ℓ < 30) and Commander (TT for 2 ≤ ℓ < 30) likelihoods with the prior imposed on the optical depth parameter τ_reio based on Eq. (4) of Planck Collaboration VI (2020, P20b). Even though this approximation was obtained for the ΛCDM scenario, the cosmological parameters recovered from CMB after one-body decay is included are very close to ΛCDM; see Tables B.1 and B.2. This allows for this approximation. Furthermore, Abellán et al. (2021) compared the results of plik_lite and the full Plik likelihoods for the more general scenario of two-body decays (which includes our model as a limiting case) and reported that the retrieved parameters agreed well. The model parameters along with their prior ranges are listed in Table 1. To evaluate the likelihood, we used all 215 data points for TT (30 ≤ ℓ ≤ 2508) and 199 data points for TE and EE (30 ≤ ℓ ≤ 1996). We tested our inference pipeline for the ΛCDM model and obtain an agreement of ∼0.1σ compared to the findings of the Planck Collaboration (see Appendix A and Fig. A.2 for more details).

3.2. Weak-lensing modelling

To model WL cosmic shear observables, we followed the approach of Schneider et al. (2022, Sch22), with some changes as specified below. Most notably, we used the Pycosmo package (Refregier et al. 2018; Tarsitano et al. 2020) combined with Class to calculate the WL shear power spectra. For the nonlinear power spectrum, we relied on the revised halo model of Takahashi et al. (2012). We included massive neutrinos with a fixed mass of 0.06 eV following the recipe from P20b. For the intrinsic alignment component, we used the nonlinear alignment model (NLA) introduced by Bridle & King (2007) and described in Hildebrandt et al. (2016).

In the following, we describe some other aspects of the modelling pipeline. We specifically focus on the implementation of DM decay, the handling of baryonic effects, and the connection to the band power data from KiDS.

3.2.1. Decaying dark matter

To include the effects of one-body decay on the nonlinear matter power spectrum, we used the fitting function of Hubert et al. (2021), which corresponds to a modified version of the fit from E15. The function is defined by the ratio P_DDM(k, z)/P_ΛCDM(k, z) = 1 − ε_nonlin(k, z), where

$$\begin{array}{c}\hfill \frac{{\epsilon }_{\mathrm{nonlin}}(k,z)}{{\epsilon }_{\mathrm{lin}}}=\frac{1+a{(k/{\text{Mpc}}^{-1})}^p}{1+b{(k/{\text{Mpc}}^{-1})}^q}f,\end{array}$$(5)

with the factors a, b, p, and q given by

$$\begin{array}{cc}& a(\tau ,z)=0.7208+2.027\left(\frac{\mathrm{Gyr}}{\tau }\right)+3.031\left(\frac{1}{1+1.1z}\right)-0.18,\hfill \\ \hfill & b(\tau ,z)=0.0120+2.786\left(\frac{\mathrm{Gyr}}{\tau }\right)+0.6699\left(\frac{1}{1+1.1z}\right)-0.09,\hfill \\ \hfill & p(\tau ,z)=1.045+1.225\left(\frac{\mathrm{Gyr}}{\tau }\right)+0.2207\left(\frac{1}{1+1.1z}\right)-0.099,\hfill \\ \hfill & q(\tau ,z)=0.992+1.735\left(\frac{\mathrm{Gyr}}{\tau }\right)+0.2154\left(\frac{1}{1+1.1z}\right)-0.056.\hfill \end{array}$$

The remaining function ε_lin describes the redshift evolution of the suppression and is given by

$$\begin{array}{cc}& {\epsilon }_{\mathrm{lin}}(\tau ,z)=\alpha {\left(\frac{\mathrm{Gyr}}{\tau }\right)}^{\beta }{\left(\frac{1}{(0.105z)+1}\right)}^{\gamma },\hfill \end{array}$$(6)

where α, β, γ are functions of ω_b, h, and ω_m = ω_b + ω_dm, that is,

$$\begin{array}{cc}& \alpha =(5.323-1.4644u-1.391v)+(-2.055+1.329u\hfill \\ \hfill & +0.8672v)w+(0.2682-0.3509u)w^2,\hfill \\ \hfill & \beta =0.9260+(0.05735-0.02690v)w\hfill \\ \hfill & +(-0.01373+0.006713v)w^2,\hfill \\ \hfill & \gamma =(9.553-0.7860v)+(0.4884+0.1754v)w\hfill \\ \hfill & +(-0.2512+0.07558v)w^2.\hfill \end{array}$$

We defined u = ω_b/0.02216, v = h/0.6776 and w = ω_m/0.14116. The fitting function is able to reproduce results from N-body simulations with an error smaller than 1% up to k = 13 h Mpc⁻¹ (Hubert et al. 2021). In order to calculate the DDM matter power spectrum at nonlinear scales, we multiplied the term (1 − ε_nonlin) with the ΛCDM power spectrum from the revised halofit model of Takahashi et al. (2012).

In Fig. 2 we illustrate the effect of DDM on the linear (solid lines) and nonlinear (dashed lines) matter power spectrum. Different colours correspond to different decay rates (Γ) for a fixed f = 1 (left panel) and different fractions (f) for a half-life time 1/Γ = 13.5 Gyr (right panel). In general, DM decay leads to a suppression of power towards small scales. This effect is amplified by nonlinear clustering. The power suppression can be understood by the fact that the clustering in the DM model is delayed compared to ΛCDM, causing galaxy groups and clusters (which dominate the power spectrum signal) to form later.

3.2.2. Baryonic feedback

Baryonic feedback effects play an important role in the WL signal (e.g. Chisari et al. 2018; van Daalen et al. 2020; Aricó et al. 2021). They lead to suppression of the matter power spectrum, which may be of similar shape to the suppression due to DDM (Hubert et al. 2021; Amon & Efstathiou 2022). In order to account for potential degeneracies between the DM and the baryonic sector, it is therefore particularly important to model baryonic effects in the DDM cases.

We used the emulator BCemu (Giri & Schneider 2021), which includes the effects of baryonic feedback on the matter power spectrum. BCemu is based on the baryonification model described in Schneider & Teyssier (2015) and Schneider et al. (2019). It has seven free model parameters describing the specifics of the gas and stellar distributions around haloes, as well as one cosmological parameter that is the baryon (f_b = Ω_b/Ω_m). We fixed four of the seven parameters and only varied the gas parameters log₁₀M_c, and θ_ej, as well as the stellar parameter η_δ. Furthermore, the baryon fraction f_b was varied in accordance with the cosmological parameters. This three-parameter model has been shown in Giri & Schneider (2021) to match the power spectra from hydrodynamical simulations at the percent level for k ≲ 12.5 h Mpc⁻¹.

3.2.3. Cosmic shear angular power spectrum with KiDS-1000

The latest catalogue released by the Kilo-Degree Survey (KiDS-1000) contains shear information of over 20 million galaxies distributed inside five tomographic bins between z ∼ 0.1 and z ∼ 1.2 (Kuijken et al. 2019). We used the band power spectrum published in A21 using the auto and cross spectra of all five tomographic bins. The corresponding covariance matrix is from Joachimi et al. (2021).

To model the cosmic shear power spectrum components from gravitational lensing (G) and the intrinsic alignment of galaxies (I), we used the modified Limber approximation (LoVerde & Afshordi 2008; Kilbinger et al. 2017), that is,

$$\begin{array}{c}\hfill C_{\mathrm{AB}}^{(ij)}(\ell )={\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}\mathrm{d}\chi \frac{W_{\mathrm{A}}^{(i)}(\chi )W_{\mathrm{B}}^{(j)}(\chi )}{f_{\mathrm{K}}^2(\chi )}{P}_{\mathrm{m}}^{\mathrm{nonlin}}(\frac{\ell +1/2}{f_{\mathrm{K}}(\chi )},z(\chi )),}\end{array}$$(7)

where A, B ∈ {G, I}. χ is the comoving radial distance, and f_K(χ) is the comoving angular diameter distance. The window functions of the gravitational and intrinsic alignment components are given by

$$\begin{array}{cc}& W_{\mathrm{G}}^{(i)}(\chi )=\frac{3H_0^2{\mathrm{\Omega }}_{\mathrm{m}}}{2c^2}\frac{f_{\mathrm{K}}(\chi )}{a(\chi )}{\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{h}}}\mathrm{d}\chi \prime n_{\mathrm{S}}^{(i)}(\chi \prime )\frac{f_{\mathrm{K}}(\chi \prime -\chi )}{f_{\mathrm{K}}(\chi \prime )},}\hfill \end{array}$$(8)

$$\begin{array}{cc}& W_{\mathrm{I}}^{(i)}(\chi )=-A_{\mathrm{IA}}{\left(\frac{1+z(\chi )}{1+z_{\mathrm{pivot}}}\right)}^{{\eta }_{\mathrm{IA}}}\frac{C_1{\rho }_{\mathrm{cr}}{\mathrm{\Omega }}_{\mathrm{m}}}{D(a(\chi ))}{n}_{\mathrm{S}}^{(i)}(\chi ),\hfill \end{array}$$(9)

where D(a) is the linear growth factor, and the $n_s^{(i)}$ terms correspond to the redshift distribution of source galaxies for each tomographic bin (i). The term C₁ρ_cr was fixed to 0.0139, and z_pivot was set to 0.3 (see Joachimi et al. 2011).

From the angular shear power shown in Eq. (7), we calculated the band power spectrum following Joachimi et al. (2021). We refer to Sch22 for more details about this procedure. The prescription for cosmic shear modelling above does not strictly rely on ΛCDM. In our case, all relevant changes to the modelling enter via modifications of the nonlinear matter power spectrum.

4. Model inference

We used the emcee package (Foreman-Mackey et al. 2013) with the stretch move ensemble method in our MCMC analyses. For the WL and the CMB setup, we assumed multivariate Gaussian likelihoods. The convergence of the chains was checked with the Gelman-Rubin criterion assuming R_c < 1.1 (Gelman & Rubin 1992). In the case of the CMB analysis, we used the covariance matrix provided alongside the Plik_lite likelihood. For the WL analysis, we relied on the band power covariance matrix published by the KiDS collaboration (Joachimi et al. 2021).

In Table 1 we provide a summary of all model parameters, including information about their priors. For the CMB analysis and the WL analysis, we sampled over 9 and 12 parameters, respectively. The combined chains contain 13 free parameters. We used flat priors for all cosmological parameters except for the optical depth τ_reio, for which we assumed a Gaussian prior with a mean τ_reio = 0.0506 and standard deviation σ_τ = 0.0086, as explained in Sect. 3.1. For cold DM abundance ω_c and primordial power spectrum amplitude A_s, we used a prior wide enough to be uninformative. In the DDM scenario, ω_c stands for the initial cold DM abundance. In terms of CLASS input variables, we set omega_cdm = (1 − f)ω_c and omega_ini_dcdm = fω_c. In the case of parameters for which WL alone is not sensitive enough (ω_b, h₀, n_s, log₁₀M_c, θ_ej, and η_δ), we defined wide prior ranges following the analyses in A21 and Sch22. Regarding the baryonic parameters, the prior ranges are limited by the range of the emulator. They comfortably include all results from hydrodynamical simulations, however (Schneider et al. 2020a,b; Giri & Schneider 2021). For the Planck absolute calibration A_planck we followed the suggestion of the Planck Collaboration⁶ and choose Gaussian prior 𝒩(1.0, 0.0025). The adopted intrinsic alignment model (NLA) assumes two free parameters A_IA and η_IA entering via Eq. (9) of Sect. 3.2.3. Following A21, for example, we set η_IA = 0, and kept only A_IA as a free parameter.

We ran six chains in total, three assuming a ΛCDM cosmology, and three including the possibility of DM decay. The three runs refer to the CMB alone, the WL alone, and the combined setup. The main results from these chains in terms of DM constraints and cosmology are shown in the next section. Further details are provided in Appendix B, where we list the best-fit values and errors for all the parameters involved in the MCMC analysis.

5. Results

The main goal of this paper is to constrain DM decays with Planck and KiDS-1000 data. However, before showing the obtained limits on the decay rate and the fraction of decaying to total DM, we discuss the effect of the DDM scenario on the S₈ difference.

In the left panel of Fig. 3, we show the posterior contours of the σ₈ − Ω_m plane for our different data and modelling choices. For the case of ΛCDM, the results from KiDS and Planck are shown in black and blue, respectively. The best-fit values and 68% errors of the combined S₈ parameter are given by

$$\begin{array}{cc}& S_8=0.{735}_{-0.024}^{+0.031}(\mathrm{WL},\mathrm{\Lambda }\mathrm{CDM}),\hfill \end{array}$$(10)

$$\begin{array}{cc}& S_8=0.841\pm 0.017(\mathrm{CMB},\mathrm{\Lambda }\mathrm{CDM}),\hfill \end{array}$$(11)

Fig. 3. Implications of two-body decay on clustering amplitude σ8 and on S8 parameter from our analysis of CMB and WL data. Left: Two-dimensional posterior contours (68% and 95% confidence levels) of the Ωm − σ8 plane from our MCMC runs for the ΛCDM and DDM scenario. Results from KiDS-1000 are shown in blue and green, and results from Planck 2018 are shown in black and yellow. The combined CMB+WL analysis is shown in red and purple. Right: One-dimensional constraints of the S8 parameter for the same models. The original results from the KiDS-1000 (Asgari et al. 2021, A21) and the Planck 2018 (Planck Collaboration VI 2020, P20b) analyses are added in grey and blue for comparison.

corresponding to a difference of 3.0σ, which we obtained using the same conventional method as was used in A21 (see their Eq. (16)). These findings agree well with the original results from the KiDS (A21) and Planck (P20b) collaborations, as shown in the right panel of Fig. 3 and Appendix A.

The posterior contours of the DDM case are shown in yellow and green for Planck and KiDS, respectively. They do not show any visible shift with respect to the ΛCDM case, except that the KiDS contours become broader, especially towards lower values of σ₈ and Ω_m. We assume this to be the result of degeneracies between the baryonic and DDM parameters. Regarding the combined S₈ parameter, the best-fitting values and 68% errors are given by

$$\begin{array}{cc}& S_8=0.{723}_{-0.027}^{+0.041}(\mathrm{WL},\mathrm{D}DM),\hfill \end{array}$$(12)

$$\begin{array}{cc}& S_8=0.841\pm 0.017(\mathrm{CMB},\mathrm{D}DM),\hfill \end{array}$$(13)

yielding an S₈ difference of 2.7σ. This small decrease in the difference is not due to a better concordance of the S₈ values, but rather to a general increase in the error budget in the DDM case of the constraints derived from the WL data.

The above point can be further quantified by investigating the decrease in the minimum chi-squared (${\chi }_{\text{min}}^2$) from the standard ΛCDM to the DDM model. The change in the Akaike information criterion $\mathrm{\Delta }\text{AIC}=\mathrm{\Delta }{\chi }_{\text{min}}^2+2(N_{\text{DDM}}-N_{\mathrm{\Lambda }\text{CDM}})$, which compensates for the increase in the goodness of fit due to the increased parameter space, gives

$$\begin{array}{cc}& \mathrm{\Delta }\mathrm{AIC}=4.0(\mathrm{WL}),\hfill \end{array}$$(14)

$$\begin{array}{cc}& \mathrm{\Delta }\mathrm{AIC}=3.9(\mathrm{CMB})\hfill \end{array}$$(15)

for the WL and the CMB case. In the definition of ΔAIC, $\mathrm{\Delta }{\chi }_{\text{min}}^2$ stands for the difference of ${\chi }_{\text{min}}^2$ between DDM and ΛCDM and N_DDM (N_ΛCDM) denotes the number of free parameters in the DDM (ΛCDM) model. Despite two more parameters, the decrease in $\mathrm{\Delta }{\chi }_{\text{min}}^2$ is not sufficient in the DDM case compared to ΛCDM (models for which the increased number of free parameters is compensated for by the better goodness-of-fit result in ΔAIC < 0).

Although there is a remaining difference between the ΛCDM and DDM models, we ran combined chains for both scenarios. In Fig. 3 they are given by the red and purple contours (left panel) and data points (right panel). As expected, the posteriors from the combined MCMC runs are located in between the two original contours. Based on these results, we can now quantify the general difference between the KiDS and Planck datasets. Following Raveri & Hu (2019), we defined a difference in the maximum a posterior (MAP), which takes the full multi-dimensional posterior distribution into account,

$$\begin{array}{cc}& Q_{\mathrm{MAP}}={\chi }_{\min ,\mathrm{comb}}^2-({\chi }_{\min ,Planck}^2+{\chi }_{\min ,\mathrm{KiDS}}^2).\hfill \end{array}$$(16)

With this definition, the difference between the two probes can be expressed as $T=\sqrt{Q_{\mathrm{MAP}}}$. For the ΛCDM and the DDM scenarios, we obtain

$$\begin{array}{cc}& T=3.4\sigma (\mathrm{\Lambda }\mathrm{CDM}),\hfill \end{array}$$(17)

$$\begin{array}{cc}& T=3.4\sigma (\mathrm{DDM}).\hfill \end{array}$$(18)

From the difference in MAP, we conclude that one-body decays do not lower the mutual difference between KiDS-1000 and Planck 2018 data. A summary of the minimum χ² values as well as the scores of the AIC and MAP is provided in Table 2.

We now turn our attention towards the constraints on one-body decay obtained by the CMB, WL, and combined datasets used in this paper. The two-dimensional constraints for the DDM parameters Γ and f are illustrated in Fig. 4. All limits are provided at the 95% confidence levels. The contours exhibit the expected hyperbolic shape, excluding the regime in the top right corner of high decay rates and larger fractions of decaying to total DM. The results from Planck (yellow contours) show much stronger constraining power than those from the KiDS data. This means that the ISW effect is currently more sensitive to DM decay than WL. However, this is likely to change in the near future due to new WL observations from Euclid (Hubert et al. 2021).

Fig. 4. Constraints on the decay rate (Γ) and the decaying-to-total DM fraction (f) for the one-body DDM scenario. Weak lensing (WL) results from KiDS-1000 and CMB results from Planck 2018 are shown in green and yellow. The combined CMB + WL analysis is shown in purple. For the limiting case of f = 1, we add other results from recent studies by Enqvist et al. (2015, E15; grey arrow), Enqvist et al. (2020, E20; red arrow), and Simon et al. (2022, S22; blue arrow). All results are provided at the 95% confidence level.

The combined CMB + WL constraints, shown as purple contours in Fig. 4, are comparable in strength to the CMB-only limits. The small differences between f = 0.2 − 0.9 are most likely caused by the inherent differences between the KiDS and Planck datasets. A similar behaviour has been reported by E15.

In Fig. 4 we compare our results to several recent studies from the literature. Because these studies only provide constraints for the limiting case of f = 1, they were added as arrows in the top part of the plot. We show findings from E15, E20, and S22. E15 used Planck 2013 data combined with nine-year WMAP polarization measurements and WL data from CFHTLens (Heymans et al. 2013), reporting Γ⁻¹ ≥ 140 Gyr. E20 combined Planck 2015 CMB data with Planck 2015 SZ cluster counts and KiDS 450 WL observations obtaining Γ⁻¹ ≥ 154 Gyr. S22 combined CMB data from Planck 2018 with the Pantheon dataset and BAO from BOSS, 6dFGS, and SDSS DR7. They reported Γ⁻¹ ≥ 260.4 Gyr.

For our limiting case of f = 1, we obtain a half-life time of Γ⁻¹ ≥ 288 Gyr. This limit is slightly stronger than that of S22 and is significantly stronger than those of E15 and E20. Compared to E15 and E20, we used a more recent dataset which is probably responsible for strengthening the constraints. Compared to S22, the differences are much smaller, which is expected because both studies used Planck 2018 data for the analysis. For the case of high decay rates and small decaying to total DM fractions, we obtain limits of f < 0.34, f < 0.07, and f < 0.03 for the KiDS, Planck, and the combined analysis at the 95% confidence level.

6. Conclusions

We have investigated the one-body DDM scenario and its effects on structure formation in the light of CMB TTTEEE data from Planck 2018 and the cosmic shear angular power spectra from the KiDS-1000 data release. The free parameters of the DDM model are the decay rate (Γ) and the decaying to total DM fraction (f).

We obtained new constraints on Γ and f from the CMB, from WL, and from the combined CMB + WL analysis. In agreement with previous results, we find that the CMB constraints are stronger than those from WL alone. This apparently surprising result is due to the ISW effect, which provides strong constraints on the late-time background evolution of the Universe.

For the limiting case of f = 1, we obtain Γ⁻¹ ≥ 288 Gyr, which is stronger than previous constraints from E15 and E20 and similar to the findings of S22. For high decay rates (Γ ∼ H₀), on the other hand, we find a limit on the decaying to total DM fraction of f < 0.03, which is based on the combination of CMB and WL data. The CMB alone provides weaker constraints of f < 0.07.

Along with the derivation of new constraints on the one-body DDM scenario, we also investigated the effect of decay on the S₈ difference reported for example by the KiDS collaboration (Heymans et al. 2012). At face value, we find a slight reduction of the difference from 3.0σ to 2.7σ from a ΛCDM to a DDM model. We showed, however, that this reduction is entirely caused by the increase in free parameters. Our maximum a posteriori probability analysis (MAP) yields no improvement from a ΛCDM to a DDM scenario.

We conclude that there is currently no evidence for a DM sector featuring one-body decay from matter to radiation. For most of the parameter space, current WL observations are not constraining enough to compete with the stringent limits obtained from the CMB radiation via the integrated Sachs-Wolfe effect. In the near future, however, results from stage-IV lensing surveys such as Euclid are expected to probe currently untested DDM scenarios.

Acknowledgments

This work is supported by the Swiss National Science Foundation under the grant number PCEFP2_181157. Nordita is supported in part by NordForsk.

References

Appendix A: ΛCDM benchmark

We compared the results of our pipeline in the case of ΛCDM to the original results of KiDS collaboration A21 and to Sch22 using a similar approach of band power modelling. For brevity, we only present the Ω_m − σ₈ contours shown in Fig. A.1. The difference in 1- and 2σ contours is marginal compared to the original KiDS-1000 results. The most significant differences (to our best knowledge) arise from the choice of baryonic prescription; the KiDS-1000 pipeline uses the one-parametric HMCODE baryonic feedback model (Mead et al. 2015), while this work adopted a four-parametric (three baryonic parameters plus the baryon-to-matter ratio) version of BCemu (Giri & Schneider 2021). Compared with Sch22, the 68% and 95% confidence intervals are slightly more extended. The modelling pipelines, which otherwise are very similar, employ a different number of baryonic parameters, specifically, eight (seven baryonic parameters plus the baryon-to-matter ratio) in the case of Sch22 and four in this work. This results in broader posterior contours in the case of Sch22. Quantitatively, our results agree better with those of A21 at 0.2σ for σ₈ and at the level of 0.1σ for Ω_m.

Fig. A.1. Comparison of ΛCDM posterior contours in the Ωm − σ8 - plane obtained in this work with two recent studies that modelled the cosmic shear band power signal.

In Fig. A.2 we show a comparison of our ΛCDM results with CMB Planck 2018 data and compare them to the result published in P20b (see Tab. 2, setup TT,TE,EE+lowE. We also included lowT at the top of this setup). We display all six inferred cosmological parameters centred on Planck values and normalized by Planck 1σ confidence intervals, thus displaying $(x-{\overline{x}}_{\mathrm{planck}})/{\sigma }_{\mathrm{x},\mathrm{planck}}$ for a parameter x. Most of the parameters agree to ∼0.1σ. The largest discrepancy is observed for n_s at the level of 0.6σ.

Fig. A.2. Values and 1σ confidence intervals of cosmological parameters resulting from our ΛCDM analysis related to the values obtained by Planck Collaboration VI (2020, P20b), (TT,TE,EE+lowE). We also guide the eye by depicting the 1, 2, and 3σ intervals as grey bands.

Appendix B: MCMC results

We present the detailed results of our MCMC analyses in Tab. B.1 (ΛCDM) and B.2 (DDM). In the top part of the tables, we show cosmological, baryonic, and DDM parameters directly sampled during the MCMC. The middle part displays the derived Ω_m, σ₈ and S₈ values, and the bottom part is dedicated to the details about the MCMC statistics (priors, likelihoods, and χ² values). A long dash indicates that a specific parameter is not relevant for a specific setup, and unconst indicates unconstrained parameters.

All Tables

All Figures

	Fig. 1. Redshift evolution of the DM abundance for different DM decay rates (Ωdm) compared to the corresponding ΛCDM model (Ωdm, ΛCDM). We assume that all DM is unstable (f = 1). The green and pink arrows indicate the sensitivity ranges of the WL data from KiDS and the ISW effect from Planck.
In the text

	Fig. 2. Suppression of the linear (solid) and nonlinear (dashed) matter power spectrum due to one-body DDM. The left panel shows the effects of varying the decay rate Γ (in Gyr−1) for a models in which all DM is allowed to decay (f = 1.0). In the right panel, the fraction of decaying to total DM is varied, and the decay rate is fixed to 1/Γ = 13.5 Gyr.
In the text

Fig. 3. Implications of two-body decay on clustering amplitude σ8 and on S8 parameter from our analysis of CMB and WL data. Left: Two-dimensional posterior contours (68% and 95% confidence levels) of the Ωm − σ8 plane from our MCMC runs for the ΛCDM and DDM scenario. Results from KiDS-1000 are shown in blue and green, and results from Planck 2018 are shown in black and yellow. The combined CMB+WL analysis is shown in red and purple. Right: One-dimensional constraints of the S8 parameter for the same models. The original results from the KiDS-1000 (Asgari et al. 2021, A21) and the Planck 2018 (Planck Collaboration VI 2020, P20b) analyses are added in grey and blue for comparison.

In the text

Fig. 4. Constraints on the decay rate (Γ) and the decaying-to-total DM fraction (f) for the one-body DDM scenario. Weak lensing (WL) results from KiDS-1000 and CMB results from Planck 2018 are shown in green and yellow. The combined CMB + WL analysis is shown in purple. For the limiting case of f = 1, we add other results from recent studies by Enqvist et al. (2015, E15; grey arrow), Enqvist et al. (2020, E20; red arrow), and Simon et al. (2022, S22; blue arrow). All results are provided at the 95% confidence level.

In the text

	Fig. A.1. Comparison of ΛCDM posterior contours in the Ωm − σ8 - plane obtained in this work with two recent studies that modelled the cosmic shear band power signal.
In the text

	Fig. A.2. Values and 1σ confidence intervals of cosmological parameters resulting from our ΛCDM analysis related to the values obtained by Planck Collaboration VI (2020, P20b), (TT,TE,EE+lowE). We also guide the eye by depicting the 1, 2, and 3σ intervals as grey bands.
In the text