© The Authors 2024

1. Introduction

The standard Λ-cold dark matter (ΛCDM) cosmology provides an outstanding description of the Universe, explaining a wide range of cosmic observations such as the cosmic microwave background (CMB), baryonic acoustic oscillations (BAO), large-scale structure formation, and Big Bang nucleosynthesis (BBN). Despite its tremendous success, there are still questions to which ΛCDM, as understood nowadays, cannot provide satisfying answers, including the nature of dark matter and dark energy. Moreover, with progressively precise measurements, several discrepancies within the model have emerged. An example of such a discrepancy is the mild yet persistent clustering amplitude tension between the CMB and weak lensing (WL) measurements, often expressed via the parameter $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$, with σ₈ and Ω_m describing the clustering amplitude and the total matter abundance. More specifically, CMB measurements from the Planck Collaboration have yielded S₈ = 0.834 ± 0.016 (Planck Collaboration VI 2020), while a variety of low-redshift surveys report consistently lower values. For example, the Kilo-Degree Survey (KiDS; Kuijken et al. 2019; Giblin et al. 2021; Hildebrandt et al. 2021) obtained a value of $S_8=0.{760}_{-0.038}^{+0.016}$ (Asgari et al. 2021), in agreement with (albeit slightly lower than) the results from the Hyper Supreme-Cam (HSC; Aihara et al. 2017; Hamana et al. 2020) and the Dark Energy Survey (DES; The Dark Energy Survey Collaboration 2005; Amon et al. 2022).

It remains unclear whether the S₈ tension emerges from an insufficient modelling of the non-linear clustering of matter (e.g. Tan et al. 2023; Aricò et al. 2023) or the modelling of cosmic shear (García-García et al. 2021) or whether one has to look beyond the standard ΛCDM. Resolving the S₈ tension can be achieved by suppressing the matter power spectrum at scales k ∼ 0.1–1 h Mpc⁻¹, which most substantially influence the clustering amplitude value S₈ (e.g. Amon & Efstathiou 2022). Such suppression may be obtained by a number of extensions of ΛCDM, such as cold-warm dark matter (Schneider et al. 2020a; Parimbelli et al. 2021), cannibal dark matter (Pappadopulo et al. 2016; Heimersheim et al. 2020), models involving interactions between dark matter and dark radiation at early times (Cyr-Racine et al. 2016; Vogelsberger et al. 2016; Rubira et al. 2023; Joseph et al. 2023), scenarios introducing an interaction between dark matter and dark energy (Pourtsidou et al. 2013; Poulin et al. 2023, and references therien) or baryons and dark energy (Ferlito et al. 2022), or models assuming unstable dark matter particles (Enqvist et al. 2015, 2020; Murgia et al. 2017; Abellán et al. 2021; Chen et al. 2021; Choi & Yanagida 2022; Bucko et al. 2023).

The class of decaying dark matter models includes several different scenarios. In the simplest case, a fraction of dark matter decays into a relativistic component (often assumed to be dark radiation). However, such a model is strongly constrained by CMB observations (see e.g. Hubert et al. 2021; Simon et al. 2022; Bucko et al. 2023), as it affects the cosmic background evolution. An alternative and only somewhat more complex scenario assumes the initial dark matter particles decay into a pair made up of a massless and a massive particle, the latter obtaining a velocity kick during the decay process. This scenario is referred to as a two-body decaying dark matter model and is denoted as ‘ΛDDM’ hereafter.

A direct consequence of the ΛDDM model is the free streaming process of the stable decay products, which alters the gravitational collapse of cosmic structures. This effect is relevant at scales set by the free-streaming length of the stable daughter particles and thus by the magnitude of the velocity kicks (v_k) they receive as a consequence of energy-momentum conservation. As a result, the matter power spectrum is suppressed during late times at scales above k ∼ 0.1 h Mpc⁻¹ (Wang & Zentner 2012).

In addition to the aforementioned reason, there are arguments motivated by particle physics to consider models in which dark matter is not stable over cosmic time. First of all, such a stability condition does not emerge naturally; it usually requires additional assumptions such as Z₂ symmetry (Hambye 2011). Moreover, there are numerous theoretical models involving dark matter decays, such as sterile neutrinos (Abazajian et al. 2012; Adhikari et al. 2017), R-parity violation (Berezinsky et al. 1991; Kim & Kim 2002), or super weakly interacting massive particles (Covi et al. 1999; Feng et al. 2003).

The two-body decaying dark matter model (ΛDDM) has been studied from various angles over the past decade (e.g. by using perturbation theory; Wang & Zentner 2012; Abellán et al. 2021; Simon et al. 2022 or N-body simulations ranging from individual galaxies; Peter & Benson 2010 to a large-scale structure; Cheng et al. 2015). For example, Mau et al. (2022) obtained constraints on the ΛDDM model based on Milky Way satellite counts, while Wang et al. (2013) and Fuß & Garny (2023) used Lyman-α forest data to constrain the two-body decay rate in the regime of low-mass splittings. Additionally, Abellán et al. (2021) and Simon et al. (2022) considered Planck CMB observations together with supernova type Ia (SNIa) data and BAO to derive constraints on two-body decays. After including priors from WL observations, they reported a reduction of the S₈ tension for a best-fitting ΛDDM model with τ = 120 Gyr and v_k/c ≃ 1.2% (Simon et al. 2022).

In this work, we perform the first WL analysis of the ΛDDM model using cosmic shear data from the KiDS-1000 survey (Asgari et al. 2021). The non-linear clustering predictions were thereby modelled using an emulator based on a suite of N-body simulations. Along with the WL analysis, we also performed a re-analysis of the Planck 2018 CMB temperature and polarisation data as well as a combined WL plus CMB modelling investigating the potential of ΛDDM to solve the S₈ clustering tension.

Our paper is organised as follows: in Sect. 2, we describe basic physics and implications of the ΛDDM model, while Sect. 3 provides a detailed overview of our modelling of WL and CMB observables. In Sect. 4, we comment on choices made in relation to model inference, and in Sect. 5, we describe our results and compare them to recent studies before concluding in Sect. 6. In Appendix A, we compare the results of our N-body simulations to previous studies. Appendix B presents a discussion on the cosmology dependence of the ΛDDM effects, and Appendix C details the effects two-body decays and baryons have on WL and CMB observables. In Appendix D, we study more closely the tension between WL and CMB data in the context of the ΛDDM model, and finally, Appendix E provides more detailed information about the parameters we obtained from model inference.

2. Two-body decaying dark matter

In the ΛDDM model, an original (‘mother’) dark matter particle decays into a slightly lighter, stable particle and a relativistic massless relic, while the energy released during the decay process is split between two product species. We describe the basic theoretical properties of the ΛDDM model in Sect. 2.1. In Sect. 2.2, we discuss our ΛDDM N-body simulations, and finally, we introduce a new emulator of the ΛDDM non-linear matter power spectra in Sect. 2.3.

2.1. Theory

The ΛDDM model is simple enough to be described by two phenomenological parameters. The first parameter is the decay rate Γ controlling the frequency of the decay processes. The second parameter corresponds to the velocity kick v_k of the massive daughter particle, which is directly linked to the mass ratio between the mother and daughter particles. Some authors replace the velocity kick magnitude v_k with the ratio ε of the rest-mass energy

$$\begin{array}{c}\hfill \epsilon =\frac{1}{2}(1-\frac{m_{\mathrm{wdm}}^2}{M_{\mathrm{dcdm}}^2}),\end{array}$$(1)

