Testing the Collisionless Nature of Dark Matter with the Radial Acceleration Relation in Galaxy Clusters

In this work, we use simulations run with four different SIDM models, as well as CDM, that were presented in Robertson et al. (2019). Three of the SIDM models have velocity-independent cross sections with isotropic scattering of σ/m = 0.1, 0.3, and 1 cm² g⁻¹, which we refer to as SIDM0.1, SIDM0.3, and SIDM1.0, respectively. The other SIDM model (hereafter vdSIDM) has a velocity-dependent and anisotropic cross section.

The vdSIDM differential cross section is (Robertson et al. 2017, 2021)

$$\begin{eqnarray}&&\displaystyle \frac{d\sigma }{d{\rm{\Omega }}}=\displaystyle \frac{{\sigma }_{0}}{4\pi {\left[1+({v}^{2}/{w}^{2}){\sin }^{2}(\theta /2)\right]}^{2}},\end{eqnarray} \tag{ 3 }$$

where w is a characteristic velocity below which the scattering is approximately isotropic with σ = σ₀. For collision velocities greater than w, scattering becomes anisotropic (favoring scattering by small angles) and the cross section decreases. Our vdSIDM model has σ₀/m = 3.04 cm² g⁻¹ and w = 560 km s⁻¹, which was chosen to reproduce the best-fit cross section in Kaplinghat et al. (2016).

To understand the macroscopic behavior of anisotropic particle interactions, we can introduce the concept of momentum-transfer cross section for vdSIDM (e.g., Robertson et al. 2017, 2021): ⁴

$$\begin{eqnarray}&&{\sigma }_{T}=2\int (1-| \cos \theta | )\displaystyle \frac{d\sigma }{d{\rm{\Omega }}}d{\rm{\Omega }}.\end{eqnarray} \tag{ 4 }$$

During collisions, the amount of momentum transferred in the direction of the collision is given as ${\rm{\Delta }}p=p(1-\cos \theta )$ with a scattering angle of θ. When the velocity increases, the vdSIDM cross section becomes more anisotropic, favoring scatter by small angles. Therefore, to better describe the effects of an anisotropic cross section, σ_T gives more weight to larger angle scattering, which contributes a larger amount of momentum transfer, while it downweights the small-angle scatter. After integrating over the solid angle, we obtain the following expressions for the total cross section (σ_tot) and the momentum-transfer cross section (σ_T ):

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{\mathrm{tot}}(v) & = & \displaystyle \frac{{\sigma }_{0}}{1+{v}^{2}/{w}^{2}},\\ {\sigma }_{T}(v) & = & {\sigma }_{0}\displaystyle \frac{4{w}^{4}}{{v}^{4}}\left[2\mathrm{ln}\left(1+\displaystyle \frac{{v}^{2}}{2{w}^{2}}\right)-\mathrm{ln}\left(1+\displaystyle \frac{{v}^{2}}{{w}^{2}}\right)\right].\end{array}\end{eqnarray} \tag{ 5 }$$

For a massive cluster with M₂₀₀ = 10¹⁵ M_⊙ at z = 0 and for a typical relative velocity between particle pairs of $\langle {v}_{\mathrm{rel}}\rangle \sim \sqrt{{{GM}}_{200}/{r}_{200}}$, we obtain σ_tot(〈v_rel〉) = 0.40 cm² g⁻¹ and σ_T (〈v_rel〉) = 0.25 cm² g⁻¹ for our vdSIDM model. We refer the reader to Robertson et al. (2021) for more details about the effective (velocity-averaged) cross section of vdSIDM halos.

