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

The radial acceleration relation (RAR) represents a tight empirical relation between the inferred total and baryonic centripetal accelerations, $g_{\rm{tot}}=GM_{\rm{tot}}(<r)/r^2$ and $g_{\rm{bar}}=GM_{\rm{bar}}(<r)/r^2$, observed in galaxies and galaxy clusters. The tight correlation between these two quantities can provide insight into the nature of dark matter. Here we use BAHAMAS, a state-of-the-art suite of cosmological hydrodynamical simulations, to characterize the RAR in cluster-scale halos for both cold and collisionless dark matter (CDM) and self-interacting dark matter (SIDM) models. SIDM halos generally have reduced central dark matter densities, which reduces the total acceleration in the central region when compared with CDM. We compare the RARs in galaxy clusters simulated with different dark matter models to the RAR inferred from CLASH observations. Our comparison shows that the cluster-scale RAR in the CDM model provides an excellent match to the CLASH RAR obtained by Tian et al. including the high-acceleration regime probed by the brightest cluster galaxies (BCGs). By contrast, models with a larger SIDM cross-section yield increasingly poorer matches to the CLASH RAR. Excluding the BCG regions results in a weaker but still competitive constraint on the SIDM cross-section. Using the RAR data outside the central $r<100$kpc region, an SIDM model with $\sigma/m=0.3$cm$^{2}$g$^{-1}$ is disfavored at the $3.8\sigma$ level with respect to the CDM model. This study demonstrates the power of the cluster-scale RAR for testing the collisionless nature of dark matter.


