The Radial Acceleration Relation in CLASH Galaxy Clusters

Combining the strong lensing, weak lensing shear and magnification effects, Umetsu et al. (2016) reconstructed the surface mass density profile of each individual cluster over a wide range of the cluster-centric radius. Umetsu et al. (2016) found that the ensemble-averaged total mass distribution of the CLASH sample is well described by a family of cuspy, outward-steepening density profiles, namely, the Navarro–Frenk–White (NFW, hereafater; Navarro et al., 1997), Einasto, and DARKexp (Hjorth & Williams, 2010) models. Of these, the NFW model best describes the CLASH lensing data (e.g., Umetsu et al., 2011b; Umetsu & Diemer, 2017). On the other hand, the single power-law, cored isothermal, and Burkert profiles were statistically disfavored by the observed CLASH lensing profile having a pronounced radial curvature (Umetsu et al., 2016).

Here we use the CLASH lensing constraints on the total mass profile $M_{\mathrm{tot}}(<r)$ of each individual CLASH cluster assuming a spherical NFW profile. The total mass $M_{\mathrm{tot}}(<r)$ of an NFW halo as a function of spherical radius $r$ is written as

where $r_{\mathrm{s}}$ and $\rho_{\mathrm{s}}$ represent the characteristic scale radius and density of the NFW profile, respectively, and $\rho_{\mathrm{s}}$ is given by

with $c_{200}\equiv r_{200}/r_{\mathrm{s}}$ the NFW concentration parameter.

For each cluster, Umetsu et al. (2016) extracted the posterior probability distributions of $(M_{200},c_{200})$ from the observed surface mass density profile assuming a spherical NFW halo, by accounting for all relevant sources of uncertainty (e.g., Umetsu et al., 2016, 2020; Miyatake et al., 2019): (i) measurement errors, (ii) cosmic noise due to projected large-scale structure uncorrelated with the cluster, (iii) statistical fluctuations of the projected cluster lensing signal due to halo triaxiality and correlated substructures. According to cosmological hydrodynamical simulations of Meneghetti et al. (2014), the CLASH sample selection is expected to be largely free from orientation bias. In fact, three-dimensional full-triaxial analyses of CLASH lensing observations found no statistical evidence for orientation bias in the CLASH sample (see Sereno et al., 2018; Chiu et al., 2018). Therefore, assuming spherical NFW halos is not expected to cause any significant bias in lensing mass estimates of the CLASH sample. The mass and concentration parameters $(M_{200},c_{200})$ of each individual CLASH cluster are summarized in Table 2 of Umetsu et al. (2016).

In our analysis, we use these posterior distributions of the NFW parameters to obtain well-characterized inference of $M_{\mathrm{tot}}(<r|M_{200},c_{200})$ for each individual cluster. We compute the total mass profile $M_{\mathrm{tot}}(<r|M_{200},c_{200})$ and its uncertainty of each cluster in the radial range at $r\geq r_{\mathrm{min}}\simeq 14$ kpc.

The X-ray emitting hot gas dominates the baryonic mass in galaxy clusters. In high-mass clusters, more than $80\%$ of the intra-cluster baryons are in the X-ray emitting hot phase (e.g., Umetsu et al., 2009; Donahue et al., 2014; Okabe et al., 2014; Chiu et al., 2018). Donahue et al. (2014) derived enclosed gas mass profiles $M_{\mathrm{gas}}(<r)$ for all CLASH clusters, finding that the $M_{\mathrm{gas}}$ profiles measured from the Chandra and XMM observations are in excellent agreement where the data overlap. Since XMM data are not available for all the clusters, we only use Chandra $M_{\mathrm{gas}}$ measurements of Donahue et al. (2014) as primary constraints on the baryonic mass content in the CLASH sample. For each cluster, we have measured $M_{\mathrm{gas}}(<r)$ values in several radial bins, where the radial range is different for each cluster, depending on the redshift and the data quality (see Figure 1).