We use N-body particle data from the Baryons and Haloes of Massive Systems (BAHAMAS) suite of cosmological hydrodynamical simulations (McCarthy et al. 2017, 2018) with WMAP 9 yr (Hinshaw et al. 2013) cosmology. BAHAMAS implements subgrid models for star formation and stellar and black hole feedback, and produces a good match to the observed stellar mass function, as well as the X-ray luminosities and gas mass fractions of galaxy groups/clusters. The simulations occupy large periodic boxes, 400h⁻¹ Mpc on a side. For the SIDM simulations, we use the BAHAMAS-SIDM suite (Robertson et al. 2019), which used the same initial conditions and subgrid models as BAHAMAS, but included an implementation of dark matter scattering. The parameters associated with the galaxy formation physics used in BAHAMAS-SIDM were kept the same as for the original BAHAMAS-CDM simulation. The friends-of-friends algorithm (Davis et al. 1985) with a linking length of 0.2 times the mean interparticle separation was run on each z = 0.375 simulation output. From each simulation we extract the 10,000 most massive friends-of-friends groups, which have spherical-overdensity masses in the range $12.5\lt {\mathrm{log}}_{10}[E(z){M}_{200}/{M}_{\odot }]\lt 15.3$.

For each halo we calculate the total enclosed mass profile, as well as the enclosed mass profile of the baryons. The center of the halo is defined by the location of the most gravitationally bound particle, and the enclosed masses are calculated at 101 different radii, logarithmically spaced between proper (as opposed to comoving) lengths of 0.1 kpc and 4 Mpc. In addition, the total density as a function of radius is calculated by taking the difference in total enclosed mass at two successive radii and dividing by the volume of the associated spherical shell. We consider the geometric mean of the inner and outer shell radii to be the radius at which this density is calculated. This density profile for each halo was used to make Figure 1.

Figure 2 shows the joint distribution of baryonic and total centripetal accelerations (g_bar, g_tot) derived from a subsample of group- and cluster-scale halos at z = 0.375, with masses E(z)M₂₀₀ > 5 × 10¹³ M_⊙. The results are shown separately for the CDM and four different SIDM models. The halo centripetal accelerations are logarithmically sampled at scales from r = 15 kpc to r = 4000 kpc.

Zoom In Zoom Out Reset image size

For each dark matter run, the magenta solid line represents the mean g_tot as a function of g_bar, or the halo RAR. The yellow dashed line represents the expectation corresponding to the cosmic mean ratio of total to baryonic mass densities, g_tot = (Ω_m/Ω_b)g_bar. In all cases, the RARs of simulated halos converge toward the cosmic mean, g_tot/g_bar = Ω_m/Ω_b, in the low-acceleration limit of g_bar ≲ 10⁻¹³ m s⁻².

The red dashed line shows the McGaugh et al. (2016) relation (Equation (2)) observed in spiral galaxies over the acceleration range $-12\lesssim {\mathrm{log}}_{10}({g}_{\mathrm{bar}}/{\rm{m}}\,{{\rm{s}}}^{-2})\lesssim -9$. Overall, the RARs of our simulated sample have a normalization that is higher than that observed at galaxy scales (McGaugh et al. 2016), suggesting a higher contribution from dark matter at a given baryonic acceleration.

The full sample of 10,000 BAHAMAS simulated halos spans a wide range of halo mass. We therefore split our sample into three mass bins: $12.5\lt {\mathrm{log}}_{10}[E(z){M}_{200}/{M}_{\odot }]\leqslant 13.7$, $13.7\lt {\mathrm{log}}_{10}[E(z){M}_{200}/{M}_{\odot }]\leqslant 14.7$, and ${\mathrm{log}}_{10}[E(z){M}_{200}/{M}_{\odot }]\gt 14.7$, to investigate the halo mass dependence of the RAR.

Figure 3 shows the resulting RARs in the three mass bins, for the five different dark matter models. For the lowest-mass bin, the RARs for different dark matter models largely overlap with each other, while for the highest-mass bin (corresponding to massive cluster halos) the slope of the RAR at high g_bar decreases with increasing SIDM cross section. The flattening feature in g_tot at high acceleration corresponds to the dark matter "cores" of approximately constant density at the center of SIDM halos. This distinguishing feature, being more significant in more massive halos, is consistent with the fact that the scattering rate is proportional to the local dark matter density and velocity dispersion (Equation (1)). This result indicates that the high-acceleration cluster-scale RAR can be used to probe the nature of dark matter. In the following analyses, we will focus on the 48 cluster-scale halos in the highest mass bin, which are more sensitive to the SIDM cross section.