INTRODUCTION
Based on modern cosmological studies (e.g. Hinshaw et al. 2013), dark matter is known to be the dominant matter component in the Universe, while its nature is still a mystery. In the current concordance cosmological model, Λ Cold Dark Matter (CDM), large-scale structures formed hierarchically, with dark matter halos growing through a series of mergers of smaller halos as well as accretion. This standard model provides a good description of the observed large-scale structure.
However, there are issues on smaller scales that are potentially challenging for the CDM model. For example, collisionless dark matter particles in CDM produce cuspy dark matter halos, where the density rises towards the halo center, which is inconsistent with the cored density profiles inferred for the halos hosting some dwarf galaxies. This is the so-where m is the dark matter particle mass. Therefore, with the highest ρ dm and v, massive galaxy clusters are crucial laboratories to search for dark matter self-interactions.
Galaxy clusters are the most massive gravitationally-bound structures resulting from the hierarchical formation process. About 85% of their mass content is invisible dark matter, with the remainder being baryons that are mostly in the form of X-ray emitting hot gas. Since cluster properties depend on the growth of structure, they contain an abundance of cosmological and astrophysical information. Several characteristic features of galaxy clusters have been used to test dark matter models, such as dark matter halo shapes (e.g. Peter et al. 2013;Umetsu et al. 2018), offsets between dark matter and galaxies in merging systems (e.g. Harvey et al. 2015;Wittman et al. 2018;Massey et al. 2018), the wobbling of brightest central galaxies ) and the amount and lensing efficiency of dark matter substructures (e.g. Jauzac et al. 2016;Tam et al. 2020;Meneghetti et al. 2020).
In the ΛCDM model, dark matter has negligible interaction with baryons, except for gravity. However, tight relations between the distribution of dark matter and of baryonic matter have been discovered. At the scale of spiral galaxies, the ratio of dynamical to baryonic masses, M tot (< r)/M bar (< r), is found to be tightly coupled with gravitational acceleration, whereas no clear correlation with other physical quantities, such as galaxy size, has been found to date (McGaugh 2004).
By analyzing rotation curves of 153 spiral galaxies, Mc-Gaugh et al. (2016) found a tight correlation between two independent observables, namely the centripetal acceleration g tot (r) = V 2 /r = GM tot (< r)/r 2 and the baryonic contribution to this acceleration g bar (r) = GM bar (< r)/r 2 : characterized by a characteristic acceleration scale, g † = 1.20 ± 0.24 × 10 −10 m s −2 . This empirical relation between the total and baryonic centripetal accelerations is referred to as the radial acceleration relation (RAR). Since then, much effort have been invested to study the RAR in various galaxy samples (e.g., Lelli et al. 2017;Rong et al. 2018;Chae et al. 2019;Brouwer et al. 2021;Oman et al. 2020). Hydrodynamical simulations in the ΛCDM framework have succeeded in reproducing the observed RAR of galaxies (e.g., Keller & Wadsley 2017;Ludlow et al. 2017;Garaldi et al. 2018;Dutton et al. 2019;Paranjape & Sheth 2021).
Recently, observational studies of the RAR have been extended to cluster-scale objects (e.g., Tian et al. 2020;Chan & Del Popolo 2020;Eckert et al. 2022a). Tian et al. (2020) studied the RAR for a subsample of 20 high-mass galaxy clusters targeted by the CLASH program (Postman et al. 2012). In their analysis, the total mass of each cluster is inferred from a combined analysis of strong and weak lensing data (Umetsu et al. 2016) and the baryonic mass from X-ray gas mass and stellar mass estimates (Donahue et al. 2014). Chan & Del Popolo (2020) analyzed X-ray data for a sample of 52 non-cool-core clusters. They obtained the cluster RAR using X-ray hydrostatic estimates for the total mass and X-ray gas mass estimates for the baryonic mass, ignoring the stellar mass contribution to the baryonic component. Eckert et al. (2022a) studied the total and baryonic mass distributions for a sample of 12 X-COP clusters with X-ray and Sunyaev-Zel'dovich (SZ) effect observations, accounting for the stellar mass contribution. They found a complex shape of the RAR that strongly departs from the RAR in galaxies.
All of these studies found that the characteristic acceleration scale g † in clusters is about an order of magnitude larger than that obtained from galaxy-scale objects. Observational RAR studies have also been extended to group-scale objects. Gopika & Desai (2021) found that g † of group-scale halos falls in between that found for galaxies and galaxy clusters. These observations suggest that there is no universal RAR that holds at all scales from galaxies to galaxy clusters.
Alternatively, the RAR observed at galaxy scales has been attributed to modified Newtonian dynamics (MOND; Milgrom 1983), which introduces a characteristic acceleration scale, g † , and modifies the dynamical law. However, MOND falls short in accounting for the total observed gravitational mass in galaxy clusters. For completeness it is worth mentioning other approaches to generalize the RAR to galaxy clusters, especially the eMOND framework of Zhao & Famaey (2012), where the parameter g † is allowed to be a monotonic increasing function of the system's escape velocity such that g † is several times greater in clusters than in field galaxies. Hodson & Zhao (2017) followed up this idea and presented fits to nearby relaxed clusters in the Chandra Sample (Vikhlinin et al. 2006). They succeed in this by tailoring g † as a specific function of the Newtonian potential of the gas and the BCG galaxies (see their Eq. 20 and their Figs A.1-A.11 for 20 galaxy clusters).
In this study, we use cosmological hydrodynamical simulations to study the RAR for simulated halos in both CDM and SIDM scenarios. We aim to explore a new method for constraining the collisionless nature of dark matter using the cluster-scale RAR, as well as to compare the RARs derived from numerical simulations with multiwavelength cluster observations.
This paper is organized as follows. Section 2 introduces the simulation data sets we use in this work. Section 3 shows the results of the halo RAR obtained from the simulations. Section 4 compares the theoretical predictions from the simulations with observational data from the CLASH program.
In Section 5, we discuss the results and implications of our findings. Finally a summary is given in Section 6.
Throughout this paper, we assume a Wilkinson Microwave Anisotropy Probe (WMAP) 9-year ΛCDM cosmology (Hinshaw et al. 2013) with Ω m = 0.287, Ω Λ = 0.713, and a Hubble constant of H 0 = 100 h km s −1 Mpc −1 with h = 0.693. We denote the critical density of the universe at a particular redshift z as ρ c (z) = 3H 2 (z)/(8πG) with H(z) the Hubble function. We also define the dimensionless expansion function as E(z) = H(z)/H 0 . We adopt the standard notation M ∆ to denote the total mass enclosed within a sphere of radius r ∆ within which the mean overdensity is ∆ × ρ c (z). We use "ln" to denote the natural logarithm.

