# Can decaying dark matter models be a solution to the S_8 tension from the thermal Sunyaev-Zel'dovich effect?

Hideki Tanimura, Marian Douspis, Nabila Aghanim, Joseph Kuruvilla

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

## To cite this version:

Hideki Tanimura, Marian Douspis, Nabila Aghanim, Joseph Kuruvilla. Can decaying dark matter models be a solution to the S_8 tension from the thermal Sunyaev-Zel'dovich effect?. Astronomy and Astrophysics - A&A, 2023, 674, pp.A222. 10.1051/0004-6361/202345882. hal-03952450

HAL Id: hal-03952450

https://hal.science/hal-03952450

Submitted on 26 Oct 2023

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>
# Testing decaying dark matter models as a solution to the $S_8$ tension with the thermal Sunyaev-Zel'dovich effect

Hideki Tanimura<fnref n="1" />,<fnref n="2" />✉, Marian Douspis<fnref n="1" />, Nabila Aghanim<fnref n="1" />, and Joseph Kuruvilla<fnref n="1" />

Received 10 January 2023 / Accepted 5 May 2023

## ABSTRACT

Considering possible solutions to the $S_8$ tension between the Planck cosmic microwave background (CMB) measurement and low-redshift probes, we extended the standard $\Lambda$CDM cosmological model by including decay of dark matter (DDM). We first tested the DDM model in which dark matter decays into a form of noninteracting dark radiation. Under this DDM model, we investigated the impacts of DDM on the Sunyaev Zel'dovich (SZ) effect by varying the decay lifetime, $\Gamma^{-1}$, including the background evolution in cosmology and the nonlinear prescription in the halo mass function. We performed a cosmological analysis under the assumption of this extended cosmological model by combining the latest high-redshift Planck CMB measurement and low-redshift measurements of the SZ power spectrum as well as the baryonic acoustic oscillations (BAO) and luminosity distances to type Ia supernovae (SNIa). Our result shows a preference for $\Gamma^{-1} \sim 220$ Gyr with a lower bound on the decay lifetime of $\sim 38$ Gyr at 95% confidence level. Additionally, we tested the other DDM model in which dark matter decays into warm dark matter and dark radiation. This model supports $\Gamma^{-1} \sim 137$ Gyr to resolve the $S_8$ tension with a lower bound on the decay lifetime of $\sim 24$ Gyr at 95% confidence level. Comparing these two models, we find that the second leads to slightly better reconciliation of the $S_8$ tension.

**Key words.** cosmic background radiation – large-scale structure of Universe – dark matter – galaxies: clusters: general – galaxies: clusters: intracluster medium

## 1. Introduction

The current standard cosmological model, called the $\Lambda$Cold dark matter (ΛCDM) model, has been supported by many observational results at different epochs, such as cosmic microwave background (CMB), big bang nucleosynthesis (BBN), and luminosity distances to type Ia supernovae (SNIa). However, as measurement precision has increased over the past few years, the validity of the $\Lambda$CDM model is under reexamination, particularly due to a discrepancy in the $S_8(\equiv \sigma_8(\Omega_m/0.3))^{0.5}$ parameter, which is the amplitude of matter density fluctuations, $\sigma_8$, scaled by the square root of the matter density, $\Omega_m$. For example, Planck measured $S_8 = 0.830 \pm 0.013$ with the CMB anisotropies ($z \sim 1100$; Planck Collaboration VI 2020), but low-redshift ($z \sim 0-1$) cosmological probes, such as the abundance of galaxy clusters by Planck (Planck Collaboration XX 2014; Planck Collaboration XXIV 2016) and gravitational lensing by the Kilo Degree Survey (KiDS; Heymans et al. 2021) and Dark Energy Survey (DES; DES Collaboration 2022), show a preference for a lower value of $S_8 \sim 0.76-0.78$, representing a tension of up to $\sim 2-3\sigma$ with the CMB constraints, referred to as the $S_8$ tension. This tension indicates that the growth rate of the cosmic structure is less than expected from the CMB measurements and may demand modifications to the standard model.

Many attempts have been made to solve the $S_8$ tension by breaking assumptions in the standard cosmological model, mainly on the nature of neutrinos (e.g., Salvati et al. 2018; Tröster et al. 2021), dark energy (e.g., Bolliet et al. 2018; Salvati et al. 2018; Lambiase et al. 2019; Di Valentino et al. 2020a,b; Tröster et al. 2021), modified gravity (e.g., Ilić et al. 2019), and others (e.g., Nunes & Vagnozzi 2021). However, none of the proposed solutions have gained wide acceptance so far. One possible solution is the decaying dark matter (DDM) model (e.g., Doroshkevich et al. 1989; Aoyama et al. 2014; Audren et al. 2014; Enqvist et al. 2015, 2020; Poulin et al. 2016; Xiao et al. 2020; Abellán et al. 2021; Chen et al. 2021; Hubert et al. 2021; Holm et al. 2023; Fuß & Garny 2022; Abellán et al. 2022; Mau et al. 2022; McCarthy & Hill 2022; Simon et al. 2022; Bucko et al. 2023). The DDM model has two features that reduce the cosmic structure growth and could naturally explain the $S_8$ tension. First, the cosmic structure growth is reduced by the decay of massive dark matter (DM) into lighter or massless daughter particles. Second, the DDM model can only change the structure formation scenario at low redshift, which is probed by clusters of galaxies, lensing, and so on, without modifying the scenario at high redshift – which can be probed by CMB – before dark matter decays.

While several scenarios can be considered for the decay modes, we consider two scenarios: one scenario in which DM decays into a form of noninteracting dark radiation (DR), hereafter referred to as the $\Lambda$DDM1 model (e.g., Audren et al. 2014; Enqvist et al. 2015, 2020; Poulin et al. 2016) and another scenario, where DM decays into one massless DR component and one massive warm DM (WDM) particle (hereafter $\Lambda$DDM2 model) that interacts only through gravity with the standard model particles (e.g., Aoyama et al. 2014; Abellán et al. 2021, 2022). The $\Lambda$DDM1 model can be parameterized with the decay rate, $\Gamma$ ($\Gamma^{-1}$ represents the lifetime of DM), in addition to the cosmological parameters under the standard $\Lambda$CDM model. The DDM2 model can be parameterized with two parameters:

<footnote n="1">¹ Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, Bâtiment 121, 91405 Orsay, France</footnote>

<footnote n="2">² Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
e-mail: hideki.tanimura@ipmu.jp</footnote>
Γ, and ε, which represents the mass-energy fraction transferred to the massless component of DR (Aoyama et al. 2014; Abellán et al. 2021, 2022). Many studies have been performed for the simplest scenario, that is, the ADDM1 model, and tight constraints have been put on the lifetime of the DM particle of Γ⁻¹ > 160 Gyr (e.g., Audren et al. 2014; Enqvist et al. 2015, 2020; Poulin et al. 2016), which is much greater than the current age of the Universe. Fewer studies have considered the extended DDM models, such as the ADDM2 scenario. For example, Abellán et al. (2022) showed that the ADDM2 can fully explain the low-S₈ measurement with a shorter lifetime of Γ⁻¹ ≈ 55 (ε/0.007)¹.⁴ Gyr using the Planck CMB measurement (Planck Collaboration VI 2020) combined with the baryonic acoustic oscillations (BAO) measurements from the 6-degree Field Galaxy Survey (6dF; Beutler et al. 2011), the Sloan Digital Sky Survey (SDSS; Ross et al. 2015), the Baryon Oscillation Spectroscopic Survey (BOSS; Alam et al. 2017), and the Extended Baryon Oscillation Spectroscopic Survey (eBOSS; Blomqvist 2019; de Sainte Agathe & Balland 2019), including the SNIa measurement from Pantheon (Scolnic et al. 2018). A further study was performed by Abellán et al. (2021), who included KiDS and DES weak-lensing data (Abbott et al. 2018; Joudaki et al. 2020; Heymans et al. 2021), with the authors concluding that there is a preference for the ADDM2 model over the standard ΛCDM model when KiDS and DES weak-lensing data are included, without degrading the fit to other cosmological datasets such as CMB, BAO, and SNIa measurements.