where m_wdm and M_dcdm denote the rest mass of the warm daughter and the decaying cold mother dark matter particles, respectively (see e.g. Blackadder & Koushiappas 2014). The momentum of the daughter particle in the centre-of-momentum frame is $p_{\mathrm{wdm}}=m_{\mathrm{wdm}}c\epsilon /\sqrt{1-2\epsilon }$ (Abellán et al. 2021), where c is the speed of light. We note that in the non-relativistic limit, we obtained a simple relation v_k = εc.

The background evolution of the energy densities of both the cold (mother) and warm (daughter) dark matter species as well as the massless, dark radiation (daughter) component can be written as (Wang & Zentner 2012)

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

$$\begin{array}{cc}& {\dot{\rho }}_{\mathrm{wdm}}+3\mathcal{H}{\rho }_{\mathrm{wdm}}=a\mathrm{\Gamma }\frac{M_{\mathrm{dcdm}}^2+m_{\mathrm{wdm}}^2}{2M_{\mathrm{dcdm}}^2}{\rho }_{\mathrm{dcdm}},\hfill \end{array}$$(3)

$$\begin{array}{cc}& {\dot{\rho }}_{\mathrm{dr}}+4\mathcal{H}{\rho }_{\mathrm{dr}}=a\mathrm{\Gamma }\frac{M_{\mathrm{dcdm}}^2-m_{\mathrm{wdm}}^2}{2M_{\mathrm{dcdm}}^2}{\rho }_{\mathrm{dcdm}},\hfill \end{array}$$(4)

or using ε instead of particle masses, as

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

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

$$\begin{array}{cc}& {\dot{\rho }}_{\mathrm{dr}}+4\mathcal{H}{\rho }_{\mathrm{dr}}=a\mathrm{\Gamma }\epsilon {\rho }_{\mathrm{dcdm}}.\hfill \end{array}$$(7)

In the above equations, ℋ is the conformal Hubble parameter, and ρ_i is the energy density of species i. The subscripts ‘dcdm’, ‘wdm’, and ‘dr’ refer to the cold, warm, and massless species. Dots denote derivatives with respect to conformal time, and a stands for the scale factor.

Next to the two model parameters Γ and v_k, we include a third parameter

$$\begin{array}{c}\hfill f=\frac{{\mathrm{\Omega }}_{\mathrm{dcdm}}}{{\mathrm{\Omega }}_{\mathrm{dcdm}}+{\mathrm{\Omega }}_{\mathrm{cdm}}},\end{array}$$(8)

that allows for a scenario where only a faction f of the total initial dark matter fluid is unstable (while the remaining dark matter corresponds to a stable CDM particle). Here, we have introduced the abundances of the stable (Ω_cdm) and unstable (Ω_dcdm) dark matter species, respectively. With the above description, one can in principle study the full parameter space of two-body decays by taking arbitrary Γ, ε, and f values, with limiting cases Γ → 0, ε → 0 or f → 0 approaching the ΛCDM cosmology and ε → 1/2 approaching one-body decays. However, as recent studies have demonstrated (Wang & Zentner 2012; Mau et al. 2022, and references therein), ΛDDM models with a very large decay rate and velocity kicks are ruled out by observations, as they lead to a strong power suppression at k < 0.1 h Mpc⁻¹. We therefore focus on the regime with late-time decays (Γ ≲ H₀) and non-relativistic velocity kicks (v_k ≪ c) throughout this paper.

2.2. Simulations

Considering only non-relativistic decays, we implemented the ΛDDM model into the N-body code PKDGRAV3 (Potter et al. 2017), a tree-based gravity solver based on fast multi-pole expansion and adaptive time stepping. Following the theoretical descriptions (2)–(4), we found that at first order, the background equations remain unmodified. With this approximation, there is no energy transfer between the radiation and dark matter components caused by the decay process. Therefore, we kept the background cosmology implementation of PKDGRAV3 unchanged and implemented only the non-relativistic velocity kicks received by the WDM particles. We note that this differs from the one-body DDM model studied in Bucko et al. (2023), where the background evolution had to be modified.

The two-body decays are implemented into PKDGRAV3 via a function pkdDecayMass, which is revisited at each global integration time step (separated by a time interval Δt). The decay probability of a given (not yet decayed) particle at time step i is P = ΓΔt. Thus, a number of $\mathrm{\Delta }{N}_{\mathrm{w}\mathrm{d}\mathrm{m}}^i=\mathrm{\Gamma }\mathrm{\Delta }t{N}_{\mathrm{d}\mathrm{c}\mathrm{d}\mathrm{m}}^i$ simulation particles undergo the decay process, where $N_{\mathrm{d}\mathrm{c}\mathrm{d}\mathrm{m}}^i$ denotes the number of unstable DCDM particles at time step i. The particles that are about to decay are chosen randomly from all CDM particles. Immediately after the decay, they obtain a uniform velocity kick in a random direction. Importantly, these particles are flagged and added to a set of already decayed particles that are excluded from the decay process occurring in future time steps.

In Fig. 1, we plot the ratios of simulated ΛDDM to ΛCDM power spectra for varying values of f, v_k, and 1/Γ (see panel a, b, and c). In general, the two-body decaying dark matter model leads to a suppression of the total matter power spectrum at small scales and leaves the large scales unchanged. The amplitude of the suppression, as well as the scale of the downturn, depends on the values of the ΛDDM parameters. The fraction of decaying dark matter as well as the decay rate both affect the amplitude of the suppression, while the value of the velocity kick primarily influences the position of the downturn along the k-axis. The latter can be understood by the fact that larger streaming velocities are able to affect the formation of structures at larger scales.

Fig. 1. Effects of two-body decays on the non-linear matter spectrum from N-body simulations. The fiducial (ΛCDM) power spectrum corresponds to a constant line equal to one, while the power suppressions for varying f, vk, Γ, and z are shown in the panels a, b, c, and d.

The redshift dependence of the power suppression is shown in panel d of Fig. 1. Not surprisingly, the amplitude of the suppression increases towards lower redshifts. This behaviour is caused by the fact that more particles decay with time, causing a reduction in the clustering process compared to ΛCDM.

We ran a suite of N-body simulations for decay rates Γ < 1/13.5 Gyr⁻¹ and velocity kicks v_k/c < 0.02. All of our simulations were run assuming a fiducial cosmology with parameters h₀ = 0.678,  Ω_m, 0 = 0.307,  Ω_Λ, 0 = 0.693,  Ω_b, 0 = 0.048,  n_s = 0.966, and σ₈ = 0.883. We obtained converged results at the scales k ∼ 0.01 − 10 h Mpc⁻¹ for box sizes of L_box = 125, 250, 512 Mpc/h and particle numbers of N = 256³, 512³, 1024³ (depending on the specific ΛDDM configuration). We compared the output of our simulations to results from previous works and found a good level of agreement (see Appendix A). As we are primarily interested in the ratio of the non-linear power spectra between the ΛDDM and ΛCDM models, one can expect that the cosmology dependence is factored out to a large extent. We tested to see whether verified that the impact of the cosmology is much smaller than the suppression due to the two-body decays (see Appendix B).

2.3. Emulating the impact of dark matter decays

In order to carry out a Bayesian inference analysis (Sect. 4), we needed a fast modelling framework to explore the vast parameter space of astrophysics, cosmology, and dark matter models. As N-body simulations are not fast enough for this purpose, we built an emulator to account for the different ΛDDM parameters. The basic characteristics of our emulator-building procedure are shown in the flowchart of Fig. 2. First, we ran ∼100 gravity-only simulations for different dark matter parameters (Γ, v_k, and f) and measured the non-linear matter power spectra up to k ∼ 6 h Mpc⁻¹ between z = 2.35 and z = 0. In the next step, we performed a principal component analysis (PCA, see e.g. Aurélien 2019) on the ratios of ΛDDM and ΛCDM matter power spectra ${\mathcal{S}}_{\mathrm{\Lambda }\mathrm{DDM}}^{\mathrm{\Gamma },v_k,f}(k,z)$. We found that five PCA components are sufficient to describe the ratio of spectra with a reconstruction error of approximately 0.1%.

