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

Jozef Bucko, Sambit K Giri, Fabian Hervas Peters, Aurel Schneider

► To cite this version:

Jozef Bucko, Sambit K Giri, Fabian Hervas Peters, Aurel Schneider. Probing the two-body decaying dark matter scenario with weak lensing and the cosmic microwave background. Astron.Astrophys., 2024, 683, pp.A152. 10.1051/0004-6361/202347844. hal-04170695

HAL Id: hal-04170695

https://hal.science/hal-04170695

Submitted on 19 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

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

Jozef Bucko, Sambit K Giri, Fabian Hervas Peters, Aurel Schneider

► To cite this version:

Jozef Bucko, Sambit K Giri, Fabian Hervas Peters, Aurel Schneider. Probing the two-body decaying dark matter scenario with weak lensing and the cosmic microwave background. Astronomy & Astrophysics - A&A, 2024, 683, pp.A152. 10.1051/0004-6361/202347844. hal-04170695

HAL Id: hal-04170695

https://hal.science/hal-04170695v1

Submitted on 19 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Probing the two-body decaying dark matter scenario with weak lensing and the cosmic microwave background

Jozef Bucko<fnref n="1" />, Sambit K. Giri<fnref n="1" />,<fnref n="2" />, Fabian Hervas Peters<fnref n="1" />,<fnref n="3" />, and Aurel Schneider<fnref n="1" />

Received 29 August 2023 / Accepted 5 October 2023

## ABSTRACT

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 (ΛDDM) scenario, including intrinsic alignment and baryonic feedback processes. We used the cosmic shear band power spectra from KiDS-1000 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 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 ΛDDM model is able to reduce the $S_8$ tension without convincingly solving it.

**Key words.** cosmological parameters – dark matter – large-scale structure of Universe

## 1. Introduction

The standard Λ-cold dark matter (ΛCDM) cosmology provides an outstanding description of the Universe, explaining a wide range of cosmic observations such as the cosmic microwave background (CMB), baryonic acoustic oscillations (BAO), large-scale structure formation, and Big Bang nucleosynthesis (BBN). Despite its tremendous success, there are still questions to which ΛCDM, as understood nowadays, cannot provide satisfying answers, including the nature of dark matter and dark energy. Moreover, with progressively precise measurements, several discrepancies within the model have emerged. An example of such a discrepancy is the mild yet persistent clustering amplitude tension between the CMB and weak lensing (WL) measurements, often expressed via the parameter $S_8 = \sigma_8 \sqrt{\Omega_m/0.3}$, with $\sigma_8$ and $\Omega_m$ describing the clustering amplitude and the total matter abundance. More specifically, CMB measurements from the Planck Collaboration have yielded $S_8 = 0.834 \pm 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 \sim 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, 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.

<footnote n="1">¹ Department of Astrophysics, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland
e-mail: jozef.bucko@uzh.ch</footnote>

<footnote n="2">² Nordita, KTH Royal Institute of Technology and Stockholm University, Hannes Alfvénsväg 12, 106 91 Stockholm, Sweden</footnote>

<footnote n="3">³ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, Astrophysique, Instrumentation et Modélisation Paris-Saclay,
91191 Gif-sur-Yvette, France</footnote>
# Probing the two-body decaying dark matter scenario with weak lensing and the cosmic microwave background

Jozef Bucko<fnref n="1" />, Sambit K. Giri<fnref n="1" />,<fnref n="2" />, Fabian Hervas Peters<fnref n="1" />,<fnref n="3" />, and Aurel Schneider<fnref n="1" />

Received 29 August 2023 / Accepted 5 October 2023

## ABSTRACT

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 (ΛDDM) scenario, including intrinsic alignment and baryonic feedback processes. We used the cosmic shear band power spectra from KiDS-1000 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 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 ΛDDM model is able to reduce the $S_8$ tension without convincingly solving it.

**Key words.** cosmological parameters – dark matter – large-scale structure of Universe

## 1. Introduction

The standard Λ-cold dark matter (ΛCDM) cosmology provides an outstanding description of the Universe, explaining a wide range of cosmic observations such as the cosmic microwave background (CMB), baryonic acoustic oscillations (BAO), large-scale structure formation, and Big Bang nucleosynthesis (BBN). Despite its tremendous success, there are still questions to which ΛCDM, as understood nowadays, cannot provide satisfying answers, including the nature of dark matter and dark energy. Moreover, with progressively precise measurements, several discrepancies within the model have emerged. An example of such a discrepancy is the mild yet persistent clustering amplitude tension between the CMB and weak lensing (WL) measurements, often expressed via the parameter $S_8 = \sigma_8 \sqrt{\Omega_m/0.3}$, with $\sigma_8$ and $\Omega_m$ describing the clustering amplitude and the total matter abundance. More specifically, CMB measurements from the Planck Collaboration have yielded $S_8 = 0.834 \pm 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 \sim 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, 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.

<footnote n="1">¹ Department of Astrophysics, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland
e-mail: jozef.bucko@uzh.ch</footnote>

<footnote n="2">² Nordita, KTH Royal Institute of Technology and Stockholm University, Hannes Alfvénsväg 12, 106 91 Stockholm, Sweden</footnote>

<footnote n="3">³ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, Astrophysique, Instrumentation et Modélisation Paris-Saclay,
91191 Gif-sur-Yvette, France</footnote>
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 \sim 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 $\tau = 120$ Gyr and $v_k/c \approx 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.

## 2. Two-body decaying dark matter

In the ΛDDM model, an original ('mother') dark matter particle decays into a slightly lighter, stable particle and a relativistic massless relic, while the energy released during the decay

process is split between two product species. We describe the basic theoretical properties of the ΛDDM model in Sect. 2.1. In Sect. 2.2, we discuss our ΛDDM N-body simulations, and finally, we introduce a new emulator of the ΛDDM non-linear matter power spectra in Sect. 2.3.

### 2.1. Theory

The ΛDDM model is simple enough to be described by two phenomenological parameters. The first parameter is the decay rate $\Gamma$ 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 $\epsilon$ of the rest-mass energy

$$ \epsilon = \frac{1}{2} \left( 1 - \frac{m_{\text{wdm}}^2}{M_{\text{dcdm}}^2} \right), \qquad (1) $$