Zoom In Zoom Out Reset image size

The velocity dependence of the vdSIDM cross section is apparent in Figure 3. For high-mass halos, vdSIDM behaves most similarly to SIDM0.3, while for the group-scale halos, vdSIDM halos are most like SIDM1.0 halos. This behavior is consistent with Robertson et al. (2019; see their Figure 1) and reflects the fact that the vdSIDM cross section decreases with increasing relative velocity between dark matter particles, and relative velocities are larger in more massive halos.

In the above analyses, a minimum radius of r = 15 kpc is applied. From cluster observations, the enclosed mass in the inner regions r ∈ [15, 100] kpc is difficult to measure. We therefore study the RARs for the same sets of halos discussed above, but with three larger minimum radii.

In Figure 4, we plot the RARs for the high-mass objects with radial cuts of r_cut = 15 kpc, 58 kpc, 98 kpc, and 206 kpc. The radial cut of 98 kpc corresponds approximately to the regime of combined strong- and weak-lensing analyses, while a radial cut of 206 kpc represents the regime of weak-lensing-only analyses for clusters at z ≳ 0.1. We recall that beyond r ∼ 100 kpc, it is challenging to distinguish the CDM and SIDM models by measuring the mass density profile of galaxy clusters, as shown in Figure 1. The bottom left panel of Figure 4 shows that even beyond r ∼ 100 kpc, the slope of the RAR for SIDM1.0 is shallower than that for CDM. Deviations in the RAR between CDM and SIDM are thus more significant than the conventional cusp–core features in the density profiles at larger radii. For the radial range of weak-lensing-only measurements, the discrepancy becomes tiny and would be almost impossible to detect. Hence, a combined strong- and weak-lensing analysis is typically required to distinguish SIDM and CDM in terms of the RAR, because this enables the measurement of the total enclosed mass down to and below approximately 100 kpc.

Zoom In Zoom Out Reset image size

The key features used to distinguish different dark matter models are largely driven by information in the inner region of cluster halos. In Appendix A, we further investigate whether the correlation between the acceleration ratio g_tot/g_bar and the clustercentric distance r, namely the g_tot/g_bar–r relation, can provide sufficient information to distinguish between different dark matter models. We find that the radial g_tot(r)/g_bar(r) (or M_tot(r)/M_bar(r)) profiles of cluster halos derived from different dark matter runs of the BAHAMAS simulation significantly overlap with each other, even in their innermost region, and hence the g_tot/g_bar–r relation has a much weaker sensitivity to the SIDM cross section than the g_tot–g_bar relation.

We characterize the halo RARs obtained in Section 3.2, focusing on cluster halos in the highest-mass bin with E(z)M₂₀₀ > 5 × 10¹⁴ M_⊙ at z = 0.375. To this end, we assume a power-law function of the form

$$\begin{eqnarray}&&{\mathrm{log}}_{10}({g}_{\mathrm{tot}}/{\rm{m}}\,{{\rm{s}}}^{-2})={b}_{1}{\mathrm{log}}_{10}({g}_{\mathrm{bar}}/{\rm{m}}\,{{\rm{s}}}^{-2})+{b}_{0},\end{eqnarray} \tag{ 7 }$$

with b ₁ the logarithmic slope and b ₀ the intercept. At g_bar ∼ 10⁻¹¹ m s⁻², the mean cluster RARs for all the dark matter runs converge to a power law with b ₁ ≈ 1.15 and b ₀ ≈ 2.68 (see the top panel of Figure 3). Then, the logarithmic slope begins to flatten gradually at g_bar ≳ 10⁻¹¹ m s⁻². At g_bar ≳ 10^−10.6 m s⁻², the cluster RAR is found to be highly sensitive to the SIDM cross section.