Fig. 2. Flowchart describing the emulation process (from left to right). The N-body simulations first need to be run to provide the matter density distribution for the given DDM parameters and redshift (f, vk, Γ, z). Next, we extracted power spectra and calculated the resulting power suppression PDDM/PΛCDM. Additionally, we performed PCA on PDDM/PΛCDM data, and using five principal components $(c_1^{\mathrm{p}\mathrm{c}\mathrm{a}}, c_2^{\mathrm{p}\mathrm{c}\mathrm{a}}, c_3^{\mathrm{p}\mathrm{c}\mathrm{a}}, c_4^{\mathrm{p}\mathrm{c}\mathrm{a}}, c_5^{\mathrm{p}\mathrm{c}\mathrm{a}})$, we recovered the original power suppression. Once the training data set was prepared, we trained the SIREN-like neural network. We used (f, vk, Γ, z) as the network input and output five PCA components. The PCA components obtained were used to reconstruct the power spectrum suppressions, which were then compared to the input ones. During the training process, the MSE loss between the input and output power suppression curves was minimised.

Next, we trained a neural network to model these five PCA components of the simulated power spectra ratios for a given parameter vector (Γ, v_k, f, z). The network output was then transformed back from the PCA representation to the power spectra ratios before being compared to the original (simulated) ratios used for the training and testing of the emulator. During the network training process, we minimised the differences between the predicted and true power spectra ratios in the training set. We considered the mean squared error (MSE) metric to quantify the differences.

We chose the architecture of the sinusoidal representation networks (SIRENs; Sitzmann et al. 2020) for building our ΛDDM emulator. The SIRENs have been successfully shown to have good interpolation and signal reconstruction properties (Sitzmann et al. 2020). The main difference compared to standard feed-forward architectures involves replacing the commonly assumed ReLU activation function with a sine function. We used the architecture with two hidden layers that each had 1024 neurons to perform the emulation task. During the training, we optimised the MSE with the Adam optimiser (Kingma & Ba 2014). To further fine-tune the SIREN architecture, we performed a hyperparameter optimisation for the network’s learning rate l_r and regularisation strength λ using the Bayesian Optimization and Hyperband (BOHB) method (Falkner et al. 2018). The entire training was performed using the PyTorch (Paszke et al. 2019) machine learning framework.

After training, we tested the emulator on separate data (i.e. the test set) and monitored the prediction mismatch in each of the 30 k-bins. We show the emulation performance in Fig. 3 and demonstrate that both 1σ and 2σ errors stay below 1%. Our emulator thus efficiently predicts the response of decays on the non-linear matter power spectrum ${\mathcal{S}}_{\mathrm{\Lambda }\mathrm{DDM}}^{\mathrm{\Gamma },v_k,f}(k,z)$. The prediction time of our emulator is a few milliseconds. With this tool, we could model the non-linear power spectrum in the presence of dark matter decays as

$$\begin{array}{c}\hfill P_{\mathrm{\Lambda }\mathrm{D}DM}^{\mathrm{nonlin}}(k,z)={\mathcal{S}}_{\mathrm{\Lambda }\mathrm{D}DM}^{\mathrm{\Gamma },v_k,f}(k,z)\times P_{\mathrm{\Lambda }\mathrm{CDM}}^{\mathrm{nonlin}}(k,z),\end{array}$$(9)

Fig. 3. Performance of the ΛDDM emulator at different scales. The precision of the two-body decaying dark matter emulation procedure was assessed by testing the emulator on the data that was not used during training. The average prediction error was calculated for each of the 30 k-bins, and the 1σ (grey) and 2σ (golden) error bars of the prediction are highlighted.

where $P_{\mathrm{\Lambda }\mathrm{CDM}}^{\mathrm{nonlin}}(k,z)$ is obtained by multiplying the non-linear power spectrum from the revised_halofit method (Takahashi et al. 2012). Our emulator to model the non-linear ΛDDM power spectrum is published as part of the publicly available code DMemu¹.

3. Data sets and modelling framework

Here we first describe the observational data that we use to study the ΛDDM model. Later, in Sects. 3.2 and 3.3, we present our modelling framework.

3.1. Data sets

Galaxy WL is a particularly promising observable to probe decaying dark matter models, as such scenarios tend to affect structure formation at late times and small cosmological scales. However, while primarily focusing on the WL analysis, we also include CMB data in our analysis.

In our study, we used the WL cosmic shear data of the KiDS-1000 data release, obtained in five redshift bins between 0.1 ≲ z ≲ 1.2 (Asgari et al. 2021). We used the band power angular spectra measured at scales 118 ≤ l ≤ 1266. Additionally, we included the Planck 2018 high-ℓ power spectra (ℓ ≥ 30) of temperature (TT), polarisation (EE), and their cross (TE) obtained from Planck Collaboration V (2020). In the following subsections, we discuss how we predicted computed these observational quantities in our modelling framework.

3.2. Cosmic shear modelling for KiDS-1000

The KiDS-1000 catalogue provides information about the shear of over 20 million galaxies, divided into five tomographic bins between z ∼ 0.1 and z ∼ 1.2. The catalogue provides the basis of the auto- and cross-correlation band power of all five redshift bins. With eight data points for each spectrum, the KiDS-1000 band power contains a total of 120 observational data points with correlated errors from the corresponding covariance matrix.

3.2.1. Non-linear matter power spectrum

We used the public version of the Boltzmann Solver CLASS (Blas et al. 2011) to calculate the ΛCDM matter power spectrum for any given set of cosmological parameters. For the non-linear modelling, we relied on the revised halofit method (Takahashi et al. 2012) implemented in CLASS. Following Planck Collaboration VI (2020), we assumed a single massive neutrino species with a fixed mass m_ν = 0.06 eV throughout this work.

The process of baryonic feedback causes gas to be expelled out of galaxies and clusters, leading to a suppression of the matter power spectrum at small cosmological scales (e.g. van Daalen et al. 2011; Schneider & Teyssier 2015; Chisari et al. 2019). We note that this suppression is similar in shape to the one caused by dark matter decays (Hubert et al. 2021), making it particularly important to include baryonic feedback in our modelling pipeline. We used BCemu² (Giri & Schneider 2021), an emulation tool providing the suppression 𝒮_b(k, z) due to baryonic feedback described in Schneider et al. (2019). The term 𝒮_b is a function of seven baryonic parameters and one cosmological parameter, namely, the ratio of the baryonic and total matter abundance. However, we used the reduced three-parameter model presented in Giri & Schneider (2021), where four parameters are fixed based on the results from hydrodynamical simulations. The three-parameter model consists of two parameters describing the gas distribution (log₁₀M_c, θ_ej) and one parameter describing the stellar mass (η_δ) around galaxies and clusters. We refer to Giri & Schneider (2021) for detailed description of the model.

Finally, the response function from the two-body dark matter decays was multiplied with the non-linear, baryonified power spectrum, as shown in Eq. (9). It The former was obtained from the emulator described in Sect. 2.3. This modular method assumes that all dependence on cosmological parameters is captured by the power spectrum from the revised halofit (see Appendix B). At the same time, the responses from the baryonification and the two-body decays remain independent of cosmology. Regarding the baryonification method, this assumption has been validated before (Schneider et al. 2020a,b). In Appendix B, we validated this assumption for the ΛDDM model.