We then account for the stellar contribution to the baryonic mass using the results of Chiu et al. (2018), who established the stellar-to-gas mass relation $f_{\mathrm{c}}(r)=M_{\mathrm{star}}/(M_{\mathrm{gas}}+M_{\mathrm{star}})$ (see their Figure 11; $f_{\mathrm{c}}$ referred to as the cold collapsed baryonic fraction), for a sample of 91 Sunyaev–Zel’dovich effect (SZE) clusters with $M_{500}>2.5\times 10^{14}M_{\odot}$ selected from the South Pole Telescope (SPT) survey (Carlstrom et al., 2011; Bleem et al., 2015). The $f_{\mathrm{c}}(r)$ relation is insensitive to the cluster redshift over a broad range out to $z\sim 1.3$ as probed by the SPT sample. In their study, the total masses were estimated from the SZE observable, the gas masses $M_{\mathrm{gas}}$ from Chandra X-ray data, and the stellar masses $M_{\mathrm{star}}$ from combined optical/near-infrared multi-band photometry. The $f_{\mathrm{c}}(r)$ relation of Chiu et al. (2018) includes the stellar mass contributions from the BCG and cluster member galaxies inside the $r_{500}$ overdensity radius.

With the mean $f_{\mathrm{c}}(r)$ relation, we estimate the total baryonic mass as $M_{\mathrm{bar}}(<r)=M_{\mathrm{gas}}(<r)/[1-f_{\mathrm{c}}(r)]$, ignoring the cluster-to-cluster scatter. The level of scatter around the mean $f_{\mathrm{c}}(r)$ relation is about $12\%$.

In the central cluster region, the baryonic mass in clusters is dominated by the stellar mass in the BCG. In this study, we model the stellar mass distribution of each BCG with the Hernquist model (Hernquist, 1990), which gives an analytical approximation to the deprojected form of de Vaucouleurs’ profile (or a Sérsic profile with index $n=4$). Then, the stellar mass $m_{\mathrm{star}}(<r)$ inside the spherical radius $r$ and the stellar gravitational acceleration $g_{\mathrm{star}}(<r)$ are expressed as

where ${\cal M}_{\mathrm{star}}$ is the total stellar mass of the BCG, and $r_{\mathrm{h}}\approx 0.551R_{\mathrm{e}}$ is a characteristic scale length of the Hernquist model, with $R_{\mathrm{e}}$ the half-light or effective radius of the de Vaucoleurs’ brightness profile.

In this study, we adopt as ${\cal M}_{\mathrm{star}}$ (see Table 1) the stellar mass estimates of CLASH BCGs from Cooke et al. (2016, their Table 1), who performed a multiwavelength analsyis on a large sample of BCGs by combining UV, optical, near-infrared, and far-infrared data sets. Since the measurement errors on BCG stellar masses were not provided in Cooke et al. (2016), we assume a fractional uncertainty of $10\%$ on ${\cal M}_{\mathrm{star}}$.

We measure the BCG effective radius $R_{\mathrm{e}}$ from the CLASH HST imaging using the GALFIT package (Peng et al., 2010). We choose to measure $R_{\mathrm{e}}$ of each BCG in the HST band corresponding to the rest-frame wavelength of $1\mu$m. The corresponding HST bands for our sample ($0.187\leq z\leq 0.686$) are all in the WFC3/IR coverage (F110W to F160W), as summarized in Table 1. We fit a single Sérsic profile to the surface brightness distribution of the BCG in the CLASH HST imaging data. The initial guess of the BCG position is based on the rest-frame UV (280 nm) measurements of Donahue et al. (2015). In Table 1, we also list the final source position ($\mathrm{R.A.},\mathrm{Decl.}$) in J2000 coordinates, Sérsic index ($n$), and effective radius ($R_{\mathrm{e}}$) from our GALFIT modeling.

We define the hot gas fraction as $f_{\mathrm{gas}}(r)=M_{\mathrm{gas}}(<r)/M_{\mathrm{tot}}(<r)$, the ratio of the gas mass $M_{\mathrm{gas}}(<r)$ to the total mass $M_{\mathrm{tot}}(<r)$ as a function of the cluster-centric radius $r$. For each cluster, we evaluate $M_{\mathrm{tot}}(<r)$ and $f_{\mathrm{gas}}(r)$ where the Chandra gas mass measurements $M_{\mathrm{gas}}(<r)$ of Donahue et al. (2014) are available.