In this paper, we use the thermal Sunyaev-Zel'dovich (SZ) effect (Sunyaev & Zeldovich 1970, 1972) as a probe for the S8 tension. To investigate the S8 parameter, the SZ effect, which is caused by the inverse Compton scattering of CMB photons by hot electrons along the line of sight, is a useful probe because the SZ power spectrum, Cℓ SZ, is sensitive to the σ8 and Ωm cosmological parameters as Cℓ SZ ∝ σ8 Ωm 3 (Planck Collaboration XXII 2016; Salvati et al. 2018). The S8 tension was first revealed by cosmological analysis of the SZ cluster number counts from Planck (Planck Collaboration XX 2014). Recent SZ results (Planck Collaboration XXII 2016; Tanimura et al. 2022) showed a similar S8 value of S8 = 0.764 +0.015 -0.018 to the weak-lensing and cluster-count observations, confirming a tension with the Planck CMB result. As a solution for the S8 tension, we assume that the origin of the tension is caused by the DDM and constrain the cosmological parameters and the DDM models using the latest SZ power spectrum measurement in Tanimura et al. (2022; hereafter T22). In our cosmological analyses, this SZ measurement is combined with other low-redshift probes that constrain the expansion history of the Universe, such as the BAO measurements from 6dFGS at z = 0.106 (Beutler et al. 2011), SDSS DR7 at z = 0.15 (Ross et al. 2015), and BOSS DR12 at z = 0.38, 0.51, 0.61 Alam et al. (2017) as well as the Pantheon SNIa catalog (Scolnic et al. 2018) (BAO+SNIa). The present paper is structured as follows. Section 2 describes the model of the SZ power spectrum, including the DDM model. Section 3 presents the constraint we place on the DDM model as a result of our cosmological analysis. We finally end with the conclusions of this study in Sect. 4.

## 2. Theoretical background

In this section, we first describe the theoretical prescription of the SZ effect, in particular the SZ angular power spectrum under the ΛCDM model by following T22, and then modify it to include the DDM models.

### 2.1. Compton y parameter

The Compton y parameter is proportional to the line-of-sight integral of electron pressure, P_e = n_e k_B T_e, where n_e is the physical electron number density, k_B is the Boltzmann constant, and T_e is the electron temperature. In an angular direction, $\hat{n}$, this Compton parameter is expressed as

$$y(\hat{n}) = \frac{\sigma_T}{m_e c^2} \int P_e(\hat{n}) d\ln(\hat{n}), \quad (1)$$

where σ_T is the Thomson cross section, m_e is the mass of one electron, c is the speed of light, and l is the physical distance. The change to the CMB temperature by the SZ effect, ΔT, at frequency v is given by

$$\frac{\Delta T}{T_{\text{CMB}}}(v, \hat{n}) = f(x) y(\hat{n}), \quad (2)$$

where T_CMB is the CMB temperature. The frequency dependence of the SZ effect is included in the pre-factor f(x) as

$$f(x) = x \coth\left(\frac{x}{2}\right) - 4 \left( x = \frac{hv}{k_B T_{\text{CMB}}} \right) \quad (3)$$

in the thermodynamic temperature unit, where h is the Planck constant.

### 2.2. SZ angular power spectrum

The SZ power spectrum can be modeled with a halo model (Cooray & Sheth 2002). The halo model consists of a “one-halo term”, which accounts for the correlation arising within an individual halo, and a “two-halo term”, which accounts for the correlation arising due to the environment surrounding a halo (Komatsu et al. 2002; Cooray & Sheth 2002). In this paper, we consider only the one-halo term because the contribution from the two-halo term to the total SZ power spectrum is minor (Komatsu & Kitayama 1999) at the scales we consider (ℓ > 60). The SZ power spectrum is then given by

$$C_{\ell}^{\text{SZ}} = \int dz \frac{d^2 V}{dz d\Omega} \int dM \frac{dn(M, z)}{dM} |\bar{y}_{\ell}(M, z)|^2, \quad (4)$$

where d²V/dzdΩ is the comoving volume element per redshift per steradian, n(M, z) is the comoving number density of halos of mass M and redshift z, called the halo mass function (HMF), and $\bar{y}_\ell(M, z)$ is the 2D Fourier transform of the y-profile of a halo, given by

$$\bar{y}_{\ell}(M, z) = \frac{\sigma_{\text{T}}}{m_e c^2} \frac{4\pi r_s}{\ell_s^2} \int dx_r x_r^2 \frac{\sin(\ell x_r / \ell_s)}{\ell x_r / \ell_s} P_e(x, M, z), \quad (5)$$

where r_s is the characteristic scale radius of the pressure profile and x_r = r/r_s is the dimensionless radial scale. d_A is the angular diameter distance and ℓ_s = d_A/r_s is the associated multipole moment. We integrate the contribution of halos in the redshift range from 0 to 3 and the mass range from 10¹³ M☉ to 5×10¹⁵ M☉ as well as the scaled radial distance, x_r, in the range from 0 to 5, following T22.

### 2.3. Halo mass function

We use the mass function from Tinker et al. (2008), as in T22, in which the number of halos per unit volume is given by

$$\frac{dn}{dM} = f(\sigma) \frac{\rho_{m,0}}{M} \frac{d\ln \sigma^{-1}}{dM}, \quad (6)$$
<figure><img src="image_3.png" /><figcaption>Fig. 1. SZ power spectra with the ΛDDM1 model (DM → DR) with different decay lifetimes of Γ⁻<fnref n="1" /> = 30 (dotted line), 100 (dash-dot line), and 300 (dashed line) Gyr. “ΛCDM” represents the ΛCDM model with stable cold dark matter (solid line). Here, the Planck 2018 cosmological parameters are assumed in all cases.</figcaption></figure>

where $\rho_{m,0}$ is the matter density at $z=0$ and $f(\sigma)$ is given by

$$f(\sigma) = A \left[ 1 + \left( \frac{\sigma}{b} \right)^{-a} \right] \exp \left( - \frac{c}{\sigma^2} \right), \qquad (7)$$

where A, a, b, c are constants calibrated with simulations in Tinker et al. (2008). Here, $\sigma$ is the standard deviation of density perturbations in a sphere of radius $R = (3M/4\pi\rho_{m,0})^{1/3}$, and is given by

$$\sigma^2 = \frac{1}{2\pi^2} \int dk k^2 P_m(k, z) |W(kR)|^2, \qquad (8)$$

where $W(kR)$ is the window function of a spherical top hat of radius $R$.

## 2.4. Universal pressure profile