3.2.2. Intrinsic alignment

The effects from intrinsic galaxy alignments were modelled assuming the non-linear alignment (NLA) model as described in Hildebrandt et al. (2016) and first published in Hirata & Seljak (2004). The intrinsic alignment enters the band power modelling via a window function (see Eq. (9) in Bucko et al. 2023) that goes into galaxy-intrinsic and intrinsic-intrinsic terms of the cosmic shear angular power spectrum. Among the two intrinsic alignment parameters A_IA and η_IA, we fixed the latter one to zero following the approach used in the standard KiDS-1000 analysis (Asgari et al. 2021).

3.2.3. Angular power spectrum

We investigated the impact of decaying dark matter on structure formation through the analysis of the cosmic shear angular power spectrum, including both auto-correlation and cross-correlation power spectra between different galaxy populations (tomographic bins). Our modelling approach follows the methodology presented in Bucko et al. (2023), and we refer the reader to that work for a more comprehensive description.

We used the multi-purpose cosmology calculation tool PyCosmo (Tarsitano et al. 2021) to model the cosmic shear angular power spectrum. The angular power was then converted into band powers following Joachimi et al. (2021). In Fig. C.1, we show all band power coefficients together with the respective error bars. The auto- and cross-band power measurements are illustrated in five KiDS-1000 redshift bins, always for multipoles between l ≃ 118 to l ≃ 1266.

3.3. CMB modelling

The Boltzmann Solver CLASS can also be used for theoretical modelling of the CMB temperature and polarisation data. Since the work of Abellán et al. (2021), CLASS comes with an implementation of the ΛDDM model in which the authors introduced a fluid approximation that reduces the computational costs significantly, and we use this model throughout our work³. To investigate the effects of ΛDDM cosmology on CMB data, we modelled the temperature and polarisation power spectra from the Planck 2018 data release (Planck Collaboration V 2020). We adopted the same methodology as the one introduced in Bucko et al. (2023) and refer the reader to that work for more details. We note, however, that our pipeline is tested for the ΛCDM cosmology reaching an excellent agreement with the results from Planck Collaboration VI (2020). A comparison can be found in Appendix A of Bucko et al. (2023).

4. Model inference

We performed a number of Markov chain Monte Carlo (MCMC) samplings in order to infer the posterior probability distribution of cosmological, baryonic, intrinsic alignment, and two-body DDM parameters based on the WL data from KiDS and the CMB observations from Planck. We employed the stretch-move ensemble method implemented within the emcee package (Foreman-Mackey et al. 2013) to sample from the posterior distribution.

An overview of the sampled parameters and the prior choices are listed in Table 1. We used flat priors for all cosmological parameters except for the optical depth τ, where we assumed a Gaussian prior 𝒩(0.0506, 0.0086), as explained in Sect. 3.1 of Bucko et al. (2023). For the dark matter abundance (ω_dm) and the amplitude of the primordial power spectrum (A_s), we used priors that are wide enough to comfortably include the WL posteriors found by Asgari et al. (2021). We note that ω_dm represents the initial dark matter abundance, and as we assume only late-time decays, this parameter describes the total initial dark matter budget in both the ΛCDM and ΛDDM scenarios. For the baryon abundance ω_b, Hubble constant h₀, and spectral index n_s, which cannot be well constrained by WL alone, we choose the prior that is as wide as possible to span the values found by surveys (e.g. CMB data) that are more sensitive to these parameters. The prior range for the intrinsic alignment parameter A_IA is wide enough to include the posterior distribution of this parameter found by Asgari et al. (2021). The Planck absolute calibration A_Planck was probed under the Gaussian prior 𝒩(1.0, 0.0025)⁴. We further imposed flat priors on the baryonic parameters log₁₀M_c, θ_ej, and η_δ covering the full range of the BCemu parameters. For the ΛDDM decay rate (Γ) and the velocity kick magnitude (v_k), we assumed flat priors for both the log₁₀Γ and log₁₀v_k spanning from ΛCDM values up to the upper boundary limited by the range of the ΛDDM emulator (see Table 1).

In both the WL and CMB setups, we assumed Gaussian likelihoods. In the case of the CMB, we used a marginalised light-weight version of the full Planck likelihood called Plik_lite (Prince & Dunkley 2019; Planck Collaboration V 2020), as it is an affordable approximation in the case of ΛDDM. Abellán et al. (2021) demonstrated that using Plik_lite for ΛDDM produces only a negligible difference on the recovered posteriors compared to a full Planck analysis. Also, Plik_lite comes with the marginalised version of the covariance matrix, which we used throughout this work. In the case of WL, we used the band power covariance matrix as published in Asgari et al. (2021). To assess the convergence of our chains, we applied the Gelman-Rubin criterion (Gelman & Rubin 1992), assuming the chains to be converged at R_c < 1.1.

We used two different metrics to assess the tension between the CMB and WL observations within the assumed cosmological model. The first one is the standard Gaussian metric (see e.g. Asgari et al. 2021), which is defined as

$$\begin{array}{c}\hfill {\tau }_{S_8}=\frac{S_8^{\mathrm{CMB}}-S_8^{\mathrm{WL}}}{\sqrt{\mathrm{Var}\left[{\mathrm{S}}_8^{\mathrm{CMB}}\right]+\mathrm{Var}\left[{\mathrm{S}}_8^{\mathrm{WL}}\right]}},\end{array}$$(10)

where Var [] stands for the variance of the S₈ measurements from the two datasets. We note that the Gaussian metric assumes a Gaussian posterior distribution of the parameter of interest and, importantly, is agnostic about how well the underlying data are being fitted. Due to these shortcomings, we also employed the difference in maximum a posteriori Q_DMAP criterion (Raveri & Hu 2019). This metric is defined as

$$\begin{array}{c}\hfill Q_{\mathrm{DMAP}}={\chi }_{\min ,\mathrm{WL}+\mathrm{CMB}}^2-({\chi }_{\min ,\mathrm{WL}}^2+{\chi }_{\min ,\mathrm{CMB}}^2),\end{array}$$(11)

where ${\chi }_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{W}\mathrm{L}}^2$, ${\chi }_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{C}\mathrm{M}\mathrm{B}}^2$, and ${\chi }_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{W}\mathrm{L}+\mathrm{C}\mathrm{M}\mathrm{B}}^2$ correspond to the minimum chi-squared values from the WL, the CMB, and the combined analysis. The Q_DMAP criterion evaluates the incapability of a combined analysis to approach the goodness of fit of the individual analyses, that is, when either of the datasets is fitted separately. In terms of σ, the tension is quantified as $\sqrt{Q_{\mathrm{DMAP}}}$.

Although the emulator presented in Sect. 2.3 accounts for the three parameters Γ, v_k, and f, we fixed the fraction of decaying dark matter to stable dark matter to unity (f = 1). This means that we restricted our analysis to the case of a universe with one initial (unstable) dark matter fluid, leaving the investigation of a multi-fluid dark matter sector to future work.

5. Results

In the first part of this section, we describe the constraints on two-body decays from WL and CMB observations and report our findings regarding baryonic physics. In the second part, we revisit the impact of the ΛDDM model on the S₈ tension before concluding with a discussion on how well the observational data can be fitted within the ΛDDM model.

5.1. Derived constraints on two-body decays