In the upper panel of Figure 1, we show the hot gas fractions $f_{\mathrm{gas}}(r)$ of all individual clusters in our CLASH sample, along with the mean profile $\langle f_{\mathrm{gas}}(r)\rangle$ of the sample. The mean $\langle f_{\mathrm{gas}}(r)\rangle$ increases with increasing cluster-centric radius $r$, approaching the cosmic baryon fraction $\Omega_{\mathrm{b}}/\Omega_{\mathrm{m}}=(15.7\pm 0.4)\%$ (Planck Collaboration et al., 2016) at $r\lower 2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}700\,\mathrm{kpc}\sim 0.5r_{500}$.

We note that Donahue et al. (2014) derived the hot gas fractions of CLASH clusters by combining their X-ray gas mass measurements with earlier CLASH weak-lensing results from Umetsu et al. (2014) or Merten et al. (2015). Our results improve upon those of Donahue et al. (2014) by using the full-lensing constraints of Umetsu et al. (2016) based on CLASH strong-lensing, weak-lensing shear and magnification data.

We compute for each cluster the baryon fraction $f_{\mathrm{bar}}(r)=M_{\mathrm{bar}}(<r)/M_{\mathrm{tot}}(<r)$ as a function of $r$. The baryonic cluster mass $M_{\mathrm{bar}}(<r)$ consists of the X-ray gas mass $M_{\mathrm{gas}}(<r)$ from Donahue et al. (2014), the stellar mass estimated as $M_{\mathrm{star}}(<r)=M_{\mathrm{gas}}(<r)\times f_{\mathrm{c}}(r)/[1-f_{\mathrm{c}}(r)]$ (Section 2.2), and the stellar mass of the BCG in the innermost cluster region (Section 2.3). As in Section 3.1, we evaluate for each cluster $M_{\mathrm{tot}}(<r)$ and $f_{\mathrm{bar}}(r)$ where the Chandra $M_{\mathrm{gas}}(<r)$ values of Donahue et al. (2014) are available.

For each cluster, we also include a single constraint on the baryon fraction $f_{\mathrm{bar}}(<r)$ in the central BCG region. The stellar mass distribution of each BCG is modeled by Equation (5), as described in Section 2.3. For 13 clusters in our sample, we have Chandra $M_{\mathrm{gas}}$ measurements (Donahue et al., 2014) lying in the central BCG region at $r\ <\langle R_{\mathrm{e}}\rangle\sim 30$ kpc (see Table 1). For these clusters, we calculate the BCG stellar mass $m_{\mathrm{star}}(<r)$ at this innermost radius of the Chandra $M_{\mathrm{gas}}$ measurements (Table 1). For the other clusters, we calculate $m_{\mathrm{star}}(<r)$ at $r_{\mathrm{min}}\simeq 14$ kpc and ignore the gas mass contribution to $M_{\mathrm{bar}}(<r_{\mathrm{min}})$. We note that, typically, the hot gas contribution in the innermost region $r<R_{\mathrm{e}}$ is subdominant compared to the BCG stellar mass (e.g., Sartoris et al., 2020).

In the lower panel of Figure 1, we show the baryon fraction profiles $f_{\mathrm{bar}}(r)$ outside the BCG region for all individual clusters in our sample. The mean $\langle f_{\mathrm{bar}}(r)\rangle$ profile of the sample is nearly constant $\sim 1/8$ (Donahue et al., 2014) with $r$. The cluster-to-cluster scatter around the mean $f_{\mathrm{bar}}(r)$ profile is 0.041 in terms of the standard deviation, and the total baryon fractions in some clusters reach the cosmic mean value. We do not find any clear radial trend in the baryon fraction profiles $f_{\mathrm{bar}}(r)$ for the CLASH sample, as in the case of spiral galaxies (McGaugh, 2004; Famaey & McGaugh, 2012).