For the electron pressure profile, we use the model from Planck Collaboration V (2013), which is used as a fiducial model in T22, and adopts the “universal” pressure profile (UPP; Nagai et al. 2007), which is a form of the generalized Navarro et al. (NFW; 1997) profile (gNFW),

$$\mathbb{P}(x_{500}) = \frac{P_0}{(c_{500}x_{500})^{\gamma}[1 + (c_{500}x_{500})^{\alpha}]^{[\beta-\alpha]/\gamma}}. \qquad (9)$$

Here, $x_{500} = r/R_{500}$ and we remind the reader that $R_{500}$ denotes 500 times the critical density. The model is defined by the following parameters: $P_0$, a normalization; $c_{500}$, a concentration parameter defined at a characteristic radius $R_{500}$; and the slopes in the central ($x_{500} \ll 1/c_{500}$), intermediate ($x_{500} \sim 1/c_{500}$), and outer regions ($x_{500} \gg 1/c_{500}$), which are given by $\gamma$, $\alpha$, and $\beta$, respectively. The scaled pressure profile for a halo with $M_{500}$ and $z$ is

$$\frac{P(r)}{P_{500}} = \mathbb{P}(x_{500}), \qquad (10)$$

with

$$P_{500} = 1.65 \times 10^{-3} \left[ \frac{H(z)}{H_0} \right]^{8/3} \times \left[ \frac{(1-b) M_{500}}{3 \times 10^{14} (h/0.7)^{-1} M_{\odot}} \right]^{2/3+\alpha_p} \left( \frac{h}{0.7} \right)^2 \text{keV cm}^{-3}, \qquad (11)$$

where $H(z)$ is the Hubble parameter at redshift $z$ and $H_0 = 100$ h km s⁻<fnref n="1" /> Mpc⁻<fnref n="1" /> is the present value. $P_{500}$ is the characteristic pressure reflecting the mass variation expected in a self-similar model of pressure evolution when $\alpha_p = 0$, purely based on gravitation (Arnaud et al. 2010). Deviation from the self-similarity appears in a variation of the scaled pressure profile, given by

$$\frac{P(r)}{P_{500}} = \mathbb{P}(x) \left[ \frac{(1-b) M_{500}}{3 \times 10^{14} (h/0.7)^{-1} M_{\odot}} \right]^{\alpha_p}, \qquad (12)$$

expressed as a function of $M_{500}$. Here, b is the hydrostatic mass bias. We note that we use the pressure model in Planck Collaboration V (2013) in which $M_{500}$ alternatively corresponds to the hydrostatic mass, $(1-b)M_{500}$. For the mass bias, we adopt $(1-b) = 0.780 \pm 0.092$ derived from the Canadian Cluster Comparison Project (Hoekstra et al. 2015; CCCP), which is consistent with most of the results from hydrodynamic simulations (Gianfagna et al. 2021). For the parameters of the generalized NFW electron pressure profile, we adopt the best-fit values of $[P_0, c_{500}, \gamma, \alpha, \beta] = [6.41, 1.81, 0.31, 1.33, 4.13]$, estimated from 62 massive nearby clusters ($10^{14.4} < M_{500} < 10^{15.3} M_{\odot}$) using the Planck SZ and XMM-Newton X-ray data in Planck Collaboration V (2013). We also adopt $\alpha_p = 0.12$ from the measurement in Arnaud et al. (2010).

## 2.5. Decaying dark matter model implementation

We include the DDM model based on the public code provided by Abellán et al. (2021; hereafter G21)<fnref n="1" />, in which the DDM model is implemented in the modified version of the Boltzmann code, CLASS (Blas et al. 2011; Lesgourgues & Tram 2011). In this code, the 1-body (ΛDDM1 model) and 2-body (ΛDDM2 model) decaying DM models are implemented by including two additional parameters with respect to ΛCDM; one is the decay rate, $\Gamma$, and the other is the mass-energy fraction transferred to DR, $\varepsilon = (1/2)(1 - m_{\text{WDM}}^2/m_{\text{DDM}}^2)$, where $0 \le \varepsilon \le 1/2$. In the two extreme cases, $\varepsilon = 0$ corresponds to the standard ΛCDM case with no DM decay and $\varepsilon = 1/2$ to the ΛDDM1 case. The values in between correspond to the ΛDDM2 case. Abellán et al. (2021) introduced a new approximation scheme that allows one to accurately and quickly compute the dynamics of the WDM linear perturbations by treating the WDM species as a viscous fluid and included the background evolution of WDM and, in particular, density perturbations, which were not included in previous studies (Vattis et al. 2019; Haridasu & Viel 2020; Clark et al. 2021). In our analyses, we use the DDM model including this new scheme.

In the ΛDDM1 model (DM → DR), the DM decay reduces the matter content in the Universe, causing a reduction in the expansion rate compared to the ΛCDM model, as shown in Fig. 1 of G21. This reduces the growth of matter-density fluctuations at small scales compared to the ΛCDM model, as shown in Fig. 4 of G21. The amount of matter power suppression increases as the DM lifetime becomes shorter.

In the ΛDDM2 model (DM → WDM+DR), the WDM component partially contributes to the matter density, leaving the expansion rate almost unchanged compared to the ΛCDM model, as shown in Fig. 1 of G21. However, as in the ΛDDM1 model, this component suppresses the growth of matter-density fluctuations due to the free-streaming of WDM at small scales compared to the ΛCDM case, as shown in Fig. 4 of G21. Also similarly to the ΛDDM1 model, the amount of matter power suppression increases as the DM lifetime becomes shorter, and

<footnote n="1">¹ https://github.com/PoulinV/class_majoron</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** SZ power spectra with the ΛDDM2 model (DM → WDM+DR) with different decay lifetimes of Γ⁻¹ = 30 (dotted line), 100 (dash-dot line), and 300 (dashed line) Gyr. “ΛCDM” represents the ΛCDM model with stable cold dark matter (solid line). The mass-energy fractions transferred to DR are ε = 0.1 (blue) in the left panel and 0.01 (green) in the right panel. Here, the Planck 2018 cosmological parameters are assumed in all cases.</figcaption></figure>

<table><caption>**Table 1.** Sampling parameters and priors.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Prior</th></tr></thead><tbody><tr><td>Baryon density</td><td>Ω<sub>b</sub>h<sup>2</sup></td><td>[0.0212–0.0234]</td></tr><tr><td>Primordial density perturbation amplitude</td><td>ln(10<sup>10</sup>A<sub>s</sub>)</td><td>[2.94–3.16]</td></tr><tr><td>Primordial density perturbation spectral index</td><td>n<sub>s</sub></td><td>[0.92–1.00]</td></tr><tr><td>Reionization optical depth</td><td>τ<sub>reio</sub></td><td>[0.01–0.10]</td></tr><tr><td>Hubble constant</td><td>H<sub>0</sub></td><td>[64–71]</td></tr><tr><td>Initial matter density of decaying dark matter</td><td>Ω<sup>ini</sup><sub>dcdm</sub> h<sup>2</sup></td><td>[0.11–0.13]</td></tr><tr><td>Inverse of decay lifetime</td><td>log<sub>10</sub> Γ</td><td>[-4 to 1]</td></tr><tr><td>Fractional mass energy transferred to DR</td><td>log<sub>10</sub> ε</td><td>[-4 to log<sub>10</sub>(0.5)]</td></tr><tr><td>CIB contamination</td><td>A<sub>CIB</sub></td><td>N(1, 0.5)</td></tr><tr><td>IR-source contamination</td><td>A<sub>IR</sub></td><td>N(1, 0.5)</td></tr><tr><td>Radio-source contamination</td><td>A<sub>rad</sub></td><td>N(1, 0.5)</td></tr><tr><td>Noise</td><td>A<sub>CN</sub></td><td>N(1, 0.5)</td></tr><tr><td>Mass bias</td><td>1 − b</td><td>N(0.78, 0.092)</td></tr></tbody></table>