SIDM Models
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 2 g −1 , 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(Robertson et al. , 2021 where w is a characteristic velocity below which the scattering is approximately isotropic with σ = σ 0 . For collision velocities greater than w, scattering becomes anisotropic (favoring scattering by small angles) and the cross-section decreases. Our vdSIDM model has σ 0 /m = 3.04 cm 2 g −1 and w = 560 km s −1 , 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. 2017Robertson et al. , 2021: 1 During collisions, the amount of momentum transferred along the collision direction is given as ∆p = p(1 − cos θ) 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 1 In Robertson et al. (2017), σ T is referred to as the modified momentumtransfer cross-section and denoted as σT .
effects of an anisotropic cross-section, σ T gives more weight to larger angle scattering which contributes 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 ): (5) For a massive cluster with M 200 = 10 15 M ⊙ at z = 0 and for a typical relative velocity between particle pairs of ⟨v rel ⟩ ∼ GM 200 /r 200 , we obtain σ tot (⟨v rel ⟩) = 0.40 cm 2 g −1 and σ T (⟨v rel ⟩) = 0.25 cm 2 g −1 for our vdSIDM model. We refer the reader to Robertson et al. (2021) for more details about the effective (velocity averaged) cross-section of vd-SIDM halos.

BAHAMAS Simulations
We use N -body particle data from the BAryons And HAloes of MAssive Systems (BAHAMAS) suite of cosmological hydrodynamical simulations (McCarthy et al. 2017(McCarthy et al. , 2018 with WMAP 9-year (Hinshaw et al. 2013) cosmology. BAHAMAS implements sub-grid 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, 400 h −1 Mpc on a side. For the SIDM simulations, we use the BAHAMAS-SIDM suite ) which used the same initial conditions and sub-grid 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 inter-particle separation was run on each z = 0.375 simulation output. From each simulation we extract the 10,000 most massive friendsof-friends groups, which have spherical-overdensity masses in the range of 12.5 < log 10 [E(z)M 200 /M ⊙ ] < 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 geomet- Figure 1. Mass density profiles for different matter components of cluster-scale halos with masses E(z)M200 > 5 × 10 14 M⊙ at z = 0.375, shown for simulations with five different dark matter models. For each component, the dashed line represents the mean density profile and the shaded region shows the standard deviation around the mean profile. The large scatter in the innermost gas density is due to the AGN feedback in central galaxies, which removes most of the surrounding gas thus leading to low central gas densities for some halos. ric 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.

CHARACTERIZATION OF THE CLUSTER-SCALE RAR IN THE BAHAMAS SIMULATIONS
With the enclosed total and baryonic mass profiles M tot (< r) = M dm (< r) + M bar (< r) and M bar (< r) measured for each individual halo (Section 2.2), we calculate their total and baryonic centripetal acceleration profiles as Equation (6) should be regarded as the definition of g tot (r) and g bar (r), not the result of assuming spherical symmetry. In this section, we aim to characterize the relationship between g tot and g bar for samples of halos selected from the BAHAMAS-CDM and -SIDM runs, focusing on massive cluster-scale objects. 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 200 > 5 × 10 13 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.

RARs in CDM and SIDM Halos
For each dark matter run, the magenta solid line represents the mean g tot as a function of g bar relation, 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 towards the cosmic mean, g tot /g bar = Ω m /Ω b , in the low-acceleration limit of g bar < ∼ 10 −13 m s −2 . The red dashed line shows the McGaugh et al. (2016) relation (Equation (2)) observed in spiral galaxies over the acceleration range of −12 < ∼ log 10 (g bar /m s −2 ) < ∼ − 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.

Mass Dependence of the Halo RAR
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 < log 10 [E(z)M 200 /M ⊙ ] ≤ 13.7, 13.7 < log 10 [E(z)M 200 /M ⊙ ] ≤ 14.7, and log 10 [E(z)M 200 /M ⊙ ] > 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 approximately constant-density dark matter "cores" 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.
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 behaviour 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, respectively. The radial cut of 98 kpc corresponds approximately to the regime of combined strong and weaklensing 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 is 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.
The key features used to distinguish different dark matter models are largely driven by information in the inner  (2)). The yellow dashed line shows the linear relationship gtot = (Ωm/Ω b )g bar corresponding to the cosmic ratio between total and baryonic mass densities. The green dashed line is the best-fit RAR observed for 20 high-mass galaxy clusters in the CLASH sample (Tian et al. 2020). The black dashed line shows the one-to-one relation.  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 compared to the g tot -g bar relation, the g tot /g bar -r relation has a much weaker sensitivity to the SIDM cross-section.