In Table 1, we summarize the best-fit values of b ₁ and b ₀ for the CDM and four SIDM runs characterized in the high-acceleration regime g_bar > 10^−10.6 m s⁻². It should be noted that we also fitted the mean RARs at g_bar > 10^−10.6 m s⁻² using the McGaugh et al. (2016) relation (Equation (2)), finding that only the CDM case can be well described by this functional form with an acceleration scale of g_† = (1.42 ± 0.06) × 10⁻⁹ m s⁻², which is much higher than the characteristic acceleration scale g_† ≈ 1.2 × 10⁻¹⁰ m s⁻² observed at galaxy scales (Section 1).

Table 1 also lists the levels of intrinsic scatter Δ_int around the mean relations obtained for the five different dark matter runs. In all cases, we find a remarkably tight distribution in $\mathrm{log}{g}_{\mathrm{tot}}$–$\mathrm{log}{g}_{\mathrm{bar}}$ space, with a slight increase in Δ_int with increasing cross section. For the CDM, SIDM0.1, and SIDM0.3 models, the values of Δ_int agree within the errors with ${{\rm{\Delta }}}_{\mathrm{int}}={0.064}_{-0.012}^{+0.013}$ (in dex) determined for the CLASH sample (Tian et al. 2020).

Here we investigate the evolution of the RAR by analyzing simulated halos at two redshifts, z = 0 and z = 0.375. For this purpose, we focus on massive cluster halos with E(z)M₂₀₀ > 5 × 10¹⁴ M_⊙ (Sections 3.2 and 3.3). At z = 0, we have a total of 82 cluster halos in this subsample.

Figure 5 shows the comparison of the cluster-scale RARs at z = 0 and z = 0.375. Compared to z = 0.375, we find a larger discrepancy at z = 0 between the CDM and SIDM results: the larger the scattering cross section, the lower the total acceleration at high g_bar (see also Section 3.2). The SIDM dependence increasing toward lower redshift is expected, because the radius out to which SIDM significantly affects density profiles is well described by the radius where there has been one scattering event per particle over the age of the halo (Robertson et al. 2021). This suggests that cluster RAR measurements for lower-z samples should provide a more sensitive test of the SIDM cross section. We find that the cluster RAR derived from vdSIDM at z = 0 resembles well that from SIDM0.3 at z = 0.

Zoom In Zoom Out Reset image size

With the aim of precisely determining the mass profiles of galaxy clusters using deep 16-band imaging with the Hubble Space Telescope (HST; Postman et al. 2012) and ground-based weak-lensing observations (Umetsu et al. 2014), a subsample of 20 CLASH clusters was X-ray-selected to be massive (>5 keV), with nearly concentric X-ray isophotes and a well-defined X-ray peak located close to the BCG. For this subsample, no lensing information was used a priori to avoid a biased sample selection. Cosmological hydrodynamical simulations suggest that the CLASH X-ray-selected subsample is mostly composed of relaxed systems (∼70%) and largely free of orientation bias (Meneghetti et al. 2014). Another subsample of five clusters was selected by their exceptional lensing strength to magnify galaxies at high redshift. These clusters often turn out to be dynamically disturbed massive systems (Umetsu 2020).

The CLASH sample spans nearly an order of magnitude in mass, 5 ≲ M₂₀₀/10¹⁴ M_⊙ ≲ 30. For each of the 25 clusters, HST weak- and strong-lensing data products are available in their central regions (Zitrin et al. 2015). Umetsu et al. (2016) combined wide-field weak-lensing data obtained primarily with Suprime-Cam on the Subaru telescope (Umetsu et al. 2014) and the HST weak- and strong-lensing constraints of Zitrin et al. (2015). For an observational determination of the RAR, Tian et al. (2020) combined the X-ray data products from Donahue et al. (2014) and the lensing data products from Umetsu et al. (2016), yielding a subsample of 20 CLASH clusters composed of 16 X-ray-selected and four lensing-selected systems. We note that five clusters of the CLASH sample were not included in the joint weak- and strong-lensing analysis performed by Umetsu et al. (2016) because of the lack of usable wide-field weak-lensing data. Consequently, they were also excluded from our analysis.