We obtained the constraints on the decay rate Γ and velocity kick magnitude v_k by marginalising over cosmology, baryons, and nuisance parameters. We display the outcome in Fig. 4, showing the ΛDDM posteriors at the 95% credible intervals as obtained from the WL (green) and CMB (orange) analyses. In both cases, the obtained constraints are in agreement with ΛCDM, showing no hint of decay in the dark matter sector. Figure 4 also includes a comparison of our findings with several recent studies: Mau et al. (2022; black) focusing on Milky Way satellites, Fuß & Garny (2023; light blue) examining the impact on Lyman-α forest, and Abellán et al. (2021) and Simon et al. (2022; pink and blue, respectively) using Planck 2018, SNIa, and BAO data. We note that our CMB posteriors are consistent with previous CMB studies, in particular with the results of Simon et al. (2022). Our limits are slightly stronger than those of Abellán et al. (2021) but exhibit similar hyperbolic contour trends. Milky Way satellites as probed by the DES collaboration (Mau et al. 2022) are sensitive to decay rate Γ and rule out half-life times of τ < 30 Gyr, while the Lyman-α appears to be less sensitive to two-body decays than Milky Way satellites. Finally, the WL data alone provides the most stringent constraints, excluding regions where the following conditions hold at 95% credible intervals:

$$\begin{array}{cc}& \tau ={\mathrm{\Gamma }}^{-1}\lesssim 125\mathrm{G}yr,and\hfill \end{array}$$(12)

$$\begin{array}{cc}& v_k\gtrsim 300\mathrm{k}m{s}^{-1}.\hfill \end{array}$$(13)

Fig. 4. Two-dimensional constraints on ΛDDM parameters Γ and vk as obtained from WL cosmic shear (green) and CMB (orange) data. We also display the results from several recent studies constraining ΛDDM parameters using different types of observables (see Sect. 5.1 for more details).

Hence, within the parameter space of ΛDDM and the datasets analysed in this study, WL data impose much stronger limits compared to CMB observations⁵. This is mainly due to the fact that, unlike for the case of one-body decays (Bucko et al. 2023), the two-body ΛDDM model does not significantly affect the background evolution of the Universe, leaving the signal from the late-time integrated Sachs-Wolfe effect unchanged. More discussions about the effects of two-body decays on the CMB signal can be found in Appendix C.

5.2. Revisiting the S₈ tension between WL and the CMB

As a next step, we investigated the extent the ΛDDM model is able to alleviate the S₈ tension between lensing and CMB data. In contrast to previous findings (e.g. Abellán et al. 2021, 2022; Tanimura et al. 2023), we argue that the ΛDDM model is unable to significantly reduce the systematic shift between the clustering signal predicted from the CMB and the KiDS-1000 WL survey.

In Fig. 5, we illustrate the Ω_m − S₈ posterior contours as obtained from the individual WL and CMB data analyses for both the ΛCDM (dark blue and black) and the ΛDDM model (green and orange). We observed only a marginal change in the WL and CMB contours when going from the ΛCDM to the ΛDDM case; both are at the 68% and 95% confidence level. However, at the 99% confidence level, the CMB contours reveal a prominent feature pointing towards low S₈-values in apparent agreement with the WL analysis.

Fig. 5. Two-dimensional posterior distributions of Ωm and S8 as resulting from MCMC analysis highlighting the 68%, 95%, and 99% confidence intervals. The WL results are shown in dark blue (ΛCDM) and green (ΛDDM), while the findings based on CMB data are displayed as black (ΛCDM) and orange (ΛDDM) contours.

Applying the Q_DMAP criterion defined in Eq. (11), we obtained a mutual tension of 3.4σ between WL and CMB data for the case of ΛCDM. When including two-body decays, the tension is reduced to 1.9σ. Using the Gaussian metric defined in Eq. (10), we obtained an S₈ tension of 3.0σ for the ΛCDM and 2.7σ for the ΛDDM. The difference between the two measures confirms the non-Gaussian shape of the CMB posteriors.

The origin of the non-Gaussian tail in the Ω − S₈ plane (visible in Fig. 5) can be better understood by looking at Fig. 6, in which we plot the posterior distribution of S₈ against the product of the ΛDDM parameters (Γ × v_k). One can see that for large values of Γ × v_k, the posterior obtained from the CMB analysis bends downwards, preferring smaller values of S₈. However, this downturn happens at a part of the parameter space that is not favoured by the WL analysis. Instead of fully overlapping with the WL posteriors (green), the CMB contours (orange) rather wrap around them only, leading to a slight overlap of the 95% confidence regions.

Fig. 6. Effects of two-body decays on S8 in the case of the WL (green contours), CMB (orange contours), and combined WL plus CMB data (purple contours). We show the 2D posterior contours of Ωm against the product of the DDM parameters at the 68% and 95% confidence level. We also display the best-fit model in each scenario as coloured stars.

Based on the analysis above, we concluded that although the two-body ΛDDM alleviates the S₈ tension to some degree, it does not provide a convincing solution to it. Furthermore, the model does not yield improved fits to either the individual CMB or WL data (see Table D.1) despite the two additional model parameters Γ and v_k. For further tests of the ΛDDM model, we refer to Appendix D.

Fitting the combined WL and CMB data, we obtained a preference for a model with non-zero DDM parameters. The best fitting values are ${log}_{10}\mathrm{\Gamma }=-2.{25}_{-0.23}^{+0.74}$ and log₁₀v_k > 2.80, compatible with the overlap region of CMB and WL posteriors. See Appendix E for the full results from our MCMC analysis. These findings align with the conclusions of Simon et al. (2022), who added a Gaussian prior with the S₈ value from KiDS to their CMB analysis. However, given the original tension between WL and CMB data, it is not surprising to obtain a preference for non-vanishing ΛDDM parameters in the combined data analysis. We want to stress that this is by no means a signal of a departure from ΛCDM but rather a natural consequence of the internal tension between the two datasets. This interpretation is confirmed by the fact that the posteriors from the individual analysis shown in Fig. 5 do not show significant overlap in the ΛCDM nor in the ΛDDM case.

The constraints inferred for all sampled parameters in our MCMC runs in the cases of CMB-only, WL-only, and the combined analysis can be found in Tables E.1 and E.2. In particular, the values of S₈ from our CMB and WL analyses are compared to the results of Asgari et al. (2021), Schneider et al. (2022) and Planck Collaboration VI (2020) in Fig. 7. Details about tests of our MCMC pipeline and the related discussion can be found in Appendix A of Bucko et al. (2023).

Fig. 7. One-dimensional constraints of the S8 parameter as inferred from the CMB, WL, and combined analyses. The original results from the KiDS-1000 (Asgari et al. 2021) and the Planck 2018 (Planck Collaboration VI 2020) analyses are included in grey and blue for comparison. We also show the results of Schneider et al. (2022) using a similar WL modelling recipe.

6. Summary and conclusion

A dark matter scenario including particle decay forms a natural extension to the minimal model of a cold, stable, and collisionless dark matter particle. In this work, we investigated the case of two-body dark matter decays where particles decay into a massless particle and a massive daughter particle, the latter obtaining a velocity kick as a result of the decay process. This model has been studied in different contexts in the past (e.g. Peter et al. 2010; Wang & Zentner 2012; Cheng et al. 2015; Fuß & Garny 2023) and has been proposed as a potential solution to the S₈ tension (Abellán et al. 2021, 2022; Simon et al. 2022).

In this paper, we performed the first fully non-linear analysis of WL and CMB data for the two-body decaying dark matter (ΛDDM) scenario. Based on a suite of N-body simulations, we constructed a neural network-based emulator to obtain fast predictions of non-linear matter power spectra for arbitrary values of the decay rate (Γ), the decay-induced velocity kick (v_k), and the fraction of decaying to total dark matter (f). We then included the emulator in our pipeline predicting WL observations and performed an MCMC analysis with WL data from KiDS-1000 and CMB data from Planck 2018.