Power-law Characterization of the Cluster-scale RAR
We characterize the halo RARs obtained in Section 3.2, focusing on cluster halos in the highest-mass bin with E(z)M 200 > 5 × 10 14 M ⊙ at z = 0.375. To this end, we assume a power-law function of the form: with b 1 the logarithmic slope and b 0 the intercept. At g bar ∼ 10 −11 m s −2 , the mean cluster RARs for all the dark matter runs converge to a power-law with b 1 ≈ 1.15 and b 0 ≈ 2.68 (see the top panel of Figure 3). Then, the logarithmic slope begins to flatten gradually at g bar > ∼ 10 −11 m s −2 . At g bar > ∼ 10 −10.6 m s −2 , 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 1 and b 0 for the CDM and four SIDM runs characterized in the high-acceleration regime g bar > 10 −10.6 m s −2 . It should NOTE-For each model, the RAR is derived for a subsample of cluster halos at z = 0.375 with masses E(z)M200 > 5 × 10 14 M⊙. A cut-off radius of rcut = 15 kcp is used. The quantities b1 and b0 represent the slope and intercept of the power-law fit (Equation (7)) in the high-acceleration region of log 10 (g bar /m s −2 ) > −10.6. The ∆int parameter is the intrinsic scatter in dex.
be noted that we also fitted the mean RARs at g bar > 10 −10.6 m s −2 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 −9 m s −2 , which is much higher than the characteristic acceleration scale g † ≈ 1.2 × 10 −10 m s −2 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 log g tot -log g 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 ∆ int = 0.064 +0.013 −0.012 (in dex) determined for the CLASH sample (Tian et al. 2020).

Evolution of the Cluster-scale RAR
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 200 > 5 × 10 14 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 r > 15 kpc r > 58 kpc r > 98 kpc r > 206 kpc 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.