The CLASH subsample analyzed by Tian et al. (2020) has a median redshift of $\overline{z}=0.377$, which closely matches our simulation snapshot at z = 0.375. The typical resolution limit of the mass reconstruction set by the HST lensing data is 10'', which corresponds to ≈50 kpc at $\overline{z}=0.377$ (Umetsu et al. 2016). It was found by Umetsu et al. (2016) that the stacked lensing signal of the CLASH X-ray-selected subsample is well described by a family of cuspy, sharply steepening density profiles, such as the Navarro–Frenk–White (hereafter NFW; Navarro et al. 1996, 1997), Einasto (Einasto 1965), and DARKexp (Hjorth & Williams 2010) profiles. Of these, the NFW model best describes the CLASH lensing data (see also Umetsu & Diemer 2017). In contrast, the single power-law, cored isothermal, and Burkert models are statistically disfavored by the averaged lensing profile having a pronounced radial curvature.

For each of the 20 clusters, Umetsu et al. (2016) performed a spherical NFW fit to the reconstructed projected mass density profile by accounting for all relevant sources of uncertainty, including measurement errors, cosmic noise due to the projection of large-scale structure uncorrelated with the cluster, statistical fluctuations of the projected cluster lensing signal due to halo triaxiality, and correlated substructures.

In this analysis, we use the CLASH RAR data set (Tian et al. 2020) published in Tian et al. (2020), which contains N_data = 84 data points in $\mathrm{log}{g}_{\mathrm{tot}}$–$\mathrm{log}{g}_{\mathrm{bar}}$ space. Tian et al. (2020) extracted total and baryonic mass estimates where possible at r = 100, 200, 400, and 600 kpc. For each cluster, they also included a single constraint at r ≲ 30 kpc in the central BCG region. These data points are sufficiently well separated from each other so as to avoid oversampling and reduce correlations between adjacent data points.

The CLASH measurements of centripetal accelerations (g_bar, g_tot) are shown in Figure 6, along with our BAHAMAS predictions of the RAR for massive cluster halos with E(z)M₂₀₀ > 5 × 10¹⁴ M_⊙, derived from five different dark matter runs at z = 0.375. The best-fit CLASH RAR obtained by Tian et al. (2020) has ${b}_{1}={0.51}_{-0.05}^{+0.04}$, ${b}_{0}=-{4.26}_{-0.47}^{+0.46}$, and ${{\rm{\Delta }}}_{\mathrm{int}}={0.064}_{-0.012}^{+0.013}$ dex, which are in excellent agreement with the best-fit RAR for the BAHAMAS-CDM run characterized in the high-acceleration region of g_bar > 10^−10.6 m s⁻² (Table 1).

Zoom In Zoom Out Reset image size

To make a quantitative comparison between the BAHAMAS simulations and the CLASH observations, we define the χ² function as

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i=1}^{{N}_{\mathrm{data}}}\displaystyle \frac{{\left({\mathrm{log}}_{10}{g}_{\mathrm{tot},i}-{\mathrm{log}}_{10}{\widehat{g}}_{\mathrm{tot}.i}\right)}^{2}}{{{\rm{\Delta }}}_{i}^{2}+{{\rm{\Delta }}}_{\mathrm{int}}^{2}},\end{eqnarray} \tag{ 8 }$$

where i runs over all data points from the CLASH data set, g_tot,i is the total acceleration at the ith data point, Δ_i is the measurement uncertainty, Δ_int ≈ 0.064 is the intrinsic scatter for the CLASH sample determined by Tian et al. (2020), and ${\widehat{g}}_{\mathrm{tot},i}$ is the theoretical prediction for the total acceleration at g_bar = g_bar,i from the BAHAMAS simulations. The resulting χ² values evaluated for different dark matter runs are listed in Table 2.