At $r=800\,\mathrm{kpc}\sim 0.6r_{500}$ (the maximum radius of our ensemble measurements), the mean baryon fraction of the CLASH sample is $\langle f_{\mathrm{bar}}\rangle=(14.2\pm 1.4)\%$, which corresponds to a depletion factor of ${\cal D}\equiv 1-f_{\mathrm{bar}}/(\Omega_{\mathrm{b}}/\Omega_{\mathrm{m}})=(10\pm 9)\%$ with respect to the cosmic mean value, $\Omega_{\mathrm{b}}/\Omega_{\mathrm{m}}=(15.7\pm 0.4)\%$. This level of depletion is not statistically significant, and it is in agreement with ${\cal D}=(18\pm 2)\%$ at $r=r_{500}$ found from an independent constraint on the SPT sample by Chiu et al. (2018), as well as with the results from numerical simulations, ${\cal D}=10\%$–$20\%$ (e.g., McCarthy et al., 2011; Barnes et al., 2017). Moreover, since the hot gas is more extended than DM, the gas fraction $f_{\mathrm{gas}}(r)$ increases with cluster-centric radius $r$, so that the depletion factor of the CLASH sample at $r_{500}$ is expected to be much less significant.

Here we quantify and characterize the relationship between the total acceleration $g_{\mathrm{tot}}(r)$ and the baryonic acceleration $g_{\mathrm{bar}}(r)$ for the CLASH sample (Figure 1). To this end, for each cluster, we extract data points where possible at $r=100$ kpc (18 clusters), $200$ kpc (20 clusters), $400$  kpc (15 clusters), and $600$  kpc (11 clusters), by using linear interpolation. These data points are sufficiently well separated from each other. Hence, for simplicity, we ignore covariances between different radial bins in our fitting procedure. Altogether, we have a total of 84 data points for our sample of 20 CLASH clusters, including 20 data points in the central BCG region.

Since the mass ratio $M_{\mathrm{tot}}(<r)/M_{\mathrm{bar}}(<r)$ is equivalent to the acceleration ratio $g_{\mathrm{tot}}(r)/g_{\mathrm{bar}}(r)$, the MDAR can be expressed as a relation between $g_{\mathrm{tot}}(r)/g_{\mathrm{bar}}(r)$ and $g_{\mathrm{bar}}$. In spiral galaxies, there exists a tight empirical relationship between the observed total acceleration $g_{\mathrm{tot}}(r)$ and the baryonic acceleration $g_{\mathrm{bar}}(r)$, namely the RAR (McGaugh et al., 2016). In fact, the MDAR and RAR are mathematically equivalent However, we point out that in the RAR, values of each of the two axes come from independent measurements.

We present our results in the form of the RAR. We model the CLASH data distribution in log-acceleration space by performing a linear regression on the relation $y=m\,x+b$ with $y=\ln(g_{\mathrm{tot}}/g_{0})$ and $x=\ln(g_{\mathrm{bar}}/g_{0})$ with a normalization scale of $g_{0}=1\,\mathrm{m}\,\mathrm{s}^{-2}$.¹¹1We note that the resulting RAR is presented in decimal logarithmic units ($\log_{10}$) in Figures 1 and 3, whereas our regression analysis uses the natural logarithm ($\ln$) of acceleration. Here we account for the uncertainties in the determinations of lensing mass $M_{\mathrm{tot}}(<r)$, gas mass $M_{\mathrm{gas}}(<r)$, and stellar mass $m_{\mathrm{star}}(<r)$ of the BCG (see Section 2). The uncertainties in $x$ and $y$ are expressed as $\sigma_{x}=\sigma(M_{\mathrm{tot}})/M_{\mathrm{tot}}$ and $\sigma_{y}=\sigma(M_{\mathrm{bar}})/M_{\mathrm{bar}}$, where $\sigma(M_{\mathrm{tot}})$ and $\sigma(M_{\mathrm{bar}})$ are the total uncertainties for the lensing and baryonic mass estimates $M_{\mathrm{tot}}$ and $M_{\mathrm{bar}}$, respectively.

The log-likelihood function is written as

where $i$ runs over all clusters and data points, and $\sigma_{i}$ includes the observational uncertainties $(\sigma_{x_{i}},\sigma_{y_{i}})$ and lognormal intrinsic scatter $\sigma_{\mathrm{int}}$ (e.g., see Umetsu et al., 2016; Okabe & Smith, 2016),

The $\sigma_{\mathrm{int}}$ parameter accounts for the intrinsic scatter around the mean RAR due to unaccounted astrophysics associated with the RAR.