We present improved constraints on the two-body decaying dark matter parameters based on the WL data from KiDS-1000. Our constraints are significantly stronger compared to previous results. Specifically, we exclude models with τ = Γ⁻¹ ≲ 125 Gyr and v_k ≳ 300 km s⁻¹. Figure 4 provides a summary of our constraints from the WL and the CMB, along with previous results from the literature.

When considering the clustering (or S₈) tension between KiDS-1000 and Planck 2018, we observed an improvement of 1.5σ with ΛDDM as compared to the original 3.4σ tension measured in the ΛCDM model using Q_DMAP criterion, which can account for the non-Gaussianity of the posteriors compared to Gaussian tension better. However, the WL and CMB posteriors overlap at a very small confidence, and therefore, we conclude that the two-body ΛDDM scenario is unable to convincingly resolve the clustering tension between WL and CMB observations. We note that previous works obtaining different conclusions (e.g. Abellán et al. 2021, 2022; Simon et al. 2022; Tanimura et al. 2023) did not include a full, self-consistent modelling of the WL signal and were therefore unable to directly test the S₈ tension in the case of ΛDDM.

A further step forward with respect to the current analysis could involve the analysis of the ΛDDM model using DES observations along with the possible addition of galaxy clustering and galaxy-galaxy lensing. Furthermore, including additional low-redshift datasets, such as eBOSS and SNIa, could lead to stronger constraints on parameters such as h₀ and Ω_b, allowing for a more precise determination of the dark matter parameters. Finally, we expect data from Euclid and the Vera C. Rubin Observatory to significantly improve current limits on dark matter decays.

The emulator of two-body decays developed in this study is now publicly available as the DMemu Python package. We welcome researchers to incorporate this package into their data analysis pipelines and further test the ΛDDM model.

Acknowledgments

We thank Douglas Potter and Jonathan Hubert for helpful discussions and technical support. This work is supported by the Swiss National Science Foundation under grant number PCEFP2_181157. Nordita is supported in part by NordForsk.

References

Appendix A: Comparing two-body decay simulations to previous work

In Fig. A.1, we present a comparison of our N-body simulations with three other studies that have investigated the same ΛDDM model (Wang & Zentner 2012; Cheng et al. 2015; Abellán et al. 2021). In the top panel of Fig.A.1, we compare linear recipes based on solving the Boltzmann hierarchy, as developed by Wang & Zentner (2012) (dashed lines) and Abellán et al. (2021) (dotted lines). The solid lines represent the results of our ΛDDM N-body simulations for f = 1.0. For the chosen values of v_k and Γ, we observed good agreement in terms of the downturn scales. However, at non-linear scales, there are more pronounced differences between the results, which is expected due to the limitations of linear calculations at higher values of k. We note that for v_k = 30000 km/s, which corresponds to 10% of the speed of light, we have concerns about the accuracy of our N-body implementation. Therefore, we did not perform model inferences for such extreme values of v_k in this work.

Fig. A.1. Comparison of our N-body simulations to previous works studying the ΛDDM model. The top panel shows how our simulations (solid lines) compare to linear (Boltzmann) codes of Wang & Zentner (2012) (dashed lines) and Abellán et al. (2021) (dotted lines) for three different decaying dark matter configurations. The bottom panel shows simulations of two-body decays published in Cheng et al. (2015) compared to our results for different decay rates Γ and velocity kick magnitudes vk.

In the bottom panel of Fig. A.1, we provide a benchmark of the ΛDDM model (for f = 1.0) by comparing it with the N-body study conducted by Cheng et al. (2015) for six different sets of dark matter parameters. We found that scenarios with smaller velocity kicks (v_k ≤ 500 km/s) exhibit agreement at the percent level across all scales. At the non-linear regime (k ≳ 1 h/Mpc) and for larger velocity kicks (v_k ≥ 1000 km/s), we observed larger deviations in the predicted power suppression. These deviations can be attributed to the different dark matter implementations employed by Cheng et al. (2015). However, it is important to note that scenarios with such a significant power suppression in the non-linear regime are ruled out by observations, which favour a much weaker decrease of power at k ≈ 1 − 10 h/Mpc (Wang & Zentner 2012). Therefore, the observed differences between our results and those of Cheng et al. (2015) are not a cause for concern.

Appendix B: Cosmology dependence of suppression of the matter power spectrum in ΛDDM with respect to ΛCDM

The effects of ΛDDM in our study are accounted for by applying a cosmology-independent boost to the ΛCDM non-linear matter power spectrum, calculated using the revised_halofit method within the CLASS code. To demonstrate the validity of this approach, we conducted a test suite of N-body simulations where we kept the ΛDDM parameters fixed at Γ⁻¹ = 26.20 Gyr and v_k = 500 km/s and varied one ΛCDM parameter at a time.

We show in Fig. B.1 that this approach provides a good approximation, as neglecting the cosmology dependence in the applied boost results only in a second-order effect. In panels (a) to (e) of this figure, we show the impact of the Hubble parameter (h₀), clustering amplitude (σ₈), spectral index (n_s), matter abundance (Ω_m), and baryon abundance (Ω_b), respectively. We display the ΛDDM boost for the fiducial cosmology of our N-body simulations with a solid salmon line, while we plot the same quantity for the best-fit cosmology of the KiDS-450 survey (Hildebrandt et al. 2016) with a dashed blue line. We also added a few additional models, represented by the dotted grey lines.

Fig. B.1. Cosmology dependence of the two-body effects on the matter power spectrum ratio between the ΛDDM and ΛCDM models. From panels (a) to (e), respectively, we show the effects of Hubble parameters (h0), clustering amplitude (σ8), spectral index (ns), matter abundance (Ωm), and baryon abundance (Ωb). With a solid salmon line, we show the fiducial cosmology of our N-body simulations, while the cosmology of the KiDS-450 survey is shown with a dashed blue line. We have also added a few additional scenarios, represented by dotted grey lines.

In this analysis, we found that the choice of cosmology has a noticeable impact on the ΛDDM boost only for the parameters σ₈ (panel b) and Ω_m (panel d). For the remaining parameters, the choice of cosmology does not significantly affect the boost. We find it is important to note that the models exhibiting substantial differences assume cosmologies that are quite distinct from the fiducial one. However, even in the cases of σ₈ and Ω_m, the observed effects are relatively small compared to the overall amplitude of the observed boost.

Appendix C: Effect of model and baryonic feedback parameters on the observables

C.1. Weak lensing

To correctly interpret our results, it is important to understand the role of model parameters in relation to observables. Figure C.1 shows the original KiDS-1000 measurements of cosmic shear band powers with corresponding error bars (black points), and we illustrate how well the WL, CMB, and WL plus CMB scenarios can fit these data. The solid green lines show the best-fit configuration in the WL-only setup in the case of ΛDDM (however, it performs very similar to the ΛCDM case). Once we calculated the band powers resulting from the CMB-only best-fit cosmology (orange lines), we recovered a significantly stronger WL signal compared to the KiDS observations. This is the consequence of larger clustering (S₈ value) being preferred by the CMB data. In the combined WL plus CMB scenario, tighter error bars on the cosmological parameters from Planck result in a best-fit cosmology that is close to the one from the CMB-only scenario, but baryons and, namely, two-body decays on small scales provide a much better fit to the cosmic shear data compared to the CMB-only scenario, as seen from the solid purple lines in the figure. To also showcase that two-body decays have a significant influence on the shear signal, we display the dashed and dotted purple lines. The former represents the same cosmology and baryonic parameters as in the solid purple case but switching off the two-body decays completely, while the latter displays the largest possible impact of the two-body decays in the case of cosmology and baryons fixed to the WL plus CMB ΛDDM best-fit scenario (solid purple lines). This impact is large enough to completely over- or underestimate the measured shear signal, as can be seen from the autocorrelation spectra of higher-redshift bins, for example. This also explains why in the case of the combined ΛDDM analysis, we recovered the constraints on the two-body parameters that are detached from ΛCDM cosmology.