To statistically characterize the level of agreement between data and simulations, we use frequentist measures of statistical significance. Specifically, we use the chance probability of exceeding the observed χ² value to quantify the significance of the match for each dark matter run. For each case, we calculate the probability to exceed (PTE), or the right-tailed p value, for a given value of χ² assuming the standard χ² probability distribution function. We adopt a significance threshold of α = 0.05 as the dividing line between satisfactory (PTE > 0.05) and unsatisfactory (PTE < 0.05) matches to the CLASH RAR data set. Since no optimization (or model fitting) is performed in our χ² evaluations, the number of degrees of freedom is N_data = 84 in all cases. The resulting values of PTE for each dark matter run are listed in Table 2. We find that the CDM run (PTE = 0.264) gives a satisfactory match, whereas all the SIDM runs give unacceptable matches to the CLASH data. Among the SIDM runs, SIDM0.1 has a PTE of 0.036, which is close to but slightly below the adopted threshold.

As a consistency check, we perform Monte Carlo simulations to derive the χ² distribution expected from the measurement errors $\{{{\rm{\Delta }}}_{i}\}{}_{i=1}^{{N}_{\mathrm{data}}}$ and the intrinsic scatter Δ_int for the CLASH RAR. In each simulation, we construct a synthetic data set $\{{g}_{\mathrm{bar},i}^{(\mathrm{MC})},{g}_{\mathrm{tot},i}^{(\mathrm{MC})}\}{}_{i=1}^{{N}_{\mathrm{data}}}$ by creating a Monte Carlo realization of random Gaussian noise ${n}_{i}\equiv {\mathrm{log}}_{10}{g}_{\mathrm{tot},i}^{(\mathrm{MC})}-{\mathrm{log}}_{10}{\widehat{g}}_{\mathrm{tot},i}=\pm \sqrt{{{{\rm{\Delta }}}_{i}}^{2}+{{\rm{\Delta }}}_{\mathrm{int}}^{2}}$ at ${\mathrm{log}}_{10}{g}_{\mathrm{bar},i}^{(\mathrm{MC})}={\mathrm{log}}_{10}{g}_{\mathrm{bar},i}^{(\mathrm{CLASH})}$ and then calculate the value of ${\chi }^{2}={\sum }_{i}{n}_{i}^{2}/({{\rm{\Delta }}}_{i}^{2}+{{\rm{\Delta }}}_{\mathrm{int}}^{2})$. We repeat this procedure 10⁶ times to generate a large set of Monte Carlo realizations and obtain the distribution of χ² values. In Table 2, we list the fraction of Monte Carlo realizations exceeding the observed value of χ² for each dark matter run. In all cases, the Monte Carlo fraction exceeding the χ² value is precisely consistent with the PTE calculated with the standard χ² distribution function. In the upper panel of Figure 7, we compare the χ² distribution of the CLASH data set constructed from our Monte Carlo simulations with the observed χ² values for the CDM and four SIDM models. This figure gives a visual summary of the χ² test (Table 2).

Zoom In Zoom Out Reset image size

We perform a likelihood-ratio test of SIDM models to quantify whether the inclusion of collisional features of dark matter is statistically warranted by the data. Velocity-independent SIDM models have one additional degree of freedom relative to the CDM model. For vdSIDM, there are two additional degrees of freedom relative to the CDM model. These models are reduced to the CDM model in the limit of a vanishing cross section. Table 3 lists for each SIDM model the differences in the χ² value ${\rm{\Delta }}{\chi }^{2}={\chi }^{2}-{\chi }_{\mathrm{CDM}}^{2}$ relative to the CDM model and the corresponding significance level. Compared with the fiducial model of CDM, SIDM0.1 is disfavored at a significance level of 4.1σ.