We perform a Markov chain Monte Carlo (MCMC) analysis to constrain the regression parameters using the emcee python package (Foreman-Mackey et al., 2013, 2019) based on an affine-invariant sampler (Goodman & Weare, 2010). We use non-informative uniform priors on $b$ and $m$ of $b\in[-100,100]$ and $m\in[-100,100]$. For the intrinsic scatter, we assume a prior that is uniformed in $\ln{\sigma_{\mathrm{int}}}$ in the range $\ln{\sigma_{\mathrm{int}}}\in[-5,1]$. Using the MCMC technique, we sample the posterior probability distributions of the regression parameters $(b,m,\sigma_{\mathrm{int}})$ over the full parameter space allowed by the priors.

From the regression analysis, we find a tight RAR for the CLASH sample in the BCG–cluster regime. Figure 2 summarizes our results. In the left panel, we show the distribution of CLASH clusters in $\log_{10}{g_{\mathrm{bar}}}$–$\log_{10}{g_{\mathrm{tot}}}$ space along with the best-fit relation (black solid line). The spread of the best-fit residuals is about 0.11 dex, as shown in the inset plot of Figure 2. The resulting constraints on the regression parameters are $m=0.51^{+0.04}_{-0.05}$, $b=-9.80^{+1.07}_{-1.08}$, and $\sigma_{\mathrm{int}}=14.7^{+2.9}_{-2.8}\%$ in terms of the MCMC-sampled posterior mean and standard deviation. The RAR for the CLASH sample is summarized as

Figure 2 also compares our results with the RAR in sprial galaxies from McGaugh et al. (2016, dashed line). The slope we obtained $m=0.51^{+0.04}_{-0.05}$ is consistent with the low acceleration limit $1/2$ of the RAR from McGaugh et al. (2016), $g_{\mathrm{tot}}=\sqrt{g_{\dagger}\,g_{\mathrm{bar}}}$ with $g_{\dagger}\simeq 1.2\times 10^{-10}$ m s^-2. On the other hand, the intercept $b$ at fixed $g_{\mathrm{bar}}$ is found to be significantly higher than that of McGaugh et al. (2016).

Here we repeat our regression analysis by fixing the slope to $m=1/2$ and rewriting the scaling relation as $g_{\mathrm{tot}}(r)=\sqrt{g_{\ddagger}\,g_{\mathrm{bar}}(r)}$ with $g_{\ddagger}$ a constant acceleration that corresponds to a certain characteristic acceleration scale. Then, the regression parameters are constrained as $g_{\ddagger}=(2.02\pm 0.11)\times 10^{-9}$ m s^-2. and $\sigma_{\mathrm{int}}=14.5^{+2.9}_{-2.8}\%$.

Using a semi-analytical model with abundance matching, Navarro et al. (2017) provided a possible explanation of the RAR in spiral galaxies within the $\Lambda$CDM framework. Here we test if the observed RAR for the CLASH sample can be explained by a semi-analytical description of cluster-scale halos in the standard $\Lambda$CDM model. Unlike the RAR in spiral galaxies based on stellar kinematics, the total acceleration $g_{\mathrm{tot}}$ in the CLASH sample has been derived assuming the NFW density profile. Our data points and semi-analytical model are thus not entirely independent in terms of the profile shape of DM.

To this end, we employ the semi-analytical model of Olamaie et al. (2012), which describes the distributions of DM and hot gas with an ideal gas equation of state in a spherical cluster halo. First, this model assumes that DM follows the NFW profile (Umetsu et al., 2011a, 2016; Niikura et al., 2015; Okabe & Smith, 2016) and the gas pressure is described by a generalized NFW profile (Nagai et al., 2007). Following Olamaie et al. (2012), we fix the values of the gas concentration and slope parameters of the generalized NFW profile to those found by Arnaud et al. (2010). Next, the system is assumed to be in hydrostatic equilibrium. We then assume a gas mass fraction of $f_{\mathrm{gas}}(r_{500})=13\%$ (e.g., Planck Collaboration et al., 2013; Donahue et al., 2014) at $r=r_{500}$ to fix the normalization of $\rho_{\mathrm{gas}}(r)$. Finally, in these calculations, we assume that the gas density $\rho_{\mathrm{gas}}$ is much smaller than the DM density $\rho_{\mathrm{DM}}$, $\rho_{\mathrm{tot}}(r)=\rho_{\mathrm{DM}}(r)+\rho_{\mathrm{gas}}(r)\approx\rho_{\mathrm{DM}}(r)$. With these assumptions, $\rho_{\mathrm{gas}}(r)$ can be fully specified by two parameters that describe the NFW density profile. We refer to Olamaie et al. (2012) for full details of the model.