where $m_{\text{wdm}}$ and $M_{\text{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_{\text{wdm}} = m_{\text{wdm}}c\epsilon/\sqrt{1-2\epsilon}$ (Abellán et al. 2021), where c is the speed of light. We note that in the non-relativistic limit, we obtained a simple relation $v_k = \epsilon 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)

$$ \dot{\rho}_{\text{dcdm}} + 3H\rho_{\text{dcdm}} = -a\Gamma\rho_{\text{dcdm}}, \qquad (2) $$

$$ \dot{\rho}_{\text{wdm}} + 3H\rho_{\text{wdm}} = a\Gamma \frac{M_{\text{dcdm}}^2 + m_{\text{wdm}}^2}{2M_{\text{dcdm}}^2} \rho_{\text{dcdm}}, \qquad (3) $$

$$ \dot{\rho}_{\text{dr}} + 4H\rho_{\text{dr}} = a\Gamma \frac{M_{\text{dcdm}}^2 - m_{\text{wdm}}^2}{2M_{\text{dcdm}}^2} \rho_{\text{dcdm}}, \qquad (4) $$

or using $\epsilon$ instead of particle masses, as

$$ \dot{\rho}_{\text{dcdm}} + 3H\rho_{\text{dcdm}} = -a\Gamma\rho_{\text{dcdm}}, \qquad (5) $$

$$ \dot{\rho}_{\text{wdm}} + 3H\rho_{\text{wdm}} = a\Gamma(1 - \epsilon)\rho_{\text{dcdm}}, \qquad (6) $$

$$ \dot{\rho}_{\text{dr}} + 4H\rho_{\text{dr}} = a\Gamma\epsilon\rho_{\text{dcdm}}. \qquad (7) $$

In the above equations, $H$ is the conformal Hubble parameter, and $\rho_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 $\Gamma$ and $v_k$, we include a third parameter

$$ f = \frac{\Omega_{\text{dcdm}}}{\Omega_{\text{dcdm}} + \Omega_{\text{cdm}}}, \qquad (8) $$

that allows for a scenario where only a fraction $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 ($\Omega_{\text{dcdm}}$) and unstable ($\Omega_{\text{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 $\Gamma$, $\epsilon$, and $f$ values, with limiting cases $\Gamma \to 0$, $\epsilon \to 0$ or $f \to 0$ approaching the ΛCDM cosmology and
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 \sim 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 $\tau = 120$ Gyr and $v_k/c \approx 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.

## 2. Two-body decaying dark matter

In the ΛDDM model, an original ('mother') dark matter particle decays into a slightly lighter, stable particle and a relativistic massless relic, while the energy released during the decay

process is split between two product species. We describe the basic theoretical properties of the ΛDDM model in Sect. 2.1. In Sect. 2.2, we discuss our ΛDDM N-body simulations, and finally, we introduce a new emulator of the ΛDDM non-linear matter power spectra in Sect. 2.3.

### 2.1. Theory

The ΛDDM model is simple enough to be described by two phenomenological parameters. The first parameter is the decay rate $\Gamma$ 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 $\epsilon$ of the rest-mass energy

$$ \epsilon = \frac{1}{2} \left( 1 - \frac{m_{\text{wdm}}^2}{M_{\text{dcdm}}^2} \right), \qquad (1) $$

where $m_{\text{wdm}}$ and $M_{\text{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_{\text{wdm}} = m_{\text{wdm}}c\epsilon/\sqrt{1-2\epsilon}$ (Abellán et al. 2021), where c is the speed of light. We note that in the non-relativistic limit, we obtained a simple relation $v_k = \epsilon 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)

$$ \dot{\rho}_{\text{dcdm}} + 3\mathcal{H}\rho_{\text{dcdm}} = -a\Gamma\rho_{\text{dcdm}}, \qquad (2) $$

$$ \dot{\rho}_{\text{wdm}} + 3\mathcal{H}\rho_{\text{wdm}} = a\Gamma \frac{M_{\text{dcdm}}^2 + m_{\text{wdm}}^2}{2M_{\text{dcdm}}^2} \rho_{\text{dcdm}}, \qquad (3) $$

$$ \dot{\rho}_{\text{dr}} + 4\mathcal{H}\rho_{\text{dr}} = a\Gamma \frac{M_{\text{dcdm}}^2 - m_{\text{wdm}}^2}{2M_{\text{dcdm}}^2} \rho_{\text{dcdm}}, \qquad (4) $$

or using $\epsilon$ instead of particle masses, as

$$ \dot{\rho}_{\text{dcdm}} + 3\mathcal{H}\rho_{\text{dcdm}} = -a\Gamma\rho_{\text{dcdm}}, \qquad (5) $$

$$ \dot{\rho}_{\text{wdm}} + 3\mathcal{H}\rho_{\text{wdm}} = a\Gamma(1 - \epsilon)\rho_{\text{dcdm}}, \qquad (6) $$

$$ \dot{\rho}_{\text{dr}} + 4\mathcal{H}\rho_{\text{dr}} = a\Gamma\epsilon\rho_{\text{dcdm}}. \qquad (7) $$

In the above equations, $\mathcal{H}$ is the conformal Hubble parameter, and $\rho_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 $\Gamma$ and $v_k$, we include a third parameter

$$ f = \frac{\Omega_{\text{dcdm}}}{\Omega_{\text{dcdm}} + \Omega_{\text{cdm}}}, \qquad (8) $$

that allows for a scenario where only a fraction $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 ($\Omega_{\text{dcdm}}$) and unstable ($\Omega_{\text{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 $\Gamma$, $\epsilon$, and $f$ values, with limiting cases $\Gamma \to 0$, $\epsilon \to 0$ or $f \to 0$ approaching the ΛCDM cosmology and
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Effects of two-body decays on the non-linear matter spectrum from N-body simulations. The fiducial (ΛCDM) power spectrum corresponds to a constant line equal to one, while the power suppressions for varying *f*, *v*<sub>*k*</sub>, Γ, and *z* are shown in the panels a, b, c, and d.</figcaption></figure>

$\varepsilon \to 1/2$ approaching one-body decays. However, as recent studies have demonstrated (Wang & Zentner 2012; Mau et al. 2022, and references therein), ΛDDM models with a very large decay rate and velocity kicks are ruled out by observations, as they lead to a strong power suppression at $k < 0.1$ h Mpc$^{-1}$. We therefore focus on the regime with late-time decays ($\Gamma \lesssim H_0$) and non-relativistic velocity kicks ($v_k \ll c$) throughout this paper.

## 2.2. Simulations

Considering only non-relativistic decays, we implemented the ΛDDM model into the N-body code PKDGRAV3 (Potter et al. 2017), a tree-based gravity solver based on fast multi-pole expansion and adaptive time stepping. Following the theoretical descriptions (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 $\Delta t$). The decay probability of a given (not yet decayed) particle at time step *i* is $P = \Gamma\Delta t$. Thus, a number of $\Delta N_{wdm}^i = \Gamma\Delta t N_{dcdm}^i$ simulation particles undergo the decay process, where $N_{dcdm}^i$ denotes the number of unstable DCDM particles at time step *i*. The particles that are about to decay are chosen randomly from all CDM particles. Immediately after the decay, they obtain a uniform velocity kick in a random direction. Importantly, these particles are flagged and added to a set of already decayed particles that are excluded from the decay process occurring in future time steps.

In Fig. 1, we plot the ratios of simulated ΛDDM to ΛCDM power spectra for varying values of *f*, *v*<sub>*k*</sub>, and $1/\Gamma$ (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 $\Gamma < 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$, $\Omega_{m,0} = 0.307$, $\Omega_{\Lambda,0} = 0.693$, $\Omega_{b,0} = 0.048$, $n_s = 0.966$, and $\sigma_8 = 0.883$. We obtained converged results at the scales $k \sim 0.01-10$ h Mpc$^{-1}$ for box sizes of $L_{\text{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).

## 2.3. Emulating the impact of dark matter decays

In order to carry out a Bayesian inference analysis (Sect. 4), we needed a fast modelling framework to explore the vast parameter space of astrophysics, cosmology, and dark matter models. As N-body simulations are not fast enough for this purpose, we built an emulator to account for the different ΛDDM parameters. The basic characteristics of our emulator-building procedure are shown in the flowchart of Fig. 2. First, we ran ~100 gravity-only simulations for different dark matter parameters ($\Gamma$, $v_k$, and *f*) and measured the non-linear matter power spectra up to $k \sim 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_{\Lambda\text{DDM}}^{\Gamma, v_k, f}(k, z)$. We found that five PCA components are sufficient to describe the ratio of spectra with a reconstruction error of approximately 0.1%.
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Effects of two-body decays on the non-linear matter spectrum from $N$-body simulations. The fiducial (ΛCDM) power spectrum corresponds to a constant line equal to one, while the power suppressions for varying $f$, $v_k$, $\Gamma$, and $z$ are shown in the panels a, b, c, and d.</figcaption></figure>

$\varepsilon \to 1/2$ approaching one-body decays. However, as recent studies have demonstrated (Wang & Zentner 2012; Mau et al. 2022, and references therein), ΛDDM models with a very large decay rate and velocity kicks are ruled out by observations, as they lead to a strong power suppression at $k < 0.1$ h Mpc$^{-1}$. We therefore focus on the regime with late-time decays ($\Gamma \lesssim H_0$) and non-relativistic velocity kicks ($v_k \ll c$) throughout this paper.

## 2.2. Simulations

Considering only non-relativistic decays, we implemented the ΛDDM model into the $N$-body code PKDGRAV3 (Potter et al. 2017), a tree-based gravity solver based on fast multi-pole expansion and adaptive time stepping. Following the theoretical descriptions (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 $\Delta t$). The decay probability of a given (not yet decayed) particle at time step $i$ is $P = \Gamma\Delta t$. Thus, a number of $\Delta N_{wdm}^i = \Gamma\Delta t N_{dcdm}^i$ simulation particles undergo the decay process, where $N_{dcdm}^i$ denotes the number of unstable DCDM particles at time step $i$. The particles that are about to decay are chosen randomly from all CDM particles. Immediately after the decay, they obtain a uniform velocity kick in a random direction. Importantly, these particles are flagged and added to a set of already decayed particles that are excluded from the decay process occurring in future time steps.

In Fig. 1, we plot the ratios of simulated ΛDDM to ΛCDM power spectra for varying values of $f$, $v_k$, and $1/\Gamma$ (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 $\Gamma < 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$, $\Omega_{m,0} = 0.307$, $\Omega_{\Lambda,0} = 0.693$, $\Omega_{b,0} = 0.048$, $n_s = 0.966$, and $\sigma_8 = 0.883$. We obtained converged results at the scales $k \sim 0.01–10$ h Mpc$^{-1}$ for box sizes of $L_{\text{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).

## 2.3. Emulating the impact of dark matter decays

In order to carry out a Bayesian inference analysis (Sect. 4), we needed a fast modelling framework to explore the vast parameter space of astrophysics, cosmology, and dark matter models. As $N$-body simulations are not fast enough for this purpose, we built an emulator to account for the different ΛDDM parameters. The basic characteristics of our emulator-building procedure are shown in the flowchart of Fig. 2. First, we ran ~100 gravity-only simulations for different dark matter parameters ($\Gamma$, $v_k$, and $f$) and measured the non-linear matter power spectra up to $k \sim 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_{\Lambda\text{DDM}}^{\Gamma, v_k, f}(k, z)$. We found that five PCA components are sufficient to describe the ratio of spectra with a reconstruction error of approximately 0.1%.
<figure><img src="image_3.png" /><figcaption>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, \Gamma, z$). Next, we extracted power spectra and calculated the resulting power suppression $P_{\text{DDM}}/P_{\text{ACDM}}$. Additionally, we performed PCA on $P_{\text{DDM}}/P_{\text{ACDM}}$ data, and using five principal components ($c_1^{\text{PCA}}, c_2^{\text{PCA}}, c_3^{\text{PCA}}, c_4^{\text{PCA}}, c_5^{\text{PCA}}$), 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, \Gamma, 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.</figcaption></figure>

Next, we trained a neural network to model these five PCA components of the simulated power spectra ratios for a given parameter vector ($\Gamma, 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 $\lambda$ 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\sigma$ and $2\sigma$ errors stay below 1%. Our emulator thus efficiently predicts the response of decays on the non-linear matter power spectrum $S_{\text{ΛDDM}}^{\Gamma, v_k, f}(k, z)$. The prediction time of our emulator is a few milliseconds. With this tool, we

<figure><img src="image_4.png" /><figcaption>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\sigma$ (grey) and $2\sigma$ (golden) error bars of the prediction are highlighted.</figcaption></figure>

could model the non-linear power spectrum in the presence of dark matter decays as

$$P_{\text{ΛDDM}}^{\text{nonlin}}(k, z) = S_{\text{ΛDDM}}^{\Gamma, v_k, f}(k, z) \times P_{\text{ΛCDM}}^{\text{nonlin}}(k, z), \quad (9)$$

where $P_{\text{ΛCDM}}^{\text{nonlin}}(k, z)$ is obtained by multiplying the non-linear power spectrum from the revised_halofit method (Takahashi et al. 2012). Our emulator to model the non-linear ΛDDM power spectrum is published as part of the publicly available code DMemu<fnref n="1" />.

<footnote n="1">¹ https://github.com/jbucko/DMemu</footnote>
<figure><img src="image_4.png" /><figcaption>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, \Gamma, z$). Next, we extracted power spectra and calculated the resulting power suppression $P_{\text{DDM}}/P_{\text{ACDM}}$. Additionally, we performed PCA on $P_{\text{DDM}}/P_{\text{ACDM}}$ data, and using five principal components ($c_1^{\text{PCA}}, c_2^{\text{PCA}}, c_3^{\text{PCA}}, c_4^{\text{PCA}}, c_5^{\text{PCA}}$), 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, \Gamma, 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.</figcaption></figure>

Next, we trained a neural network to model these five PCA components of the simulated power spectra ratios for a given parameter vector ($\Gamma, 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 $\lambda$ 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\sigma$ and $2\sigma$ errors stay below 1%. Our emulator thus efficiently predicts the response of decays on the non-linear matter power spectrum $S_{\text{ADDM}}^{\Gamma, v_k, f}(k, z)$. The prediction time of our emulator is a few milliseconds. With this tool, we

<figure><img src="image_5.png" /><figcaption>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\sigma$ (grey) and $2\sigma$ (golden) error bars of the prediction are highlighted.</figcaption></figure>

could model the non-linear power spectrum in the presence of dark matter decays as

$$P_{\Lambda\text{DDM}}^{\text{nonlin}}(k, z) = S_{\Lambda\text{DDM}}^{\Gamma, v_k, f}(k, z) \times P_{\Lambda\text{CDM}}^{\text{nonlin}}(k, z), \quad (9)$$

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

<footnote n="1">¹ https://github.com/jbucko/DMemu</footnote>
# 3. Data sets and modelling framework

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

## 3.1. Data sets

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

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

## 3.2. Cosmic shear modelling for KiDS-1000

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

### 3.2.1. Non-linear matter power spectrum

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

The process of baryonic feedback causes gas to be expelled out of galaxies and clusters, leading to a suppression of the matter power spectrum at small cosmological scales (e.g. van Daalen et al. 2011; Schneider & Teyssier 2015; Chisari et al. 2019). We note that this suppression is similar in shape to the one caused by dark matter decays (Hubert et al. 2021), making it particularly important to include baryonic feedback in our modelling pipeline. We used BCemu<fnref n="2" /> (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, \theta_{ej}$) and one parameter describing

the stellar mass ($\eta_\delta$) 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 was then obtained from the emulator described in Sect. 2.3. This modular method assumes that all dependence on cosmological parameters is captured by the power spectrum from the revised halofit (see Appendix B). At the same time, the responses from the baryonification and the two-body decays remain independent of cosmology. Regarding the baryonification method, this assumption has been validated before (Schneider et al. 2020a,b). In Appendix B, we validated this assumption for the ΛDDM model.

### 3.2.2. Intrinsic alignment

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

### 3.2.3. Angular power spectrum

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

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

<footnote n="2">² The BCemu code can be found at https://github.com/sambit-giri/BCemu</footnote>

<footnote n="3">³ https://github.com/PoulinV/class_decays</footnote>
# 3. Data sets and modelling framework

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

## 3.1. Data sets

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

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

## 3.2. Cosmic shear modelling for KiDS-1000

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

### 3.2.1. Non-linear matter power spectrum

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

The process of baryonic feedback causes gas to be expelled out of galaxies and clusters, leading to a suppression of the matter power spectrum at small cosmological scales (e.g. van Daalen et al. 2011; Schneider & Teyssier 2015; Chisari et al. 2019). We note that this suppression is similar in shape to the one caused by dark matter decays (Hubert et al. 2021), making it particularly important to include baryonic feedback in our modelling pipeline. We used BCemu<fnref n="2" /> (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, \theta_{ej}$) and one parameter describing

the stellar mass ($\eta_\delta$) around galaxies and clusters. We refer to Giri & Schneider (2021) for detailed description of the model.

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

### 3.2.2. Intrinsic alignment

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

### 3.2.3. Angular power spectrum

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

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

<footnote n="2">² The BCemu code can be found at https://github.com/sambit-giri/BCemu</footnote>

<footnote n="3">³ https://github.com/PoulinV/class_decays</footnote>
<table><caption>**Table 1.** Parameter samples in MCMC chains.</caption><thead><tr><th>Parameter name</th><th>Acronym</th><th>Prior</th><th>Range</th></tr></thead><tbody><tr><td>(Initial) DM abundance</td><td>ω<sub>dm</sub></td><td>Flat</td><td>[0.051, 0.255]</td></tr><tr><td>Baryon abundance</td><td>ω<sub>b</sub></td><td>Flat</td><td>[0.019, 0.026]</td></tr><tr><td>Scalar amplitude</td><td>ln(10<sup>10</sup>A<sub>s</sub>)</td><td>Flat</td><td>[1.0, 5.0]</td></tr><tr><td>Hubble constant</td><td>h<sub>0</sub></td><td>Flat</td><td>[0.6, 0.8]</td></tr><tr><td>Spectral index</td><td>n<sub>s</sub></td><td>Flat</td><td>[0.9, 1.03]</td></tr><tr><td>Optical depth</td><td>τ<sub>reio</sub></td><td>Normal</td><td>N(0.0506, 0.0086)</td></tr><tr><td>Intrinsic alignment amplitude</td><td>A<sub>IA</sub></td><td>Flat</td><td>[0.0, 2.0]</td></tr><tr><td>Planck calibration parameter</td><td>A<sub>Planck</sub></td><td>Normal</td><td>N(1.0, 0.0025)</td></tr><tr><td>First gas parameter (BCemu)</td><td>log<sub>10</sub> M<sub>c</sub></td><td>Flat</td><td>[11.0, 15.0]</td></tr><tr><td>Second gas parameter (BCemu)</td><td>θ<sub>ej</sub></td><td>Flat</td><td>[2.0, 8.0]</td></tr><tr><td>Stellar parameter (BCemu)</td><td>η<sub>δ</sub></td><td>Flat</td><td>[0.05, 0.40]</td></tr><tr><td>Decay rate (Gyr<sup>-1</sup>)</td><td>log<sub>10</sub> Γ</td><td>Flat</td><td>[-4.00, -1.13]</td></tr><tr><td>Velocity kick magnitude (km s<sup>-1</sup>)</td><td>log<sub>10</sub> v<sub>k</sub></td><td>Flat</td><td>[0.00, 3.69]</td></tr></tbody></table>

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

## 4. Model inference

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

An overview of the sampled parameters and the prior choices are listed in Table 1. We used flat priors for all cosmological parameters except for the optical depth $\tau$, 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 ($\omega_{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 $\omega_{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 $\Lambda$CDM and $\Lambda$DDM scenarios. For the baryon abundance $\omega_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$, $\theta_{ej}$, and $\eta_{\delta}$ covering the full range of the BCemu parameters. For the $\Lambda$DDM decay rate ($\Gamma$) and the velocity kick magnitude ($v_k$), we assumed flat priors for both the $log_{10} \Gamma$ and $log_{10} v_k$ spanning from $\Lambda$CDM values up to the upper boundary limited by the range of the $\Lambda$DDM emulator (see Table 1).

In both the WL and CMB setups, we assumed Gaussian likelihoods. In the case of the CMB, we used a marginalised light-weight version of the full Planck likelihood called Plikann

(Prince & Dunkley 2019; Planck Collaboration V 2020), as it is an affordable approximation in the case of $\Lambda$DDM. Abellán et al. (2021) demonstrated that using Plikann for $\Lambda$DDM produces only a negligible difference on the recovered posteriors compared to a full Planck analysis. Also, Plikann 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

$$ \tau_{S_8} = \frac{S_8^{\text{CMB}} - S_8^{\text{WL}}}{\sqrt{\text{Var}[S_8^{\text{CMB}}] + \text{Var}[S_8^{\text{WL}}]}}, \quad (10) $$

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

$$ Q_{\text{DMAP}} = \chi^2_{\text{min,WL+CMB}} - (\chi^2_{\text{min,WL}} + \chi^2_{\text{min,CMB}}), \quad (11) $$

where $\chi^2_{\text{min,WL}}$, $\chi^2_{\text{min,CMB}}$, and $\chi^2_{\text{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 $\sigma$, the tension is quantified as $\sqrt{Q_{DMAP}}$.

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

## 5. Results

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

### 5.1. Derived constraints on two-body decays

We obtained the constraints on the decay rate $\Gamma$ and velocity kick magnitude $v_k$ by marginalising over cosmology, baryons, and nuisance parameters. We display the outcome in Fig. 4, showing the $\Lambda$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 $\Lambda$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-$\alpha$ forest, and Abellán et al. (2021) and

<footnote n="4">⁴ https://wiki.cosmos.esa.int/planck-legacy-archive/index.php/CMB_spectrum_%26_Likelihood_Code</footnote>
<table><caption>**Table 1.** Parameter samples in MCMC chains.</caption><thead><tr><th>Parameter name</th><th>Acronym</th><th>Prior</th><th>Range</th></tr></thead><tbody><tr><td>(Initial) DM abundance</td><td>ω<sub>dm</sub></td><td>Flat</td><td>[0.051, 0.255]</td></tr><tr><td>Baryon abundance</td><td>ω<sub>b</sub></td><td>Flat</td><td>[0.019, 0.026]</td></tr><tr><td>Scalar amplitude</td><td>ln(10<sup>10</sup>A<sub>s</sub>)</td><td>Flat</td><td>[1.0, 5.0]</td></tr><tr><td>Hubble constant</td><td>h<sub>0</sub></td><td>Flat</td><td>[0.6, 0.8]</td></tr><tr><td>Spectral index</td><td>n<sub>s</sub></td><td>Flat</td><td>[0.9, 1.03]</td></tr><tr><td>Optical depth</td><td>τ<sub>reio</sub></td><td>Normal</td><td>N(0.0506, 0.0086)</td></tr><tr><td>Intrinsic alignment amplitude</td><td>A<sub>IA</sub></td><td>Flat</td><td>[0.0, 2.0]</td></tr><tr><td>Planck calibration parameter</td><td>A<sub>Planck</sub></td><td>Normal</td><td>N(1.0, 0.0025)</td></tr><tr><td>First gas parameter (BCemu)</td><td>log<sub>10</sub> M<sub>c</sub></td><td>Flat</td><td>[11.0, 15.0]</td></tr><tr><td>Second gas parameter (BCemu)</td><td>θ<sub>ej</sub></td><td>Flat</td><td>[2.0, 8.0]</td></tr><tr><td>Stellar parameter (BCemu)</td><td>η<sub>δ</sub></td><td>Flat</td><td>[0.05, 0.40]</td></tr><tr><td>Decay rate (Gyr<sup>-1</sup>)</td><td>log<sub>10</sub> Γ</td><td>Flat</td><td>[-4.00, -1.13]</td></tr><tr><td>Velocity kick magnitude (km s<sup>-1</sup>)</td><td>log<sub>10</sub> v<sub>k</sub></td><td>Flat</td><td>[0.00, 3.69]</td></tr></tbody></table>

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

## 4. Model inference

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

An overview of the sampled parameters and the prior choices are listed in Table 1. We used flat priors for all cosmological parameters except for the optical depth $\tau$, 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 ($\omega_{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 $\omega_{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 $\Lambda$CDM and $\Lambda$DDM scenarios. For the baryon abundance $\omega_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$, $\theta_{ej}$, and $\eta_{\delta}$ covering the full range of the BCemu parameters. For the $\Lambda$DDM decay rate ($\Gamma$) and the velocity kick magnitude ($v_k$), we assumed flat priors for both the $log_{10} \Gamma$ and $log_{10} v_k$ spanning from $\Lambda$CDM values up to the upper boundary limited by the range of the $\Lambda$DDM emulator (see Table 1).

In both the WL and CMB setups, we assumed Gaussian likelihoods. In the case of the CMB, we used a marginalised light-weight version of the full Planck likelihood called Plikann

(Prince & Dunkley 2019; Planck Collaboration V 2020), as it is an affordable approximation in the case of $\Lambda$DDM. Abellán et al. (2021) demonstrated that using Plikann for $\Lambda$DDM produces only a negligible difference on the recovered posteriors compared to a full Planck analysis. Also, Plikann 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

$$ \tau_{S_8} = \frac{S_8^{\text{CMB}} - S_8^{\text{WL}}}{\sqrt{\text{Var}[S_8^{\text{CMB}}] + \text{Var}[S_8^{\text{WL}}]}}, \quad (10) $$

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

$$ Q_{\text{DMAP}} = \chi^2_{\text{min,WL+CMB}} - (\chi^2_{\text{min,WL}} + \chi^2_{\text{min,CMB}}), \quad (11) $$

where $\chi^2_{\text{min,WL}}$, $\chi^2_{\text{min,CMB}}$, and $\chi^2_{\text{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 $\sigma$, the tension is quantified as $\sqrt{Q_{DMAP}}$.

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

## 5. Results

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

### 5.1. Derived constraints on two-body decays

We obtained the constraints on the decay rate $\Gamma$ and velocity kick magnitude $v_k$ by marginalising over cosmology, baryons, and nuisance parameters. We display the outcome in Fig. 4, showing the $\Lambda$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 $\Lambda$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-$\alpha$ forest, and Abellán et al. (2021) and

<footnote n="4">⁴ https://wiki.cosmos.esa.int/planck-legacy-archive/index.php/CMB_spectrum_%26_Likelihood_Code</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Two-dimensional constraints on $\Lambda$DDM parameters $\Gamma$ and $v_k$ as obtained from WL cosmic shear (green) and CMB (orange) data. We also display the results from several recent studies constraining $\Lambda$DDM parameters using different types of observables (see Sect. 5.1 for more details).</figcaption></figure>

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 $\Gamma$ and rule out half-life times of $\tau < 30$ Gyr, while the Lyman-$\alpha$ 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:

$$ \tau = \Gamma^{-1} \le 125 \text{ Gyr, and} \qquad (12) $$

$$ v_k \ge 300 \text{ km s}^{-1}. \qquad (13) $$

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

## 5.2. Revisiting the $S_8$ tension between WL and the CMB

As a next step, we investigated the extent the $\Lambda$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 $\Lambda$DDM model is unable to significantly reduce the systematic shift between the clustering signal predicted from the CMB and the KiDS-1000 WL survey.

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Two-dimensional posterior distributions of $\Omega_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 ($\Lambda$CDM) and green ($\Lambda$DDM), while the findings based on CMB data are displayed as black ($\Lambda$CDM) and orange ($\Lambda$DDM) contours.</figcaption></figure>

In Fig. 5, we illustrate the $\Omega_m - S_8$ posterior contours as obtained from the individual WL and CMB data analyses for both the $\Lambda$CDM (dark blue and black) and the $\Lambda$DDM model (green and orange). We observed only a marginal change in the WL and CMB contours when going from the $\Lambda$CDM to the $\Lambda$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.

<footnote n="5">⁵ We find that it is important to note that in our analysis, we do not include the CMB lensing effect.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Two-dimensional constraints on $\Lambda$DDM parameters $\Gamma$ and $v_k$ as obtained from WL cosmic shear (green) and CMB (orange) data. We also display the results from several recent studies constraining $\Lambda$DDM parameters using different types of observables (see Sect. 5.1 for more details).</figcaption></figure>

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 $\Gamma$ and rule out half-life times of $\tau < 30$ Gyr, while the Lyman-$\alpha$ 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:

$$ \tau = \Gamma^{-1} \le 125 \text{ Gyr, and} \qquad (12) $$

$$ v_k \ge 300 \text{ km s}^{-1}. \qquad (13) $$

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

## 5.2. Revisiting the $S_8$ tension between WL and the CMB

As a next step, we investigated the extent the $\Lambda$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 $\Lambda$DDM model is unable to significantly reduce the systematic shift between the clustering signal predicted from the CMB and the KiDS-1000 WL survey.

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Two-dimensional posterior distributions of $\Omega_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 ($\Lambda$CDM) and green ($\Lambda$DDM), while the findings based on CMB data are displayed as black ($\Lambda$CDM) and orange ($\Lambda$DDM) contours.</figcaption></figure>

In Fig. 5, we illustrate the $\Omega_m - S_8$ posterior contours as obtained from the individual WL and CMB data analyses for both the $\Lambda$CDM (dark blue and black) and the $\Lambda$DDM model (green and orange). We observed only a marginal change in the WL and CMB contours when going from the $\Lambda$CDM to the $\Lambda$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.

<footnote n="5">⁵ We find that it is important to note that in our analysis, we do not include the CMB lensing effect.</footnote>
<figure><img src="image_7.png" /><figcaption>**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 $\Omega_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.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

Applying the $Q_{DMAP}$ criterion defined in Eq. (11), we obtained a mutual tension of $3.4\sigma$ between WL and CMB data for the case of $\Lambda$CDM. When including two-body decays, the tension is reduced to $1.9\sigma$. Using the Gaussian metric defined in Eq. (10), we obtained an $S_8$ tension of $3.0\sigma$ for the $\Lambda$CDM and $2.7\sigma$ for the $\Lambda$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 $\Omega - 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 $\Lambda$DDM parameters ($\Gamma \times v_k$). One can see that for large values of $\Gamma \times 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 $\Lambda$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 $\Gamma$ and $v_k$. For further tests of the $\Lambda$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} \Gamma = -2.25_{-0.23}^{+0.74}$ 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 non-vanishing $\Lambda$DDM parameters in the combined data analysis. We

want to stress that this is by no means a signal of a departure from $\Lambda$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 $\Lambda$CDM nor in the $\Lambda$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).

## 6. Summary and conclusion

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

In this paper, we performed the first fully non-linear analysis of WL and CMB data for the two-body decaying dark matter ($\Lambda$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 ($\Gamma$), 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
<figure><img src="image_8.png" /><figcaption>**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 $\Omega_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.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>

Applying the $Q_{DMAP}$ criterion defined in Eq. (11), we obtained a mutual tension of $3.4\sigma$ between WL and CMB data for the case of $\Lambda$CDM. When including two-body decays, the tension is reduced to $1.9\sigma$. Using the Gaussian metric defined in Eq. (10), we obtained an $S_8$ tension of $3.0\sigma$ for the $\Lambda$CDM and $2.7\sigma$ for the $\Lambda$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 $\Omega - 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 $\Lambda$DDM parameters ($\Gamma \times v_k$). One can see that for large values of $\Gamma \times 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 $\Lambda$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 $\Gamma$ and $v_k$. For further tests of the $\Lambda$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} \Gamma = -2.25_{-0.23}^{+0.74}$ 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 non-vanishing $\Lambda$DDM parameters in the combined data analysis. We

want to stress that this is by no means a signal of a departure from $\Lambda$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 $\Lambda$CDM nor in the $\Lambda$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).

## 6. Summary and conclusion

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

In this paper, we performed the first fully non-linear analysis of WL and CMB data for the two-body decaying dark matter ($\Lambda$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 ($\Gamma$), the decay-induced velocity kick ($v_k$), and the fraction of decaying to total dark matter ($f$). We then included the emulator in our pipeline predicting WL
observations and performed an MCMC analysis with WL data from KiDS-1000 and CMB data from Planck 2018.

We present improved constraints on the two-body decaying dark matter parameters based on the WL data from KiDS-1000. Our constraints are significantly stronger compared to previous results. Specifically, we exclude models with $\tau = \Gamma^{-1} \le 125$ Gyr and $v_k \ge 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\sigma$ with $\Lambda$DDM as compared to the original $3.4\sigma$ tension measured in the $\Lambda$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 $\Lambda$DDM scenario is unable to convincingly resolve the clustering tension between WL and CMB observations. We note that previous works obtaining different conclusions (e.g. Abellán et al. 2021, 2022; Simon et al. 2022; Tanimura et al. 2023) did not include a full, self-consistent modelling of the WL signal and were therefore unable to directly test the $S_8$ tension in the case of $\Lambda$DDM.

A further step forward with respect to the current analysis could involve the analysis of the $\Lambda$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 $\Omega_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 $\Lambda$DDM model.

*Acknowledgements.* 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 e-prints [arXiv:1204.5379]

Abellán, G. F., Murgia, R., & Poulin, V. 2021, Phys. Rev. D, 104, 123533

Abellán, G. F., Murgia, R., Poulin, V., & Lavalle, J. 2022, Phys. Rev. D, 105, 063525

Adhikari, R., Agostini, M., Ky, N. A., et al. 2017, JCAP, 2017, 025

Aihara, H., Arimoto, N., Armstrong, R., et al. 2017, PASJ, 70

Amon, A., & Efstathiou, G. 2022, MNRAS, 516, 5355

Amon, A., Gruen, D., Troxel, M. A., et al. 2022, Phys. Rev. D, 105, 023514

Aricò, G., Angulo, R. E., Zennaro, M., et al. 2023, A&A, 678, A109

Asgari, M., Lin, C.-A., Joachimi, B., et al. 2021, A&A, 645, A104

Aurélien, G. 2019, Hands-on Machine Learning with Scikit-Learn, Keras, and TensorFlow: Concepts, Tools, and Techniques to Build Intelligent Systems (OReilly Media, Inc.)

Berezinsky, V., Masiero, A., & Valle, J. 1991, Phys. Lett. B, 266, 382

Blackadder, G., & Koushiappas, S. M. 2014, Phys. Rev. D, 90, 103527

Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 2011, 034

Bucko, J., Giri, S. K., & Schneider, A. 2023, A&A, 672, A157

Chen, A., Huterer, D., Lee, S., et al. 2021, Phys. Rev. D, 103, 123528

Cheng, D., Chu, M.-C., & Tang, J. 2015, JCAP, 2015, 009

Chisari, N. E., Mead, A. J., Joudaki, S., et al. 2019, Open J. Astrophys., 2, 4

Choi, G., & Yanagida, T. T. 2022, Phys. Lett. B, 827, 136954

Covi, L., Kim, J. E., & Roszkowski, L. 1999, Phys. Rev. Lett., 82, 4180

Cyr-Racine, F.-Y., Sigurdson, K., Zavala, J., et al. 2016, Phys. Rev. D, 93, 123527

Enqvist, K., Nadathur, S., Sekiguchi, T., & Takahashi, T. 2015, JCAP, 09, 067

Enqvist, K., Nadathur, S., Sekiguchi, T., & Takahashi, T. 2020, JCAP, 04, 015

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

Feng, J. L., Rajaraman, A., & Takayama, F. 2003, Phys. Rev. Lett., 91, 011302

Ferlito, F., Vagnozzi, S., Mota, D. F., & Baldi, M. 2022, MNRAS, 512, 1885

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fuß, L., & Garny, M. 2023, JCAP, 2023, 020

García-García, C., Ruiz-Zapatero, J., Alonso, D., et al. 2021, JCAP, 2021, 030

Gelman, A., & Rubin, D. B. 1992, Stat. Sci., 7, 457

Giblin, B., Heymans, C., Asgari, M., et al. 2021, A&A, 645, A105

Giri, S. K., & Schneider, A. 2021, JCAP, 2021, 046

Hamana, T., Shirasaki, M., Miyazaki, S., et al. 2020, PASJ, 72, 16

Hambye, T. 2011, in Proceedings of Identification of Dark Matter 2010 - PoS(IDM2010), 110, 098

Heimersheim, S., Schöneberg, N., Hooper, D. C., & Lesgourgues, J. 2020, JCAP, 2020, 016

Hildebrandt, H., Viola, M., Heymans, C., et al. 2016, MNRAS, 465, 1454

Hildebrandt, H., van den Busch, J. L., Wright, A. H., et al. 2021, A&A, 647, A124

Hirata, C. M., Seljak, U., & c. v., 2004, Phys. Rev. D, 70, 063526

Hubert, J., Schneider, A., Potter, D., Stadel, J., & Giri, S. K. 2021, JCAP, 2021, 040

Jeffreys, H. 1961, Theory of Probability (Clarendon)

Joachimi, B., Lin, C.-A., Asgari, M., et al. 2021, A&A, 646, A129

Joseph, M., Aloni, D., Schmaltz, M., Sivarajan, E. N., & Weiner, N. 2023, Phys. Rev. D, 108, 023520

Kim, H. B., & Kim, J. E. 2002, Phys. Lett. B, 527, 18

Kingma, D. P., & Ba, J. 2014, arXiv e-prints [arXiv:1412.6980]

Kuijken, K., Heymans, C., Dvornik, A., et al. 2019, A&A, 625, A2

Mau, S., Nadler, E. O., Wechsler, R. H., et al. 2022, ApJ, 932, 128

Murgia, R., Merle, A., Viel, M., Totzauer, M., & Schneider, A. 2017, JCAP, 2017, 046

Pappadopulo, D., Ruderman, J. T., & Trevisan, G. 2016, Phys. Rev. D, 94, 035005

Parimbelli, G., Scelfo, G., Giri, S., et al. 2021, JCAP, 2021, 044

Paszke, A., Gross, S., Massa, F., et al. 2019, Advances in Neural Information Processing Systems 32 (Curran Associates, Inc.), 8024

Peter, A. H. G., & Benson, A. J. 2010, Phys. Rev. D, 80, 123521

Peter, A. H. G., Moody, C. E., & Kamionkowski, M. 2010, Phys. Rev. D, 81, 103501

Planck Collaboration V. 2020, A&A, 641, A5

Planck Collaboration VI. 2020, A&A, 641, A6

Potter, D., Stadel, J., & Teyssier, R. 2017, Comput. Astrophys. Cosmol., 4, 2

Poulin, V., Bernal, J. L., Kovetz, E., & Kamionkowski, M. 2023, Phys. Rev. D, 107, 123538

Poursidou, A., Skordis, C., & Copeland, E. J. 2013, Phys. Rev. D, 88, 083505

Prince, H., & Dunkley, J. 2019, Phys. Rev. D, 100, 083502

Raveri, M., & Hu, W. 2019, Phys. Rev. D, 99, 043506

Rubira, H., Mazoun, A., & Garny, M. 2023, JCAP, 2023, 034

Schöneberg, N., Abellán, G. F., Sánchez, A. P., et al. 2022, Phys. Rep., 984, 1

Schneider, A., & Teyssier, R. 2015, JCAP, 1512, 049

Schneider, A., Teyssier, R., Stadel, J., et al. 2019, JCAP, 03, 020

Schneider, A., Refregier, A., Grandis, S., et al. 2020a, JCAP, 04, 020

Schneider, A., Stoira, N., Refregier, A., et al. 2020b, JCAP, 04, 019

Schneider, A., Giri, S. K., Amodeo, S., & Refregier, A. 2022, MNRAS, 514, 3802

Simon, T., Abellán, G. F., Du, P., Poulin, V., & Tsai, Y. 2022, Phys. Rev. D, 106, 023516

Sitzmann, V., Martel, J. N. P., Bergman, A. W., Lindell, D. B., & Wetzstein, G. 2020, arXiv e-prints [arXiv:2006.09661]

Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152

Tan, T., Zuercher, D., Fluri, J., et al. 2023, MNRAS, 522, 3766

Tanimura, H., Douspis, M., Aghanim, N., & Kuruvilla, J. 2023, A&A, 674, A222

Tarsitano, F., Schmitt, U., Refregier, A., et al. 2021, Astron. Comput., 36, 100484

The Dark Energy Survey Collaboration 2005, arXiv e-prints [arXiv:astro-ph/0510346]

van Daalen, M. P., Schaye, J., Booth, C. M., & Vecchia, C. D. 2011, MNRAS, 415, 3649

Vogelsberger, M., Zavala, J., Cyr-Racine, F.-Y., et al. 2016, MNRAS, 460, 1399

Wang, M.-Y., Croft, R. A. C., Peter, A. H. G., Zentner, A. R., & Purcell, C. W. 2013, Phys. Rev. D, 88, 123515

Wang, M.-Y., & Zentner, A. R. 2012, Phys. Rev. D, 85, 043514
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 $\tau = \Gamma^{-1} \lesssim 125$ Gyr and $v_k \gtrsim 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\sigma$ with $\Lambda$DDM as compared to the original $3.4\sigma$ tension measured in the $\Lambda$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 $\Lambda$DDM scenario is unable to convincingly resolve the clustering tension between WL and CMB observations. We note that previous works obtaining different conclusions (e.g. Abellán et al. 2021, 2022; Simon et al. 2022; Tanimura et al. 2023) did not include a full, self-consistent modelling of the WL signal and were therefore unable to directly test the $S_8$ tension in the case of $\Lambda$DDM.

A further step forward with respect to the current analysis could involve the analysis of the $\Lambda$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 $\Omega_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 $\Lambda$DDM model.

Acknowledgements. 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 e-prints [arXiv:1204.5379]

Abellán, G. F., Murgia, R., & Poulin, V. 2021, Phys. Rev. D, 104, 123533

Abellán, G. F., Murgia, R., Poulin, V., & Lavalle, J. 2022, Phys. Rev. D, 105, 063525

Adhikari, R., Agostini, M., Ky, N. A., et al. 2017, JCAP, 2017, 025

Aihara, H., Arimoto, N., Armstrong, R., et al. 2017, PASJ, 70

Amon, A., & Efstathiou, G. 2022, MNRAS, 516, 5355

Amon, A., Gruen, D., Troxel, M. A., et al. 2022, Phys. Rev. D, 105, 023514

Aricò, G., Angulo, R. E., Zennaro, M., et al. 2023, A&A, 678, A109

Asgari, M., Lin, C.-A., Joachimi, B., et al. 2021, A&A, 645, A104

Aurélien, G. 2019, Hands-on Machine Learning with Scikit-Learn, Keras, and TensorFlow: Concepts, Tools, and Techniques to Build Intelligent Systems (OReilly Media, Inc.)

Berezinsky, V., Masiero, A., & Valle, J. 1991, Phys. Lett. B, 266, 382

Blackadder, G., & Koushiappas, S. M. 2014, Phys. Rev. D, 90, 103527

Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 2011, 034

Bucko, J., Giri, S. K., & Schneider, A. 2023, A&A, 672, A157

Chen, A., Huterer, D., Lee, S., et al. 2021, Phys. Rev. D, 103, 123528

Cheng, D., Chu, M.-C., & Tang, J. 2015, JCAP, 2015, 009

Chisari, N. E., Mead, A. J., Joudaki, S., et al. 2019, Open J. Astrophys., 2, 4

Choi, G., & Yanagida, T. T. 2022, Phys. Lett. B, 827, 136954

Covi, L., Kim, J. E., & Roszkowski, L. 1999, Phys. Rev. Lett., 82, 4180

Cyr-Racine, F.-Y., Sigurdson, K., Zavala, J., et al. 2016, Phys. Rev. D, 93, 123527

Enqvist, K., Nadathur, S., Sekiguchi, T., & Takahashi, T. 2015, JCAP, 09, 067

Enqvist, K., Nadathur, S., Sekiguchi, T., & Takahashi, T. 2020, JCAP, 04, 015

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

Feng, J. L., Rajaraman, A., & Takayama, F. 2003, Phys. Rev. Lett., 91, 011302

Ferlito, F., Vagnozzi, S., Mota, D. F., & Baldi, M. 2022, MNRAS, 512, 1885

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fuß, L., & Garny, M. 2023, JCAP, 2023, 020

García-García, C., Ruiz-Zapatero, J., Alonso, D., et al. 2021, JCAP, 2021, 030

Gelman, A., & Rubin, D. B. 1992, Stat. Sci., 7, 457

Giblin, B., Heymans, C., Asgari, M., et al. 2021, A&A, 645, A105

Giri, S. K., & Schneider, A. 2021, JCAP, 2021, 046

Hamana, T., Shirasaki, M., Miyazaki, S., et al. 2020, PASJ, 72, 16

Hambye, T. 2011, in Proceedings of Identification of Dark Matter 2010 - PoS(IDM2010), 110, 098

Heimersheim, S., Schöneberg, N., Hooper, D. C., & Lesgourgues, J. 2020, JCAP, 2020, 016

Hildebrandt, H., Viola, M., Heymans, C., et al. 2016, MNRAS, 465, 1454

Hildebrandt, H., van den Busch, J. L., Wright, A. H., et al. 2021, A&A, 647, A124

Hirata, C. M., Seljak, U., & c. v., 2004, Phys. Rev. D, 70, 063526

Hubert, J., Schneider, A., Potter, D., Stadel, J., & Giri, S. K. 2021, JCAP, 2021, 040

Jeffreys, H. 1961, Theory of Probability (Clarendon)

Joachimi, B., Lin, C.-A., Asgari, M., et al. 2021, A&A, 646, A129

Joseph, M., Aloni, D., Schmaltz, M., Sivarajan, E. N., & Weiner, N. 2023, Phys. Rev. D, 108, 023520

Kim, H. B., & Kim, J. E. 2002, Phys. Lett. B, 527, 18

Kingma, D. P., & Ba, J. 2014, arXiv e-prints [arXiv:1412.6980]

Kuijken, K., Heymans, C., Dvornik, A., et al. 2019, A&A, 625, A2

Mau, S., Nadler, E. O., Wechsler, R. H., et al. 2022, ApJ, 932, 128

Murgia, R., Merle, A., Viel, M., Totzauer, M., & Schneider, A. 2017, JCAP, 2017, 046

Pappadopulo, D., Ruderman, J. T., & Trevisan, G. 2016, Phys. Rev. D, 94, 035005

Parimbelli, G., Scelfo, G., Giri, S., et al. 2021, JCAP, 2021, 044

Paszke, A., Gross, S., Massa, F., et al. 2019, Advances in Neural Information Processing Systems 32 (Curran Associates, Inc.), 8024

Peter, A. H. G., & Benson, A. J. 2010, Phys. Rev. D, 80, 123521

Peter, A. H. G., Moody, C. E., & Kamionkowski, M. 2010, Phys. Rev. D, 81, 103501

Planck Collaboration V. 2020, A&A, 641, A5

Planck Collaboration VI. 2020, A&A, 641, A6

Potter, D., Stadel, J., & Teyssier, R. 2017, Comput. Astrophys. Cosmol., 4, 2

Poulin, V., Bernal, J. L., Kovetz, E., & Kamionkowski, M. 2023, Phys. Rev. D, 107, 123538

Poursidou, A., Skordis, C., & Copeland, E. J. 2013, Phys. Rev. D, 88, 083505

Prince, H., & Dunkley, J. 2019, Phys. Rev. D, 100, 083502

Raveri, M., & Hu, W. 2019, Phys. Rev. D, 99, 043506

Rubira, H., Mazoun, A., & Garny, M. 2023, JCAP, 2023, 034

Schöneberg, N., Abellán, G. F., Sánchez, A. P., et al. 2022, Phys. Rep., 984, 1

Schneider, A., & Teyssier, R. 2015, JCAP, 1512, 049

Schneider, A., Teyssier, R., Stadel, J., et al. 2019, JCAP, 03, 020

Schneider, A., Refregier, A., Grandis, S., et al. 2020a, JCAP, 04, 020

Schneider, A., Stoira, N., Refregier, A., et al. 2020b, JCAP, 04, 019

Schneider, A., Giri, S. K., Amodeo, S., & Refregier, A. 2022, MNRAS, 514, 3802

Simon, T., Abellán, G. F., Du, P., Poulin, V., & Tsai, Y. 2022, Phys. Rev. D, 106, 023516

Sitzmann, V., Martel, J. N. P., Bergman, A. W., Lindell, D. B., & Wetzstein, G. 2020, arXiv e-prints [arXiv:2006.09661]

Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152

Tan, T., Zuercher, D., Fluri, J., et al. 2023, MNRAS, 522, 3766

Tanimura, H., Douspis, M., Aghanim, N., & Kuruvilla, J. 2023, A&A, 674, A222

Tarsitano, F., Schmitt, U., Refregier, A., et al. 2021, Astron. Comput., 36, 100484

The Dark Energy Survey Collaboration 2005, arXiv e-prints [arXiv:astro-ph/0510346]

van Daalen, M. P., Schaye, J., Booth, C. M., & Vecchia, C. D. 2011, MNRAS, 415, 3649

Vogelsberger, M., Zavala, J., Cyr-Racine, F.-Y., et al. 2016, MNRAS, 460, 1399

Wang, M.-Y., Croft, R. A. C., Peter, A. H. G., Zentner, A. R., & Purcell, C. W. 2013, Phys. Rev. D, 88, 123515

Wang, M.-Y., & Zentner, A. R. 2012, Phys. Rev. D, 85, 043514
# Appendix A: Comparing two-body decay simulations to previous work

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

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

mentation. Therefore, we did not perform model inferences for such extreme values of $v_k$ in this work.

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

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

The effects of ΛDDM in our study are accounted for by applying a cosmology-independent boost to the ΛCDM non-linear matter power spectrum, calculated using the revised halofit method within the CLASS code. To demonstrate the validity of this approach, we conducted a test suite of *N*-body simulations where we kept the ΛDDM parameters fixed at $\Gamma^{-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 second-order effect. In panels (a) to (e) of this figure, we show the impact of the Hubble parameter ($h_0$), clustering amplitude ($\sigma_8$), spectral index ($n_s$), matter abundance ($\Omega_m$), and baryon abundance ($\Omega_b$), respectively. We display the ΛDDM boost for the fiducial cosmology of our *N*-body simulations with a solid salmon line, while we plot the same quantity for the best-fit cosmology of the KiDS-450 survey (Hildebrandt et al. 2016) with a dashed blue line. We also added a few additional models, represented by the dotted grey lines.

In this analysis, we found that the choice of cosmology has a noticeable impact on the ΛDDM boost only for the parameters $\sigma_8$ (panel b) and $\Omega_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 $\sigma_8$ and $\Omega_m$, the observed effects are relatively small compared to the overall amplitude of the observed boost.
<figure><img src="image_10.png" /></figure>

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

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

## C.1. Weak lensing

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

<figure><img src="image_11.png" /></figure>

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

## C.2. Baryon feedback

We used three free parameters in the baryonic effects emulator ($\log_{10} M_c, \theta_{ej}, \eta_\delta$) 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 $\eta_\delta$. This is evident from Fig. C.2 and can be inferred from Tab. E.2, as varying $\eta_\delta$ does not affect the modelled WL observables. On the other hand, the gas profile parameters $\log_{10} M_c$ and $\theta_{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 WL-only analysis and $\log_{10} M_c > 13.8, (unconst)$ from the combined analysis. Similarly, for $\theta_{ej}$, we found $\theta_{ej} < 5.45, (5.57)$ in the case of ΛCDM (ΛDDM) from the WL-only analysis and $\theta_{ej} > 5.88, (unconst)$ from the combined analysis. These results indicate that the values of these baryonic parameters are mutually exclusive when comparing the WL-only and combined (ΛCDM) scenarios. The underlying reasons for this behaviour are discussed in the following sections.
<figure><img src="image_10.png" /></figure>

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

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

## C.1. Weak lensing

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

<figure><img src="image_11.png" /></figure>

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

## C.2. Baryon feedback

We used three free parameters in the baryonic effects emulator ($\log_{10} M_c, \theta_{ej}, \eta_\delta$) 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 $\eta_\delta$. This is evident from Fig. C.2 and can be inferred from Tab. E.2, as varying $\eta_\delta$ does not affect the modelled WL observables. On the other hand, the gas profile parameters $\log_{10} M_c$ and $\theta_{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 WL-only analysis and $\log_{10} M_c > 13.8, (unconst)$ from the combined analysis. Similarly, for $\theta_{ej}$, we found $\theta_{ej} < 5.45, (5.57)$ in the case of ΛCDM (ΛDDM) from the WL-only analysis and $\theta_{ej} > 5.88, (unconst)$ from the combined analysis. These results indicate that the values of these baryonic parameters are mutually exclusive when comparing the WL-only and combined (ΛCDM) scenarios. The underlying reasons for this behaviour are discussed in the following sections.
<figure><img src="image_12.png" /></figure>

**Fig. C.2.** Effects of baryons on band power spectra illustrated using different baryonic parameter configurations. Band power auto- and cross-correlations for five tomographic bins with respective error bars as measured by the KiDS-1000 survey are shown as black dots. Solid dark blue (green) lines represent the best fit in the case of ΛCDM (ΛDDM), and the cyan, magenta, and orange lines depict the model predictions for different modifications of the ΛCDM best fit, varying baryonic parameters. In the legend, we include the minimum $\chi^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$, $\theta_{ej}$, and $\eta_\delta$ 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 $\eta_\delta$ does not alter the shear signal. However, adjustments to the gas profile parameters log$_{10} M_c$ and $\theta_{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 \sim 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$.

<figure><img src="image_13.png" /></figure>

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

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 $\theta_{ej}$ in the ΛDDM model. In both panels, it is evident that fitting the WL-only data favours a weak suppression of the power spectrum, indicated by the green contours located in the bottom-left region. Conversely, when combining WL and CMB, strong suppression is required either from baryons or from DDM, as indicated by the purple contours. It is conceivable that having even stronger baryonic effects at our disposal could entirely replace the need for the DDM and provide a satisfactory fit to both datasets individually, though the physical justification for such high baryonic feedback remains an open question.

### 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 $\Gamma$ and $v_k$ (dashed lines). In the smaller bottom panels, we show the difference between all these predictions with the ΛCDM best-fit cosmology. We observed that the two-body decays do not change the CMB signal significantly, meaning that the scatter of the CMB data is larger than the observable effects of two-body decays.
<figure><img src="image_12.png" /></figure>

**Fig. C.2.** Effects of baryons on band power spectra illustrated using different baryonic parameter configurations. Band power auto- and cross-correlations for five tomographic bins with respective error bars as measured by the KiDS-1000 survey are shown as black dots. Solid dark blue (green) lines represent the best fit in the case of ΛCDM (ΛDDM), and the cyan, magenta, and orange lines depict the model predictions for different modifications of the ΛCDM best fit, varying baryonic parameters. In the legend, we include the minimum $\chi^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$, $\theta_{ej}$, and $\eta_\delta$ 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 $\eta_\delta$ does not alter the shear signal. However, adjustments to the gas profile parameters log$_{10} M_c$ and $\theta_{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 \sim 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$.

<figure><img src="image_13.png" /></figure>

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

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 $\theta_{ej}$ in the ΛDDM model. In both panels, it is evident that fitting the WL-only data favours a weak suppression of the power spectrum, indicated by the green contours located in the bottom-left region. Conversely, when combining WL and CMB, strong suppression is required either from baryons or from DDM, as indicated by the purple contours. It is conceivable that having even stronger baryonic effects at our disposal could entirely replace the need for the DDM and provide a satisfactory fit to both datasets individually, though the physical justification for such high baryonic feedback remains an open question.

### 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 $\Gamma$ and $v_k$ (dashed lines). In the smaller bottom panels, we show the difference between all these predictions with the ΛCDM best-fit cosmology. We observed that the two-body decays do not change the CMB signal significantly, meaning that the scatter of the CMB data is larger than the observable effects of two-body decays.
<figure><img src="image_14.png" /><figcaption>**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 χ² value that each configuration yields.</figcaption></figure>

## 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

$$ \Delta \text{AIC} = \Delta \chi_{\text{min}}^2 + 2(N_{\text{\Lambda DDM}} - N_{\text{\Lambda CDM}}), \quad (\text{D.1}) $$

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(-\Delta\text{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

$$ \Delta \text{AIC} < \Delta \text{AIC}_0 = -2.3. \quad (\text{D.2}) $$

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

<table><caption>**Table D.1.** Minimum χ² 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.</caption><thead><tr><th></th><th>WL</th><th>CMB</th><th>Combined</th><th>√Q<sub>DMAP</sub></th></tr></thead><tbody><tr><td>χ<sup>2</sup><sub>min</sub> (ΛCDM)</td><td>158.7</td><td>580.2</td><td>750.5</td><td>3.4σ</td></tr><tr><td>χ<sup>2</sup><sub>min</sub> (ΛDDM)</td><td>158.6</td><td>580.2</td><td>742.6</td><td>1.9σ</td></tr><tr><td>Δχ<sup>2</sup><sub>min</sub></td><td>-0.1</td><td>0.0</td><td>-7.9</td><td>-</td></tr><tr><td>ΔAIC</td><td>3.9</td><td>4.0</td><td>-3.9</td><td>-</td></tr></tbody></table>

fit to the data. In the combined scenario, ΔAIC = -3.9 < ΔAIC₀. This means that the two additional parameters in the combined analysis are efficient. However, the combined ΛDDM fit is still worse by Δχ² = 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 χ² as a minimum χ² among all sampled configurations in our MCMC chain. For a summary of combined MCMC analysis results, we refer to Tab. D.1.
<figure><img src="image_14.png" /><figcaption>**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 χ² value that each configuration yields.</figcaption></figure>

## 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

$$ \Delta \text{AIC} = \Delta \chi_{\text{min}}^2 + 2(N_{\text{\Lambda DDM}} - N_{\text{\Lambda CDM}}), \quad (\text{D.1}) $$

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(-\Delta\text{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

$$ \Delta \text{AIC} < \Delta \text{AIC}_0 = -2.3. \quad (\text{D.2}) $$

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

<table><caption>**Table D.1.** Minimum χ² 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.</caption><thead><tr><th></th><th>WL</th><th>CMB</th><th>Combined</th><th>√Q<sub>DMAP</sub></th></tr></thead><tbody><tr><td>χ<sup>2</sup><sub>min</sub> (ΛCDM)</td><td>158.7</td><td>580.2</td><td>750.5</td><td>3.4σ</td></tr><tr><td>χ<sup>2</sup><sub>min</sub> (ΛDDM)</td><td>158.6</td><td>580.2</td><td>742.6</td><td>1.9σ</td></tr><tr><td>Δχ<sup>2</sup><sub>min</sub></td><td>-0.1</td><td>0.0</td><td>-7.9</td><td>-</td></tr><tr><td>ΔAIC</td><td>3.9</td><td>4.0</td><td>-3.9</td><td>-</td></tr></tbody></table>

fit to the data. In the combined scenario, ΔAIC = -3.9 < ΔAIC₀. This means that the two additional parameters in the combined analysis are efficient. However, the combined ΛDDM fit is still worse by Δχ² = 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 χ² as a minimum χ² among all sampled configurations in our MCMC chain. For a summary of combined MCMC analysis results, we refer to Tab. D.1.
## Appendix E: MCMC results

In this section, we provide more details about the parameters resulting from our MCMC analysis. In Tab. E.1, we report the findings from our ΛCDM analysis (taken from Bucko et al. 2023), while in Tab. E.2, we summarise the actual ΛDDM model results. In both cases, we provide separate constraints from the WL, CMB, and combined analyses. In the upper part of the

tables, we show the posteriors of the parameters sampled by in our MCMC runs, showing mean (best fit) values with corresponding upper and lower deviations. The dashed fields in the tables indicate the parameters were not present in a given MCMC setup, 'unconst' is used when parameter constraints could not be obtained from MCMC analysis, and we state no uncertainties whenever referring to a value constant throughout the inference.<fnref n="6" />

<table><caption>**Table E.1.** Summary of our MCMC analysis assuming the ΛCDM model. We report individual results separately based on WL (KiDS-1000) 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 $\chi^2$ values (bottom).</caption><thead><tr><th>Parameter</th><th>WL ΛCDM<br>68% limits</th><th>CMB ΛCDM<br>68% limits</th><th>WL + CMB ΛCDM<br>68% limits</th></tr></thead><tbody><tr><td>ω<sub>dm</sub></td><td>0.146(0.081)<sup>+0.034</sup><sub>-0.056</sub></td><td>0.1208(0.1200) ± 0.0014</td><td>0.1182(0.1188) ± 0.0012</td></tr><tr><td>ω<sub>b</sub></td><td>unconst(0.02455)</td><td>0.02231(0.02236) ± 0.00015</td><td>0.02248(0.02253) ± 0.00013</td></tr><tr><td>ln(10<sup>10</sup>A<sub>s</sub>)</td><td>2.56(3.50)<sup>+0.61</sup><sub>-0.80</sub></td><td>3.050(3.110) ± 0.017</td><td>3.039(3.069) ± 0.017</td></tr><tr><td>h<sub>0</sub></td><td>unconst(0.6041)</td><td>0.6700(0.6734) ± 0.0061</td><td>0.6815(0.6799)<sup>+0.0050</sup><sub>-0.0057</sub></td></tr><tr><td>n<sub>s</sub></td><td>unconst(0.9377)</td><td>0.9622(0.9640) ± 0.0043</td><td>0.9678(0.9668) ± 0.0041</td></tr><tr><td>τ<sub>reio</sub></td><td>—</td><td>0.0566(0.0841) ± 0.0083</td><td>0.0540(0.0677) ± 0.0077</td></tr><tr><td>A<sub>IA</sub></td><td>0.75(0.89)<sup>+0.33</sup><sub>-0.38</sub></td><td>—</td><td>0.63(0.44)<sup>+0.25</sup><sub>-0.31</sub></td></tr><tr><td>A<sub>planck</sub></td><td>—</td><td>1.0005(1.0031) ± 0.0025</td><td>1.0003(1.0005) ± 0.0025</td></tr><tr><td>log<sub>10</sub> M<sub>c</sub></td><td>&lt; 13.1(12.6)</td><td>—</td><td>&gt; 13.8(15.0)</td></tr><tr><td>θ<sub>ej</sub></td><td>&lt; 5.45(2.23)</td><td>—</td><td>&gt; 5.88(7.74)</td></tr><tr><td>η<sub>δ</sub></td><td>unconst(0.21)</td><td>—</td><td>unconst(0.13)</td></tr><tr><td>log<sub>10</sub>(Γ × Gyr)</td><td>—</td><td>—</td><td>—</td></tr><tr><td>log<sub>10</sub>(v<sub>k</sub> × s/km)</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Ω<sub>m</sub></td><td>0.347<sup>+0.066</sup><sub>-0.110</sub></td><td>0.3189 ± 0.0087</td><td>0.3043 ± 0.0070</td></tr><tr><td>σ<sub>8</sub></td><td>0.70<sup>+0.11</sup><sub>-0.13</sub></td><td>0.8154 ± 0.0080</td><td>0.8029 ± 0.0073</td></tr><tr><td>S<sub>8</sub></td><td>0.735<sup>+0.031</sup><sub>-0.024</sub></td><td>0.841 ± 0.017</td><td>0.809 ± 0.014</td></tr><tr><td>ln(prior)</td><td>5.99</td><td>7.9(0.8)<sup>+1.2</sup><sub>-0.4</sub></td><td>13.9(12.9)<sup>+1.0</sup><sub>-0.3</sub></td></tr><tr><td>ln(lik<sub>WL</sub>)</td><td>-81.3(-79.3)<sup>+1.2</sup><sub>-0.5</sub></td><td>—</td><td>-83.5(-83.3)<sup>+1.7</sup><sub>-1.0</sub></td></tr><tr><td>ln(lik<sub>CMB</sub>)</td><td>—</td><td>-294.4(-290.1)<sup>+2.0</sup><sub>-1.1</sub></td><td>-295.7(-291.9)<sup>+2.8</sup><sub>-1.6</sub></td></tr><tr><td>χ<sup>2</sup><sub>min</sub></td><td>158.7</td><td>580.2</td><td>750.5</td></tr></tbody></table>

<footnote n="6">⁶ adapted from Bucko et al. (2023)</footnote>
## 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.<fnref n="6" />

<table><caption>**Table E.1.** Summary of our MCMC analysis assuming the ΛCDM model. We report individual results separately based on WL (KiDS-1000) 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 $\chi^2$ values (bottom).</caption><thead><tr><th>Parameter</th><th>WL ΛCDM<br>68% limits</th><th>CMB ΛCDM<br>68% limits</th><th>WL + CMB ΛCDM<br>68% limits</th></tr></thead><tbody><tr><td>ω<sub>dm</sub></td><td>0.146(0.081)<sup>+0.034</sup><sub>-0.056</sub></td><td>0.1208(0.1200) ± 0.0014</td><td>0.1182(0.1188) ± 0.0012</td></tr><tr><td>ω<sub>b</sub></td><td>unconst(0.02455)</td><td>0.02231(0.02236) ± 0.00015</td><td>0.02248(0.02253) ± 0.00013</td></tr><tr><td>ln(10<sup>10</sup>A<sub>s</sub>)</td><td>2.56(3.50)<sup>+0.61</sup><sub>-0.80</sub></td><td>3.050(3.110) ± 0.017</td><td>3.039(3.069) ± 0.017</td></tr><tr><td>h<sub>0</sub></td><td>unconst(0.6041)</td><td>0.6700(0.6734) ± 0.0061</td><td>0.6815(0.6799)<sup>+0.0050</sup><sub>-0.0057</sub></td></tr><tr><td>n<sub>s</sub></td><td>unconst(0.9377)</td><td>0.9622(0.9640) ± 0.0043</td><td>0.9678(0.9668) ± 0.0041</td></tr><tr><td>τ<sub>reio</sub></td><td>—</td><td>0.0566(0.0841) ± 0.0083</td><td>0.0540(0.0677) ± 0.0077</td></tr><tr><td>A<sub>IA</sub></td><td>0.75(0.89)<sup>+0.33</sup><sub>-0.38</sub></td><td>—</td><td>0.63(0.44)<sup>+0.25</sup><sub>-0.31</sub></td></tr><tr><td>A<sub>planck</sub></td><td>—</td><td>1.0005(1.0031) ± 0.0025</td><td>1.0003(1.0005) ± 0.0025</td></tr><tr><td>log<sub>10</sub> M<sub>c</sub></td><td>&lt; 13.1(12.6)</td><td>—</td><td>&gt; 13.8(15.0)</td></tr><tr><td>θ<sub>ej</sub></td><td>&lt; 5.45(2.23)</td><td>—</td><td>&gt; 5.88(7.74)</td></tr><tr><td>η<sub>δ</sub></td><td>unconst(0.21)</td><td>—</td><td>unconst(0.13)</td></tr><tr><td>log<sub>10</sub>(Γ × Gyr)</td><td>—</td><td>—</td><td>—</td></tr><tr><td>log<sub>10</sub>(v<sub>k</sub> × s/km)</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Ω<sub>m</sub></td><td>0.347<sup>+0.066</sup><sub>-0.110</sub></td><td>0.3189 ± 0.0087</td><td>0.3043 ± 0.0070</td></tr><tr><td>σ<sub>8</sub></td><td>0.70<sup>+0.11</sup><sub>-0.13</sub></td><td>0.8154 ± 0.0080</td><td>0.8029 ± 0.0073</td></tr><tr><td>S<sub>8</sub></td><td>0.735<sup>+0.031</sup><sub>-0.024</sub></td><td>0.841 ± 0.017</td><td>0.809 ± 0.014</td></tr><tr><td>ln(prior)</td><td>5.99</td><td>7.9(0.8)<sup>+1.2</sup><sub>-0.4</sub></td><td>13.9(12.9)<sup>+1.0</sup><sub>-0.3</sub></td></tr><tr><td>ln(lik<sub>WL</sub>)</td><td>-81.3(-79.3)<sup>+1.2</sup><sub>-0.5</sub></td><td>—</td><td>-83.5(-83.3)<sup>+1.7</sup><sub>-1.0</sub></td></tr><tr><td>ln(lik<sub>CMB</sub>)</td><td>—</td><td>-294.4(-290.1)<sup>+2.0</sup><sub>-1.1</sub></td><td>-295.7(-291.9)<sup>+2.8</sup><sub>-1.6</sub></td></tr><tr><td>χ<sup>2</sup><sub>min</sub></td><td>158.7</td><td>580.2</td><td>750.5</td></tr></tbody></table>

<footnote n="6">⁶ adapted from Bucko et al. (2023)</footnote>
<table><caption>**Table E.2.** Same as Tab. E.1 but assuming the ΛDDM model.</caption><thead><tr><th>Parameter</th><th>WL ΛDDM<br>68% limits</th><th>CMB ΛDDM<br>68% limits</th><th>WL + CMB ΛDDM<br>68% limits</th></tr></thead><tbody><tr><td>ω<sub>dm</sub></td><td>0.139(0.106)<sup>+0.039</sup><sub>-0.055</sub></td><td>0.1208(0.1190) ± 0.0014</td><td>0.1191(0.1193) ± 0.0015</td></tr><tr><td>ω<sub>b</sub></td><td>unconst(0.02103)</td><td>0.02231(0.02242) ± 0.00015</td><td>0.02241(0.02245) ± 0.00016</td></tr><tr><td>ln(10<sup>10</sup>A<sub>s</sub>)</td><td>2.69(2.93)<sup>+0.56</sup><sub>-0.94</sub></td><td>3.052(3.116) ± 0.017</td><td>3.045(3.104) ± 0.017</td></tr><tr><td>h<sub>0</sub></td><td>unconst(0.6093)</td><td>0.6703(0.6785) ± 0.0062</td><td>0.6772(0.6766) ± 0.0067</td></tr><tr><td>n<sub>s</sub></td><td>unconst(0.9295)</td><td>0.9620(0.9667) ± 0.0044</td><td>0.9653(0.9665) ± 0.0046</td></tr><tr><td>τ</td><td>—</td><td>0.0571(0.0913) ± 0.0082</td><td>0.0557(0.0823) ± 0.0085</td></tr><tr><td>A<sub>IA</sub></td><td>0.77(0.86) ± 0.33</td><td>—</td><td>0.68(0.90)<sup>+0.29</sup><sub>-0.33</sub></td></tr><tr><td>A<sub>planck</sub></td><td>—</td><td>1.0004(1.0007) ± 0.0025</td><td>1.0005(1.0027) ± 0.0024</td></tr><tr><td>log<sub>10</sub> M<sub>c</sub></td><td>&lt; 13.2(13.1)</td><td>—</td><td>unconst(12.4)</td></tr><tr><td>θ<sub>ej</sub></td><td>&lt; 5.57(2.08)</td><td>—</td><td>unconst(2.67)</td></tr><tr><td>η<sub>δ</sub></td><td>unconst(0.11)</td><td>—</td><td>&lt; 0.27(0.27)</td></tr><tr><td>log<sub>10</sub>(Γ × Gyr)</td><td>unconst(−2.12)</td><td>unconst(−2.26)</td><td>−2.25(−1.77)<sup>+0.74</sup><sub>−0.23</sub></td></tr><tr><td>log<sub>10</sub>(v<sub>k</sub> × s/km)</td><td>unconst(0.74)</td><td>unconst(2.51)</td><td>&gt; 2.80(3.47)</td></tr><tr><td>Ω<sub>m</sub></td><td>0.338(0.341)<sup>+0.077</sup><sub>-0.11</sub></td><td>0.3200(0.3084) ± 0.0087</td><td>0.3100(0.3101)<sup>+0.0084</sup><sub>-0.0095</sub></td></tr><tr><td>σ<sub>8</sub></td><td>0.720(0.693)<sup>+0.094</sup><sub>-0.150</sub></td><td>0.812(0.837)<sup>+0.012</sup><sub>-0.0064</sub></td><td>0.776(0.752)<sup>+0.037</sup><sub>-0.031</sub></td></tr><tr><td>S<sub>8</sub></td><td>0.740(0.738) ± 0.033</td><td>0.839(0.849)<sup>+0.021</sup><sub>-0.016</sub></td><td>0.788(0.765)<sup>+0.028</sup><sub>-0.024</sub></td></tr><tr><td>ln(prior)</td><td>3.6</td><td>14.1(4.1)<sup>+1.3</sup><sub>-0.4</sub></td><td>11.4(5.2)<sup>+1.2</sup><sub>-0.4</sub></td></tr><tr><td>ln(lik<sub>WL</sub>)</td><td>−81.3(−79.3)<sup>+1.3</sup><sub>-0.5</sub></td><td>—</td><td>−82.6(−80.5)<sup>+2.0</sup><sub>-1.1</sub></td></tr><tr><td>ln(lik<sub>CMB</sub>)</td><td>—</td><td>−294.3(−290.1)<sup>+2.0</sup><sub>-1.1</sub></td><td>−295.5(−290.8)<sup>+2.5</sup><sub>-1.3</sub></td></tr><tr><td>χ<sup>2</sup><sub>min</sub></td><td>158.6</td><td>580.2</td><td>742.6</td></tr></tbody></table>
<table><caption>**Table E.2.** Same as Tab. E.1 but assuming the ΛDDM model.</caption><thead><tr><th>Parameter</th><th>WL ΛDDM<br>68% limits</th><th>CMB ΛDDM<br>68% limits</th><th>WL + CMB ΛDDM<br>68% limits</th></tr></thead><tbody><tr><td>ω<sub>dm</sub></td><td>0.139(0.106)<sup>+0.039</sup><sub>-0.055</sub></td><td>0.1208(0.1190) ± 0.0014</td><td>0.1191(0.1193) ± 0.0015</td></tr><tr><td>ω<sub>b</sub></td><td>unconst(0.02103)</td><td>0.02231(0.02242) ± 0.00015</td><td>0.02241(0.02245) ± 0.00016</td></tr><tr><td>ln(10<sup>10</sup>A<sub>s</sub>)</td><td>2.69(2.93)<sup>+0.56</sup><sub>-0.94</sub></td><td>3.052(3.116) ± 0.017</td><td>3.045(3.104) ± 0.017</td></tr><tr><td>h<sub>0</sub></td><td>unconst(0.6093)</td><td>0.6703(0.6785) ± 0.0062</td><td>0.6772(0.6766) ± 0.0067</td></tr><tr><td>n<sub>s</sub></td><td>unconst(0.9295)</td><td>0.9620(0.9667) ± 0.0044</td><td>0.9653(0.9665) ± 0.0046</td></tr><tr><td>τ</td><td>—</td><td>0.0571(0.0913) ± 0.0082</td><td>0.0557(0.0823) ± 0.0085</td></tr><tr><td>A<sub>IA</sub></td><td>0.77(0.86) ± 0.33</td><td>—</td><td>0.68(0.90)<sup>+0.29</sup><sub>-0.33</sub></td></tr><tr><td>A<sub>planck</sub></td><td>—</td><td>1.0004(1.0007) ± 0.0025</td><td>1.0005(1.0027) ± 0.0024</td></tr><tr><td>log<sub>10</sub> M<sub>c</sub></td><td>&lt; 13.2(13.1)</td><td>—</td><td>unconst(12.4)</td></tr><tr><td>θ<sub>ej</sub></td><td>&lt; 5.57(2.08)</td><td>—</td><td>unconst(2.67)</td></tr><tr><td>η<sub>δ</sub></td><td>unconst(0.11)</td><td>—</td><td>&lt; 0.27(0.27)</td></tr><tr><td>log<sub>10</sub>(Γ × Gyr)</td><td>unconst(-2.12)</td><td>unconst(-2.26)</td><td>-2.25(-1.77)<sup>+0.74</sup><sub>-0.23</sub></td></tr><tr><td>log<sub>10</sub>(v<sub>k</sub> × s/km)</td><td>unconst(0.74)</td><td>unconst(2.51)</td><td>&gt; 2.80(3.47)</td></tr><tr><td>Ω<sub>m</sub></td><td>0.338(0.341)<sup>+0.077</sup><sub>-0.11</sub></td><td>0.3200(0.3084) ± 0.0087</td><td>0.3100(0.3101)<sup>+0.0084</sup><sub>-0.0095</sub></td></tr><tr><td>σ<sub>8</sub></td><td>0.720(0.693)<sup>+0.094</sup><sub>-0.150</sub></td><td>0.812(0.837)<sup>+0.012</sup><sub>-0.0064</sub></td><td>0.776(0.752)<sup>+0.037</sup><sub>-0.031</sub></td></tr><tr><td>S<sub>8</sub></td><td>0.740(0.738) ± 0.033</td><td>0.839(0.849)<sup>+0.021</sup><sub>-0.016</sub></td><td>0.788(0.765)<sup>+0.028</sup><sub>-0.024</sub></td></tr><tr><td>ln(prior)</td><td>3.6</td><td>14.1(4.1)<sup>+1.3</sup><sub>-0.4</sub></td><td>11.4(5.2)<sup>+1.2</sup><sub>-0.4</sub></td></tr><tr><td>ln(lik<sub>WL</sub>)</td><td>-81.3(-79.3)<sup>+1.3</sup><sub>-0.5</sub></td><td>—</td><td>-82.6(-80.5)<sup>+2.0</sup><sub>-1.1</sub></td></tr><tr><td>ln(lik<sub>CMB</sub>)</td><td>—</td><td>-294.3(-290.1)<sup>+2.0</sup><sub>-1.1</sub></td><td>-295.5(-290.8)<sup>+2.5</sup><sub>-1.3</sub></td></tr><tr><td>χ<sup>2</sup><sub>min</sub></td><td>158.6</td><td>580.2</td><td>742.6</td></tr></tbody></table>