These results are based on the CLASH RAR at BCG–cluster scales, which includes, for each cluster, a single constraint in the central BCG region at r ≲ 30 kpc. However, since the typical resolution of CLASH mass reconstructions is Δr ≈ 50 kpc (Section 4.1), the mass distribution is not resolved in the BCG region and thus the CLASH constraints on g_tot at the BCG scale are model dependent to some extent. Moreover, the distribution of baryonic and dark matter in cluster cores is sensitive to baryonic physics (e.g., Cui et al. 2018).

We therefore repeat the tests described above using core-excised RAR data. Excluding the central r < 100 kpc region from the CLASH data set, we have 64 data points at r ∈ [100, 600] kpc. The results of the χ² test with the core-excised data set are summarized in Table 2. Of the five BAHAMAS dark matter runs, CDM and SIDM0.1 provide satisfactory matches to the core-excised CLASH RAR, at a significance level of α = 0.05. Excluding the BCG regions results in a weaker but still competitive constraint on the SIDM cross section (Table 3). With a likelihood-ratio test, we find that the core-excised CLASH RAR data disfavor the SIDM0.3 model at the 3.8σ level with respect to the CDM model.

We note that in the above analysis, the measurement uncertainty of ${\mathrm{log}}_{10}{g}_{\mathrm{bar}}$ (which is typically ∼7% of that of ${\mathrm{log}}_{10}{g}_{\mathrm{tot}}$ and thus negligible in the analysis) is not taken into account. In Appendix B, we perform an alternative analysis in terms of the g_tot/g_bar–g_bar relation, which is referred to as the mass discrepancy–acceleration relation (MDAR; see Famaey & McGaugh 2012). In this MDAR analysis, we can explicitly account for the uncertainty in the g_bar measurements. We find that the inclusion of the measurement uncertainty in g_bar has only a minor impact on the statistical inference and does not change our conclusions regarding the acceptance of dark matter models.

In this work, we analyzed the CLASH RAR data set of Tian et al. (2020), which consists of N_data = 84 data points (g_bar, g_tot) inferred from the multiwavelength CLASH observations of 20 high-mass galaxy clusters (Section 4.1). In Tian et al. (2020), the measurements of centripetal accelerations were sparsely sampled over a sufficiently wide range of clustercentric distances so as to reduce the covariance between adjacent data for each cluster. Thus, our comparison of the BAHAMAS simulations and CLASH measurements (Equation (8)) involves a two-step procedure, which ignores the covariance and does not fully exploit all the information contained in the data. As a result, our analysis is likely to overestimate the significance of our constraints. In principle, these limitations can be overcome by using a forward-modeling method. In particular, likelihood-free approaches based on forward simulations allow us to bypass the need for a direct evaluation of the likelihood function assuming Gaussian statistics, which avoids the complex derivation of the covariance matrix in an inherently complex problem (e.g., Tam et al. 2022).

Another potential source of systematic uncertainty is the smoothing of the inner density profile due to cluster miscentering (e.g., Johnston et al. 2007). Because of the CLASH selection, our cluster sample exhibits, on average, a small positional offset between the BCG and the X-ray peak, characterized by an rms offset of ∼40 kpc (Umetsu et al. 2014, 2016). This level of offset is comparable to the typical effective radius R_e of CLASH BCGs (〈R_e〉 ∼ 30 kpc; Section 4) but sufficiently small compared to the range of cluster radii of interest (say, r ≳ 100 kpc). Moreover, since the RAR method uses the same aperture to compare total and baryonic accelerations, the miscentering effect is not expected to significantly affect the SIDM constraint, although it could potentially contribute to the scatter in the RAR inferred from cluster observations.