Following the procedure outlined above, we can describe the average properties of our cluster sample, which includes 16 X-ray-selected and 4 high-magnification-selected CLASH clusters of Umetsu et al. (2016). Here we adopt $M_{200}=1.55\times 10^{15}M_{\odot}$ and $c_{200}=3.28$ to describe the DM distribution $\rho_{\mathrm{DM}}(r)$ of our sample with a median redshift of $z=0.377$. Given the NFW parameters, we compute the gas density profile $\rho_{\mathrm{gas}}(r)$ for the CLASH sample. We use $f_{\mathrm{c}}(r)$ of Chiu et al. (2018) (Section 2.2) to account for the stellar mass contribution to the baryonic mass. With these average profiles $\rho_{\mathrm{DM}}(r)$, $\rho_{\mathrm{gas}}(r)$, and $f_{\mathrm{c}}(r)$, we can predict the total and baryonic gravitational acceleration profiles, $g_{\mathrm{tot}}(r)$ and $g_{\mathrm{bar}}(r)$, for the CLASH sample. Here we compute the acceleration profiles at $r=100,200,400$, and $600$ kpc, as done in our CLASH analysis.

We also model the intrinsic scatter around the average profiles $g_{\mathrm{tot}}(r)$ and $g_{\mathrm{bar}}(r)$ due to cluster-to-cluster variations in the DM and baryonic distributions. For the DM distribution, we assign intrinsic scatter in $c_{200}$ with a lognormal intrinsic dispersion of $30\%$ (e.g., Bhattacharya et al., 2013). For the baryonic distribution, we assign intrinsic scatter in $f_{\mathrm{c}}$ with a Gaussian dispersion of $12\%$ (Chiu et al., 2018). We employ Monte-Carlo simulations to evaluate the intrinsic dispersions around the average acceleration profiles $g_{\mathrm{tot}}(r)$ and $g_{\mathrm{bar}}(r)$ at each cluster-centric radius. The result is shown in Figure 3.

The 20 galaxy clusters in our sample are all high-mass systems selected for the CLASH survey (Postman et al., 2012). It is thus reasonable to adopt a set of average properties in the semi-analytical model to study the new RAR. As shown in Figure 3, all the model points (color-coded squares) match the distribution of data points (color open circles) and the mean relation (black solid line) well.

The inferred level of $1\sigma$ intrinsic scatter estimated by Monte-Carlo simulations appears to be larger than the data distribution. This is similar to the findings in spiral galaxies (see, e.g., Di Cintio & Lelli, 2016; Desmond, 2017; Lelli et al., 2017). However, it should be noted that the CLASH sample is dominated by relaxed systems (Meneghetti et al., 2014) because of the CLASH selection based on X-ray morphology regularity (Postman et al., 2012). As a result, the CLASH sample is predicted to have a much smaller level of intrinsic scatter in the $c_{200}$–$M_{200}$ relation ($16\%$; Meneghetti et al., 2014; Umetsu et al., 2016).

The CLASH RAR (Equation (8)) expresses $g_{\mathrm{tot}}(r)$ as a function of $g_{\mathrm{bar}}(r)$. Thus, one can obtain the baryon fraction $f_{\mathrm{bar}}(r)$ if $g_{\mathrm{bar}}(r)$ is known, because $f_{\mathrm{bar}}(r)=M_{\mathrm{bar}}(<r)/M_{\mathrm{tot}}(<r)=g_{\mathrm{bar}}(r)/g_{\mathrm{tot}}(r)$. If we approximate the CLASH RAR as $g_{\mathrm{tot}}\approx\sqrt{g_{\ddagger}\,g_{\mathrm{bar}}}$, then the baryon fraction has the simple form,