Fig. C.1. Effects of two-body decays on band power spectra illustrated using different forward-modelled ΛDDM configurations. Band power auto- and cross-correlations for five tomographic bins with respective error bars as measured by the KiDS-1000 survey are shown as black dots. Solid green (orange) lines represent the best-fit configurations of the WL (CMB) analysis projected on the band powers. The purple lines represent the different configurations of the WL plus CMB setup in ΛDDM cosmology. The solid lines show the best fit, while the dashed and dotted lines respectively illustrate the impact of underestimating and overestimating the two-body decays, keeping all the remaining parameters fixed to those of the WL plus CMB best fit. In the legend, we include the minimal χ2 value for every plotted configuration.

C.2. Baryon feedback

We used three free parameters in the baryonic effects emulator (log₁₀ M_c, θ_ej, η_δ) as proposed in Giri & Schneider (2021). From our analysis, we found that the WL data does not provide any information about the stellar population parameter η_δ. This is evident from Fig. C.2 and can be inferred from Tab. E.2, as varying η_δ does not affect the modelled WL observables. On the other hand, the gas profile parameters log₁₀M_c and θ_ej have an impact on the WL observables. For log₁₀M_c, we obtained log₁₀M_c < 13.1, (13.2) for ΛCDM (ΛDDM) from the WL-only analysis and log₁₀M_c > 13.8, (unconst) from the combined analysis. Similarly, for θ_ej, we found θ_ej < 5.45, (5.57) in the case of ΛCDM (ΛDDM) from the WL-only analysis and θ_ej > 5.88, (unconst) from the combined analysis. These results indicate that the values of these baryonic parameters are mutually exclusive when comparing the WL-only and combined (ΛCDM) scenarios. The underlying reasons for this behaviour are discussed in the following sections.

Fig. C.2. Effects of baryons on band power spectra illustrated using different baryonic parameter configurations. Band power auto- and cross-correlations for five tomographic bins with respective error bars as measured by the KiDS-1000 survey are shown as black dots. Solid dark blue (green) lines represent the best fit in the case of ΛCDM (ΛDDM), and the cyan, magenta, and orange lines depict the model predictions for different modifications of the ΛCDM best fit, varying baryonic parameters. In the legend, we include the minimum χ2 value for every plotted configuration.

Fig. C.2 demonstrates the impact of baryonic feedback on the KiDS band powers by varying the log M_c, θ_ej, and η_δ parameters within the BCemu framework. The best-fit configurations for ΛCDM (ΛDDM) are depicted in dark blue (green), while the coloured lines represent the variation of baryonic feedback strength in the ΛCDM best-fit case. Notably, we observe that modifying the stellar population parameter η_δ does not alter the shear signal. However, adjustments to the gas profile parameters log₁₀M_c and θ_ej affect the model predictions at scales l > 300.

Fig. C.1 and C.2 provide insights into how the baryonic and ΛDDM parameters impact the band powers. Both sets of parameters exhibit a similar qualitative effect, resulting in a suppression of the signal at small scales. However, the influence of baryons is relatively subtle compared to the effects of DDM, as indicated by the available priors. Also, we note that baryonic effects primarily manifest on the smallest scales, with no discernible impact below l ∼ 300. In contrast, ΛDDM can influence the signal on larger scales, particularly when considering scenarios with large velocity kicks.

The analysis thus reveals that when considering KiDS data alone, the preferred values for baryonic feedback parameters (and in the ΛDDM case, the two-body DDM parameters) tend towards the lower end. However, when combined with CMB data, these parameters drive the cosmology towards higher S₈ values, resulting in an excessive boost to the WL signal, as discussed earlier. The baryonic parameters can partially counteract this effect in the ΛCDM scenario by decreasing the signal on small scales, which explains the extreme values observed for the gas profile parameters in the combined analysis with CMB. Nevertheless, the limited suppression capability of baryons is inadequate to fully accommodate the WL signal, leading to deteriorated fits for both WL and CMB.

In contrast, in the ΛDDM scenario, the DDM parameters can more effectively suppress the band power signal, resulting in improved fits for both WL and CMB compared to the ΛCDM setup. This can be observed in Fig. C.3, where the left panel depicts the combination of DDM parameters with log₁₀M_c and the right panel illustrates their combination with θ_ej in the ΛDDM model. In both panels, it is evident that fitting the WL-only data favours a weak suppression of the power spectrum, indicated by the green contours located in the bottom-left region. Conversely, when combining WL and CMB, strong suppression is required either from baryons or from DDM, as indicated by the purple contours. It is conceivable that having even stronger baryonic effects at our disposal could entirely replace the need for the DDM and provide a satisfactory fit to both datasets individually, though the physical justification for such high baryonic feedback remains an open question.

Fig. C.3. Effects of baryons and two-body decays in the ΛDDM universe as obtained from probing the WL (green) and WL plus CMB (purple) scenarios. Both baryons and two-body decays affect the matter power spectrum in qualitatively similar ways, and we observed a preference for weak suppression in the WL setup while requiring much stronger impact in the combined WL plus CMB analysis.

C.3. Cosmic microwave background

In Fig. C.4, we show how two-body decays influence the predictions for CMB observables. We display temperature power spectra in the left panel and polarisation spectra in the right panel and add their cross-correlations in the middle panel. In the top row, we show the CMB spectra for ΛCDM (black) and ΛDDM (green) best fits, while we plot the predictions for the ΛCDM best-fit configuration extended by different combinations of Γ and v_k (dashed lines). In the smaller bottom panels, we show the difference between all these predictions with the ΛCDM best-fit cosmology. We observed that the two-body decays do not change the CMB signal significantly, meaning that the scatter of the CMB data is larger than the observable effects of two-body decays.

Fig. C.4. Planck TTTEEE angular power spectra as affected by different decaying dark matter parameter configurations. Black (green) solid lines represent the best fit in the case of ΛCDM (ΛDDM), and dashed curves depict the model predictions for different modifications of the ΛCDM best fit varying DDM parameters. In the legend, we include the χ2 value that each configuration yields.

Appendix D: Additional quantification of the tension

In the main text, we introduced two different metrics to estimate tension within a model at hand, namely, the Gaussian tension, defined in Eq. (10), and the difference in maximum a posteriori, following Eq. (11). The latter criterion is not subject to the assumption of Gaussian posteriors but cannot deal with overfitting (i.e. the number of additional degrees of freedom our new model possesses). The last criterion we used to assess the efficiency of the ΛDDM model is the change in the Akaike information criterion (AIC) defined as

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{AIC}=\mathrm{\Delta }{\chi }_{\min }^2+2(N_{\mathrm{\Lambda }\mathrm{D}DM}-N_{\mathrm{\Lambda }\mathrm{CDM}}),\end{array}$$(D.1)

where N_ℳ is a number of free parameters in model ℳ. To determine whether a new model is preferred over the ΛCDM, substantial evidence against it is required based on Jeffreys’ scale (Jeffreys 1961) using p < exp(−ΔAIC/2) (Schöneberg et al. 2022), where p = 10^0.5. This implies that the difference in AIC (ΔAIC) between the new model and ΛCDM should satisfy the condition

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{AIC}<\mathrm{\Delta }{\mathrm{AIC}}_0=-2.3.\end{array}$$(D.2)