**Notes.** Sampling cosmological parameters and nuisance parameters are listed. Nuisance parameters include A<sub>CIB</sub>, A<sub>IR</sub>, A<sub>rad</sub>, A<sub>CN</sub>, and 1 − b. [min-max] corresponds to a flat prior with minimum and maximum values. N(mean, variance) corresponds to a Gaussian prior with a given mean and variance.

the scale of the power suppression is determined by the free-streaming length of WDM, similarly to that induced by massive neutrinos.

In order to include these DDM models in the model of SZ power spectrum, we modified the comoving volume, HMF, and the pressure profile of the halo in Eq. (4). First, the comoving volume at a given redshift depends on the background evolution of assumed cosmological models, and the model including the DDM can be directly computed with G21.

Second, we consider the modification of the HMF. The theoretical prescription of the SZ power spectrum, including the DDM model, was studied in Takahashi et al. (2004). The authors included the impact of DDM on the HMF by modifying the Press–Schechter formalism (Press & Schechter 1974). We instead use the HMF from Tinker et al. (2008) but replace the matter power spectrum in Eq. (8) and ρ<sub>m,0</sub> in Eq. (6) with ones including the DDM model. The HMF in the ΛDDM1 model was tested using N-body simulations in Enqvist et al. (2020, E20). E20 checked the deviation between their HMF in the ΛDDM1

model and the Tinker HMF form in the Planck ΛCDM model and provided its fitting formula with an accuracy of ~20% for halos with 10<sup>14</sup>–10<sup>15</sup> h⁻¹ M<sub>☉</sub> at 0 < z < 1. We checked our modified Tinker HMF in the ΛDDM1 model by comparing with the results of E20. Our HMF shows a similar trend to that found by E20: the relative discrepancy between the HMFs in the ΛDDM1 and ΛCDM models increases as the halo mass increases and the rate of DDM increases (or the DDM lifetime decreases). In addition, we find that our modified HMF is consistent with that of E20 to within ~5% below 10<sup>14</sup> h⁻¹ M<sub>☉</sub> and to within ~16% in the range of 10<sup>14</sup>–10<sup>15</sup> h⁻¹ M<sub>☉</sub> around our best-fit value of Γ⁻¹ ~ 220 Gyr, which is within the accuracy of the HMF of E20. With the ΛDDM2 model, a similar result using numerical simulations has not yet been published; however, our best-fit value of ε is small, namely of ~0.002, as shown later in Table 2 and the impact of the ε parameter on the SZ power spectrum is minor, as shown in Fig. 3. Therefore, our modified HMF in the ΛDDM2 model would have a similar accuracy to the one in the case of the ΛDDM1 model.
<table><caption>**Table 2.** Best cosmological parameter estimates under the ΛDDM1 and ΛDDM2 models using the Planck CMB data without or with SZ data.</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="2">ΛDDM1</th><th colspan="2">ΛDDM2</th></tr><tr><th>Planck CMB</th><th>+BAO+SNIa+SZ</th><th>Planck CMB</th><th>+BAO+SNIa+SZ</th></tr></thead><tbody><tr><td>100 Ω<sub>b</sub>h<sup>2</sup></td><td>0.0223 ± 0.0002</td><td>0.0224 ± 0.0001</td><td>0.0223 ± 0.0002</td><td>0.0226 ± 0.0002</td></tr><tr><td>ln (10<sup>10</sup> A<sub>s</sub>)</td><td>3.045<sup>+0.017</sup><sub>-0.016</sub></td><td>3.045<sup>+0.021</sup><sub>-0.018</sub></td><td>3.045<sup>+0.016</sup><sub>-0.015</sub></td><td>3.042<sup>+0.020</sup><sub>-0.016</sub></td></tr><tr><td>n<sub>s</sub></td><td>0.9642 ± 0.0045</td><td>0.9683<sup>+0.0053</sup><sub>-0.0041</sub></td><td>0.9641<sup>+0.0044</sup><sub>-0.0043</sub></td><td>0.9726<sup>+0.0041</sup><sub>-0.0069</sub></td></tr><tr><td>τ<sub>reio</sub></td><td>0.0542<sup>+0.0081</sup><sub>-0.0078</sub></td><td>0.0566<sup>+0.0099</sup><sub>-0.0084</sub></td><td>0.0560<sup>+0.0096</sup><sub>-0.0083</sub></td><td>0.0560<sup>+0.0082</sup><sub>-0.0085</sub></td></tr><tr><td>H<sub>0</sub></td><td>67.2 ± 0.6</td><td>68.3<sup>+0.6</sup><sub>-0.5</sub></td><td>67.2 ± 0.6</td><td>68.6<sup>+0.7</sup><sub>-0.6</sub></td></tr><tr><td>Ω<sup>ini</sup><sub>dcdm</sub> h<sup>2</sup></td><td>0.1202 ± 0.0014</td><td>0.1187 ± 0.0012</td><td>0.1202 ± 0.0014</td><td>0.1172<sup>+0.0016</sup><sub>-0.0014</sub></td></tr><tr><td>log<sub>10</sub> Γ</td><td>Not preferred</td><td>-2.34 ± 0.23</td><td>Not constrained</td><td>-2.14<sup>+0.23</sup><sub>-0.36</sub></td></tr><tr><td>log<sub>10</sub> ε</td><td>–</td><td>–</td><td>Not constrained</td><td>-2.68<sup>+1.07</sup><sub>-0.87</sub></td></tr><tr><td>Ω<sub>m</sub></td><td>0.3168<sup>+0.0089</sup><sub>-0.0087</sub></td><td>0.288<sup>+0.009</sup><sub>-0.012</sub></td><td>0.3167<sup>+0.0087</sup><sub>-0.0084</sub></td><td>0.298 ± 0.009</td></tr><tr><td>σ<sub>8</sub></td><td>0.8116<sup>+0.0078</sup><sub>-0.0076</sub></td><td>0.8013<sup>+0.0091</sup><sub>-0.0082</sub></td><td>0.8108<sup>+0.0080</sup><sub>-0.0085</sub></td><td>0.7932<sup>+0.0139</sup><sub>-0.0198</sub></td></tr><tr><td>S<sub>8</sub></td><td>0.834 ± 0.017</td><td>0.7852<sup>+0.016</sup><sub>-0.023</sub></td><td>0.833 ± 0.017</td><td>0.7901<sup>+0.019</sup><sub>-0.024</sub></td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** SZ power spectra with the ΛDDM2 model (DM → WDM+DR) for the lifetime of Γ⁻¹ = 200 Gyr with different values of ε, namely, ε = 0.1 (dotted line), 0.01 (dash-dot line), and 0.001 (dashed line). “ΛCDM” represents the ΛCDM model with stable cold dark matter (solid line). Here, the Planck 2018 cosmological parameters are assumed in all cases.</figcaption></figure>