Since the 20 CLASH clusters in our sample are of similar size and acceleration profile, we can infer an average relation between the baryonic acceleration $g_{\mathrm{bar}}$ and the cluster-centric radius $r$. Making use of the X-ray gas mass measurements from Donahue et al. (2014) and the stellar mass correction from Chiu et al. (2018), we obtain such an average relation, i.e., the average baryonic acceleration as a function of $r$, $\langle g_{\mathrm{bar}}(r)\rangle$. Then, the best-fit CLASH RAR (Equation (8)) together with $\langle g_{\mathrm{bar}}(r)\rangle$ gives an empirical relation for the baryon fraction as a function of $r$.

In Figure 4, we compare the $f_{\mathrm{bar}}(r)$ profile inferred from the best-fit CLASH RAR (Equation (8)) with the observed distribution of CLASH baryon fractions shown in Figure 1. The gray circles are the inferred baryon fractions at different cluster-centric radii, $r$. The thick black line shows the observed mean $\langle f_{\mathrm{bar}}(r)\rangle$ profile averaged over the CLASH sample (see Figure 1). The $f_{\mathrm{bar}}(r)$ profile inferred from the best-fit RAR and the mean $\langle f_{\mathrm{bar}}(r)\rangle$ profile for the CLASH sample are consistent within $r=400$ kpc, but deviate from each other beyond $400$ kpc. This is because the best-fit RAR model predicts a radially decreasing $f_{\mathrm{bar}}(r)$ profile, whereas the observed CLASH $\langle f_{\mathrm{bar}}(r)\rangle$ profile is nearly constant ($\sim 1/8$; see Figure 1) in the intracluster regime.

We also compare these baryon fraction profiles with the corresponding values predicted by our semi-analytical model, which are denoted by the orange, green, blue and purple squares for $r=100,200,400$, and $600$ kpc, respectively. The values of all model points agree well with the best-fit RAR and the observed mean values for the CLASH sample. Our semi-analytical model predicts a nearly constant $f_{\mathrm{bar}}(r)$ profile, in good agreement with the mean relation of the CLASH sample (black solid line).

To assess how the data deviate from the best-fit CLASH RAR, we show in Figure 5 the distribution of best-fit residuals as a function of cluster radius $r$. We note again that both the semi-analytical model and the CLASH RAR assume the NFW density profile for the total matter distribution. The residuals are defined as the difference of $\log_{10}(g_{\mathrm{tot}})$ between the observational data and the best-fit RAR (Equation (8)). As shown in the figure, the residuals at small $r$ ($\approx 14$, 100, and 200 kpc) distribute symmetrically across the best-fit RAR (the zero residual line). However, at larger cluster radii $r\lower 2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}400$ kpc, the residuals deviate systemically toward negative values. The semi-analytical model (color-coded squares) exhibits a similar trend as seen for the observational data.

As discussed earlier, Equation (1) in the MOND framework is consistent with the observed RAR on galaxy scales with $g_{\mathrm{bar}}\equiv g_{\mathrm{M}}$, without introducing DM. Here we test this hypothesis with our CLASH RAR results and infer the level of residual missing mass on the BCG-cluster scale within the MOND framework.

Since the total acceleration $g_{\mathrm{tot}}$ in the CLASH RAR implies $g_{\mathrm{tot}}\approx\sqrt{g_{\mathrm{bar}}\,g_{\ddagger}}$, we can relate $g_{\mathrm{bar}}$ and $g_{\mathrm{M}}$ by

The mass ratio between $M_{\mathrm{M}}$ and $M_{\mathrm{bar}}$ is $M_{\mathrm{M}}/M_{\mathrm{bar}}=g_{\mathrm{M}}/g_{\mathrm{bar}}$. Our formulation implies that the level of residual missing mass in MOND depends on the baryonic acceleration $g_{\mathrm{bar}}$. In our CLASH data, the largest baryonic acceleration in the BCG regime is $2.1\times 10^{-10}$ m s^-2 and the smallest one in the intracluster regime is $1.3\times 10^{-11}$ m s^-2. The corresponding mass ratio $M_{\mathrm{M}}/M_{\mathrm{bar}}$ thus ranges from $2.7$ to $7.3$, increasing with decreasing baryonic acceleration $g_{\mathrm{bar}}$ (or increasing cluster-centric radius).

