EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 666 (October 2022) A&A, 666 (2022) A41 Full HTML
Open Access
Issue
A&A
Volume 666, October 2022
Article Number A41
Number of page(s) 10
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202243205
Published online 03 October 2022

© D. Eckert et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1. Introduction

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 ρ(r)∝r−1 in the core, gradually steepening toward the outskirts asympotically as ρ(r)∝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 dlnρ/dlnr ∝ r−α, 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 α regulates the curvature of the profile, with the model profiles getting progressively more curved with increasing α. In the CDM framework, the value of α is set primarily by the slope ns of the primordial power spectrum (Ludlow & Angulo 2017; Brown et al. 2020), corresponding to α ∼ 0.18 for ns = 0.96 (Planck Collaboration XIII 2016). The value of α 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 σ/m ∼ 1 − 5 cm2 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 (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 σ/m ≲ 0.3 − 1 cm2 g−1.

In this paper, we provide high-precision estimates of the Einasto shape parameter α 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 Ωm = 0.308, ΩΛ = 0.692, and H0 = 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 × R500.

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 available1. 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,

d P gas d r = − ρ gas G M ( < r ) r 2 . $$ \begin{aligned} \frac{\mathrm{d}P_{\rm gas}}{\mathrm{d}r} = - \rho _{\rm gas}\frac{GM(<r)}{r^2}. \end{aligned} $$(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 ≲20% level out to R200.

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 hydromass2. 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,

ρ Einasto ( r ) = ρ s exp [ − 2 α ( ( r r s ) α − 1 ) ] . $$ \begin{aligned} \rho _{\rm Einasto}(r) = \rho _s \exp \left[-\frac{2}{\alpha } \left( \left( \frac{r}{ r_s}\right)^\alpha - 1\right) \right]. \end{aligned} $$(2)

The density profile was integrated numerically over the volume to determine the cumulative model mass as a function of radius, Mmod(< 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 3 D ( r ) = P 0 + ∫ r r 0 ρ gas G M mod ( < r ′ ) r ′ 2 d r ′ , $$ \begin{aligned} P_{3\mathrm D}(r) = P_0 + \int _{r}^{r_0} \rho _{\rm gas} \frac{G M_{\rm mod}(<r^\prime )}{r^{\prime 2}}\,dr^\prime , \end{aligned} $$(3)

with r0 being the outermost radius of the SZ pressure profile and P0 ≡ P(r0) 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 ρs, we optimized for the unitless normalization c, which is related to ρs as

ρ s = Δ 3 c 3 log ( 1 + c ) − c / ( 1 + c ) ρ c ( z ) , $$ \begin{aligned} \rho _s = \frac{\Delta }{3} \frac{c^3}{\log (1 + c) - c / (1 + c)} \rho _c(z), \end{aligned} $$(4)

with Δ = 200 being the chosen fit overdensity and ρc(z)=3H(z)2/8πG being the critical density of the Universe. In addition, to enhance the stability of the procedure, we optimized for μ ≡ 1/α (Mamon & Łokas 2005). All together, our model includes the following four free parameters: the scale radius rs, the normalization c, the inverse μ of the Einasto shape parameter α, and the integration constant P0. Weak Gaussian priors were set to the Einasto model parameters (see Table 1), whereas for P0 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. 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.

Table 1.

Normal priors on the Einasto fit parameters. Here Pm and dPm denote the outermost SZ pressure value and its error.

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

thumbnail 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 R500. In the right-hand panel, we show the posterior distributions of the fitted Einasto model parameters and the covariance between them.

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 α (or rather its inverse μ) with good precision. For each system, we used the posterior distribution of μ to determine the corresponding value of α and its uncertainty. In Fig. 2 we show the estimated values of α for all the systems. In the majority of cases, the measured values of α are in the range 0.15–0.3, with three systems (A2319, A644, and A2255) exhibiting larger values of α (≳0.5).

Table 2.

Result of Einasto fitting procedure for the 12 X-COP clusters.

thumbnail Fig. 2.

Measurements of the Einasto shape parameter α 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.

To determine the average value of α 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 α and the intrinsic scatter. Through this procedure, we found a mean value ⟨α⟩=0.22 ± 0.04, with a substantial intrinsic scatter σ α = 0 . 14 − 0.05 + 0.10 $ \sigma_{\alpha}=0.14_{-0.05}^{+0.10} $. 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 Ri, with R1 = R500, and we determined the X-ray centroid Δi within each of them. The centroid shift is then defined as

w = 1 R 500 1 N − 1 ∑ i = 1 N ( Δ i − Δ 500 ) 2 , $$ \begin{aligned} { w} = \frac{1}{R_{500}} \sqrt{\frac{1}{N-1}\sum _{i=1}^N (\Delta _i - \Delta _{500})^2} , \end{aligned} $$(5)

with Δi − Δ500 being the difference in the centroid position between aperture Ri and R500. 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 α toward high values.

To avoid the potential impact ofmiscentering 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 α may be unreliable. We again fit the mean value of α and its intrinsic scatter in the regular subsample, and found that the mean value slightly decreases, ⟨αrel⟩=0.19 ± 0.03. The scatter in the population substantially decreases compared to the full population to σ α , rel = 0 . 06 − 0.03 + 0.04 $ \sigma_{\alpha, \mathrm{rel}}=0.06_{-0.03}^{+0.04} $, showing that the systems for which reliable measurements of α can be made exhibit very similar values for the Einasto shape parameter.

The average value of α in the regular subsample is significantly lower than the value of 0.29  ±  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 ∼R2500. 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  ±  0.07 measured by Umetsu et al. (2014) on the stacked weak lensing shear profile of the CLASH sample, although our constraints are substantially tighter.

3.2. Dependence of αEinasto on the SIDM cross section

To find how α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 M200 >  3 × 1014 M⊙ by redshift z = 0. Four separate simulations start from identical initial conditions, then they model SIDM with velocity-independent cross section σ/m = 0 cm2 g−1 (CDM), 0.1 cm2 g−1 (SIDM0.1), 0.3 cm2 g−1 (SIDM0.3), and 1.0 cm2 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 M200 >  3 × 1014 M⊙ 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ρ over 40 radial bins, logarithmically spaced between 0.01 × R200 and R200. 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 αEinasto for each simulated cluster, as well as the median and 16th to 84th percentiles for different SIDM cross sections. The median value of α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 M200 >  5 × 1014 M⊙ does not significantly change the average value of α.

thumbnail Fig. 3.

Simulations of the structural properties of massive DM halos in the CDM and SIDM frameworks. Left-hand panel: Einasto shape parameter αEinasto for individual simulated halos with a varying SIDM cross section (0.1 cm2 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σ percentiles of αEinasto values for the various simulations runs.

Since we are fitting the total density profile (including baryons), it is important to check whether the dependence of α 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 α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 α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 ΛCDM model, we compared our observed data against predictions of SIDM simulations. We used the values of α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 xoff, 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 xoff can be viewed as a DM analog of the centroid shift, such that a selection of halos with low xoff should correspond to our selection of X-ray regular systems. We selected a subsample of 114 relaxed systems with xoff <  0.05 to be compared with the regular X-COP subsample. We found that αEinasto slightly depends on xoff, with the mean value of the CDM sample increasing from 0.16 in the full population to 0.18 in the low xoff 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 α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 α. 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 xoff sample. In the left-hand panel of Fig. 4, we show the mean value of α and the intrinsic scatter in the low xoff subsample for the four simulations with varying σ/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 σ/m by comparing the mean value of α to the value expected in the various simulation sets. To interpolate between the discrete values of σ/m, we describe the relation between α and σ/m as

thumbnail Fig. 4.

Relation between the Einasto shape parameter αEinasto and the DM self-interaction cross section σ/m for several sets of numerical simulations with varying σ/m (Robertson et al. 2021). The left-hand panel shows the fitted mean values of α 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 α on σ/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 αEinasto and the 1σ uncertainty in the regular X-COP subsample, with the corresponding constraints on σ/m.

α Einasto = α 0 + α 1 ( σ / m 1 cm 2 g − 1 ) γ , $$ \begin{aligned} \alpha _{\rm Einasto} = \alpha _0 + \alpha _1 \left( \frac{\sigma /m}{1\, \mathrm{cm}^2\,\mathrm{g}^{-1}}\right)^\gamma , \end{aligned} $$(6)

that is α0 represents the mean value of α in CDM and the dependence of α on σ/m is described as a power law. For the low xoff subsample, we found α0 = 0.178, α1 = 0.20, and γ = 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 α into upper limits on σ/m, we converted each value of ⟨αrel⟩ in the output chain into a value of σ/m through Eq. (6) and constructed a posterior distribution for σ/m. As expected, the posterior distribution is consistent with zero. We then computed a one-tailed 95% upper limit as the 95th percentile of the posterior distribution,

σ m < 0.19 cm 2 g − 1 , at v DM ≈ 1000 ( 95 % confidence level ) . $$ \begin{aligned} \frac{\sigma }{m} < 0.19 \, \mathrm{cm}^2\,\mathrm{g}^{-1},\ \mathrm{at\ } {v}_{\mathrm{DM} } \approx 1000\ \mathrm{(95\%\,confidence\, level)}. \end{aligned} $$(7)

Similar constraints on σ/m are obtained if we use a different method to relate α to σ/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 α 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 σ/m ≲ 1 cm2 g−1 (e.g., Randall et al. 2008) or σ/m <  0.47 cm2 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 σ/m ≲ 0.39 cm2 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 σ/m ≲ 0.13 cm2 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 cm2 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 cm2 g−1 at vDM  ∼  10 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, ∝ v DM − 4 $ {\propto} {v}_{\mathrm{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 α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.

  • Most X-COP clusters have an Einasto index of 0.15 <  α <  0.3, with mean ⟨α⟩=0.22 ± 0.04. The substantial intrinsic scatter, σ α = 0 . 14 − 0.05 + 0.10 $ \sigma_{\alpha}=0.14_{-0.05}^{+0.10} $, is dominated by a few outliers with α ≳ 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 ⟨α⟩=0.19 ± 0.03. Their intrinsic scatter, σ α = 0 . 06 − 0.03 + 0.05 $ \sigma_\alpha=0.06_{-0.03}^{+0.05} $, is so small that α appears close to universal.

  • In CDM and SIDM simulations, we find that cluster mass profiles become more curved as the SIDM cross section σ/m increases, with α ≈ 0.18 + 0.20((σ/m)/(1 cm2 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 σ/m <  0.19 cm2 g−1 (95% c.l.) on the cross section for DM interactions at a collision velocity vDM ≈ 1000 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

  1. Andrade, K. E., Fuson, J., Gad-Nasr, S., et al. 2022, MNRAS, 510, 54 [Google Scholar]
  2. Angelinelli, M., Vazza, F., Giocoli, C., et al. 2020, MNRAS, 495, 864 [NASA ADS] [CrossRef] [Google Scholar]
  3. Biffi, V., Borgani, S., Murante, G., et al. 2016, ApJ, 827, 112 [NASA ADS] [CrossRef] [Google Scholar]
  4. Booth, C. M., & Schaye, J. 2009, MNRAS, 398, 53 [Google Scholar]
  5. Bradač, M., Allen, S. W., Treu, T., et al. 2008, ApJ, 687, 959 [CrossRef] [Google Scholar]
  6. Brown, S. T., McCarthy, I. G., Diemer, B., et al. 2020, MNRAS, 495, 4994 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bullock, J. S., & Boylan-Kolchin, M. 2017, ARA&A, 55, 343 [Google Scholar]
  8. Buote, D. A., & Humphrey, P. J. 2012, MNRAS, 421, 1399 [NASA ADS] [CrossRef] [Google Scholar]
  9. Cassano, R., Ettori, S., Giacintucci, S., et al. 2010, ApJ, 721, L82 [Google Scholar]
  10. Davé, R., Spergel, D. N., Steinhardt, P. J., & Wandelt, B. D. 2001, ApJ, 547, 574 [CrossRef] [Google Scholar]
  11. Eckert, D., Ettori, S., Pointecouteau, E., et al. 2017, Astron. Nachr., 338, 293 [Google Scholar]
  12. Eckert, D., Ghirardini, V., Ettori, S., et al. 2019, A&A, 621, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Eckert, D., Finoguenov, A., Ghirardini, V., et al. 2020, Open J. Astrophys., 3, 12 [NASA ADS] [CrossRef] [Google Scholar]
  14. Eckert, D., Ettori, S., Pointecouteau, E., van der Burg, R. F. J., & Loubser, S. I. 2022, A&A, 662, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Einasto, J. 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87 [NASA ADS] [Google Scholar]
  16. Ettori, S., Ghirardini, V., Eckert, D., et al. 2019, A&A, 621, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Euclid Collaboration (Adam, R., et al.) 2019, A&A, 627, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Gastaldello, F., Limousin, M., Foex, G., et al. 2014, MNRAS, 442, L76 [CrossRef] [Google Scholar]
  19. Ghirardini, V., Ettori, S., Eckert, D., et al. 2018, A&A, 614, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Ghirardini, V., Eckert, D., Ettori, S., et al. 2019, A&A, 621, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Harvey, D., Massey, R., Kitching, T., Taylor, A., & Tittley, E. 2015, Science, 347, 1462 [Google Scholar]
  22. Harvey, D., Robertson, A., Massey, R., & McCarthy, I. G. 2019, MNRAS, 488, 1572 [NASA ADS] [CrossRef] [Google Scholar]
  23. Hurier, G., Macías-Pérez, J. F., & Hildebrandt, S. 2013, A&A, 558, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Kahlhoefer, F., Schmidt-Hoberg, K., Frandsen, M. T., & Sarkar, S. 2014, MNRAS, 437, 2865 [CrossRef] [Google Scholar]
  25. Kaplinghat, M., Keeley, R. E., Linden, T., & Yu, H.-B. 2014, Phys. Rev. Lett., 113, 021302 [NASA ADS] [CrossRef] [Google Scholar]
  26. Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340 [Google Scholar]
  27. Lau, E. T., Kravtsov, A. V., & Nagai, D. 2009, ApJ, 705, 1129 [NASA ADS] [CrossRef] [Google Scholar]
  28. Le Brun, A. M. C., McCarthy, I. G., Schaye, J., & Ponman, T. J. 2014, MNRAS, 441, 1270 [NASA ADS] [CrossRef] [Google Scholar]
  29. Loubser, S. I., Babul, A., Hoekstra, H., et al. 2020, MNRAS, 496, 1857 [Google Scholar]
  30. Ludlow, A. D., & Angulo, R. E. 2017, MNRAS, 465, L84 [Google Scholar]
  31. Ludlow, A. D., Bose, S., Angulo, R. E., et al. 2016, MNRAS, 460, 1214 [Google Scholar]
  32. Macciò, A. V., Dutton, A. A., van den Bosch, F. C., et al. 2007, MNRAS, 378, 55 [Google Scholar]
  33. Mamon, G. A., & Łokas, E. L. 2005, MNRAS, 362, 95 [Google Scholar]
  34. Mantz, A. B., Allen, S. W., & Morris, R. G. 2016, MNRAS, 462, 681 [CrossRef] [Google Scholar]
  35. McCarthy, I. G., Schaye, J., Bird, S., & Le Brun, A. M. C. 2017, MNRAS, 465, 2936 [Google Scholar]
  36. McCarthy, I. G., Bird, S., Schaye, J., et al. 2018, MNRAS, 476, 2999 [NASA ADS] [CrossRef] [Google Scholar]
  37. Meneghetti, M., Yoshida, N., Bartelmann, M., et al. 2001, MNRAS, 325, 435 [NASA ADS] [CrossRef] [Google Scholar]
  38. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [Google Scholar]
  39. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [Google Scholar]
  40. Navarro, J. F., Hayashi, E., Power, C., et al. 2004, MNRAS, 349, 1039 [Google Scholar]
  41. Nelson, K., Lau, E. T., & Nagai, D. 2014, ApJ, 792, 25 [NASA ADS] [CrossRef] [Google Scholar]
  42. O’Hara, T. B., Mohr, J. J., Bialek, J. J., & Evrard, A. E. 2006, ApJ, 639, 64 [Google Scholar]
  43. Planck Collaboration V. 2013, A&A, 550, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Planck Collaboration XXIX. 2014, A&A, 571, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Predehl, P., Andritschke, R., Arefiev, V., et al. 2021, A&A, 647, A1 [EDP Sciences] [Google Scholar]
  47. Quynh Lan, N., Mathews, G. J., Arielle Phillips, L., et al. 2021, Mod. Phys. Lett. A, 36, 2130001 [NASA ADS] [CrossRef] [Google Scholar]
  48. Ragagnin, A., Saro, A., Singh, P., & Dolag, K. 2021, MNRAS, 500, 5056 [Google Scholar]
  49. Randall, S. W., Markevitch, M., Clowe, D., Gonzalez, A. H., & Bradač, M. 2008, ApJ, 679, 1173 [NASA ADS] [CrossRef] [Google Scholar]
  50. Rasia, E., Ettori, S., Moscardini, L., et al. 2006, MNRAS, 369, 2013 [CrossRef] [Google Scholar]
  51. Rasia, E., Meneghetti, M., & Ettori, S. 2013, Astron. Rev., 8, 40 [Google Scholar]
  52. Ren, T., Kwa, A., Kaplinghat, M., & Yu, H.-B. 2019, Phys. Rev. X, 9, 031020 [NASA ADS] [Google Scholar]
  53. Robertson, A., Massey, R., & Eke, V. 2017, MNRAS, 467, 4719 [NASA ADS] [CrossRef] [Google Scholar]
  54. Robertson, A., Harvey, D., Massey, R., et al. 2019, MNRAS, 488, 3646 [NASA ADS] [CrossRef] [Google Scholar]
  55. Robertson, A., Massey, R., Eke, V., Schaye, J., & Theuns, T. 2021, MNRAS, 501, 4610 [NASA ADS] [CrossRef] [Google Scholar]
  56. Rocha, M., Peter, A. H. G., Bullock, J. S., et al. 2013, MNRAS, 430, 81 [NASA ADS] [CrossRef] [Google Scholar]
  57. Sagunski, L., Gad-Nasr, S., Colquhoun, B., Robertson, A., & Tulin, S. 2021, JCAP, 2021, 024 [CrossRef] [Google Scholar]
  58. Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJ Comput. Sci., 2, e55 [Google Scholar]
  59. Spergel, D. N., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760 [NASA ADS] [CrossRef] [Google Scholar]
  60. Tulin, S., & Yu, H.-B. 2018, Phys. Rep., 730, 1 [Google Scholar]
  61. Umetsu, K., Medezinski, E., Nonino, M., et al. 2014, ApJ, 795, 163 [NASA ADS] [CrossRef] [Google Scholar]
  62. Valli, M., & Yu, H.-B. 2018, Nat. Astron., 2, 907 [NASA ADS] [CrossRef] [Google Scholar]
  63. Weißmann, A., Böhringer, H., Šuhada, R., & Ameglio, S. 2013, A&A, 549, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Yoshida, N., Springel, V., White, S. D. M., & Tormen, G. 2000, ApJ, 544, L87 [NASA ADS] [CrossRef] [Google Scholar]

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 Ptot by numerically integrating Eq. 1, and we assumed that the measured thermal pressure Pth is a fraction of the total pressure,

P th ( r ) = P tot ( r ) − P NT ( r ) , $$ \begin{aligned} P_{\rm th}(r) = P_{\rm tot}(r) - P_{\rm NT}(r), \end{aligned} $$(A.1)

with PNT(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 Pth profile with an Einasto model in the radial range [0.01 − 1]R200 and compared the retrieved value of α with the input one. In the left-hand panel of Fig. A.1, we show the ratio of the measured valued of α in the HSE assumption (αHSE) to the true input value (αTrue) for a range of input values of α. We can see that the HSE assumption typically overestimates α by 2 − 8% depending on the assumed nonthermal pressure profile and the true value of α, corresponding to differences of 0.01 for αTrue = 0.18. Therefore, we conclude that the hydrostatic assumption has little impact on the results presented here.

thumbnail Fig. A.1.

Impact of systematic uncertainties. Left: Ratio between the value of α derived under the assumption of hydrostatic equilibrium (αHSE) and the true value of α 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 α. 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).

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 iLoubser 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]R200 range with a single Einasto profile and quantified the changes in α with respect to the input value. We found that α is anticorrelated with the BCG stellar mass, with larger values for the stellar mass corresponding to lower values of α. For an input value α = 0.18, the retrieved value of α changes by ±0.02, which is smaller than the statistical uncertainties in our measurements.

All Tables

Table 1.

Normal priors on the Einasto fit parameters. Here Pm and dPm denote the outermost SZ pressure value and its error.

In the text
Table 2.

Result of Einasto fitting procedure for the 12 X-COP clusters.

In the text

All Figures

thumbnail 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 R500. In the right-hand panel, we show the posterior distributions of the fitted Einasto model parameters and the covariance between them.
In the text
thumbnail Fig. 2.

Measurements of the Einasto shape parameter α 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.
In the text
thumbnail Fig. 3.

Simulations of the structural properties of massive DM halos in the CDM and SIDM frameworks. Left-hand panel: Einasto shape parameter αEinasto for individual simulated halos with a varying SIDM cross section (0.1 cm2 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σ percentiles of αEinasto values for the various simulations runs.
In the text
thumbnail Fig. 4.

Relation between the Einasto shape parameter αEinasto and the DM self-interaction cross section σ/m for several sets of numerical simulations with varying σ/m (Robertson et al. 2021). The left-hand panel shows the fitted mean values of α 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 α on σ/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 αEinasto and the 1σ uncertainty in the regular X-COP subsample, with the corresponding constraints on σ/m.
In the text
thumbnail Fig. A.1.

Impact of systematic uncertainties. Left: Ratio between the value of α derived under the assumption of hydrostatic equilibrium (αHSE) and the true value of α 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 α. 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).
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.