# Constraints on dark matter self-interaction from the internal density profiles of X-COP galaxy clusters

D. Eckert, S. Ettori, A. Robertson, R. Massey, E. Pointecouteau, D. Harvey,

I.G. McCarthy

► To cite this version:

D. Eckert, S. Ettori, A. Robertson, R. Massey, E. Pointecouteau, et al., Constraints on dark matter self-interaction from the internal density profiles of X-COP galaxy clusters. Astronomy and Astrophysics - A&A, 2022, 666, pp.A41. 10.1051/0004-6361/202243205. [hal-03669286](https://hal.science/hal-03669286)

HAL Id: hal-03669286

https://hal.science/hal-03669286

Submitted on 11 Apr 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Constraints on dark matter self-interaction from the internal density profiles of X-COP galaxy clusters

D. Eckert<fnref n="1" />, S. Ettori<fnref n="2" />,<fnref n="3" />, A. Robertson<fnref n="4" />,<fnref n="5" />, R. Massey<fnref n="4" />, E. Pointecouteau<fnref n="6" />, D. Harvey<fnref n="7" />, and I. G. McCarthy<fnref n="8" />

Received 27 January 2022 / Accepted 7 June 2022

## ABSTRACT

The fundamental properties of the postulated dark matter (DM) affect the internal structure of gravitationally bound structures. In the cold dark matter paradigm, DM particles interact only via gravity. Their distribution is well represented by an Einasto profile with shape parameter $\alpha \approx 0.18$ in the smallest dwarf galaxies or the most massive galaxy clusters alike. Conversely, if DM particles self-interact via additional forces, we expect the mass density profiles of DM halos to flatten in their central regions, thereby increasing the Einasto shape parameter. We measured the structural properties of 12 massive galaxy clusters from observations of their hot gaseous atmosphere, using the X-ray observatory *XMM-Newton*, and of the Sunyaev-Zel'dovich effect using the *Planck* all-sky survey. After removing morphologically disturbed systems, we measured Einasto shape parameters with mean $\langle\alpha\rangle = 0.19 \pm 0.03$ and intrinsic scatter $\sigma_{\alpha} = 0.06$, which is in close agreement with the prediction of the cold dark matter paradigm. We used cosmological hydrodynamical simulations of cluster formation with self-interacting DM (BAHAMAS-SIDM) to determine how the Einasto shape parameter depends on the self-interaction cross section. We used the fitted relation to turn our measurements of $\alpha$ into constraints on the self-interaction cross section, which imply $\sigma/m < 0.19 \text{ cm}^2 \text{ g}^{-1}$ (95% confidence level) at collision velocity $v_{\text{DM-DM}} \sim 1000 \text{ km s}^{-1}$. This is lower than the interaction cross section required for DM self-interactions to solve the core-cusp problem in dwarf spheroidal galaxies, unless the cross section is a strong function of velocity.

**Key words.** dark matter – X-rays: galaxies: clusters – galaxies: clusters: general – galaxies: clusters: intracluster medium

## 1. Introduction

In the currently favored Λ cold dark matter (ΛCDM) cosmological paradigm, gravitationally bound structures form hierarchically through the merging of smaller structures and accretion of material from the large-scale structure. The universality of the structure formation process should leave its imprint on the internal structure of DM halos, such that the density profiles of all halos from dwarf galaxies to galaxy clusters should follow a universal Navarro-Frenk-White (NFW) shape (Navarro et al. 1996, 1997) in which DM forms a central cusp with $\rho(r) \propto r^{-1}$ in the core, gradually steepening toward the outskirts asymptotically as $\rho(r) \propto r^{-3}$. Recent N-body simulations have shown that the Einasto (1965) functional form, in which the density profile is represented as a rolling power-law index $d \ln \rho/d \ln r \propto r^{-\alpha}$, provides a better representation of the shape of simulated halos (Navarro et al. 2004; Klypin et al. 2016; Ludlow et al. 2016; Brown et al. 2020; Ragagnin et al. 2021). The Einasto index $\alpha$ regulates the curvature of the profile, with the model profiles getting progressively more curved with increasing $\alpha$. In the CDM framework, the value of $\alpha$ is set primarily by the slope $n_s$ of the primordial power spectrum (Ludlow & Angulo 2017; Brown et al. 2020), corresponding to $\alpha \sim 0.18$ for $n_s = 0.96$

(Planck Collaboration XIII 2016). The value of $\alpha$ should be close to universal in all halos, with the exception of a dependence on the “peak height”, that is to say the amplitude of a local density fluctuation relative to its surroundings (Ludlow & Angulo 2017).

Given the universality of the structural properties of DM halos in ΛCDM, deviations from the expected universal shape can be used to set constraints on the nature of DM. In particular, a nonvanishing DM self-interaction probability could impact the internal structure of halos since collisional processes would transport heat across high-density regions, thereby homogenizing the mass distribution and flattening the observed profiles in the cores of halos (Yoshida et al. 2000; Rocha et al. 2013; Robertson et al. 2019, 2021). Self-interacting dark matter (SIDM) has been proposed as a possible solution to the core-cusp problem in dwarf spheroidal galaxies (see Tulin & Yu 2018, for a review). SIDM cross sections in the range $\sigma/m \sim 1–5 \text{ cm}^2 \text{ g}^{-1}$ would produce central densities in broad agreement with the values observed in dwarf galaxies (e.g., Spergel & Steinhardt 2000; Davé et al. 2001; Valli & Yu 2018; Ren et al. 2019).

On galaxy cluster scales, constraints on the SIDM cross section were obtained using dissociative cluster mergers

<footnote n="1">¹ Department of Astronomy, University of Geneva, Ch. d'Ecogia 16, 1290 Versoix, Switzerland
e-mail: Dominique.Eckert@unige.ch</footnote>

<footnote n="2">² INAF, Osservatorio di Astrofisica e Scienza dello Spazio, via Piero Gobetti 93/3, 40129 Bologna, Italy</footnote>

<footnote n="3">³ INFN, Sezione di Bologna, viale Berti Pichat 6/2, 40127 Bologna, Italy</footnote>

<footnote n="4">⁴ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="5">⁵ Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="6">⁶ IRAP, Université de Toulouse, CNRS, CNES, UPS, Toulouse, France</footnote>

<footnote n="7">⁷ Lorentz Institute, Leiden University, Niels Bohrweg 2, Leiden 2333, CA, The Netherlands</footnote>

<footnote n="8">⁸ Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK</footnote>
(Randall et al. 2008; Bradač et al. 2008; Kahlhoefer et al. 2014; Gastaldello et al. 2014; Harvey et al. 2015; Robertson et al. 2017), strong gravitational lensing (Meneghetti et al. 2001; Sagunski et al. 2021; Andrade et al. 2022), and the wobbling of the brightest central galaxies (Harvey et al. 2019). All the aforementioned studies have concluded that the behavior of DM in galaxy clusters is consistent with the collisionless scenario, with typical upper limits on the SIDM cross section at the level of $\sigma/m \le 0.3-1$ cm$^2$ g$^{-1}$.

In this paper, we provide high-precision estimates of the Einasto shape parameter $\alpha$ using the data of the XMM-Newton Cluster Outskirts Project (X-COP, Eckert et al. 2017), a very large program on XMM-Newton providing a deep X-ray mapping of a sample of 12 massive galaxy clusters selected from the Planck Sunyaev-Zel'dovich (SZ) survey. In a companion paper (Eckert et al. 2022, hereafter Paper I), we present a general framework to recover the mass profiles of galaxy clusters from X-ray and SZ data under the assumption that the gas is in hydrostatic equilibrium within the DM potential well. Here we apply our framework to recover the parameters of the Einasto functional form by fitting the X-COP data over a wide range in radius, which we subsequently use to search for the effect of DM self-interaction. Throughout the paper, we assume a Planck 2015 ΛCDM cosmology with $\Omega_m = 0.308$, $\Omega_\Lambda = 0.692$, and $H_0 = 67.8$ km s$^{-1}$ Mpc$^{-1}$ (Planck Collaboration XIII 2016).

## 2. Data analysis

### 2.1. The sample

The X-COP sample (Eckert et al. 2017) consists of a set of 12 massive galaxy clusters in the redshift range $0.04 < z < 0.1$. The sample contains the most significant SZ detections from the Planck PSZ1 catalog (Planck Collaboration XXIX 2014) with the exception of three highly perturbed systems for which the validity of a radially averaged analysis cannot be guaranteed. For each cluster, a minimum of five XMM-Newton pointings are available to cover the gas distribution uniformly, at least out to $2 \times R_{500}$.

The XMM-Newton data were reduced in a homogeneous way (Ghirardini et al. 2019; Eckert et al. 2019; Ettori et al. 2019) and the reduced data and thermodynamic profiles (gas density and spectroscopic temperature) are publicly available<fnref n="1" />. In addition, Planck SZ pressure profiles were derived for all systems following Planck Collaboration V (2013) from the SZ signal reconstructed using the MILCA algorithm (Hurier et al. 2013). For more details on the X-COP program and on the XMM-Newton and Planck data reduction, we refer readers to Ghirardini et al. (2019).

### 2.2. Hydrostatic mass reconstruction

To determine the total mass profile, we assumed that the gas is fully thermalized within the gravitational well set by the DM, such that the pressure gradient locally balances the gravitational force,

$$ \frac{dP_{\text{gas}}}{dr} = -\rho_{\text{gas}} \frac{GM(<r)}{r^2}. \qquad (1) $$

The validity of the assumption of hydrostatic equilibrium (HSE) in the case of the X-COP clusters is discussed in detail in Ettori et al. (2019), Eckert et al. (2019), and in Paper I. With

the exception of one peculiar system for which there is clear evidence of a breakdown of hydrostatic equilibrium (A2319, Ghirardini et al. 2018), the comparison of HSE masses with other techniques (Ettori et al. 2019) and the fitted gas and baryon fractions (Eckert et al. 2019) indicate that the HSE assumption is valid at the $\le 20\%$ level out to $R_{200}$.

The methodology adopted to recover the mass profiles from X-ray surface brightness and spectroscopic temperature profiles and SZ pressure profiles is presented in detail in Paper I. The code used to reconstruct the mass profiles from the X-ray and SZ data includes a wide range of effects (nonparametric gas density deprojection, PSF convolution, spectroscopic temperature weights, and profile covariance matrices) into an efficient Bayesian optimization framework. The code is publicly available in the form of the Python package hydromass<fnref n="2" />. The accuracy of the reconstruction technique was tested on mock X-ray observations of a synthetic NFW cluster in hydrostatic equilibrium, and the code was found to reproduce the input mass profile with an accuracy of better than 3%.

In this paper, we focus on the reconstruction of Einasto mass model parameters. Namely, we describe the mass density profile using the Einasto (1965) parametric form,

$$ \rho_{\text{Einasto}}(r) = \rho_s \exp \left[ -\frac{2}{\alpha} \left( \left( \frac{r}{r_s} \right)^{\alpha} - 1 \right) \right]. \qquad (2) $$

The density profile was integrated numerically over the volume to determine the cumulative model mass as a function of radius, $M_{\text{mod}}(<r)$. The gas density profile is described as a linear combination of King functions (Eckert et al. 2020) and fitted jointly to the X-ray surface brightness profile. The hydrostatic equilibrium equation Eq. (1) was then integrated to predict the 3D pressure profile,

$$ P_{3D}(r) = P_0 + \int_{r}^{r_0} \rho_{\text{gas}} \frac{GM_{\text{mod}}(<r')}{r'^2} dr', \qquad (3) $$

with $r_0$ being the outermost radius of the SZ pressure profile and $P_0 = P(r_0)$ being the integration constant corresponding to the pressure at the outer boundary. Finally, the model pressure was projected along the line of sight and convolved with the instrumental PSF to predict the expected spectroscopic temperature profile. For more details on the mass reconstruction technique and a thorough discussion of the associated systematic uncertainties, we refer readers to Paper I.

Following Ettori et al. (2019), rather than fitting for the characteristic density $\rho_s$, we optimized for the unitless normalization $c$, which is related to $\rho_s$ as

$$ \rho_s = \frac{\Delta}{3} \frac{c^3}{\log(1+c) - c/(1+c)} \rho_c(z), \qquad (4) $$

with $\Delta = 200$ being the chosen fit overdensity and $\rho_c(z) = 3H(z)^2/8\pi G$ being the critical density of the Universe. In addition, to enhance the stability of the procedure, we optimized for $\mu \equiv 1/\alpha$ (Mamon & Łokas 2005). All together, our model includes the following four free parameters: the scale radius $r_s$, the normalization $c$, the inverse $\mu$ of the Einasto shape parameter $\alpha$, and the integration constant $P_0$. Weak Gaussian priors were set to the Einasto model parameters (see Table 1), whereas for $P_0$ we set a Gaussian prior with a mean and standard deviation equal to the outermost SZ pressure value and its uncertainty. The adopted priors are described in detail in Table 1.

<footnote n="1">¹ https://dominiqueeckert.wixsite.com/xcop</footnote>

<footnote n="2">² https://github.com/domeckert/hydromass</footnote>
<table><caption>**Table 1.** Normal priors on the Einasto fit parameters. Here $P_m$ and $dP_m$ denote the outermost SZ pressure value and its error.</caption><thead><tr><th>Parameter</th><th>Mean</th><th>σ</th><th>Min</th><th>Max</th></tr></thead><tbody><tr><td>r<sub>s</sub> [kpc]</td><td>700</td><td>300</td><td>100</td><td>3000</td></tr><tr><td>c</td><td>1.8</td><td>1.5</td><td>0</td><td>10</td></tr><tr><td>μ</td><td>5</td><td>3</td><td>0.2</td><td>20</td></tr><tr><td>P<sub>0</sub></td><td>P<sub>m</sub></td><td>dP<sub>m</sub></td><td>P<sub>m</sub> - 2dP<sub>m</sub></td><td>P<sub>m</sub> + 2dP<sub>m</sub></td></tr></tbody></table>

In every case, the posterior distributions are much narrower than the prior and they clearly pull away from it; therefore, the results of this paper weakly depend on the choice of the prior. Using uniform priors over the range given in Table 1 also returns similar results and uncertainties for all parameters.

The optimization was performed using the No U-Turn Sampler (NUTS) as implemented in PyMC3 (Salvatier et al. 2016). Four independent chains were run in parallel, each with 1000 tuning steps and 1000 output samples. In Paper I we show that the Einasto model provides a good representation of the data at hand and closely follows the results obtained with a nonparametric technique that makes no assumption on the shape of the mass profile. As an example, in Fig. 1 we show the result of the fitting procedure for A1795. The left-hand panel shows the X-ray spectroscopic temperature profile and the combined X and SZ temperature profile obtained by dividing the SZ pressure by the X-ray gas density. The best-fit 3D and projected model and the corresponding $1\sigma$ error are displayed as well. The data are very well represented by the model throughout the entire radial range. We note that constraints over a wide radial range are required to break the degeneracy between the model parameters and determine the Einasto shape parameter with good accuracy. In Appendix A, we quantify the impact of a hydrostatic bias on the recovered Einasto indices and show that the resulting values of $\alpha$ are mildly affected by deviations from HSE. Moreover, we note that while most galaxy clusters should exhibit an elliptical shape, our results are mildly affected by the assumption of spherical symmetry (Buote & Humphrey 2012).

# 3. Results

## 3.1. Einasto shape parameter

We applied our mass profile reconstruction technique to all 12 X-COP systems. The resulting parameters are provided in Table 2. As can be seen in Table 2, despite the substantial covariance between the parameters, the data at hand are of a sufficient quality to determine the Einasto shape parameter $\alpha$ (or rather its inverse $\mu$) with good precision. For each system, we used the posterior distribution of $\mu$ to determine the corresponding value of $\alpha$ and its uncertainty. In Fig. 2 we show the estimated values of $\alpha$ for all the systems. In the majority of cases, the measured values of $\alpha$ are in the range 0.15–0.3, with three systems (A2319, A644, and A2255) exhibiting larger values of $\alpha$ ($\ge 0.5$).

To determine the average value of $\alpha$ and its intrinsic scatter across the population, we describe the population with a constant mean value and a normal intrinsic scatter. The total observed scatter is assumed to be the quadratic sum of the intrinsic and the statistical scatter. We then fit the data in PyMC3 using this simple model and recovered output chains for the mean value of $\alpha$ and the intrinsic scatter. Through this procedure, we found a mean value $\langle\alpha\rangle = 0.22 \pm 0.04$, with a substantial intrinsic scatter $\sigma_{\alpha} = 0.14^{+0.10}_{-0.05}$. The majority of the scatter can be

attributed to the three outliers quoted above (A2319, A644, and A2255). Incidentally, in Ettori et al. (2019), the mass profiles of these same systems were found to be better represented by models exhibiting a central core (Burkert, nonsingular isothermal sphere). A closer look at these systems shows that they are all unrelaxed systems with large deviations from spherical symmetry. To investigate this point, we used the [0.7–1.2] keV mosaic images of our systems to measure their centroid shift $w$, which quantifies the change in the position of the X-ray centroid as a function of the used aperture. Namely, we define $N$ apertures of decreasing radius $R_i$, with $R_1 = R_{500}$, and we determined the X-ray centroid $\Delta_i$ within each of them. The centroid shift is then defined as

$$w = \frac{1}{R_{500}} \sqrt{\frac{1}{N-1} \sum_{i=1}^{N} (\Delta_i - \Delta_{500})^2}, \quad (5)$$

with $\Delta_i - \Delta_{500}$ being the difference in the centroid position between aperture $R_i$ and $R_{500}$. The centroid shift is known to be an accurate indicator of a cluster’s dynamical state (Rasia et al. 2013), that is to say the more relaxed a cluster is, the smaller the value of $w$. The values of $w$ for all X-COP clusters are given in Table 2.

Most studies (e.g., O’Hara et al. 2006; Cassano et al. 2010; Weißmann et al. 2013) place the boundary between relaxed and unrelaxed systems in the range $w = 0.01–0.02$. We can see that all the quoted outliers exhibit values of $w$ greater than 0.02, which shows that their X-ray morphology is unrelaxed. A number of possible issues may render our mass reconstruction unreliable in unrelaxed systems. First, the validity of the hydrostatic equilibrium assumption in unrelaxed systems is unclear, as these systems are likely undergoing merging events at the present time. Second, the position of the X-ray centroid may not coincide with the bottom of the potential well, since the gas experiences hydrodynamical effects such as ram pressure that may displace it from the underlying DM distribution. In case our profiles are miscentered with respect to the bottom of the potential well, we expect the reconstructed mass profiles to be substantially flatter than the true profiles, which would bias the reconstructed values of $\alpha$ toward high values.

To avoid the potential impact of miscentering and hydrostatic bias, we selected a subsample of X-ray regular systems with $w < 0.015$. This subsample excludes four objects with high centroid shifts, for which the offset between the X-ray peak and the bottom of the gravitational potential may be large and thus the value of $\alpha$ may be unreliable. We again fit the mean value of $\alpha$ and its intrinsic scatter in the regular subsample, and found that the mean value slightly decreases, $\langle\alpha_{rel}\rangle = 0.19 \pm 0.03$. The scatter in the population substantially decreases compared to the full population to $\sigma_{\alpha,rel} = 0.06^{+0.04}_{-0.03}$, showing that the systems for which reliable measurements of $\alpha$ can be made exhibit very similar values for the Einasto shape parameter.

The average value of $\alpha$ in the regular subsample is significantly lower than the value of $0.29 \pm 0.04$ reported by Mantz et al. (2016) using Chandra data in a sample of 40 relaxed clusters. The difference may arise from the radial range probed in the two studies, as in most cases the data used by Mantz et al. (2016) do not extend beyond $\sim R_{2500}$. Alternatively, the systems used by Mantz et al. (2016) may have larger peak heights than the X-COP clusters since the Einasto parameter should be a strong function of the peak height (Klypin et al. 2016). Conversely, our results agree with the value of $0.19 \pm 0.07$ measured by Umetsu et al. (2014) on the stacked weak lensing shear
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Example Einasto model fit for A1795. The left-hand panel shows the spectroscopic X-ray temperatures (red points) and combined X and SZ temperatures obtained by dividing the SZ pressure by the gas density (orange stars). The projected spectroscopic-like model fitted to the X-ray temperatures is shown by the blue curve, whereas the 3D model temperature profile is shown as the green curve. The dotted vertical line indicates the fitted location of $R_{500}$. In the right-hand panel, we show the posterior distributions of the fitted Einasto model parameters and the covariance between them.</figcaption></figure>

<table><caption>**Table 2.** Result of Einasto fitting procedure for the 12 X-COP clusters.</caption><thead><tr><th>Cluster</th><th><i>M</i><sub>500,Einasto</sub><br>10<sup>14</sup> <i>M</i><sub>☉</sub></th><th><i>r</i><sub>s</sub><br>[kpc]</th><th><i>ρ</i><sub>s</sub><br>10<sup>-27</sup> g cm<sup>-3</sup></th><th>μ</th><th><i>w</i><br>10<sup>-3</sup></th><th><i>R</i><sub>min</sub><br><i>R</i><sub>200</sub></th><th><i>R</i><sub>max</sub><br><i>R</i><sub>200</sub></th></tr></thead><tbody><tr><td>A85</td><td>6.35<sup>+0.08</sup><sub>-0.08</sub></td><td>1688<sup>+160</sup><sub>-153</sub></td><td>1.08<sup>+0.21</sup><sub>-0.16</sub></td><td>10.88<sup>+0.41</sup><sub>-0.41</sub></td><td>3.85 ± 0.02</td><td>0.009</td><td>1.03</td></tr><tr><td>A644</td><td>6.06<sup>+0.30</sup><sub>-0.30</sub></td><td>334<sup>+13</sup><sub>-12</sub></td><td>35.86<sup>+3.17</sup><sub>-3.15</sub></td><td>1.82<sup>+0.21</sup><sub>-0.19</sub></td><td>20.97 ± 0.15</td><td>0.013</td><td>0.87</td></tr><tr><td>A1644</td><td>2.97<sup>+0.09</sup><sub>-0.09</sub></td><td>982<sup>+238</sup><sub>-174</sub></td><td>1.94<sup>+0.87</sup><sub>-0.65</sub></td><td>6.57<sup>+1.12</sup><sub>-1.03</sub></td><td>13.76 ± 0.08</td><td>0.010</td><td>1.30</td></tr><tr><td>A1795</td><td>4.36<sup>+0.12</sup><sub>-0.10</sub></td><td>323<sup>+16</sup><sub>-13</sub></td><td>23.82<sup>+2.06</sup><sub>-2.16</sub></td><td>4.05<sup>+0.24</sup><sub>-0.20</sub></td><td>2.36 ± 0.02</td><td>0.012</td><td>1.29</td></tr><tr><td>A2029</td><td>7.76<sup>+0.21</sup><sub>-0.22</sub></td><td>488<sup>+36</sup><sub>-39</sub></td><td>14.18<sup>+2.52</sup><sub>-1.84</sub></td><td>6.12<sup>+0.43</sup><sub>-0.44</sub></td><td>0.85 ± 0.01</td><td>0.012</td><td>1.13</td></tr><tr><td>A2142</td><td>9.41<sup>+0.20</sup><sub>-0.19</sub></td><td>1003<sup>+147</sup><sub>-103</sub></td><td>3.89<sup>+0.89</sup><sub>-0.88</sub></td><td>6.04<sup>+0.58</sup><sub>-0.49</sub></td><td>4.51 ± 0.03</td><td>0.007</td><td>1.07</td></tr><tr><td>A2255</td><td>5.58<sup>+0.29</sup><sub>-0.27</sub></td><td>699<sup>+49</sup><sub>-39</sub></td><td>6.62<sup>+1.07</sup><sub>-1.02</sub></td><td>1.33<sup>+0.41</sup><sub>-0.30</sub></td><td>31.40 ± 0.26</td><td>0.028</td><td>1.24</td></tr><tr><td>A2319</td><td>7.79<sup>+0.09</sup><sub>-0.10</sub></td><td>437<sup>+7</sup><sub>-6</sub></td><td>21.77<sup>+0.78</sup><sub>-0.82</sub></td><td>2.02<sup>+0.11</sup><sub>-0.10</sub></td><td>33.08 ± 0.17</td><td>0.020</td><td>1.12</td></tr><tr><td>A3158</td><td>4.32<sup>+0.14</sup><sub>-0.15</sub></td><td>730<sup>+176</sup><sub>-129</sub></td><td>4.34<sup>+1.99</sup><sub>-1.47</sub></td><td>5.49<sup>+1.23</sup><sub>-0.99</sub></td><td>5.87 ± 0.04</td><td>0.011</td><td>1.06</td></tr><tr><td>A3266</td><td>7.44<sup>+0.14</sup><sub>-0.13</sub></td><td>1528<sup>+81</sup><sub>-132</sub></td><td>1.53<sup>+0.28</sup><sub>-0.14</sub></td><td>6.45<sup>+0.40</sup><sub>-0.44</sub></td><td>30.84 ± 0.16</td><td>0.009</td><td>0.92</td></tr><tr><td>RXC1825</td><td>3.93<sup>+0.09</sup><sub>-0.09</sub></td><td>579<sup>+135</sup><sub>-74</sub></td><td>6.34<sup>+2.08</sup><sub>-2.18</sub></td><td>6.18<sup>+1.68</sup><sub>-1.22</sub></td><td>8.03 ± 0.06</td><td>0.013</td><td>1.19</td></tr><tr><td>ZW1215</td><td>6.65<sup>+0.30</sup><sub>-0.32</sub></td><td>689<sup>+85</sup><sub>-81</sub></td><td>6.86<sup>+1.97</sup><sub>-1.39</sub></td><td>2.96<sup>+0.50</sup><sub>-0.52</sub></td><td>3.73 ± 0.03</td><td>0.013</td><td>0.94</td></tr></tbody></table>

**Notes.** Cluster name; value of $M_{500}$ and uncertainty from the best-fit Einasto model; parameter values and uncertainties for the Einasto parameters $r_s$, $\rho_s$, and $\mu = 1/\alpha$ (see Eq. (2)); centroid shift $w$ (see Eq. (5)); and lower and upper bounds of the fitting range as a fraction of $R_{200}$.

profile of the CLASH sample, although our constraints are substantially tighter.

## 3.2. Dependence of $\alpha_{Einasto}$ on the SIDM cross section

To find how $\alpha_{Einasto}$ is expected to depend on the SIDM cross section, we used the BAHAMAS-SIDM cosmological simulations (McCarthy et al. 2017; Robertson et al. 2021). These simulate the growth of structures in a large volume of the Universe (400 h$^{-1}$ Mpc on a side), which includes approximately 230 clusters with mass $M_{200} > 3 \times 10^{14} M_{\odot}$ by redshift $z = 0$. Four separate simulations start from identical initial conditions,

then they model SIDM with velocity-independent cross section $\sigma/m = 0 \text{ cm}^2 \text{ g}^{-1}$ (CDM), $0.1 \text{ cm}^2 \text{ g}^{-1}$ (SIDM0.1), $0.3 \text{ cm}^2 \text{ g}^{-1}$ (SIDM0.3), and $1.0 \text{ cm}^2 \text{ g}^{-1}$ (SIDM1). During each time step, SIDM particles have a small chance of elastically scattering off nearby neighbors (Robertson et al. 2019). Baryonic particles are described using the BAHAMAS model (McCarthy et al. 2017), which includes a wide range of physical processes for galaxy evolution (cooling, star formation, and stellar and AGN feedback) and reproduces both the galaxy stellar mass function and X-ray scaling relations.

We fit an Einasto profile to the distribution of the total mass (DM + baryons) for all haloes with $M_{200} > 3 \times 10^{14} M_{\odot}$
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Measurements of the Einasto shape parameter $\alpha$ for the 12 X-COP clusters. The black line and gray shaded area show the fitted mean value and its error, whereas the magenta curve shows the fit to the “unperturbed” subsample with $w < 0.0015$.</figcaption></figure>

in each simulation’s $z = 0$ snapshot. The spherically averaged density profile of each simulated halo was calculated about the most bound particle, and the fit was done by minimizing the sum of squared residuals in $\log \rho$ over 40 radial bins, logarithmically spaced between $0.01 \times R_{200}$ and $R_{200}$. The adopted radial range matches well with the range over which the observational constraints were obtained (see Table 2). The Einasto model was found to provide a good representation of the mass profiles, although the isothermal Jeans model proposed by Kaplinghat et al. (2014) and tested on simulated halos by Robertson et al. (2021) describes the simulated halos more accurately. Figure 3 shows the best-fit value of $\alpha_{Einasto}$ for each simulated cluster, as well as the median and 16th to 84th percentiles for different SIDM cross sections. The median value of $\alpha_{Einasto}$ increases from 0.18 for CDM (as previously known; Navarro et al. 2004; Ludlow et al. 2016) to 0.35 for SIDM1, demonstrating the shallower central slope of SIDM clusters. Increasing the mass threshold to $M_{200} > 5 \times 10^{14} M_{\odot}$ does not significantly change the average value of $\alpha$.

Since we are fitting the total density profile (including baryons), it is important to check whether the dependence of $\alpha$ on the SIDM cross section could be degenerate with uncertainties associated with the evolution of baryons. The main source of uncertainty in the baryonic model is the implementation of AGN feedback, which affects the star formation efficiency and the gas density profile. BAHAMAS implements the Booth & Schaye (2009) model, in which the central black holes store accreted energy until it exceeds a given threshold. The stored energy is then released thermally to reheat the surrounding medium. A low energy threshold provides gentle, continuous energy injection, whereas a high value for the energy threshold leads to burstier AGN feedback. The default BAHAMAS model implements the AGN 7.8 model, which was found to reproduce the X-ray scaling relations and gas fractions in galaxy clusters and groups accurately (Le Brun et al. 2014; McCarthy et al. 2017). For comparison, Fig. 3 also shows the median value and scatter of $\alpha_{Einasto}$ in CDM simulations with a decreased (AGN 7.6) or increased (AGN 8.0) energy threshold (McCarthy et al. 2018), which bracket the range of parameters reproducing the observed scaling relations and gas fractions. The insignificant dependence

on the energy threshold demonstrates that $\alpha_{Einasto}$ provides a clean test of DM self-interactions, regardless of the adopted prescription for baryonic physics.

### 3.3. Constraints on the SIDM cross section

To test for deviations from the $\Lambda$CDM model, we compared our observed data against predictions of SIDM simulations. We used the values of $\alpha_{Einasto}$ fitted to individual halos in each simulation set and applied the same procedure as for the observed halos, that is to say we fit the data set with a constant mean value and a normal intrinsic scatter. In addition, to allow for a meaningful comparison with the regular X-COP subsample, we selected a subsample of relaxed systems according to the DM offset parameter $x_{off}$, which determines the offset between the most bound particle and the center of mass of the halo, relative to the system’s virial radius (Macciò et al. 2007). It is important to note that $x_{off}$ can be viewed as a DM analog of the centroid shift, such that a selection of halos with low $x_{off}$ should correspond to our selection of X-ray regular systems. We selected a subsample of 114 relaxed systems with $x_{off} < 0.05$ to be compared with the regular X-COP subsample. We found that $\alpha_{Einasto}$ slightly depends on $x_{off}$, with the mean value of the CDM sample increasing from 0.16 in the full population to 0.18 in the low $x_{off}$ subsample. We note that the effect goes in the opposite direction with respect to the observations, where restricting to the X-ray regular subsample slightly reduces the mean value of $\alpha_{Einasto}$. We interpret this difference as resulting from miscentering in dynamically active systems, for which finding the bottom of the potential well is difficult, which flattens the mass profiles in the innermost regions and thus raises the value of $\alpha$. Conversely, numerical simulations are always able to pinpoint the most bound particle in 3D. To make a meaningful comparison to the full sample, applying the exact same analysis procedure to mock observations extracted from numerical simulations would be needed, which is beyond the scope of this paper.

Given the observational uncertainty for the unrelaxed clusters, we focus here on the comparison of the X-ray regular subsample with the simulated low $x_{off}$ sample. In the left-hand panel of Fig. 4, we show the mean value of $\alpha$ and the intrinsic scatter in the low $x_{off}$ subsample for the four simulations with varying $\sigma/m$. We compare the results with the mean values obtained for the full X-COP sample and the X-ray regular subsample. We can see that the fitted mean value for the regular subsample agrees very well with the CDM expectation, both in terms of the mean value and the intrinsic scatter. We can thus set an upper limit on $\sigma/m$ by comparing the mean value of $\alpha$ to the value expected in the various simulation sets. To interpolate between the discrete values of $\sigma/m$, we describe the relation between $\alpha$ and $\sigma/m$ as

$$ \alpha_{Einasto} = \alpha_0 + \alpha_1 \left( \frac{\sigma/m}{1 \text{ cm}^2 \text{ g}^{-1}} \right)^{\gamma}, \qquad (6) $$

that is $\alpha_0$ represents the mean value of $\alpha$ in CDM and the dependence of $\alpha$ on $\sigma/m$ is described as a power law. For the low $x_{off}$ subsample, we found $\alpha_0 = 0.178$, $\alpha_1 = 0.20$, and $\gamma = 0.63$. This model accurately reproduces the data, as can be seen in the right-hand panel of Fig. 4, where the model curve is compared with the mean values of the simulation sets.

To translate our constraints on $\alpha$ into upper limits on $\sigma/m$, we converted each value of $\langle \alpha_{rel} \rangle$ in the output chain into a value of $\sigma/m$ through Eq. (6) and constructed a posterior distribution for $\sigma/m$. As expected, the posterior distribution is consistent with zero. We then computed a one-tailed 95% upper limit as
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Simulations of the structural properties of massive DM halos in the CDM and SIDM frameworks. **Left-hand panel:** Einasto shape parameter $\alpha_{Einasto}$ for individual simulated halos with a varying SIDM cross section (0.1 cm$^2$ g$^{-1}$, cyan; 0.3, purple; and 1.0, red; Robertson et al. 2021). For comparison, we also show the effect of varying baryonic physics (reference model, black; models with lower (yellow) and higher (green) feedback efficiency; McCarthy et al. 2017). **Right-hand panel:** median and $1\sigma$ percentiles of $\alpha_{Einasto}$ values for the various simulations runs.</figcaption></figure>

the 95th percentile of the posterior distribution,

$$ \frac{\sigma}{m} < 0.19 \text{ cm}^2 \text{ g}^{-1}, \text{ at } v_{\text{DM}} \approx 1000 \text{ (95\% confidence level).} \quad (7) $$

Similar constraints on $\sigma/m$ are obtained if we use a different method to relate $\alpha$ to $\sigma/m$ (linear interpolation, cubic spline, or second-order polynomial), which shows that our result is robust to changes in the interpolation method. Moreover, the statistical uncertainties in the mean value of $\alpha$ within each simulation set are much smaller than the observational errors (see the left-hand panel of Fig. 4) and thus they can safely be ignored.

The upper limit quoted here is one of the most stringent to date. For comparison, offsets between DM and baryons in individual merging clusters imply $\sigma/m \lesssim 1 \text{ cm}^2 \text{ g}^{-1}$ (e.g., Randall et al. 2008) or $\sigma/m < 0.47 \text{ cm}^2 \text{ g}^{-1}$ at a 95% confidence level from an ensemble of 72 merging systems (Harvey et al. 2015). The wobbling of the brightest cluster galaxies around the center of DM, which also probes the radial mass distribution, constrains $\sigma/m \lesssim 0.39 \text{ cm}^2 \text{ g}^{-1}$ at a 95% confidence level (Harvey et al. 2019). Using a similar technique to that presented in this paper, but from the gravitational lensing profiles of eight clusters, Andrade et al. (2022) found $\sigma/m \lesssim 0.13 \text{ cm}^2 \text{ g}^{-1}$ at a 95% confidence level. Our results are of similar precision and, combined with the studies mentioned above, imply that SIDM cross sections larger than $0.1 \text{ cm}^2 \text{ g}^{-1}$ would be in tension with existing data.

Removing the DM cusps in dwarf spheroidal galaxies, and solving the “core-cusp problem” (see Bullock & Boylan-Kolchin 2017, for a review), requires cross sections larger than $1 \text{ cm}^2 \text{ g}^{-1}$

at $v_{DM} \sim 10 \text{ km s}^{-1}$ (Spergel & Steinhardt 2000; Davé et al. 2001; Quynh Lan et al. 2021). This appears impossible for SIDM models in which the cross section is constant. However, many particle physics models predict a cross section that is a strong decreasing function of velocity, $\propto v_{DM}^{-4}$ above the (unknown) mass scale of the dark sector force mediator particle. In these models, SIDM may still be a viable solution to the core-cusp problem.

Generally speaking, the technique developed here opens the possibility of probing SIDM at much higher precision using future surveys such as eROSITA (Predehl et al. 2021) and Euclid (Euclid Collaboration 2019), which will detect tens of thousands of clusters and will potentially improve the constraints on $\alpha_{Einasto}$. The technique proposed here can also be applied to systems of lower mass to directly test for variation in the SIDM cross section with velocity.

## 4. Conclusion

We have measured the internal DM structure of galaxy clusters and compared our observations with clusters in CDM and SIDM simulations. Our results can be summarized as follows:

- The hydrostatic mass profiles of the 12 X-COP galaxy clusters are well represented by an Einasto model over more than two decades in radius. The data quality is sufficient to break the degeneracy between the various parameters of the Einasto model, and determine the Einasto index with high precision.
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Relation between the Einasto shape parameter $\alpha_{Einasto}$ and the DM self-interaction cross section $\sigma/m$ for several sets of numerical simulations with varying $\sigma/m$ (Robertson et al. 2021). The left-hand panel shows the fitted mean values of $\alpha$ in the relaxed subsample for the CDM and SIDM cases, compared to the average value in the full X-COP sample (ALL) and in the regular X-ray subsample. In the right-hand panel, we show the dependence of $\alpha$ on $\sigma/m$. The solid blue line shows a fit to the relation with a power law, whereas the dashed vertical line indicates the average value in CDM. The black line and gray shaded area show the mean value of $\alpha_{Einasto}$ and the $1\sigma$ uncertainty in the regular X-COP subsample, with the corresponding constraints on $\sigma/m$.</figcaption></figure>

- Most X-COP clusters have an Einasto index of $0.15 < \alpha < 0.3$, with mean $\langle\alpha\rangle = 0.22 \pm 0.04$. The substantial intrinsic scatter, $\sigma_{\alpha} = 0.14^{+0.10}_{-0.05}$, is dominated by a few outliers with $\alpha \gtrsim 0.5$, all of which appear morphologically disturbed, as indicated by a large centroid shift. The high Einasto indices are therefore likely a systematic effect in the reconstruction of the mass profiles, caused by hydrostatic bias and/or miscentering. The subset of eight X-ray regular systems with low centroid shift have mean $\langle\alpha\rangle = 0.19 \pm 0.03$. Their intrinsic scatter, $\sigma_{\alpha} = 0.06^{+0.05}_{-0.03}$, is so small that $\alpha$ appears close to universal.

- In CDM and SIDM simulations, we find that cluster mass profiles become more curved as the SIDM cross section $\sigma/m$ increases, with $\alpha \approx 0.18 + 0.20((\sigma/m)/(1\,\text{cm}^2\,\text{g}^{-1}))^{0.63}$. This predicted trend is robust to changes in the adopted prescription for subgrid baryonic processes.

- Comparing the relaxed subset of observed clusters with those in BAHAMAS-SIDM simulations allowed us to set an upper limit $\sigma/m < 0.19\,\text{cm}^2\,\text{g}^{-1}$ (95% c.l.) on the cross section for DM interactions at a collision velocity $v_{DM} \approx 1000\,\text{km s}^{-1}$. This is one of the most stringent constraints to date, and substantially lower than the value required to explain the core-cusp problem in dwarf spheroidal galaxies, unless the cross section is a strong decreasing function of collision velocity.

Future studies could apply our technique to large samples of galaxies, groups, and clusters. It has the potential to improve the constraints on the fundamental self-interaction cross section of DM substantially, or even to measure it as a function of velocity.

## References

* Andrade, K. E., Fuson, J., Gad-Nasr, S., et al. 2022, MNRAS, 510, 54

* Angelinelli, M., Vazza, F., Giocoli, C., et al. 2020, MNRAS, 495, 864

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

* Booth, C. M., & Schaye, J. 2009, MNRAS, 398, 53

* Bradač, M., Allen, S. W., Treu, T., et al. 2008, ApJ, 687, 959

* Brown, S. T., McCarthy, I. G., Diemer, B., et al. 2020, MNRAS, 495, 4994

* Bullock, J. S., & Boylan-Kolchin, M. 2017, ARA&A, 55, 343

* Buote, D. A., & Humphrey, P. J. 2012, MNRAS, 421, 1399

* Cassano, R., Ettori, S., Giacintucci, S., et al. 2010, ApJ, 721, L82

* Davé, R., Spergel, D. N., Steinhardt, P. J., & Wandelt, B. D. 2001, ApJ, 547, 574

* Eckert, D., Ettori, S., Pointecouteau, E., et al. 2017, Astron. Nachr., 338, 293

* Eckert, D., Ghirardini, V., Ettori, S., et al. 2019, A&A, 621, A40

* Eckert, D., Finoguenov, A., Ghirardini, V., et al. 2020, Open J. Astrophys., 3, 12

* Eckert, D., Ettori, S., Pointecouteau, E., van der Burg, R. F. J., & Loubser, S. I. 2022, A&A, 662, A123

* Einasto, J. 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87

* Ettori, S., Ghirardini, V., Eckert, D., et al. 2019, A&A, 621, A39

* Euclid Collaboration (Adam, R., et al.) 2019, A&A, 627, A23

* Gastaldello, F., Limousin, M., Foex, G., et al. 2014, MNRAS, 442, L76

* Ghirardini, V., Ettori, S., Eckert, D., et al. 2018, A&A, 614, A7

* Ghirardini, V., Eckert, D., Ettori, S., et al. 2019, A&A, 621, A41

* Harvey, D., Massey, R., Kitching, T., Taylor, A., & Tittley, E. 2015, Science, 347, 1462

* Harvey, D., Robertson, A., Massey, R., & McCarthy, I. G. 2019, MNRAS, 488, 1572

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

* Kahlhoefer, F., Schmidt-Hoberg, K., Frandsen, M. T., & Sarkar, S. 2014, MNRAS, 437, 2865

* Kaplinghat, M., Keeley, R. E., Linden, T., & Yu, H.-B. 2014, Phys. Rev. Lett., 113, 021302

* Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340

* Lau, E. T., Kravtsov, A. V., & Nagai, D. 2009, ApJ, 705, 1129

* Le Brun, A. M. C., McCarthy, I. G., Schaye, J., & Ponman, T. J. 2014, MNRAS, 441, 1270

* Loubser, S. I., Babul, A., Hoekstra, H., et al. 2020, MNRAS, 496, 1857

* Ludlow, A. D., & Angulo, R. E. 2017, MNRAS, 465, L84

* Ludlow, A. D., Bose, S., Angulo, R. E., et al. 2016, MNRAS, 460, 1214

* Macciò, A. V., Dutton, A. A., van den Bosch, F. C., et al. 2007, MNRAS, 378, 55

* Mamon, G. A., & Łokas, E. L. 2005, MNRAS, 362, 95

* Mantz, A. B., Allen, S. W., & Morris, R. G. 2016, MNRAS, 462, 681

* McCarthy, I. G., Schaye, J., Bird, S., & Le Brun, A. M. C. 2017, MNRAS, 465, 2936

* McCarthy, I. G., Bird, S., Schaye, J., et al. 2018, MNRAS, 476, 2999

* Meneghetti, M., Yoshida, N., Bartelmann, M., et al. 2001, MNRAS, 325, 435

* Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

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

* Navarro, J. F., Hayashi, E., Power, C., et al. 2004, MNRAS, 349, 1039

* Nelson, K., Lau, E. T., & Nagai, D. 2014, ApJ, 792, 25

* O'Hara, T. B., Mohr, J. J., Bialek, J. J., & Evrard, A. E. 2006, ApJ, 639, 64

* Planck Collaboration V. 2013, A&A, 550, A131
Planck Collaboration XIII. 2016, A&A, 594, A13

Planck Collaboration XXIX. 2014, A&A, 571, A29

Predehl, P., Andritschke, R., Arefiev, V., et al. 2021, A&A, 647, A1

Quynh Lan, N., Mathews, G. J., Arielle Phillips, L., et al. 2021, Mod. Phys. Lett. A, 36, 2130001

Ragagnin, A., Saro, A., Singh, P., & Dolag, K. 2021, MNRAS, 500, 5056

Randall, S. W., Markevitch, M., Clowe, D., Gonzalez, A. H., & Bradač, M. 2008, ApJ, 679, 1173

Rasia, E., Ettori, S., Moscardini, L., et al. 2006, MNRAS, 369, 2013

Rasia, E., Meneghetti, M., & Ettori, S. 2013, Astron. Rev., 8, 40

Ren, T., Kwa, A., Kaplinghat, M., & Yu, H.-B. 2019, Phys. Rev. X, 9, 031020

Robertson, A., Massey, R., & Eke, V. 2017, MNRAS, 467, 4719

Robertson, A., Harvey, D., Massey, R., et al. 2019, MNRAS, 488, 3646

Robertson, A., Massey, R., Eke, V., Schaye, J., & Theuns, T. 2021, MNRAS, 501, 4610

Rocha, M., Peter, A. H. G., Bullock, J. S., et al. 2013, MNRAS, 430, 81

Sagunski, L., Gad-Nasr, S., Colquhoun, B., Robertson, A., & Tulin, S. 2021, JCAP, 2021, 024

Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJ Comput. Sci., 2, e55

Spergel, D. N., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760

Tulin, S., & Yu, H.-B. 2018, Phys. Rep., 730, 1

Umetsu, K., Medezinski, E., Nonino, M., et al. 2014, ApJ, 795, 163

Valli, M., & Yu, H.-B. 2018, Nat. Astron., 2, 907

Weißmann, A., Böhringer, H., Šuhada, R., & Ameglio, S. 2013, A&A, 549, A19

Yoshida, N., Springel, V., White, S. D. M., & Tormen, G. 2000, ApJ, 544, L87
## Appendix A: Impact of hydrostatic bias on the shape parameter

The constraints on the Einasto shape parameter presented in this paper were derived under the assumption that the gas is in hydrostatic equilibrium within the gravitational potential set by the DM (see Sect. 2.2). While the comparison with weak lensing masses (Ettori et al. 2019) and the measured gas fractions (Eckert et al. 2019) indicate that the assumption of equilibrium is well justified in X-COP clusters, and that the selection of the regular subsample based on the centroid shift excludes the most problematic cases, we do expect a certain level of nonthermal pressure to be present within the systems of interest. The relative importance of nonthermal pressure is expected to be nonuniform in the cluster’s volume (e.g., Rasia et al. 2006; Lau et al. 2009; Nelson et al. 2014; Biffi et al. 2016; Angelinelli et al. 2020), which affects the shape of the recovered mass profiles, and hence the Einasto shape parameter.

To quantify the potential impact of hydrostatic bias on our results, we considered a fiducial cluster with a mass profile described by the Einasto functional form (Eq. 2) and a typical gas density profile (Ghirardini et al. 2019). We computed its total

pressure profile $P_{tot}$ by numerically integrating Eq. 1, and we assumed that the measured thermal pressure $P_{th}$ is a fraction of the total pressure,

$$P_{th}(r) = P_{tot}(r) - P_{NT}(r), \qquad (A.1)$$

with $P_{NT}(r)$ being the nonthermal pressure, which is expected to be dominated by kinetic motions. We considered two possible cases: first, the nonthermal pressure profile predicted by Nelson et al. (2014), which takes all the gas motions into account; and, second, the functional form introduced by Angelinelli et al. (2020), which only selects the random, turbulent motions as a source of outward pressure. We then fit the modified $P_{th}$ profile with an Einasto model in the radial range $[0.01-1]R_{200}$ and compared the retrieved value of $\alpha$ with the input one. In the left-hand panel of Fig. A.1, we show the ratio of the measured value of $\alpha$ in the HSE assumption ($\alpha_{HSE}$) to the true input value ($\alpha_{True}$) for a range of input values of $\alpha$. We can see that the HSE assumption typically overestimates $\alpha$ by 2–8% depending on the assumed nonthermal pressure profile and the true value of $\alpha$, corresponding to differences of 0.01 for $\alpha_{True} = 0.18$. Therefore, we conclude that the hydrostatic assumption has little impact on the results presented here.

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

<figure><img src="image_8.png" /><figcaption>Fig. A.1. Impact of systematic uncertainties. Left: Ratio between the value of $\alpha$ derived under the assumption of hydrostatic equilibrium ($\alpha_{HSE}$) and the true value of $\alpha$ when assuming two possible nonthermal pressure profiles (Angelinelli et al. (2020), blue; Nelson et al. (2014), orange) as a function of the input value of $\alpha$. Right: Stellar mass of the BCG within 100 kpc in the BAHAMAS CDM (red) and SIDM1 simulations (blue) as a function of the halo mass. The black stars show the measured BCG stellar masses within the same aperture in X-COP clusters (see Paper I).</figcaption></figure>
# Appendix B: Impact of the brightest cluster galaxy

Since measurements of the stellar mass profiles are currently only available for a subset of X-COP systems, in this paper we fit an Einasto profile to the total gravitational field and compared the results to simulations including a detailed modeling of the baryonic components (gas + stars, see McCarthy et al. 2017). In case the modeling of baryonic processes implemented in BAHAMAS simulations does not reproduce the observed baryonic mass profiles accurately, the corresponding comparison between simulations and observations may be biased. In particular, the stellar content of the brightest cluster galaxy (BCG) dominates the gravitational field in the innermost regions, thus differences in the BCG mass profiles between simulations and observations could bias the comparison of the Einasto shape parameter. To verify this point, we determined the stellar masses of the BCGs in the simulation runs and compared them with the stellar masses of the BCGs in X-COP clusters (see Paper I Loubser et al. 2020). In the right-hand panel of Fig. A.1, we show the BCG stellar masses within a 100 kpc radius for the CDM and SIDM1 cases as well as the X-COP BCG stellar masses within the same aperture.

We can see that the stellar masses of BCGs in the BAHAMAS model agree with the observations well and do not differ from the CDM to the SIDM case significantly, which shows that the comparisons presented in this paper are mildly impacted by the stellar content of the BCG.

While the analysis presented in the left-hand panel of Fig. A.1 does not show any potential issue associated with the BCG mass profile, it is worthwhile to quantify the potential impact of wrong BCG mass profiles on the Einasto shape parameter. Similar to the tests performed in Appendix A, we considered a fiducial Einasto profile, but this time fitted to the sum of BCG and DM profiles, and we modified the stellar mass of the BCG by a conservative factor of 2 in the upward and downward directions. We then refit the corresponding total mass profiles in the $[0.01-1]R_{200}$ range with a single Einasto profile and quantified the changes in $\alpha$ with respect to the input value. We found that $\alpha$ is anticorrelated with the BCG stellar mass, with larger values for the stellar mass corresponding to lower values of $\alpha$. For an input value $\alpha = 0.18$, the retrieved value of $\alpha$ changes by $\pm0.02$, which is smaller than the statistical uncertainties in our measurements.