Hence, the CLASH RAR confirms the existence of residual missing mass in MOND on the BCG-cluster scale, as found in the literature (Sanders, 1999, 2003; Pointecouteau & Silk, 2005; Famaey & McGaugh, 2012). Furthermore, it reveals a more substantial level of discrepancy or residual missing mass $M_{\mathrm{M}}/M_{\mathrm{bar}}$ in the low $g_{\mathrm{bar}}$ regime, which corresponding to the cluster outskirts. The typical level of residual missing mass found in the literature is $\langle M_{\mathrm{M}}/M_{\mathrm{bar}}\rangle\sim 2$. However, it should be note that a fair comparison of the CLASH RAR with MOND will require a relativistic extension of MOND to properly interpret the gravitational lensing data for the CLASH sample. Some possibilities in this approach for cluster lensing have been discussed in the literature (Bruneton et al., 2009; Zhao & Famaey, 2012; Famaey & McGaugh, 2012).

We recall that the RAR in galaxies has a characteristic acceleration scale of $g_{\dagger}\simeq 1.2\times 10^{-10}$ m s^-2, above which the acceleration of the system asymptotically tends to Newtonian dynamics without DM (i.e., $g_{\mathrm{tot}}=g_{\mathrm{bar}}$), while below which the acceleration asymptotically tends to $g_{\mathrm{tot}}=\sqrt{g_{\dagger}\,g_{\mathrm{bar}}}$. Our CLASH results indicate that $g_{\mathrm{tot}}\approx\sqrt{g_{\ddagger}\,g_{\mathrm{bar}}}$ (see Equation (9)) on BCG–cluster scales, with a characteristic acceleration scale of $g_{\ddagger}\gg g_{\dagger}$. However, in our CLASH sample, it is not clear whether or not the acceleration of the system will approach to $g_{\mathrm{tot}}=g_{\mathrm{bar}}$ in the high acceleration limit, $g_{\mathrm{bar}}\gg g_{\ddagger}$.

If the CLASH RAR holds in general for other galaxy clusters, then there exists a kinematic relation or law in galaxy clusters, in analogy to Kepler’s law of planetary motion that comes from Newtonian dynamics. Let us take the example of the baryonic Tully–Fisher relation (BTFR; e.g., McGaugh, 2011, 2012; Famaey & McGaugh, 2012), which comes from the low acceleration end of the RAR of spiral galaxies. In rotationally support systems such as spiral galaxies, the centripetal acceleration is provided by the total gravitational acceleration, i.e., $v^{2}/r=g_{\mathrm{tot}}$ where $v$ is the circular speed. In the small acceleration regime of $g_{\mathrm{bar}}\ll g_{\dagger}$, the RAR in galaxies gives $g_{\mathrm{tot}}\approx\sqrt{g_{\dagger}\,g_{\mathrm{bar}}}$. By expressing the baryonic acceleration as $g_{\mathrm{bar}}=GM_{\mathrm{bar}}/r^{2}$ where $M_{\mathrm{bar}}$ is the total baryonic mass inside $r$, we obtain the BTFR, $v^{4}=Gg_{{\dagger}}M_{\mathrm{bar}}$. For systems supported by random motions, such as elliptical galaxies and galaxy clusters, $g_{\mathrm{tot}}\propto\sigma^{2}/r$ with $\sigma_{v}$ the velocity dispersion. If a system that follows the RAR is in the low acceleration regime, $\sigma_{v}^{4}\propto g_{{\ddagger}}M_{\mathrm{bar}}$ is anticipated. We refer to this kinematic law as the baryonic Faber–Jackson relation (BFJR).

In the literature, the BFJR has not been confirmed in galaxy clusters. However, scaling relations between total cluster mass and X-ray mass proxies (e.g., X-ray gas temperature; Sanders, 1994; Ettori et al., 2004; Angus et al., 2008; Famaey & McGaugh, 2012), and that between X-ray luminosity and galaxy velocity dispersion (Xue & Wu, 2000; Sanders, 2010; Zhang et al., 2011; Nastasi et al., 2014), have been firmly established based on multiwavelength observations and numerical simulations. If the total baryonic mass in galaxy clusters is tightly coupled with thermodynamic properties of the hot gas in hydrostatic equilibrium, we may expect a correlation between the total baryonic mass and galaxy velocity dispersion.