The total and baryonic mass profiles of galaxy clusters, M_tot(< r) and M_bar(< r), inferred from observational data are model-dependent to some extent. Tian et al. (2020) found that the observed surface brightness distribution of CLASH BCGs is well described by a de Vaucouleurs profile, which is a special case of the Sérsic model, with a Sérsic index of n = 4. Accordingly, they modeled the 3D stellar mass distribution in the CLASH BCGs with a Hernquist profile, which closely resembles the de Vaucouleurs law in projection (Hernquist 1990; see their Figure 4), especially at ≳0.5R/R_e where Tian et al. (2020) extracted stellar mass estimates of the CLASH BCGs.

On the other hand, Tian et al. (2020) modeled the total mass distribution of CLASH clusters with an NFW profile, which best describes the stacked lensing profile of the CLASH sample (Umetsu et al. 2016; see Section 4.1). However, the choice of the cuspy NFW model implicitly assumes collisionless CDM, which could lead to a biased inference of the cluster RAR in an SIDM cosmology. In particular, this NFW assumption in our analysis is likely to underestimate the goodness of fit for the SIDM models. In future studies, it will thus be important to consider more flexible and self-consistent mass models, such as the density profiles of Einasto (1965) and Diemer & Kravtsov (2014), with an additional parameter to capture the radial curvature of the central density profile that depends on the SIDM cross section (Eckert et al. 2022b).

Self-interactions between dark matter particles are expected to make halos more spherical than triaxial halos of collisionless CDM, especially in the central region where the scattering rate is largest. Self-interactions in the optically thin regime also reduce the central dark matter densities, transforming a cusp into a core. Moreover, offsets between the galactic and dark-matter centroids in merging clusters can be used to constrain the SIDM cross section. Previous studies placed upper limits on the self-interaction cross section of dark matter particles using such observed density features in galaxy clusters (for a review, see Tulin & Yu 2018).

Peter et al. (2013) compared the halo ellipticities inferred from lensing and X-ray observations with cosmological simulations with SIDM cross sections of σ/m = 0.03, 0.1, and 1 cm² g⁻¹. They found that the strong-lensing measurement of the cluster ellipticity for MS 213723 (Miralda-Escudé 2002) is compatible with an SIDM cross section of σ/m = 1 cm² g⁻¹, whereas the X-ray shape measurement of the isolated elliptical galaxy NGC 720 (Buote et al. 2002) is consistent with σ/m = 0.1 cm² g⁻¹.

Using the galaxy–dark-matter offset measured in the moving subcluster of the Bullet Cluster, Randall et al. (2008) placed an upper limit of σ/m < 1.25 cm² g⁻¹ at the 68% confidence level (CL). Harvey et al. (2015) performed an ensemble analysis of offset measurements for 72 substructures in 30 systems, including both major and minor mergers, and set an upper limit of σ/m < 0.47 cm² g⁻¹ at the 95% CL. Wittman et al. (2018) revisited the analysis of Harvey et al. (2015) using more comprehensive data and carefully reinterpreted their refined offset measurements, finding that the SIDM constraint of Harvey et al. (2015) is relaxed to σ/m ≲ 2 cm² g⁻¹.

Analyzing observations of X-rays and the SZ effect, Eckert et al. (2022b) constrained the structural parameters of the mass density profiles for a sample of 12 massive X-COP clusters, assuming that the intracluster gas is in hydrostatic equilibrium. They used the BAHAMAS-SIDM simulations to construct an empirical scaling relation between the Einasto shape parameter and the velocity-independent SIDM cross section. With this relation and the assumption of hydrostatic equilibrium, they obtained an upper limit of σ/m < 0.19 cm² g⁻¹ at the 95% CL.

In this study, we have obtained competitive constraints on the SIDM cross section using cluster RAR measurements. By comparing the CLASH RAR data set with the mean RARs derived from the BAHAMAS simulations, we are able to reject an SIDM model with σ/m = 0.1 cm² g⁻¹ at the 4.1σ CL with respect to the CDM model. Excluding the central r < 100 kpc region, we find that an SIDM model with σ/m = 0.3 cm² g⁻¹ is disfavored at the 3.8σ level with respect to CDM.