COMPARISON WITH OBSERVATIONAL DATA
In this section, we compare the cluster-scale RARs derived from the BAHAMAS simulations with observations of galaxy clusters. An observational determination of the cluster RAR relies on accurate measurements of both total and baryonic mass profiles, M tot (< r) and M bar (< r), for a sizable sample of galaxy clusters. The baryonic mass in galaxy clusters is dominated by the hot intracluster gas, except in their central region where the BCG dominates the total bary-onic mass (e.g., Sartoris et al. 2020). 2 . Thus, baryonic mass estimates for both components are essential. Furthermore, unbiased estimates for the total mass in galaxy clusters are critical for a robust determination of the cluster RAR.
In this context, Tian et al. (2020) determined the RAR at BCG-cluster scales for a sample of 20 high-mass galaxy clusters targeted by the CLASH program, by combining weak and strong gravitational lensing data (Umetsu et al. 2014(Umetsu et al. , 2016Merten et al. 2015;Zitrin et al. 2015), X-ray gas mass measurements (Donahue et al. 2014), and BCG stellar mass estimates. The stellar mass contribution from member galaxies was statistically corrected for. To date, this is the only study that uses gravitational lensing data to directly probe the total acceleration g tot in galaxy clusters. By contrast, other studies based on hydrostatic estimates for g tot could potentially bias the true underlying RAR. In this study, we thus focus on the CLASH RAR studied by Tian et al. (2020).

CLASH Data
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 groundbased 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 200 /10 14 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 4 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 in our analysis.
The CLASH subsample analyzed by Tian et al. (2020) has a median redshift of 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 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 (Navarro et al. 1996(Navarro et al. , 1997, hereafter NFW), 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 den-sity 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. (2021), which contains N data = 84 data points in log g tot -log g bar space. Tian et al. (2020) extracted total and baryonic mass estimates where possible at r = 100 kpc, 200 kpc, 400 kpc, 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 BA-HAMAS predictions of the RAR for massive cluster halos with E(z)M 200 > 5 × 10 14 M ⊙ , derived from five different dark matter runs at z = 0.375. The best-fit CLASH RAR obtained by Tian et al. (2020) −0.47 , and ∆ int = 0.064 +0.013 −0.012 dex, which is 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 −2 (Table 1).

Statistical Comparison
To make a quantitative comparison between the BA-HAMAS simulations and the CLASH observations, we define the χ 2 function as 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 g tot,i is the theoretical prediction for the total acceleration at g bar = g bar,i from the BAHAMAS simulations. The resulting χ 2 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 χ 2 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 righttailed p-value, for a given value of χ 2 assuming the standard χ 2 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 χ 2 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 that is close to but slightly below the adopted threshold.
As a consistency check, we perform Monte-Carlo simulations to derive the χ 2 distribution expected from the measurement errors {∆ i } N data i=1 and the intrinsic scatter ∆ int for the CLASH RAR. In each simulation, we construct a synthetic data set {g by creating a Monte-Carlo realization of random Gaussian noise n i ≡ log 10 g (MC) tot,i − log 10 g tot,i = ± ∆ i 2 + ∆ 2 int at log 10 g (MC) bar,i = log 10 g (CLASH) bar,i and then calculate the value of χ 2 = i n 2 i /(∆ 2 i + ∆ 2 int ). We repeat this procedure 10 6 times to generate a large set of Monte-Carlo realizations and obtain the distribution of χ 2 values. In Table 2, we list the fraction of Monte-Carlo realizations exceeding the observed value of χ 2 for each dark matter run. In all cases, the Monte-Carlo fraction of exceeding the χ 2 value is precisely consistent with the PTE calculated with the standard χ 2 distribution function. In the upper panel of Figure 7, we compare the χ 2 distribution of the CLASH data set constructed from our Monte-Carlo simulations with the observed χ 2 values for the CDM and four SIDM models. This figure gives a visual summary of the χ 2 test (Table 2).
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. Velocityindependent 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 the vanishing cross-section. Table 3 lists for each SIDM model the differences in the χ 2 value ∆χ 2 = χ 2 − χ 2 CDM 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 BCGcluster 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).  a Observed χ 2 value between the CLASH data and each dark matter model.
b Probability to exceed the observed χ 2 value assuming the standard χ 2 probability distribution function.
c Fraction of Monte-Carlo realizations exceeding the observed value of χ 2 .
We therefore repeat the tests described above using coreexcised 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 χ 2 test with the coreexcised data set are summarized in Table 2. Of the five BA-HAMAS dark matter runs, CDM and SIDM0.1 provide satisfactory matches to the core-excised CLASH RAR, at a sig-nificance 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 log 10 g bar (which is typically ∼ 7% of that of log 10 g 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.

DISCUSSION
In this section, we first discuss current limitations and possible improvements of SIDM constraints from measurements of the cluster RAR. Then, we compare our results to previous astrophysical constraints on SIDM.

Possible Systematics and Improvements
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 forwardmodeling 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 co-variance 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 X-ray peak, characterized by an rms offset of ∼ 40 kpc (Umetsu et al. 2014(Umetsu et al. , 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 Sersic model, with a Sersic 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 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 Einasto (Einasto 1965) and Diemer-Kravtsov (Diemer & Kravtsov 2014) density profiles, 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).

Comparison with Other Work
Self-interactions between dark matter particles are expected to make halos more spherical compared to 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 2 g −1 . They found that the strong-lensing measurement of the cluster ellipticity for MS 2137-23 (Miralda-Escudé 2002) is compatible with an SIDM cross-section of σ/m = 1 cm 2 g −1 , whereas the X-ray shape measurement of the isolated elliptical galaxy NGC 720 (Buote et al. 2002) is consistent with σ/m = 0.1 cm 2 g −1 .
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 2 g −1 at the 68% 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 2 g −1 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 2 g −1 .
Analyzing X-ray and SZ effect observations, 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 2 g −1 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 2 g −1 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 2 g −1 is disfavored at the 3.8σ level with respect to CDM.

SUMMARY AND CONCLUSIONS
The RAR is a tight relation between the total and baryonic centripetal accelerations, g tot and g bar , inferred for galaxies and galaxy clusters (McGaugh et al. 2016;Tian et al. 2020). This tight empirical correlation offers a new possibility of testing the collisionless nature of dark matter at galaxy-cluster scales.
As a first step toward this goal, we studied in this paper the RAR in simulated halos for both CDM and SIDM models, using the BAHAMAS suite of cosmological hydrodynamical simulations (McCarthy et al. 2017(McCarthy et al. , 2018Robertson et al. 2019). We analyzed simulations at z = 0 and z = 0.375 run with four different SIDM models (SIDM0.1, SIDM0.3, SIDM1.0, and vdSIDM; see Section 2), as well as collisionless CDM.
For each dark matter model, we have determined the mean g tot as a function of g bar , or the halo RAR, in halos of different mass bins (Section 3; see Figure 2). We find that the slope of the halo RAR at high acceleration (g bar > ∼ 10 −10.6 m s −2 ) decreases with increasing SIDM cross-section. This flattening feature at high g bar is more significant in more massive halos at lower redshift (Figures 3 and 5), consistent with the fact that the scattering rate is proportional to the local dark matter density and velocity dispersion (Equation (1)). This suggests that the high-g bar cluster-scale RAR for lowredshift samples can be used to probe the nature of dark matter.
Focusing on massive cluster halos at z = 0.375, we have also characterized the slope (b 1 ), intercept (b 0 ), and intrinsic scatter (∆ int ) of the mean RAR at high g bar for different dark matter models (see Table 1). In all cases, we find a remarkably tight distribution in log g tot -log g bar space, with a slight increase in ∆ int with increasing SIDM cross-section. We find that only the CDM case can be well described by Equation (2) proposed by McGaugh et al. (2016), with an acceleration scale of g † = (1.42 ± 0.06) × 10 −9 m s −2 . This is much higher than the characteristic acceleration scale of g † ≈ 1.2×10 −10 m s −2 observed at galaxy scales (McGaugh et al. 2016).
We have compared the halo RARs from the BAHAMAS-CDM and -SIDM runs to the cluster RAR inferred from CLASH observations (Section 4; see Tables 2 and 3). Our comparison shows that the RAR in the CDM model provides an excellent match to the CLASH RAR (Tian et al. 2020). This comparison includes the high-g bar regime probed by the BCGs. By contrast, models with a larger SIDM cross-section (hence with a greater flattening in g tot at high g bar ) yield increasingly poorer matches to the CLASH RAR. Excluding the BCG regions, we obtain a weaker but still competitive constraint on the SIDM cross-section. Using the RAR data outside the central r < 100 kpc region, we find that an SIDM model with σ/m = 0.3 cm 2 g −1 is disfavored at the 3.8σ level with respect to the CDM model. However, it should be noted that the choice of the NFW model used for fitting the CLASH lensing data (Tian et al. 2020) implicitly assumes collisionless CDM, which could lead to a biased inference of the cluster RAR in an SIDM cosmology. As a result, this NFW assumption is likely to underestimate the goodness of fit for the SIDM models. In future work, it will be important to use more flexible mass models with an additional parameter to describe the central density slope that depends on the SIDM cross-section (Eckert et al. 2022b).
In this study, we have demonstrated the power and potential of the RAR for testing the collisionless nature of dark matter. Thus far, the cluster RAR has been determined using gravitational lensing only for the CLASH sample. To place stringent and robust constraints on the SIDM cross-section, it is necessary to increase the sample of clusters for which lensing and X-ray observations are available over a broad radial range down to r ∼ 100 kpc (Section 3.2). For nearby clusters at z < 0.1, such mass measurements can be obtained from wide-field weak-lensing observations (e.g., Okabe et al. 2014). For clusters at higher redshifts, combined strong and weak lensing is required to distinguish SIDM and CDM using the RAR. The ongoing CHEX-MATE project will provide such ideal multiwavelength data sets of high quality, for a minimally biased, signal-to-noise-limited sample of 118 Planck galaxy clusters at 0.05 < z < 0.6 detected through the SZ effect (CHEX-MATE Collaboration et al. 2021). Extending this work to the CHEX-MATE sample will thus be a substantial step toward understanding the collisionless nature of dark matter.

ACKNOWLEDGMENTS
We thank the anonymous referee for valuable comments that improved the clarity of the paper. We acknowledge fruitful discussions with Yong Tian, Teppei Okumura, and Stefano Ettori. This work is supported by the Ministry of Science and Technology of Taiwan

A. CORRELATION BETWEEN THE TOTAL AND BARYONIC MATTER DISTRIBUTIONS
The key features used to distinguish different dark matter models mainly come from the inner region of cluster halos (r < ∼ 200 kpc), as shown in Figure 4. Here, we further investigate whether the correlation between the acceleration ratio g tot /g bar (or the total-to-baryonic mass ratio, M tot /M bar ) and the clustercentric distance r can provide sufficient information to distinguish between different dark matter models.
The top panel of Figure 8 shows the mean and standard deviation profiles of g tot (r)/g bar (r) as a function of cluster radius r, derived for each of the dark matter runs of the BAHAMAS simulation at z = 0.375 using the high-mass subsample with E(z)M 200 > 5 × 10 14 M ⊙ . The values of the logarithmic scatter of g tot (r)/g bar (r) evaluated in the inner region of cluster halos (r < 400 kpc or r < 200 kpc) are listed in Table 4. The radial range of r < 400 kpc corresponds approximately to the high-acceleration region with g bar > 10 −10.6 m s −2 used to characterize the intrinsic scatter around the mean RAR, or the g bar -g tot relation (Table 1).
Overall, the magnitude of the scatter in the g tot /g bar -r relation is about a factor of ∼ 1.6 larger than that of the RAR (see Table 1). As a result, the radial g tot (r)/g bar (r) profiles of cluster halos derived from the different dark matter runs significantly overlap with each other, even in their innermost region (Figure 8). Therefore, we conclude that compared to the RAR or MDAR, the g tot /g bar -r relation (i.e., the M tot (r)/M bar (r) profile) has a much weaker sensitivity to the SIDM cross-section.
To understand the tightness of RAR/MDAR, we split the simulation data in log(g tot /g bar )-log g bar space into six dif- Figure 8. Ratio of the total acceleration gtot to the baryonic acceleration g bar as a function of cluster radius r (upper panel) and as a function of baryonic acceleration g bar (lower panel) for massive cluster-scale halos with E(z)M200 > 5 × 10 14 M⊙ at z = 0.375. The results are shown separately for five different dark matter simulation runs. For each dark matter model, the solid line represents the mean relation and the shaded area shows the standard deviation around the mean. ferent radial bins ranging from 15 kpc to 1200 kpc and show them separately in the left panels of Figure 9. In each radial bin, the MDAR sequence spans a wide range in g bar , indicating that even dark matter particles at large clustercentric distances can contribute to the high acceleration region. The distribution of g tot /g bar in a given radial bin is obtained by projecting each set of color-coded data onto the y-axis in Figure 9, which results in a large scatter around the mean value ⟨g tot /g bar ⟩. Therefore, the tight correlation in the RAR/MDAR appears to be a unique signature of the empirical coupling between the total and baryonic components, which increases the ease of distinguishing between CDM and SIDM models.  Figure 9. The gtot/g bar -g bar diagram derived for a subsample of massive cluster-scale halos with E(z)M200 > 5 × 10 14 M⊙ at z = 0.375. The results are shown separately for five different dark matter runs of the BAHAMAS simulation. In the left panels, the data are divided into six different radial bins ranging from r = 15 kpc to 1200 kpc. All points in each panel are color-coded according to their radial bin. In the right panels, the corresponding histogram distribution is shown for each dark matter run. In each panel, the red solid line represents the mean relation and the red shaded area show the standard deviation around the mean relation.
The physical origin and nature of this tight empirical coupling between the total and baryonic accelerations are still under debate. In the future, it will be useful to search for different combinations of the total and baryonic mass properties of galaxy clusters to identify other low-scatter probes of the SIDM cross-section.

B. ACCOUNTING FOR THE MEASUREMENT UNCERTAINTY OF THE BARYONIC ACCELERATION
In this appendix, we explicitly account for the measurement uncertainty of g bar in our analysis of the CLASH RAR data set. To this end, we utilize the MDAR, namely g tot /g bar as a function of g bar , and follow the analysis procedure described in Section 4.2. The χ 2 function in this analysis is defined by [log 10 (g tot,i /g bar,i ) − log 10 ( g tot,i / g bar,i )] 2 ∆ 2 i + ∆ 2 int , (B1) where i runs over all data points from the CLASH data set, ∆ 2 i = ∆ 2 (log 10 g tot,i ) + ∆ 2 (log 10 g bar,i ) represents the total uncertainty including the measurement errors in both log 10 g bar and log 10 g tot , and ∆ int is the lognormal intrinsic scatter of g tot at fixed g bar . Figure 10 shows the distribution of the CLASH-MDAR data and the BAHAMAS predictions of the mean g tot /g barg bar relation derived from five different dark matter runs. The resulting χ 2 values and corresponding PTE values evaluated for the respective dark matter runs are listed in Table 5.  a Observed χ 2 value between the CLASH data and each dark matter model.
b Probability to exceed the observed χ 2 value assuming the standard χ 2 probability distribution function.
When the central BCG constraint is included, only the CDM model gives a satisfactory match to the CLASH data set, which is consistent with the result of the RAR-based analysis given in Section 4.2. However, we note that the SIDM0.1 model now has a PTE value of 0.047, which is very close to but slightly below the significance threshold of α = 0.05. When excluding the central BCG region of r < 100 kpc, we find that both CDM and SIDM0.1 provide satisfactory matches to the CLASH data. All these results are consistent with the findings based on the RAR analysis presented in Section 4.2, suggesting that the inclusion of the measurement uncertainty of g bar has only a minor impact on the statistical inference and it does not change our conclusions of the accepted dark matter models.