Third, we consider the modification of the pressure profile of the halo. Currently, the model of this pressure profile is constrained by the pressure measurements in galaxy clusters with X-ray or the SZ effect. For example, in Planck Collaboration V (2013) and Pointecouteau et al. (2021), the SZ profiles of galaxy clusters were stacked and fitted with the gNFW model with five parameters in Eq. (9). We include the impact from the DDM model on the pressure profile model in Planck Collaboration V (2013) by modifying the time evolution in Eq. (11) and replacing the H(z) with the one including the DDM model. The pressure profile of the halo is also given as a function of mass, and this relation might be modified when the DDM model is included. However, we assume that this is not the case. This assumption would not be valid if the DM decay were to depend its density, but the DDM model we consider here does not have such a dependency.

The SZ power spectra, including these modifications, are shown in Fig. 1 for the ΛDDM1 model with different decay lifetimes of Γ⁻¹ = 30, 100, and 300 Gyr, which are compared

with the ΛCDM model with stable CDM under the Planck 2018 cosmology of Table 1 in Planck Collaboration VI (2020). As expected, the DM lifetime Γ⁻¹ determines the depth of the suppression, and the amount of power suppression increases for smaller lifetimes.

We also show the SZ power spectra for the ΛDDM2 model in Fig. 2 with different decay lifetimes of Γ⁻¹ = 30, 100, and 300 Gyr and also with different mass-energy fractions transferred to DR with ε = 0.1 and 0.01. As in the case of the DDM1 model, the lifetime Γ⁻¹ determines the depth of the suppression, which increases for smaller lifetimes. Here, ε determines the free-streaming scale of WDM and the suppression scale of the matter power spectrum, reducing the number of halos in the HMF. Because this reduction is more dramatic for massive halos, the SZ power spectrum is more suppressed at large scales, as shown in Fig. 3.

## 2.6. Emulator of SZ power spectrum

We compute the $C_\ell^{SZ}$ in our cosmological analysis with the Monte Carlo Markov chain (MCMC) in Sect. 3.2, and therefore the $C_\ell^{SZ}$ computation time has to be minimized. To achieve this, we adopted the machine learning technique that uses the Random Forest algorithm developed<fnref n="2" /> by Douspis et al. (2022) and used in Gorce et al. (2022), and modified it to build an emulator adapted to the DDM models. We built about 50 000 SZ power spectra with a random sampling of cosmological parameters within the ranges listed in Table 1 at 50 different multipoles between 2 and 3072 and used 60%, 20%, and 20% of the sample for the training, validation, and test, respectively. We optimized the hyperparameters of the Random Forest, such as numbers of trees and depths to be 40 and 30, respectively, to maximize the final score of the prediction. This resulted in the accuracy of the constructed model being better than 1% in our analysis range of 60 < ℓ < 1411, which is well below the uncertainty in the current SZ measurement of ~30% used in our analysis.

## 3. Cosmological analysis

We use the measurement of the SZ angular power spectrum in T22 to explain the $S_8$ tension and constrain the DDM model.

<footnote n="2">² and distributed at https://sz-power-spectra.osups.universite-paris-saclay.fr/</footnote>
In our cosmological analyses, this SZ measurement is also combined with other low-redshift probes of BAO+SNIa. These measurements are compared to the theoretical prescription with the DDM model described in Sect. 2.

### 3.1. Compton y-map

