Probing the two-body decaying dark matter scenario with weak lensing and the cosmic microwave background

Decaying dark matter (DDM) scenarios have recently regained attention due to their potential ability to resolve the well-known clustering (or $S_8$) tension between weak lensing (WL) and cosmic microwave background (CMB) measurements. In this paper, we investigate a well-established model where the original dark matter particle decays into a massless particle and a massive daughter particle. The latter obtains a velocity kick during the decay process that results in the suppression of the matter power spectrum at scales that are observable with WL shear observations. We perform the first fully non-linear WL analysis of this two-body decaying dark matter ($\Lambda$DDM) scenario, including intrinsic alignment and baryonic feedback processes. We used the cosmic shear band power spectra from \textit{KiDS-1000} data and combined it with temperature and polarisation data from \texttt{Planck} in order to constrain the $\Lambda$DDM model. We report new limits on the decay rate and mass splitting parameters that are significantly stronger than previous results, especially in the case of low-mass splittings. Regarding the $S_8$ tension, we found a reduction from about 3 to 2 $\sigma$, depending on which statistical measure is applied. We therefore conclude that the two-body $\Lambda$DDM model is able to reduce the $S_8$ tension without convincingly solving it.


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), largescale 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 = σ 8 √ Ω m /0.3, with σ 8 and Ω m describing the clustering amplitude and the total matter abundance.More specifically, CMB measurements from the Planck Collaboration have yielded S 8 = 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.016  −0.038 (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 8 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 8 tension can be achieved by suppressing the matter power spectrum at scales k ∼ 0.1-1 h Mpc −1 , which most substantially influence the clustering amplitude value S 8 (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(Enqvist et al. , 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 −1 (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 2 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 8 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 8 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.

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.

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 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 wdm = m wdm cε/ √ 1 − 2ε (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.
Next to the two model parameters Γ and v k , we include a third parameter 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

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 ( 5)-(7), 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 ∆N i wdm = Γ∆tN i dcdm simulation particles undergo the decay process, where N i dcdm 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.
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 −1 and velocity kicks v k /c < 0.02.All of our simulations were run assuming a fiducial cosmology with parameters h 0 = 0.678, Ω m,0 = 0.307, Ω Λ,0 = 0.693, Ω b,0 = 0.048, n s = 0.966, and σ 8 = 0.883.We obtained converged results at the scales k ∼ 0.01−10 h Mpc −1 for box sizes of L box = 125, 250, 512 Mpc/h and particle numbers of N = 256 3 , 512 3 , 1024 3 (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).

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 −1 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 S Γ,v k , f ΛDDM (k, z).We found that five PCA components are sufficient to describe the ratio of spectra with a reconstruction error of approximately 0.1%.A152, page 3 of 15 The N-body simulations first need to be run to provide the matter density distribution for the given DDM parameters and redshift ( f, v k , Γ, z).Next, we extracted power spectra and calculated the resulting power suppression P DDM /P ΛCDM .Additionally, we performed PCA on P DDM /P ΛCDM data, and using five principal components (c pca 1 , c pca 2 , c pca 3 , c pca 4 , c pca 5 ), we recovered the original power suppression.Once the training data set was prepared, we trained the SIREN-like neural network.We used ( f, v k , Γ, 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 S Γ,v k , f ΛDDM (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 where P nonlin ΛCDM (k, z) is obtained by multiplying the nonlinear 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 DMemu1 .

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.

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.

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.

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 BCemu2 (Giri & Schneider 2021), an emulation tool providing the suppression S b (k, z) due to baryonic feedback described in Schneider et al. (2019).The term S 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 10 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.

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

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 crosscorrelation 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.

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 work3 .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).Notes.Next to the parameter, we show the prior type used and the range for a flat prior or mean and standard deviation when a Gaussian prior was used.In the top part of the table, cosmological parameters are stated followed by the intrinsic alignment and Planck absolute calibration parameter, three baryonic parameters, and two free parameters characterising the decaying dark matter.

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 twobody DDM parameters based on the WL data from KiDS and the CMB observations from Planck.We employed the stretchmove 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 N(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 0 , 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 N(1.0, 0.0025) 4 .We further imposed flat priors on the baryonic parameters log 10 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 10 Γ and log 10 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 lightweight version of the full Planck likelihood called Plik_lite 4 https://wiki.cosmos.esa.int/planck-legacy-archive/index.php/CMB_spectrum_%26_Likelihood_Code(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 where Var [] stands for the variance of the S 8 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 where χ 2 min,WL , χ 2 min,CMB , and χ 2 min,WL+CMB 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 √ Q 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.

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 8 tension before concluding with a discussion on how well the observational data can be fitted within the ΛDDM model.

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)  Fig. 4. Two-dimensional constraints on ΛDDM parameters Γ and v k 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).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: v k 300 km s −1 . (13) Hence, within the parameter space of ΛDDM and the datasets analysed in this study, WL data impose much stronger limits compared to CMB observations5 .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.

Revisiting the S 8 tension between WL and the CMB
As a next step, we investigated the extent the ΛDDM model is able to alleviate the S 8 tension between lensing and CMB data.
In contrast to previous findings (e.g.Abellán et al. 2021Abellán et al. , 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 8 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 8 -values in apparent agreement with the WL analysis.
A152, page 7 of 15 Fig. 6.Effects of two-body decays on S 8 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.
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 8 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 8 plane (visible in Fig. 5) can be better understood by looking at Fig. 6, in which we plot the posterior distribution of S 8 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 8 .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.
Based on the analysis above, we concluded that although the two-body ΛDDM alleviates the S 8 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 Γ = −2.25 +0.74 −0.23 and log 10 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 8 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 nonvanishing Λ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 8 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).

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 8 tension (Abellán et al. 2021(Abellán et al. , 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 A152, page 8 of 15 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 τ = Γ −1 125 Gyr and v k 300 km s −1 .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 8 ) 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. 2021Abellán et al. , 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 8 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 0 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.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 10 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 8 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

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 bestfit 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.

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.where N M is a number of free parameters in model M. 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 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 Table D.1.Minimum χ 2 values resulting from sampling with KiDS and Planck data separately as well as from the combined analysis.In the last column, the tension between the two datasets is shown by the difference in maximum a posteriori criterion.The last two rows show differences of χ 2 min between the ΛCDM and ΛDDM cases, subtracting the former from the latter, and the difference in Akaike information criterion ∆AIC between the ΛCDM and ΛDDM scenarios.

WL
CMB Combined √ Q DMAP χ 2 min (ΛCDM) 158.7 580.2 750.5 3.4σ χ 2 min (ΛDDM) 158.6 580.2 742.6 1.9σ ∆χ 2 min -0.1 0.0 -7.9 -∆AIC 3.9 4.0 -3.9 fit to the data.In the combined scenario, ∆AIC = −3.9< ∆AIC 0 .This means that the two additional parameters in the combined analysis are efficient.However, the combined ΛDDM fit is still worse by ∆χ 2 min = 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 8 tension, in which the ΛDDM model fails.Throughout the work, we obtained the value of χ 2 min as a minimum χ 2 among all sampled configurations in our MCMC chain.For a summary of combined MCMC analysis results, we refer to Tab.D.1.A152, page 13 of 15

Fig. 2 .
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, v k , Γ, z).Next, we extracted power spectra and calculated the resulting power suppression P DDM /P ΛCDM .Additionally, we performed PCA on P DDM /P ΛCDM data, and using five principal components (c pca 1 , c pca 2 , c pca 3 , c pca 4 , c pca 5 ), we recovered the original power suppression.Once the training data set was prepared, we trained the SIREN-like neural network.We used ( f, v k , Γ, 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.

Fig. 3 .
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.

Fig. 5 .
Fig.5.Two-dimensional posterior distributions of Ω m and S 8 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.

Fig. 7 .
Fig. 7. One-dimensional constraints of the S 8 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.

Fig
Fig. C.2. Effects of baryons on band power spectra illustrated using different baryonic parameter configurations.Band power auto-and crosscorrelations 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
Fig. C.2 demonstrates the impact of baryonic feedback onthe 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 10 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 8 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.
Fig. C.3, where the left panel depicts the combination of DDM parameters with log 10 M c and the right panel illustrates their combination with θ ej in the ΛDDM model.In both panels, it is evident that fitting the WLonly 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
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.
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 ∆AIC = ∆χ 2 min + 2(N ΛDDM − N ΛCDM ), (D.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 , v k , Γ, and z are shown in the panels a, b, c, and d.

Table 1 .
Parameter samples in MCMC chains.

Table E .
2. Same as Tab.E.1 but assuming the ΛDDM model.