From the above equation, we obtained ΔAIC = 3.9 and ΔAIC = 4.0 in the WL and CMB scenarios, respectively, and this adding of two free parameters did not lead to an improved

fit to the data. In the combined scenario, ΔAIC = −3.9 < ΔAIC₀. This means that the two additional parameters in the combined analysis are efficient. However, the combined ΛDDM fit is still worse by $\mathrm{\Delta }{\chi }_{\mathrm{m}\mathrm{i}\mathrm{n}}^2=3.8$ compared to what the ΛDDM scenario yields when treating the two datasets separately. We believe that fitting the data separately provides a better indicator of resolving the S₈ tension, in which the ΛDDM model fails. Throughout the work, we obtained the value of ${\chi }_{\mathrm{m}\mathrm{i}\mathrm{n}}^2$ as a minimum χ² among all sampled configurations in our MCMC chain. For a summary of combined MCMC analysis results, we refer to Tab. D.1.

Appendix E: MCMC results

In this section, we provide more details about the parameters resulting from our MCMC analysis. In Tab. E.1, we report the findings from our ΛCDM analysis (taken from Bucko et al. 2023), while in Tab. E.2, we summarise the actual ΛDDM model results. In both cases, we provide separate constraints from the WL, CMB, and combined analyses. In the upper part of the tables, we show the posteriors of the parameters sampled by in our MCMC runs, showing mean (best fit) values with corresponding upper and lower deviations. The dashed fields in the tables indicate the parameters were not present in a given MCMC setup, ’unconst’ is used when parameter constraints could not be obtained from MCMC analysis, and we state no uncertainties whenever referring to a value constant throughout the inference.⁶

All Tables

All Figures

	Fig. 1. Effects of two-body decays on the non-linear matter spectrum from N-body simulations. The fiducial (ΛCDM) power spectrum corresponds to a constant line equal to one, while the power suppressions for varying f, vk, Γ, and z are shown in the panels a, b, c, and d.
In the text

Fig. 2. Flowchart describing the emulation process (from left to right). The N-body simulations first need to be run to provide the matter density distribution for the given DDM parameters and redshift (f, vk, Γ, z). Next, we extracted power spectra and calculated the resulting power suppression PDDM/PΛCDM. Additionally, we performed PCA on PDDM/PΛCDM data, and using five principal components $(c_1^{\mathrm{p}\mathrm{c}\mathrm{a}}, c_2^{\mathrm{p}\mathrm{c}\mathrm{a}}, c_3^{\mathrm{p}\mathrm{c}\mathrm{a}}, c_4^{\mathrm{p}\mathrm{c}\mathrm{a}}, c_5^{\mathrm{p}\mathrm{c}\mathrm{a}})$, we recovered the original power suppression. Once the training data set was prepared, we trained the SIREN-like neural network. We used (f, vk, Γ, z) as the network input and output five PCA components. The PCA components obtained were used to reconstruct the power spectrum suppressions, which were then compared to the input ones. During the training process, the MSE loss between the input and output power suppression curves was minimised.

In the text

	Fig. 3. Performance of the ΛDDM emulator at different scales. The precision of the two-body decaying dark matter emulation procedure was assessed by testing the emulator on the data that was not used during training. The average prediction error was calculated for each of the 30 k-bins, and the 1σ (grey) and 2σ (golden) error bars of the prediction are highlighted.
In the text

	Fig. 4. Two-dimensional constraints on ΛDDM parameters Γ and vk as obtained from WL cosmic shear (green) and CMB (orange) data. We also display the results from several recent studies constraining ΛDDM parameters using different types of observables (see Sect. 5.1 for more details).
In the text

	Fig. 5. Two-dimensional posterior distributions of Ωm and S8 as resulting from MCMC analysis highlighting the 68%, 95%, and 99% confidence intervals. The WL results are shown in dark blue (ΛCDM) and green (ΛDDM), while the findings based on CMB data are displayed as black (ΛCDM) and orange (ΛDDM) contours.
In the text

	Fig. 6. Effects of two-body decays on S8 in the case of the WL (green contours), CMB (orange contours), and combined WL plus CMB data (purple contours). We show the 2D posterior contours of Ωm against the product of the DDM parameters at the 68% and 95% confidence level. We also display the best-fit model in each scenario as coloured stars.
In the text

	Fig. 7. One-dimensional constraints of the S8 parameter as inferred from the CMB, WL, and combined analyses. The original results from the KiDS-1000 (Asgari et al. 2021) and the Planck 2018 (Planck Collaboration VI 2020) analyses are included in grey and blue for comparison. We also show the results of Schneider et al. (2022) using a similar WL modelling recipe.
In the text

Fig. A.1. Comparison of our N-body simulations to previous works studying the ΛDDM model. The top panel shows how our simulations (solid lines) compare to linear (Boltzmann) codes of Wang & Zentner (2012) (dashed lines) and Abellán et al. (2021) (dotted lines) for three different decaying dark matter configurations. The bottom panel shows simulations of two-body decays published in Cheng et al. (2015) compared to our results for different decay rates Γ and velocity kick magnitudes vk.

In the text

Fig. B.1. Cosmology dependence of the two-body effects on the matter power spectrum ratio between the ΛDDM and ΛCDM models. From panels (a) to (e), respectively, we show the effects of Hubble parameters (h0), clustering amplitude (σ8), spectral index (ns), matter abundance (Ωm), and baryon abundance (Ωb). With a solid salmon line, we show the fiducial cosmology of our N-body simulations, while the cosmology of the KiDS-450 survey is shown with a dashed blue line. We have also added a few additional scenarios, represented by dotted grey lines.

In the text

Fig. C.1. Effects of two-body decays on band power spectra illustrated using different forward-modelled ΛDDM configurations. Band power auto- and cross-correlations for five tomographic bins with respective error bars as measured by the KiDS-1000 survey are shown as black dots. Solid green (orange) lines represent the best-fit configurations of the WL (CMB) analysis projected on the band powers. The purple lines represent the different configurations of the WL plus CMB setup in ΛDDM cosmology. The solid lines show the best fit, while the dashed and dotted lines respectively illustrate the impact of underestimating and overestimating the two-body decays, keeping all the remaining parameters fixed to those of the WL plus CMB best fit. In the legend, we include the minimal χ2 value for every plotted configuration.

In the text

Fig. C.2. Effects of baryons on band power spectra illustrated using different baryonic parameter configurations. Band power auto- and cross-correlations for five tomographic bins with respective error bars as measured by the KiDS-1000 survey are shown as black dots. Solid dark blue (green) lines represent the best fit in the case of ΛCDM (ΛDDM), and the cyan, magenta, and orange lines depict the model predictions for different modifications of the ΛCDM best fit, varying baryonic parameters. In the legend, we include the minimum χ2 value for every plotted configuration.

In the text

	Fig. C.3. Effects of baryons and two-body decays in the ΛDDM universe as obtained from probing the WL (green) and WL plus CMB (purple) scenarios. Both baryons and two-body decays affect the matter power spectrum in qualitatively similar ways, and we observed a preference for weak suppression in the WL setup while requiring much stronger impact in the combined WL plus CMB analysis.
In the text

	Fig. C.4. Planck TTTEEE angular power spectra as affected by different decaying dark matter parameter configurations. Black (green) solid lines represent the best fit in the case of ΛCDM (ΛDDM), and dashed curves depict the model predictions for different modifications of the ΛCDM best fit varying DDM parameters. In the legend, we include the χ2 value that each configuration yields.
In the text