A new all-sky y map was reconstructed in T22 from the Planck 2020 data release 4 (Planck Collaboration Int. LVII 2020)<fnref n="3" /> in HEALpix<fnref n="4" /> format (Górski et al. 2005) with a pixel resolution of $N_{side} = 2048$ (~1.7 arcmin). The reconstruction was performed based on the modified internal linear combination algorithm (MILCA; Hurier et al. 2013 by combining six Planck HFI frequency maps between 100 and 857 GHz so that the spectral response of the Compton y parameter is unity.

We use sky masks suitable for the analysis of the y maps provided in the 2015 Planck data release (Planck Collaboration XXII 2016). These masks cover point sources and regions around the Galactic plane, excluding about 58% of the sky. However, the residual extragalactic foreground emissions may remain in the y map. We therefore model them including radio and infrared (IR) point sources with the Planck sky model (Delabrouille et al. 2013) and cosmic infrared background (CIB) fluctuations from Maniyar et al. (2021).

Our cosmological analysis uses two y-maps reconstructed from the first and last half-ring Planck frequency maps and computes their cross-power spectrum to avoid the bias induced by the noise in the auto-power spectrum. We also limit our analysis to the multipole range of $60 < l$ in order to minimize the contribution from the foreground and to $l < 1411$ in order to minimize the contribution from the noise.

### 3.2. Maximum likelihood analysis

Cosmological constraints can be obtained by fitting the SZ power spectrum measurement with the SZ and foreground models simultaneously. In our model, we consider four components: SZ, CIB, radio point sources, and IR point sources. We also include an instrumental noise, CN, in the model. Finally, the observed SZ power spectrum, $C_{\ell}^{\text{obs}}$, is modeled by

$$C_{\ell}^{\text{obs}} = C_{\ell}^{\text{SZ}}(\Theta, l-b) + A_{\text{CIB}} C_{\ell}^{\text{CIB}} + A_{\text{IR}} C_{\ell}^{\text{IR}} + A_{\text{rad}} C_{\ell}^{\text{rad}} + A_{\text{CN}} C_{\ell}^{\text{CN}}, \quad (13)$$

where $C_{\ell}^{\text{SZ}}(\Theta, l-b)$ is the SZ power spectrum, $\Theta$ is the set of free cosmological parameters shown below in Eq. (14), $l-b$ is the hydrostatic mass bias, $C_{\ell}^{\text{CIB}}$ is the CIB power spectrum, $C_{\ell}^{\text{IR}}$ and $C_{\ell}^{\text{rad}}$ are the IR and radio source power spectra, and $C_{\ell}^{\text{CN}}$ is the empirical model for the instrumental noise. For the cosmological parameters, we vary

$$\Theta_{\Lambda\text{DDM1}} = \{\Omega_b h^2, \ln(10^{10} A_s), n_s, \tau_{\text{reio}}, H_0, \Omega_{\text{dcdm}}^{\text{ini}} h^2, \log_{10} \Gamma\} \quad (14)$$

in the $\Lambda$DDM1 model and

$$\Theta_{\Lambda\text{DDM2}} = \{\Omega_b h^2, \ln(10^{10} A_s), n_s, \tau_{\text{reio}}, H_0, \Omega_{\text{dcdm}}^{\text{ini}} h^2, \log_{10} \Gamma, \log_{10} \varepsilon\} \quad (15)$$

in the $\Lambda$DDM2 model, where $\Omega_b h^2$ is the baryon density, $\ln(10^{10} A_s)$ is the primordial density perturbation amplitude, $n_s$ is the primordial density perturbation spectral index, $\tau_{\text{reio}}$ is the Reionization optical depth, $H_0$ is the Hubble constant, $\Omega_{\text{dcdm}}^{\text{ini}} h^2$

is the initial matter density of decaying dark matter, $\log_{10} \Gamma$ is the decay rate in logarithmic scale, and $\log_{10} \varepsilon$ is the fractional mass energy transferred to DR in logarithmic scale.

We performed the MCMC likelihood analysis using the MONTEPYTHON-v3 (Audren et al. 2013; Brinckmann & Lesgourgues 2019). We first performed the cosmological analysis only with the high-redshift Planck CMB data from high-$l$ TT,TE,EE lite and low-$l$ TT, EE (Planck Collaboration VI 2020) under the $\Lambda$DDM1 and $\Lambda$DDM2 models, respectively, without the low-redshift data of the SZ, BAO, and SNIa. We then combined the CMB data with the low-redshift probes and compared the results. We used flat priors for the cosmological parameters and Gaussian priors for nuisance parameters. The sampling parameters and priors used in our cosmological analysis are summarized in Table 1.

### 3.3. Results

Figure 4 shows the result of our cosmological analysis with the $\Lambda$DDM1 model. In this model, DM decays into DR and $\Omega_m$ decreases compared to that in the $\Lambda$CDM. Therefore, the lower $S_8$ value is mainly induced by the lower value of $\Omega_m$, and not by that of $\sigma_8$, as shown in Fig. 4. This figure shows the posterior distributions of the cosmological parameters, $\log_{10} \Gamma, S_8, \sigma_8$, and $\Omega_m$ with 68% and 95% confidence interval contours when only the Planck CMB data are used (gray) and the Planck CMB data are combined with the BAO + SNIa + SZ (red). The posterior distributions of other cosmological parameters are summarized in Table 2 with 68% confidence interval, and we find that they are consistent between the two cases.

The CMB data can constrain the DDM model through the integrated Sachs-Wolfe (ISW) effect and the lensing of the CMB as shown in Figs. 5 and 13 in Aoyama et al. (2014) and put a tight constraint of $\Gamma^{-1} > 160$ Gyr on the $\Lambda$DDM1 model (e.g., Audren et al. 2014; Enqvist et al. 2015, 2020; Poulin et al. 2016). In our analysis, a combination of the high-redshift CMB data and low-redshift BAO + SNIa + SZ data provides the best-fit value on the DM decay lifetime of

$$\Gamma^{-1} = 220.5^{+155.8}_{-91.5} \text{ [Gyr]} \quad (\Lambda\text{DDM1}), \quad (16)$$

showing that the value of $\Gamma^{-1} \sim 220$ Gyr is better able to resolve the $S_8$ tension between the Planck CMB measurement and the SZ measurement. Given the uncertainty on our measurement, we also constrained a lower bound on the DM decay lifetime of ~38 Gyr at 95% confidence level but did not obtain a tighter constraint than that based on the CMB data.

Figure 5 shows the result of our cosmological analysis with the $\Lambda$DDM2 model. In this model, DM decays into WDM, but $\Omega_m$ does not decrease unlike the case in the $\Lambda$DDM1 model. (We note that the scale ranges of $\Omega_m$ and $\sigma_8$ are the same in Figs. 4 and 5.) Rather, the $\sigma_8$ is lowered by the matter power suppression at small scales because of the free-streaming length of WDM, which is similar to that caused by massive neutrinos. The lower $S_8$ value is therefore mainly a result of the lower value of $\sigma_8$, and not of the lower value of $\Omega_m$, as shown in Fig. 5. This figure shows the posterior distributions of the cosmological parameters $\log_{10} \Gamma, S_8, \log_{10} \varepsilon, \sigma_8$, and $\Omega_m$ with 68% and 95% confidence interval contours when only the Planck CMB data are used (gray) and the Planck CMB data are combined with the BAO + SNIa + SZ data (red). The posterior distributions of other cosmological parameters are summarized in Table 2 with 68% confidence interval, and we find that they are consistent between the two cases.

<footnote n="3">³ https://pla.esac.esa.int</footnote>

<footnote n="4">⁴ http://healpix.sourceforge.net/</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 4. Posterior distributions of the cosmological parameters under the $\Lambda$DDM1 model, $\log_{10}\Gamma$, $S_8$, $\sigma_8$, and $\Omega_m$ with 68% and 95% confidence interval contours obtained from our cosmological analysis, when only the Planck CMB data are used (gray) and the Planck CMB data are combined with the BAO + SNIa + SZ data (red). The latest result with the SZ power spectrum from T22 is shown in blue. We note that only $\Omega_m$ and $\sigma_8$ parameters are varied in T22, while the six $\Lambda$CDM cosmological parameters are varied in our analysis as shown in Table 1.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Posterior distributions of the cosmological parameters under the $\Lambda$DDM2 model, $\log_{10}\Gamma$, $\log_{10}\varepsilon$, $S_8$, $\sigma_8$, and $\Omega_m$ with 68% and 95% confidence interval contours obtained from our cosmological analysis, when only the Planck CMB data are used (gray) and the Planck CMB data are combined with the BAO + SNIa + SZ data (red). The latest result with the SZ power spectrum from T22 is shown in blue. We note that only $\Omega_m$ and $\sigma_8$ parameters are varied in T22, while the six $\Lambda$CDM cosmological parameters are varied in our analysis as shown in Table 1.</figcaption></figure>

Similarly to the case with the $\Lambda$DDM1 model, a combination of the high-redshift CMB data and low-redshift BAO+SNIa+SZ data provides the best-fit value on the DM decay lifetime of

$$ \Gamma^{-1} = 137.3_{-56.6}^{+176.8} \text{ [Gyr]} \quad (\Lambda\text{DDM2}), \quad (17) $$

which gives a lower bound on the decay lifetime of DM of ~24 Gyr at 95% confidence level.

In summary, we find that both models give lower $S_8$ values than the value from the Planck CMB measurement and reconcile the $S_8$ tension. Therefore, both models may provide a solution for the tension. We also compare the best-fit $\chi^2$ values in our cosmological analyses. These are 1021.94 and 1020.67 with the $\Lambda$DDM1 and $\Lambda$DDM2 models, respectively, and we find a slight preference for the $\Lambda$DDM2 model compared to the $\Lambda$DDM1 model.

However, one of the DDM parameters, $\varepsilon$, is not well constrained in our analysis. To constrain this parameter, precise measurements of the shape of the SZ power spectrum are required as shown in Fig. 3. However, the uncertainty in the current SZ measurement is ~30%, which prevents us from determining the value of this parameter. For this purpose, more sensitive data at small scales will be useful, such as those provided by AdvACT Henderson et al. (2016) and SPT-3G Benson et al. (2014) and the upcoming Simons Observatory Ade et al. (2019) and CMB-S4 Abazajian et al. (2019).

### 3.4. Systematic effects

Finally, we consider systematic uncertainties in our analysis. As discussed in Planck Collaboration XXII (2016) and T22, the mass bias and the pressure profile model are not well constrained, which affects the cosmological analysis with the SZ

power spectrum. We summarize the results in Tables 3 and 4 for the $\Lambda$DDM1 and $\Lambda$DDM2 models.

First, in order to investigate the impact on our DDM constraints from the mass bias, we replaced the mass bias prior on the CCCP with two others: one from the “weighting the giants” weak lensing measurements (von der Linden et al. 2014; WtG) and the other from cosmological hydrodynamical simulations (Biffi et al. 2016; BIFFI). Table 3 shows a comparison of the effects of these two mass bias priors on $\Gamma_{-1}$ and its lower bound at 95% confidence level. The result shows that the $\Gamma^{-1}$ value increases as the mass bias increases (or $1-b$ decreases). The reason is that the amplitude of the SZ power spectrum model shifts downwards when the mass bias increases, as in Eq. (11), which reduces the $S_8$ discrepancy between the Planck CMB and SZ results and thus requires a smaller amount of DDM or, in other words, a longer decay lifetime.

In addition, we investigated the impact of the pressure profile model on our DDM constraints by replacing P13 with two other models: one based on the combination of XMM-Newton measurements and numerical simulations (Arnaud et al. 2010; A10) and the other based on the analysis with combined data from Planck and ACT (Pointecouteau et al. 2021; PACT21). Table 3 shows a comparison of the effects of these models on $\Gamma^{-1}$ and its lower band at 95% confidence level. The result shows that the $\Gamma^{-1}$ value of PACT21 is higher than the others. The reason is that the amplitude of the PACT21 pressure profile is slightly lower than the others, as shown in Fig. 5 of Pointecouteau et al. (2021). This shifts the amplitude of the SZ power spectrum model downwards, reduces the $S_8$ discrepancy between the Planck CMB and SZ results, and requires less DDM.
<table><caption>**Table 3.** $\Gamma^{-1}$ constraints obtained with different mass bias priors and pressure profile models under $\Lambda$DDM1.</caption><thead><tr><th>Mass bias prior</th><th>1 - b</th><th>Γ<sup>-1</sup> [Gyr] (lower bound)</th></tr></thead><tbody><tr><td>CCCP (fiducial)</td><td>0.780 ± 0.092</td><td>220.5<sup>+155.8</sup><sub>-91.5</sub> (>37.5)</td></tr><tr><td>WtG</td><td>0.688 ± 0.072</td><td>255.1<sup>+119.8</sup><sub>-108.8</sub> (>37.5)</td></tr><tr><td>BIFFI</td><td>0.877 ± 0.015</td><td>156.1<sup>+82.2</sup><sub>-51.2</sub> (>53.7)</td></tr><tr><th>Pressure model</th><th></th><th>Γ<sup>-1</sup> [Gyr] (lower bound)</th></tr><tr><td>P13 (fiducial)</td><td></td><td>220.5<sup>+155.8</sup><sub>-91.5</sub> (>37.5)</td></tr><tr><td>PACT21</td><td></td><td>276.1<sup>+143.4</sup><sub>-131.5</sub> (>13.1)</td></tr><tr><td>A10</td><td></td><td>209.3<sup>+110.8</sup><sub>-73.3</sub> (>62.7)</td></tr></tbody></table>

<table><caption>**Table 4.** $\Gamma^{-1}$ constraints obtained with different mass bias priors and pressure profile models under $\Lambda$DDM2.</caption><thead><tr><th>Mass bias prior</th><th>1 - b</th><th>Γ<sup>-1</sup> [Gyr] (lower bound)</th></tr></thead><tbody><tr><td>CCCP (fiducial)</td><td>0.780 ± 0.092</td><td>137.3<sup>+176.8</sup><sub>-56.6</sub> (>24.1)</td></tr><tr><td>WtG</td><td>0.688 ± 0.072</td><td>190.1<sup>+157.8</sup><sub>-81.3</sub> (>27.5)</td></tr><tr><td>BIFFI</td><td>0.877 ± 0.015</td><td>102.2<sup>+133.6</sup><sub>-37.4</sub> (>27.4)</td></tr><tr><th>Pressure model</th><th></th><th>Γ<sup>-1</sup> [Gyr] (lower bound)</th></tr><tr><td>P13 (fiducial)</td><td></td><td>137.3<sup>+176.8</sup><sub>-56.6</sub> (>24.1)</td></tr><tr><td>PACT21</td><td></td><td>181.9<sup>+198.9</sup><sub>-87.8</sub> (>6.3)</td></tr><tr><td>A10</td><td></td><td>141.1<sup>+164.6</sup><sub>-54.9</sub> (>31.3)</td></tr></tbody></table>

## 4. Summary and conclusion

To solve the $S_8$ tension between the *Planck* CMB measurement and low-redshift probes, we extend the $\Lambda$CDM model, including a decaying dark matter (DDM) model. Two DDM models are tested in our study: one DDM model ($\Lambda$DDM1) where DM decays into a form of noninteracting dark radiation (DR), which is parameterized with the decay rate, $\Gamma$, and another model ($\Lambda$DDM2), where the DM decays into warm dark matter (WDM) and DR, which is parameterized with the decay rate, $\Gamma$, and the mass-energy fraction transferred to the massless component, $\epsilon$. For the low-redshift probe, we use the Sunyaev Zel'dovich effect and compute the impact of DDM on the SZ power spectrum by varying the DDM parameters, including the background evolution in cosmology and nonlinear prescription in the halo mass function. The result shows the suppression of the SZ power spectrum relative to the $\Lambda$CDM model due to DDM. As one expects, the suppression becomes more remarkable as the decay time becomes shorter. We combine this SZ data with other low-redshift probes, namely BAO+SNIa, which constrain the expansion history of the Universe.

We performed a cosmological analysis and compared the results when only the *Planck* CMB data are used and when the *Planck* CMB data are combined with the BAO, SNIa, and SZ data from T22. Under the $\Lambda$DDM1 model, the result shows a preference for $\Gamma^{-1} = 220.5^{+155.8}_{-91.5}$ Gyr to resolve the tension between the *Planck* CMB and SZ measurements on the $S_8$ parameter, in agreement with literature supporting $\Gamma^{-1} > 160$ Gyr. Given the uncertainty on our measurement, we also set a lower bound on the DM decay lifetime of ~38 Gyr at 95% confidence level. We also performed a cosmological analysis under the $\Lambda$DDM2 model using the *Planck* CMB and BAO+SNIa+SZ measurements. The result shows a preference for $\Gamma^{-1} = 137.3^{+176.8}_{-56.6}$ Gyr to resolve the $S_8$ tension. This result provides a lower bound on the DM decay lifetime of ~24 Gyr at 95% confidence level. As a result, both models give lower $S_8$ values than the value from the *Planck* CMB measurement and reconcile the $S_8$ tension. Thus, both models may provide a solution for the tension. In addition, we compared the best-fit $\chi^2$ values in these cosmological analyses and find a slight preference for the $\Lambda$DDM2 model compared to the $\Lambda$DDM1 model.

We also checked the systematic uncertainty in our analysis with the SZ power spectrum, which includes the impact of the mass bias and pressure profile model. The lower bound changes depending on the choice of these models. This indicates that the DDM model can only be more accurately constrained if the systematic uncertainties originating from the internal structures in galaxy clusters – such as contributions of nonthermal pressure and baryonic feedback effects – are well understood and modeled.

Finally, current studies based on weak lensing and galaxy clustering as low-redshift probes indicate that a lower $S_8$ value may be caused by a lower $\sigma_8$, and not a lower $\Omega_m$ (Di Valentino et al. 2021). To investigate this using the SZ signal, additional high-resolution and high-sensitivity data will be needed, such as those of AdvACT and SPT-3G and the upcoming Simons Observatory and CMB-S4.

**Acknowledgements.** This research has been supported by the funding for the Baryon Picture of the Cosmos (ByoPiC) project from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme grant agreement ERC-2015-AdG 695561. The authors acknowledge fruitful discussions with the members of the ByoPiC project (https://byopic.eu/team). This work is supported by World Premier International Research Center Initiative (WPI), MEXT, Japan. Kavli IPMU was established by World Premier International Research Center Initiative (WPI), MEXT, Japan. Kavli IPMU is supported by World Premier International Research Center Initiative (WPI), MEXT, Japan.

## References

* Abazajian, K., Addison, G., Adshead, P., et al. 2019, arXiv e-prints [arXiv:1907.04473]

* Abbott, T. M. C., Abdalla, F. B., Alarcon, A., et al. 2018, Phys. Rev. D, 98, 043526

* 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

* Ade, P., Aguirre, J., Ahmed, Z., et al. 2019, J. Cosmol. Astropart. Phys., 2019, 056

* Alam, S., Ata, M., Bailey, S., et al. 2017, MNRAS, 470, 2617

* Aoyama, S., Sekiguchi, T., Ichiki, K., & Sugiyama, N. 2014, J. Cosmol. Astropart. Phys., 2014, 021

* Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92

* Audren, B., Lesgourgues, J., Benabed, K., & Prunet, S. 2013, J. Cosmol. Astropart. Phys., 2013, 001

* Audren, B., Lesgourgues, J., Mangano, G., Serpico, P. D., & Tram, T. 2014, J. Cosmol. Astropart. Phys., 2014, 028

* Benson, B. A., Ade, P. A. R., Ahmed, Z., et al. 2014, in Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VII, eds. W. S. Holland, & J. Zmuidzinas, SPIE Conf. Ser., 9153, 91531P

* Beutler, F., Blake, C., Colless, M., et al. 2011, MNRAS, 416, 3017

* Biffi, V., Borgani, S., Murante, G., et al. 2016, ApJ, 827, 112

* Blas, D., Lesgourgues, J., & Tram, T. 2011, J. Cosmol. Astropart. Phys., 2011, 034

* Blomqvist, M., du Mas des Bourboux, H., Busca, N. G., et al. 2019, A&A, 629, A86

* Bolliet, B., Comis, B., Komatsu, E., & Macías-Pérez, J. F. 2018, MNRAS, 477, 4957

* Brinckmann, T., & Lesgourgues, J. 2019, Phys. Dark Univ., 24

* 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, 123528P

* Clark, S. J., Vattis, K., & Koushiappas, S. M. 2021, Phys. Rev. D, 103, 043014P
Cooray, A., & Sheth, R. 2002, Phys. Rep., 372, 1

de Sainte Agathe, V., Balland, C., du Mas des Bourboux, H., et al. 2019, A&A, 629, A85

Delabrouille, J., Betoule, M., Melin, J. B., et al. 2013, A&A, 553, A96

DES Collaboration 2022, Phys. Rev. D, 105, 023520P

Di Valentino, E., Melchiorri, A., Mena, O., & Vagnozzi, S. 2020a, Phys. Dark Univ., 30, 100666

Di Valentino, E., Melchiorri, A., Mena, O., & Vagnozzi, S. 2020b, Phys. Rev. D, 101, 063502

Di Valentino, E., Anchordoqui, L. A., Akarsu, Ö., et al. 2021, Astropart. Phys., 131, 102604

Doroshkevich, A. G., Khlopov, M. I., & Klypin, A. A. 1989, MNRAS, 239, 923

Douspis, M., Salvati, L., Gorce, A., & Aghanim, N. 2022, A&A, 659, A99

Enqvist, K., Nadathur, S., Sekiguchi, T., & Takahashi, T. 2015, J. Cosmol. Astropart. Phys., 2015, 067

Enqvist, K., Nadathur, S., Sekiguchi, T., & Takahashi, T. 2020, J. Cosmol. Astropart. Phys., 2020, 015

Fuß, L., & Garny, M. 2022, arXiv e-prints [arXiv:2210.06117]

Gianfagna, G., De Petris, M., Yepes, G., et al. 2021, MNRAS, 502, 5115

Gorce, A., Douspis, M., & Salvati, L. 2022, A&A, 662, A122

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Haridasu, B. S., & Viel, M. 2020, MNRAS, 497, 1757

Henderson, S. W., Allison, R., Austermann, J., et al. 2016, J. Low Temperat. Phys., 184, 772

Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140

Hoekstra, H., Herbonnet, R., Muzzin, A., et al. 2015, MNRAS, 449, 685

Holm, E. B., Herold, L., Hannestad, S., Nygaard, A., & Tram, T. 2023, Phys. Rev. D, 107, L021303

Hubert, J., Schneider, A., Potter, D., Stadel, J., & Giri, S. K. 2021, J. Cosmol. Astropart. Phys., 2021, 040

Hurier, G., Macías-Pérez, J. F., & Hildebrandt, S. 2013, A&A, 558, A118

Ilić, S., Sakr, Z., & Blanchard, A. 2019, A&A, 631, A96

Joudaki, S., Hildebrandt, H., Traykova, D., et al. 2020, A&A, 638, L1

Komatsu, E., & Kitayama, T. 1999, ApJ, 526, L1

Komatsu, E., Kitayama, T., Refregier, A., Spergel, D. N., & Pen, U. L. 2002, in The Ninth Marcel Grossmann Meeting, eds. V. G. Gurzadyan, R. T. Jantzen, & R. Ruffini, 2189

Lambiase, G., Mohanty, S., Narang, A., & Parashari, P. 2019, Eur. Phys. J. C, 79, 141

Lesgourgues, J., & Tram, T. 2011, J. Cosmol. Astropart. Phys., 2011, 032

Maniyar, A., Béthermin, M., & Lagache, G. 2021, A&A, 645, A40

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

McCarthy, F., & Hill, J. C. 2022, arXiv e-prints [arXiv:2210.14339]

Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Nunes, R. C., & Vagnozzi, S. 2021, MNRAS, 505, 5427

Planck Collaboration V. 2013, A&A, 550, A131

Planck Collaboration XX. 2014, A&A, 571, A20

Planck Collaboration XXII. 2016, A&A, 594, A22

Planck Collaboration XXIV. 2016, A&A, 594, A24

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

Planck Collaboration Int. LVII. 2020, A&A, 643, A42

Pointecouteau, E., Santiago-Bautista, I., Douspis, M., et al. 2021, A&A, 651, A73

Poulin, V., Serpico, P. D., & Lesgourgues, J. 2016, J. Cosmol. Astropart. Phys., 2016, 036

Press, W. H., & Schechter, P. 1974, ApJ, 187, 425

Ross, A. J., Samushia, L., Howlett, C., et al. 2015, MNRAS, 449, 835

Salvati, L., Douspis, M., & Aghanim, N. 2018, A&A, 614, A13

Scolnic, D. M., Jones, D. O., Rest, A., et al. 2018, ApJ, 859, 101

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

Sunyaev, R. A., & Zeldovich, Y. B. 1970, Ap&SS, 7, 3

Sunyaev, R. A., & Zeldovich, Y. B. 1972, Comments Astrophys. Space Phys., 4, 173

Takahashi, K., Oguri, M., & Ichiki, K. 2004, MNRAS, 352, 311

Tanimura, H., Douspis, M., Aghanim, N., & Salvati, L. 2022, MNRAS, 509, 300

Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709

Tröster, T., Asgari, M., Blake, C., et al. 2021, A&A, 649, A88

Vattis, K., Koushiappas, S. M., & Loeb, A. 2019, Phys. Rev. D, 99, 121302

von der Linden, A., Mantz, A., Allen, S. W., et al. 2014, MNRAS, 443, 1973

Xiao, L., Zhang, L., An, R., Feng, C., & Wang, B. 2020, J. Cosmology Astropart. Phys., 2020, 045