Issue 
A&A
Volume 683, March 2024



Article Number  A152  
Number of page(s)  15  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202347844  
Published online  15 March 2024 
Probing the twobody decaying dark matter scenario with weak lensing and the cosmic microwave background
^{1}
Department of Astrophysics, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland
email: jozef.bucko@uzh.ch
^{2}
Nordita, KTH Royal Institute of Technology and Stockholm University, Hannes Alfens väg 12, 106 91 Stockholm, Sweden
^{3}
Université ParisSaclay, Université Paris Cité, CEA, CNRS, Astrophysique, Instrumentation et Modélisation ParisSaclay, 91191 GifsurYvette, France
Received:
29
August
2023
Accepted:
5
October
2023
Decaying dark matter (DDM) scenarios have recently regained attention due to their potential ability to resolve the wellknown clustering (or S_{8}) tension between weak lensing (WL) and cosmic microwave background (CMB) measurements. In this paper, we investigate a wellestablished 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 nonlinear WL analysis of this twobody decaying dark matter (ΛDDM) scenario, including intrinsic alignment and baryonic feedback processes. We used the cosmic shear band power spectra from KiDS1000 data and combined it with temperature and polarisation data from Planck in order to constrain the Λ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 lowmass splittings. Regarding the S_{8} tension, we found a reduction from about 3 to 2σ, depending on which statistical measure is applied. We therefore conclude that the twobody ΛDDM model is able to reduce the S_{8} tension without convincingly solving it.
Key words: cosmological parameters / dark matter / largescale structure of Universe
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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), 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 , 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 lowredshift surveys report consistently lower values. For example, the KiloDegree Survey (KiDS; Kuijken et al. 2019; Giblin et al. 2021; Hildebrandt et al. 2021) obtained a value of (Asgari et al. 2021), in agreement with (albeit slightly lower than) the results from the Hyper SupremeCam (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 nonlinear clustering of matter (e.g. Tan et al. 2023; Aricò et al. 2023) or the modelling of cosmic shear (GarcíaGarcí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 coldwarm 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 (CyrRacine 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 twobody 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 freestreaming length of the stable daughter particles and thus by the magnitude of the velocity kicks (v_{k}) they receive as a consequence of energymomentum 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), Rparity violation (Berezinsky et al. 1991; Kim & Kim 2002), or super weakly interacting massive particles (Covi et al. 1999; Feng et al. 2003).
The twobody 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 Nbody simulations ranging from individual galaxies; Peter & Benson 2010 to a largescale 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 twobody decay rate in the regime of lowmass 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 twobody decays. After including priors from WL observations, they reported a reduction of the S_{8} tension for a bestfitting Λ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 KiDS1000 survey (Asgari et al. 2021). The nonlinear clustering predictions were thereby modelled using an emulator based on a suite of Nbody simulations. Along with the WL analysis, we also performed a reanalysis 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 Nbody simulations to previous studies. Appendix B presents a discussion on the cosmology dependence of the ΛDDM effects, and Appendix C details the effects twobody 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. Twobody 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 Nbody simulations, and finally, we introduce a new emulator of the ΛDDM nonlinear 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 restmass 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 centreofmomentum frame is (Abellán et al. 2021), where c is the speed of light. We note that in the nonrelativistic 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)
or using ε instead of particle masses, as
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
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 twobody decays by taking arbitrary Γ, ε, and f values, with limiting cases Γ → 0, ε → 0 or f → 0 approaching the ΛCDM cosmology and ε → 1/2 approaching onebody 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^{−1}. We therefore focus on the regime with latetime decays (Γ ≲ H_{0}) and nonrelativistic velocity kicks (v_{k} ≪ c) throughout this paper.
2.2. Simulations
Considering only nonrelativistic decays, we implemented the ΛDDM model into the Nbody code PKDGRAV3 (Potter et al. 2017), a treebased gravity solver based on fast multipole 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 nonrelativistic velocity kicks received by the WDM particles. We note that this differs from the onebody DDM model studied in Bucko et al. (2023), where the background evolution had to be modified.
The twobody 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 simulation particles undergo the decay process, where 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 twobody 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 kaxis. 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 twobody decays on the nonlinear matter spectrum from Nbody 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. 
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 Nbody 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 nonlinear 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 twobody 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 Nbody simulations are not fast enough for this purpose, we built an emulator to account for the different ΛDDM parameters. The basic characteristics of our emulatorbuilding procedure are shown in the flowchart of Fig. 2. First, we ran ∼100 gravityonly simulations for different dark matter parameters (Γ, v_{k}, and f) and measured the nonlinear 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 . 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 Nbody 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 , we recovered the original power suppression. Once the training data set was prepared, we trained the SIRENlike 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 feedforward 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 finetune 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 kbins. 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 nonlinear matter power spectrum . The prediction time of our emulator is a few milliseconds. With this tool, we could model the nonlinear power spectrum in the presence of dark matter decays as
Fig. 3. Performance of the ΛDDM emulator at different scales. The precision of the twobody 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 kbins, and the 1σ (grey) and 2σ (golden) error bars of the prediction are highlighted. 
where is obtained by multiplying the nonlinear power spectrum from the revised_halofit method (Takahashi et al. 2012). Our emulator to model the nonlinear ΛDDM power spectrum is published as part of the publicly available code DMemu^{1}.
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 KiDS1000 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 KiDS1000
The KiDS1000 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 crosscorrelation band power of all five redshift bins. With eight data points for each spectrum, the KiDS1000 band power contains a total of 120 observational data points with correlated errors from the corresponding covariance matrix.
3.2.1. Nonlinear 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 nonlinear 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^{2} (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 threeparameter model presented in Giri & Schneider (2021), where four parameters are fixed based on the results from hydrodynamical simulations. The threeparameter 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 twobody dark matter decays was multiplied with the nonlinear, 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 twobody 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 nonlinear 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 galaxyintrinsic and intrinsicintrinsic 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 KiDS1000 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 autocorrelation 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 multipurpose 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 crossband power measurements are illustrated in five KiDS1000 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^{3}. 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 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 (ForemanMackey 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 latetime 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 𝒩(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).
Parameter samples in MCMC chains.
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 (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 GelmanRubin 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 , , and correspond to the minimum chisquared 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 .
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 multifluid dark matter sector to future work.
5. Results
In the first part of this section, we describe the constraints on twobody 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.
5.1. Derived constraints on twobody 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 halflife times of τ < 30 Gyr, while the Lymanα appears to be less sensitive to twobody 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:
Fig. 4. Twodimensional 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). 
Hence, within the parameter space of ΛDDM and the datasets analysed in this study, WL data impose much stronger limits compared to CMB observations^{5}. This is mainly due to the fact that, unlike for the case of onebody decays (Bucko et al. 2023), the twobody ΛDDM model does not significantly affect the background evolution of the Universe, leaving the signal from the latetime integrated SachsWolfe effect unchanged. More discussions about the effects of twobody decays on the CMB signal can be found in Appendix C.
5.2. 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. 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 KiDS1000 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.
Fig. 5. Twodimensional 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. 
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 twobody 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 nonGaussian shape of the CMB posteriors.
The origin of the nonGaussian 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.
Fig. 6. Effects of twobody 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 bestfit model in each scenario as coloured stars. 
Based on the analysis above, we concluded that although the twobody Λ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 nonzero DDM parameters. The best fitting values are 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 CMBonly, WLonly, 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).
Fig. 7. Onedimensional constraints of the S_{8} parameter as inferred from the CMB, WL, and combined analyses. The original results from the KiDS1000 (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 twobody 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, 2022; Simon et al. 2022).
In this paper, we performed the first fully nonlinear analysis of WL and CMB data for the twobody decaying dark matter (ΛDDM) scenario. Based on a suite of Nbody simulations, we constructed a neural networkbased emulator to obtain fast predictions of nonlinear matter power spectra for arbitrary values of the decay rate (Γ), the decayinduced 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 KiDS1000 and CMB data from Planck 2018.
We present improved constraints on the twobody decaying dark matter parameters based on the WL data from KiDS1000. 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 KiDS1000 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 nonGaussianity 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 twobody Λ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, selfconsistent 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 galaxygalaxy lensing. Furthermore, including additional lowredshift 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 twobody 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.
The BCemu code can be found at https://github.com/sambitgiri/BCemu
adapted from Bucko et al. (2023)
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
 Abazajian, K. N., Acero, M. A., Agarwalla, S. K., et al. 2012, arXiv eprints [arXiv:1204.5379] [Google Scholar]
 Abellán, G. F., Murgia, R., & Poulin, V. 2021, Phys. Rev. D, 104, 123533 [CrossRef] [Google Scholar]
 Abellán, G. F., Murgia, R., Poulin, V., & Lavalle, J. 2022, Phys. Rev. D, 105, 063525 [CrossRef] [Google Scholar]
 Adhikari, R., Agostini, M., Ky, N. A., et al. 2017, JCAP, 2017, 025 [CrossRef] [Google Scholar]
 Aihara, H., Arimoto, N., Armstrong, R., et al. 2017, PASJ, 70 [Google Scholar]
 Amon, A., & Efstathiou, G. 2022, MNRAS, 516, 5355 [CrossRef] [Google Scholar]
 Amon, A., Gruen, D., Troxel, M. A., et al. 2022, Phys. Rev. D, 105, 023514 [NASA ADS] [CrossRef] [Google Scholar]
 Aricò, G., Angulo, R. E., Zennaro, M., et al. 2023, A&A, 678, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Asgari, M., Lin, C.A., Joachimi, B., et al. 2021, A&A, 645, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Aurélien, G. 2019, Handson Machine Learning with ScikitLearn, Keras, and TensorFlow: Concepts, Tools, and Techniques to Build Intelligent Systems (OReilly Media, Inc.) [Google Scholar]
 Berezinsky, V., Masiero, A., & Valle, J. 1991, Phys. Lett. B, 266, 382 [NASA ADS] [CrossRef] [Google Scholar]
 Blackadder, G., & Koushiappas, S. M. 2014, Phys. Rev. D, 90, 103527 [NASA ADS] [CrossRef] [Google Scholar]
 Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 2011, 034 [Google Scholar]
 Bucko, J., Giri, S. K., & Schneider, A. 2023, A&A, 672, A157 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chen, A., Huterer, D., Lee, S., et al. 2021, Phys. Rev. D, 103, 123528 [NASA ADS] [CrossRef] [Google Scholar]
 Cheng, D., Chu, M.C., & Tang, J. 2015, JCAP, 2015, 009 [CrossRef] [Google Scholar]
 Chisari, N. E., Mead, A. J., Joudaki, S., et al. 2019, Open J. Astrophys., 2, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Choi, G., & Yanagida, T. T. 2022, Phys. Lett. B, 827, 136954 [NASA ADS] [CrossRef] [Google Scholar]
 Covi, L., Kim, J. E., & Roszkowski, L. 1999, Phys. Rev. Lett., 82, 4180 [NASA ADS] [CrossRef] [Google Scholar]
 CyrRacine, F.Y., Sigurdson, K., Zavala, J., et al. 2016, Phys. Rev. D, 93, 123527 [NASA ADS] [CrossRef] [Google Scholar]
 Enqvist, K., Nadathur, S., Sekiguchi, T., & Takahashi, T. 2015, JCAP, 09, 067 [CrossRef] [Google Scholar]
 Enqvist, K., Nadathur, S., Sekiguchi, T., & Takahashi, T. 2020, JCAP, 04, 015 [CrossRef] [Google Scholar]
 Falkner, S., Klein, A., & Hutter, F. 2018, in Proceedings of the 35th International Conference on Machine Learning, eds. J. Dy, & A. Krause (Stockholmsmässan, Stockholm Sweden: PMLR), Proc. Mach. Learn. Res., 80, 1437 [Google Scholar]
 Feng, J. L., Rajaraman, A., & Takayama, F. 2003, Phys. Rev. Lett., 91, 011302 [NASA ADS] [CrossRef] [Google Scholar]
 Ferlito, F., Vagnozzi, S., Mota, D. F., & Baldi, M. 2022, MNRAS, 512, 1885 [CrossRef] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
 Fuß, L., & Garny, M. 2023, JCAP, 2023, 020 [Google Scholar]
 GarcíaGarcía, C., RuizZapatero, J., Alonso, D., et al. 2021, JCAP, 2021, 030 [CrossRef] [Google Scholar]
 Gelman, A., & Rubin, D. B. 1992, Stat. Sci., 7, 457 [Google Scholar]
 Giblin, B., Heymans, C., Asgari, M., et al. 2021, A&A, 645, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Giri, S. K., & Schneider, A. 2021, JCAP, 2021, 046 [Google Scholar]
 Hamana, T., Shirasaki, M., Miyazaki, S., et al. 2020, PASJ, 72, 16 [Google Scholar]
 Hambye, T. 2011, in Proceedings of Identification of Dark Matter 2010  PoS(IDM2010), 110, 098 [CrossRef] [Google Scholar]
 Heimersheim, S., Schöneberg, N., Hooper, D. C., & Lesgourgues, J. 2020, JCAP, 2020, 016 [CrossRef] [Google Scholar]
 Hildebrandt, H., Viola, M., Heymans, C., et al. 2016, MNRAS, 465, 1454 [Google Scholar]
 Hildebrandt, H., van den Busch, J. L., Wright, A. H., et al. 2021, A&A, 647, A124 [EDP Sciences] [Google Scholar]
 Hirata, C. M., Seljak, U. C. V., 2004, Phys. Rev. D, 70, 063526 [NASA ADS] [CrossRef] [Google Scholar]
 Hubert, J., Schneider, A., Potter, D., Stadel, J., & Giri, S. K. 2021, JCAP, 2021, 040 [CrossRef] [Google Scholar]
 Jeffreys, H. 1961, Theory of Probability (Clarendon) [Google Scholar]
 Joachimi, B., Lin, C.A., Asgari, M., et al. 2021, A&A, 646, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Joseph, M., Aloni, D., Schmaltz, M., Sivarajan, E. N., & Weiner, N. 2023, Phys. Rev. D, 108, 023520 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, H. B., & Kim, J. E. 2002, Phys. Lett. B, 527, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Kingma, D. P., & Ba, J. 2014, arXiv eprints [arXiv:1412.6980] [Google Scholar]
 Kuijken, K., Heymans, C., Dvornik, A., et al. 2019, A&A, 625, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mau, S., Nadler, E. O., Wechsler, R. H., et al. 2022, ApJ, 932, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Murgia, R., Merle, A., Viel, M., Totzauer, M., & Schneider, A. 2017, JCAP, 2017, 046 [CrossRef] [Google Scholar]
 Pappadopulo, D., Ruderman, J. T., & Trevisan, G. 2016, Phys. Rev. D, 94, 035005 [NASA ADS] [CrossRef] [Google Scholar]
 Parimbelli, G., Scelfo, G., Giri, S., et al. 2021, JCAP, 2021, 044 [CrossRef] [Google Scholar]
 Paszke, A., Gross, S., Massa, F., et al. 2019, Advances in Neural Information Processing Systems 32 (Curran Associates, Inc.), 8024 [Google Scholar]
 Peter, A. H. G., & Benson, A. J. 2010, Phys. Rev. D, 82, 123521 [NASA ADS] [CrossRef] [Google Scholar]
 Peter, A. H. G., Moody, C. E., & Kamionkowski, M. 2010, Phys. Rev. D, 81, 103501 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration V. 2020, A&A, 641, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Potter, D., Stadel, J., & Teyssier, R. 2017, Comput. Astrophys. Cosmol., 4, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Poulin, V., Bernal, J. L., Kovetz, E., & Kamionkowski, M. 2023, Phys. Rev. D, 107, 123538 [NASA ADS] [CrossRef] [Google Scholar]
 Pourtsidou, A., Skordis, C., & Copeland, E. J. 2013, Phys. Rev. D, 88, 083505 [NASA ADS] [CrossRef] [Google Scholar]
 Prince, H., & Dunkley, J. 2019, Phys. Rev. D, 100, 083502 [NASA ADS] [CrossRef] [Google Scholar]
 Raveri, M., & Hu, W. 2019, Phys. Rev. D, 99, 043506 [NASA ADS] [CrossRef] [Google Scholar]
 Rubira, H., Mazoun, A., & Garny, M. 2023, JCAP, 2023, 034 [CrossRef] [Google Scholar]
 Schöneberg, N., Abellán, G. F., Sánchez, A. P., et al. 2022, Phys. Rep., 984, 1 [CrossRef] [Google Scholar]
 Schneider, A., & Teyssier, R. 2015, JCAP, 1512, 049 [Google Scholar]
 Schneider, A., Teyssier, R., Stadel, J., et al. 2019, JCAP, 03, 020 [CrossRef] [Google Scholar]
 Schneider, A., Refregier, A., Grandis, S., et al. 2020a, JCAP, 04, 020 [Google Scholar]
 Schneider, A., Stoira, N., Refregier, A., et al. 2020b, JCAP, 04, 019 [CrossRef] [Google Scholar]
 Schneider, A., Giri, S. K., Amodeo, S., & Refregier, A. 2022, MNRAS, 514, 3802 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, T., Abellán, G. F., Du, P., Poulin, V., & Tsai, Y. 2022, Phys. Rev. D, 106, 023516 [NASA ADS] [CrossRef] [Google Scholar]
 Sitzmann, V., Martel, J. N. P., Bergman, A. W., Lindell, D. B., & Wetzstein, G. 2020, arXiv eprints [arXiv:2006.09661] [Google Scholar]
 Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152 [Google Scholar]
 Tan, T., Zuercher, D., Fluri, J., et al. 2023, MNRAS, 522, 3766 [NASA ADS] [CrossRef] [Google Scholar]
 Tanimura, H., Douspis, M., Aghanim, N., & Kuruvilla, J. 2023, A&A, 674, A222 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tarsitano, F., Schmitt, U., Refregier, A., et al. 2021, Astron. Comput., 36, 100484 [NASA ADS] [CrossRef] [Google Scholar]
 The Dark Energy Survey Collaboration 2005, arXiv eprints [arXiv:astroph/0510346] [Google Scholar]
 van Daalen, M. P., Schaye, J., Booth, C. M., & Vecchia, C. D. 2011, MNRAS, 415, 3649 [Google Scholar]
 Vogelsberger, M., Zavala, J., CyrRacine, F.Y., et al. 2016, MNRAS, 460, 1399 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, M.Y., Croft, R. A. C., Peter, A. H. G., Zentner, A. R., & Purcell, C. W. 2013, Phys. Rev. D, 88, 123515 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, M.Y., & Zentner, A. R. 2012, Phys. Rev. D, 85, 043514 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Comparing twobody decay simulations to previous work
In Fig. A.1, we present a comparison of our Nbody 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 Nbody 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 nonlinear 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 Nbody implementation. Therefore, we did not perform model inferences for such extreme values of v_{k} in this work.
Fig. A.1. Comparison of our Nbody 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 twobody decays published in Cheng et al. (2015) compared to our results for different decay rates Γ and velocity kick magnitudes v_{k}. 
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 Nbody 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 nonlinear 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 nonlinear 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 cosmologyindependent boost to the ΛCDM nonlinear 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 Nbody simulations where we kept the ΛDDM parameters fixed at Γ^{−1} = 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 secondorder effect. In panels (a) to (e) of this figure, we show the impact of the Hubble parameter (h_{0}), clustering amplitude (σ_{8}), spectral index (n_{s}), matter abundance (Ω_{m}), and baryon abundance (Ω_{b}), respectively. We display the ΛDDM boost for the fiducial cosmology of our Nbody simulations with a solid salmon line, while we plot the same quantity for the bestfit cosmology of the KiDS450 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 twobody 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 (h_{0}), clustering amplitude (σ_{8}), spectral index (n_{s}), matter abundance (Ω_{m}), and baryon abundance (Ω_{b}). With a solid salmon line, we show the fiducial cosmology of our Nbody simulations, while the cosmology of the KiDS450 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 σ_{8} (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 σ_{8} 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 KiDS1000 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 bestfit configuration in the WLonly setup in the case of ΛDDM (however, it performs very similar to the ΛCDM case). Once we calculated the band powers resulting from the CMBonly bestfit cosmology (orange lines), we recovered a significantly stronger WL signal compared to the KiDS observations. This is the consequence of larger clustering (S_{8} 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 bestfit cosmology that is close to the one from the CMBonly scenario, but baryons and, namely, twobody decays on small scales provide a much better fit to the cosmic shear data compared to the CMBonly scenario, as seen from the solid purple lines in the figure. To also showcase that twobody 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 twobody decays completely, while the latter displays the largest possible impact of the twobody decays in the case of cosmology and baryons fixed to the WL plus CMB ΛDDM bestfit 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 higherredshift bins, for example. This also explains why in the case of the combined ΛDDM analysis, we recovered the constraints on the twobody parameters that are detached from ΛCDM cosmology.
Fig. C.1. Effects of twobody decays on band power spectra illustrated using different forwardmodelled ΛDDM configurations. Band power auto and crosscorrelations for five tomographic bins with respective error bars as measured by the KiDS1000 survey are shown as black dots. Solid green (orange) lines represent the bestfit 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 twobody 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_{10} 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_{10}M_{c} and θ_{ej} have an impact on the WL observables. For log_{10}M_{c}, we obtained log_{10}M_{c} < 13.1, (13.2) for ΛCDM (ΛDDM) from the WLonly analysis and log_{10}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 WLonly 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 WLonly 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 crosscorrelations for five tomographic bins with respective error bars as measured by the KiDS1000 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 bestfit configurations for ΛCDM (ΛDDM) are depicted in dark blue (green), while the coloured lines represent the variation of baryonic feedback strength in the ΛCDM bestfit 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 twobody 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 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 bottomleft 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 twobody decays in the ΛDDM universe as obtained from probing the WL (green) and WL plus CMB (purple) scenarios. Both baryons and twobody 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 twobody 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 crosscorrelations 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 bestfit 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 twobody decays do not change the CMB signal significantly, meaning that the scatter of the CMB data is larger than the observable effects of twobody 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
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
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_{0}. This means that the two additional parameters in the combined analysis are efficient. However, the combined ΛDDM fit is still worse by 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 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.
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 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.
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.^{6}
Summary of our MCMC analysis assuming the ΛCDM model. We report individual results separately based on WL (KiDS1000) and CMB (Planck 2018) data only as well as values inferred from the combined MCMC chain. We show the mean (best fit) values of the sampled (top) and derived (middle) parameters as well as the obtained prior, likelihood, and χ^{2} values (bottom).
Same as Tab. E.1 but assuming the ΛDDM model.
All Tables
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 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.
Summary of our MCMC analysis assuming the ΛCDM model. We report individual results separately based on WL (KiDS1000) and CMB (Planck 2018) data only as well as values inferred from the combined MCMC chain. We show the mean (best fit) values of the sampled (top) and derived (middle) parameters as well as the obtained prior, likelihood, and χ^{2} values (bottom).
All Figures
Fig. 1. Effects of twobody decays on the nonlinear matter spectrum from Nbody 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. 

In the text 
Fig. 2. Flowchart describing the emulation process (from left to right). The Nbody 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 , we recovered the original power suppression. Once the training data set was prepared, we trained the SIRENlike 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. 

In the text 
Fig. 3. Performance of the ΛDDM emulator at different scales. The precision of the twobody 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 kbins, and the 1σ (grey) and 2σ (golden) error bars of the prediction are highlighted. 

In the text 
Fig. 4. Twodimensional 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). 

In the text 
Fig. 5. Twodimensional 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. 

In the text 
Fig. 6. Effects of twobody 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 bestfit model in each scenario as coloured stars. 

In the text 
Fig. 7. Onedimensional constraints of the S_{8} parameter as inferred from the CMB, WL, and combined analyses. The original results from the KiDS1000 (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 Nbody 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 twobody decays published in Cheng et al. (2015) compared to our results for different decay rates Γ and velocity kick magnitudes v_{k}. 

In the text 
Fig. B.1. Cosmology dependence of the twobody 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 (h_{0}), clustering amplitude (σ_{8}), spectral index (n_{s}), matter abundance (Ω_{m}), and baryon abundance (Ω_{b}). With a solid salmon line, we show the fiducial cosmology of our Nbody simulations, while the cosmology of the KiDS450 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 twobody decays on band power spectra illustrated using different forwardmodelled ΛDDM configurations. Band power auto and crosscorrelations for five tomographic bins with respective error bars as measured by the KiDS1000 survey are shown as black dots. Solid green (orange) lines represent the bestfit 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 twobody 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 crosscorrelations for five tomographic bins with respective error bars as measured by the KiDS1000 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 twobody decays in the ΛDDM universe as obtained from probing the WL (green) and WL plus CMB (purple) scenarios. Both baryons and twobody 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.