THE DENSITY PROFILES OF MASSIVE, RELAXED GALAXY CLUSTERS. II. SEPARATING LUMINOUS AND DARK MATTER IN CLUSTER CORES

The internal structure of dark matter (DM) halos is a key prediction of the cold dark matter (CDM) paradigm. Numerical simulations following the detailed structure of collisionless CDM halos (e.g., Navarro et al. 1996; Ghigna et al. 2000; Diemand et al. 2005; Graham et al. 2006; Gao et al. 2012b) generically produce a central density cusp with ρ_DM ∼ r⁻¹, characteristic of the Navarro–Frenk–White (NFW; Navarro et al. 1996) form, probably becoming slightly shallower on very small scales (e.g., Navarro et al. 2010). However, simulations are only beginning to make predictions for DM halos that include baryons, which could profoundly reshape their host halos. The structure of real DM halos thus contains important information about galaxy formation, but there is currently no theoretical consensus on the magnitude or even sign of these baryonic effects, particularly over a wide range in mass. Additionally, the microphysics of the unknown DM particle could become important in the densest regions, and the inner structure of DM halos may therefore provide valuable indirect clues to its nature (e.g., Spergel & Steinhardt 2000; Abazajian et al. 2001; Kaplinghat 2005; Peter et al. 2010).

Given the current uncertainty, observations are clearly in a good position to guide theoretical efforts. However, measurements of DM mass profiles are extremely challenging and are usually limited by confusion with baryons, the small dynamic range of the observations, and degeneracies that are inherent to individual mass probes (e.g., velocity anisotropy). Clusters of galaxies are promising locations to make progress. Accurate mass measures are available via many independent observational probes, especially gravitational lensing and X-ray emission (see Allen et al. 2011; Kneib & Natarajan 2011 for reviews, and references in Paper I). As we have shown (Newman et al. 2009, 2011), combining stellar kinematics with strong and weak gravitational lensing yields constraints over three decades in radius. This is comparable to the best simulations and is thus suitable for detailed comparison of the DM profile shape if the baryonic mass can be constrained. On small scales in relaxed clusters, the latter is dominated by stars in the central brightest cluster galaxy (BCG).

Sand et al. (2002, 2004) demonstrated the utility of combining resolved stellar kinematics with strong lensing to constrain two-component mass models, i.e., the BCG stars and DM halo separately (see Miralda-Escude 1995; Natarajan & Kneib 1996). Sand et al. (2004) studied six clusters and inferred a mean 〈β〉 = 0.52 ± 0.05, where ρ_DM∝r^−β, significantly shallower than an NFW cusp having β = 1. They further noted a possible variation in β from cluster to cluster. Sand et al. (2008) improved on this analysis for two clusters (MS2137 and A383) by relaxing the assumption of axial symmetry in the lensing analysis, instead conducting a full two-dimensional study. They found that this did not alter their earlier findings, but noted that the inferred DM slope β is sensitive to the adopted scale radius r_s, which could only be constrained by additional mass probes at larger radii. This was implemented by Newman et al. (2009) in A611 through the addition of weak-lensing data. In Newman et al. (2011), we further extended the methodology in A383 by constraining the role of projection effects (i.e., line of sight (l.o.s.) ellipticity; Gavazzi 2005) via a comparison of X-ray and lensing data. We also presented a radially extended velocity dispersion profile measured in a very deep spectroscopic observation. In both A611 and A383, we confirmed a shallow inner DM cusp with β < 0.3 (68% CL) and β = 0.59^+0.30_−0.35, respectively.⁴

In Paper I of the present series, we presented strong and weak lensing and stellar kinematic data for a sample of seven massive (M₂₀₀ = 0.4–2 × 10¹⁵ M_☉), relaxed galaxy clusters at z = 0.19–0.31. This built upon our earlier papers by enlarging the sample of clusters with the highest-quality data: weak lensing measured using deep multi-color imaging, primarily from the Subaru telescope, extended stellar kinematic profiles in the BCG obtained primarily at the Keck telescopes, and multiply imaged sources located in Hubble Space Telescope (HST) imaging (25 strongly lensed sources in total, of which 21 have spectroscopic redshifts). We showed that the total inner density profile is remarkably well described by numerical simulations containing only CDM at radii ≳ 5–10 kpc, despite the significant contribution of stellar mass on these scales.

Here we extend Paper I by dissecting the stellar and DM contributions, using improved versions of the techniques developed in our earlier papers. We first show how the mass content of the BCGs can be constrained using information from the entire sample. We then isolate the DM density profiles and quantify their behavior on small scales. We show that the DM profiles become shallower than NFW models within ≈30 kpc, roughly the typical effective radius of the BCGs. Furthermore, the inner DM density profiles exhibit likely variation from cluster to cluster, and this variation is correlated with the properties of the BCG. Finally, we interpret our results in the context of the recent theoretical literature, focusing on the interactions between baryons and DM in galaxy clusters and the possibility that cores in galaxy clusters are imprints of DM particle physics.

Throughout this paper we adopt a ΛCDM cosmology with Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹. Error bars and upper limits encompass the 68% confidence interval. When pairs of errors are quoted, they refer to the random and systematic components, respectively.

Whereas the total density profiles were studied in Paper I, the goal of this paper is to use our two-component fits to separate the stellar and dark mass contributions in the cluster cores. All aspects of the data and modeling were discussed extensively in Paper I (Section 7). Here, we provide a summary of the features most relevant for this paper. First, the stellar mass in the BCG is modeled based on fits to the surface luminosity in HST imaging. A uniform stellar mass-to-light ratio ϒ_* is assumed within each BCG. As discussed in Paper I (Sections 5.1 and 9.3), this is justified by the mild or null color gradients observed over the relevant radial interval. Non-BCG cluster galaxies—relevant as perturbations in the strong lens model—are included via scaling relations based on the fundamental plane (Paper I, Section 7).

Second, the cluster-scale smooth DM halo is parameterized using either a generalized NFW (gNFW) model

$$\begin{equation} \rho _{{\rm DM}}(r) = \frac{\rho _s}{(r/r_s)^{\beta } (1 + r/r_s)^{3-\beta }} \end{equation} \tag{ 1 }$$

or a cored NFW (cNFW) model with

$$\begin{equation} \rho _{{\rm DM}}(r) = \frac{b \rho _s}{(1 + b r/r_s)(1+r/r_s)^2}. \end{equation} \tag{ 2 }$$

Both models have the same large-scale behavior (ρ_DM∝r⁻³), but the gNFW model contains a central power-law cusp with dlog ρ_DM/dlog r → −β as r → 0, while the cNFW model asymptotes to a constant-density core within a characteristic radius r_core = r_s/b. Both models contain the NFW profile in the limits β = 1 and r_core → 0 and, therefore, allow us to explore deviations from canonical CDM halos in the central regions. Broad, uninformative priors are placed on the halo parameters (ρ_s, r_s, and β or b) and ϒ_* (Paper I, Table 7).

Based on the close alignment between the optical centers of the BCGs and both the X-ray centroids (typically separated by ≃ 3 kpc, comparable to the measurement errors) and the lensing-derived centers of mass, we fix the center of the halo to that of the BCG (Paper I, Sections 2 and 7.3). Furthermore, the mass estimates derived from lensing agree well with independent X-ray observations, which constrains the l.o.s. ellipticity of the halo to be mild in all clusters except A383 (see Paper I, Section 8, and Newman et al. 2011).

We do not specifically distinguish the hot gas in the intracluster medium (ICM), which is thus implicitly included in the halo in our models. Since the distribution of the ICM is similar to that of the halo and comprises only a ≃ 10% mass fraction (e.g., Allen et al. 2004), subtracting the ICM to isolate the DM has very little effect on the slope of the density profile (Δdlog ρ/dlog r ≲ 0.05; Newman et al. 2009; Sommer-Larsen & Limousin 2010), which is the main focus of this paper.

These models are constrained by three data sets. First, the mass on scales of ≃ 100 kpc to 3 Mpc is constrained using gravitational shear (weak lensing) measured in deep, multi-color images primarily from the Subaru telescope. Second, the angular positions and redshifts of background galaxies that are strongly lensed by the clusters precisely constrain the mass from ≃ 20 to 100 kpc, varying from cluster to cluster. In total we located 25 multiply imaged sources, of which 21 have spectroscopic redshifts (7 were first presented in Paper I). Finally, the most unique aspect of our analysis is the inclusion of spatially resolved stellar kinematics within the BCGs. These measures are derived from long-slit spectra primarily obtained at the Keck telescopes. The stellar kinematic data typically extend to R ≈ 10–20 kpc and display a very homogenous shape that rises with radius, indicating a rising total mass-to-light ratio as expected at the centers of massive clusters. As demonstrated in Paper I (Section 9), the mass models provide good fits to the full range of data.

For the purpose of distinguishing dark and stellar mass, the most important physical assumptions are that stellar mass follows light, as justified above, and that the DM halo is adequately described by a gNFW- or cNFW-like profile. The precise parametric form is not as critical as the presumption that the DM density turns over smoothly at small radii—either to a power-law cusp in the gNFW case or to a constant density in the cNFW models—without a sharp upturn on small scales. This is reasonable: by design these profiles describe pure CDM halos in the appropriate limits, and although the effects of adding baryons are uncertain, adiabatic contraction prescriptions (Gnedin et al. 2004, 2011) predict DM profiles that are well fit by gNFW models over the relevant range of radii when applied to the halos and BCGs representative of our sample. Due to the density and radial extent of observational constraints (extended stellar kinematic profiles, strongly lensed galaxies usually at multiple redshifts, weak lensing), we emphasize that we are able to consider quite general families of DM halos in each cluster.

In individual clusters, there is a degeneracy between the stellar mass-to-light ratio ϒ_*V = M_*/L_V and the inner DM slope. This is illustrated in Figure 1, which shows results for the mass models summarized in Section 2 and derived in Paper I (Section 9). This degeneracy is expected, since stellar mass in the BCG can be traded against DM. Owing to the multiplicity of constraints described above, particularly kinematic measurements at small radii in the stellar-dominated regime, the model degeneracy is not complete, and each cluster does carry information on both ϒ_* and β or b.

Zoom In Zoom Out Reset image size

It is already evident in Figure 1 that most of the clusters in our sample prefer a DM inner slope that is shallower than an NFW profile (i.e., β < 1), consistent with our previous findings (Sand et al. 2002, 2004, 2008; Newman et al. 2009, 2011). However, it is also clear that the precision of the constraints on the inner slope could be increased if additional information regarding ϒ_* is available. Indeed, most clusters are consistent with a wide range of ϒ_* when viewed in isolation, due to the uncertainty arising from the degeneracy described above. Furthermore, the figure suggests a possible variation from cluster to cluster in the DM inner slope, but this conclusion may be contingent upon substantial variations in ϒ_* as well. We do not have strong a priori expectations about the possible variation from cluster to cluster in the DM inner slope, particularly recalling the uncertain role of baryons in theoretical predictions. There are, however, several strong reasons to believe that the true physical variation in ϒ_* within our sample is small.

First, Figure 1 shows estimates of the stellar mass-to-light ratio ϒ^SPS_*V derived by fitting stellar population synthesis (SPS) models to the broadband colors of the BCGs (see Paper I, Section 5.2). Currently, SPS models cannot predict absolute masses more accurately than a factor of ≃ 2, primarily due to the unknown stellar initial mass function (IMF), which we discuss further in Section 5.1. On the other hand, relative stellar masses are more robust, especially within a homogeneous galaxy population. As the bottom right panel of Figure 1 demonstrates, the range in ϒ^SPS_*V within our sample at a fixed IMF is small. Assuming a Chabrier (2003) IMF, the median 〈ϒ^SPS_*V〉 = 2.2; the full range is only 1.80–2.32, and the rms scatter is 9%.⁵

Second, the rms dispersion in the absolute luminosities L_V of the BCGs in our sample is only 0.1 dex. This small variation is consistent with previous studies of BCGs as "standard candles" with uniform luminosities and colors (e.g., Sandage 1972; Postman & Lauer 1995; Collins & Mann 1998; Bernardi et al. 2007). Finally, the environments of the BCGs are the same: by construction they are all central galaxies in massive clusters, and their central velocity dispersions are comparable. It would be very surprising if this uniformity in luminosity and ϒ^SPS_*V, which are thought to derive from a similar assembly history, were the result of a conspiracy that masks larger variations in stellar mass. Instead, based on these physical similarities, it is very likely that the BCGs in our sample have similar stellar masses and ϒ_*V. As we discuss in Section 5.1, this is further supported by recent, independent studies.

With the well-motivated assumption that the BCGs in our sample have a similar ϒ_*V, we can use the full sample of seven clusters to jointly constrain its value, thereby improving the precision and robustness of our measurements of the DM profile. Before embarking on this, we consider how to handle the small variations in ϒ_*V that we do anticipate, despite the overall similarity. The sample spans a redshift range z = 0.19–0.31, so some mild passive evolution is expected. Additionally, the BCGs with the lowest ϒ^SPS_*V estimates show optical emission lines and far-infrared photometry indicative of ongoing star formation (although it involves a small fraction of the stellar mass; see Paper I). These BCGs reside in the cool core clusters, consistent with earlier studies (Bildfell et al. 2008; Loubser et al. 2009; Sanderson et al. 2009).

Therefore, a more precise technique is to define ϒ_* for each cluster relative to the SPS measurement:

$$\begin{equation} \log \alpha _{{\rm SPS}} = \log \Upsilon _{*{\rm V}}/ \Upsilon _{*{\rm V}}^{{\rm SPS}}. \end{equation} \tag{ 3 }$$

We can then use the full cluster sample to constrain 〈log α_SPS〉, which parameterizes a common, systematic offset from photometrically derived stellar mass-to-light ratios. As described in Section 5.1, the most probable source for large systematic offsets is an IMF that differs from that assumed in the SPS models: in this case, that of Chabrier. However, our analysis does not depend on the physical origin of the offset, only that is it common among our BCGs. Since the variation in ϒ^SPS_*V is small compared to the range of ϒ_*V explored in our fits (25% versus a factor of 5.3), this approach is not radically different from assuming a common ϒ_*V. However, it improves on that assumption by making use of SPS models to adjust for small differences in ϒ_*V arising from age and dust, while making no assumption on the validity of their absolute mass scale.

3.1. Constraining the Stellar Mass Scale

Figure 2 shows the probability distribution for log α_SPS derived in each cluster. The uncertainty in log α_SPS arises from two sources: that in the ϒ_*V derived from dynamics and lensing, and the uncertainty in ϒ^SPS_*V arising from random photometric errors. In Paper I, we estimated the latter as σ_SPS = 0.07 dex. Thus, the probability distributions for log α_SPS are derived by broadening those for log ϒ_*V by a Gaussian with a dispersion of σ_SPS.

Zoom In Zoom Out Reset image size

We have argued that there are strong a priori reasons to expect that α_SPS is uniform across our sample of BCGs. Using the probability distributions in Figure 2, we can ask whether the lensing and kinematic data are indeed consistent with this assumption. One way to quantify this is to suppose that the true distribution of log α_SPS is Gaussian and infer its intrinsic dispersion σ_log α. The formalism for inferring the probability distribution P(σ_log α) was discussed in Paper I (Section 9, and see Bolton et al. 2012). The preference for non-zero intrinsic scatter can then be assessed by

$$\begin{equation} \Delta P = \sqrt{2 \ln [P(\sigma _{\log \alpha } = \sigma _{{\rm peak}}) / P(\sigma _{\log \alpha } = 0)]}, \end{equation} \tag{ 4 }$$

where σ_peak is the location of the maximum of P(σ_log α). For a Gaussian distribution, ΔP is the number of standard deviations from the mean. We find ΔP = 0.85, i.e., a <1σ preference for intrinsic scatter. Thus, the lensing and kinematic data are consistent with (although they alone cannot prove) our assumption that there is little intrinsic variation in α_SPS within our sample.

With the physically motivated assumption that α_SPS is the same for each BCG, we can constrain its common value simply by multiplying the seven independent probability distributions. The results are shown by the thick curves in Figure 2. Very similar values of log α_SPS = 0.28 ± 0.05 and 0.26 ± 0.05 are derived using the gNFW and cNFW models, respectively, demonstrating that these results do not strongly depend on the exact halo model. Given the closeness of these results, in the following analysis we adopt log α_SPS = 0.27 ± 0.05. Despite marginalizing over fairly general parameterizations of the DM profile, we are able to obtain informative results due to the high density of observational constraints and the sample size.

Taking the median 〈ϒ^SPS_*V〉 = 2.2, we find that log α_SPS = 0.27 corresponds to ϒ_*V = 4.1. In Section 4.3, we describe sources of systematic uncertainty leading to a final estimate log α_SPS = 0.27 ± 0.05^+0.10_−0.16. In Section 5.1, we discuss the physical implications of this result and compare to the recent literature on the stellar mass-to-light ratio and IMF in early-type galaxies.

We now turn to the inner DM density profiles. In our earlier papers, we studied the inner DM density slope β by marginalizing over the uncertainty in ϒ_* separately in each cluster. With the benefit of a larger sample with improved data, we have now combined constraints from seven clusters to arrive at a joint measurement of the stellar mass scale α_SPS (Section 3.1). Incorporating this information, we can now conduct our analysis in a more physically consistent way that recognizes the homogeneity of the BCGs, as well as further reducing the remaining degeneracies between dark and stellar mass.

Technically, we implement the joint constraint on log α_SPS via importance sampling (e.g., Lewis & Bridle 2002), reweighting the Markov chain samples to effectively convert our flat prior on log α_SPS to a Gaussian with mean 〈log α_SPS〉 = 0.27 and dispersion σ = (σ²_α + σ_SPS²)^1/2 = 0.09. Here σ_α = 0.05 dex is the uncertainty in 〈log α_SPS〉, and σ_SPS = 0.07 dex is the random error in ϒ^SPS_*V for each BCG. The latter accounts for the fact that α_SPS refers to a systematic offset from SPS-based mass estimates, but random errors due to photometric noise remain in each cluster.⁶

4.1. Dark and Stellar Mass Profiles

Figure 3 shows the resulting spherically averaged density profiles for the DM halo, BCG stars, and their sum. The results based on gNFW and cNFW models are again quite similar, showing that the choice of parameterization does not strongly affect the derived density profiles. We do not detect an overall preference for one model over the other: the ratio of the total Bayesian evidence is consistent with unity.⁷

Zoom In Zoom Out Reset image size

The black line segment in each panel spans r/r₂₀₀ = 0.003–0.03, which is the interval over which the total density slope γ_tot was defined in Paper I. Its slope r^−1.13 is the average measured in CDM-only cluster simulations from the Phoenix project (Gao et al. 2012b). As quantified in Paper I, the stars and DM sum to produce a slope very close to CDM-only simulations over this interval.

Now we can see that both stars and DM contribute significantly to the mass in this regime: stars dominate the density in the inner radius, while virtually all the mass is DM at the outer radius. This demonstrates a tight coordination between the inner DM profile and the distribution of stars: the NFW-like density slope is not a property of the DM halo or the BCG alone, but of their sum. As noted in Paper I, at yet smaller radii r ≲ 5–10 kpc where stars are dominant—well within the mean effective radius 〈R_e〉 = 30 kpc—the total density profile generally steepens.

As expected, if the total density is NFW-like, then the DM profiles become shallower only on scales where the BCG contributes significantly, roughly within R_e. As we describe in Section 5, our results thus do not conflict with other studies that claim the DM alone follows an NFW profile but are confined to r ≳ R_e. The stellar mass density in our models reaches that of the DM at a median radius of 〈r〉 = 7 kpc. In terms of enclosed mass, equality occurs at 〈r〉 = 12 kpc. Within 5 kpc the median DM fraction is 〈f_DM〉 = 25%, similar to massive field ellipticals (e.g., Auger et al. 2010a), but within their three-dimensional half-light radii r_h the BCGs are far more DM-dominated: 〈f_DM〉 = 80%.

4.2. Inner DM Density Slopes and Core Radii

Figure 4 shows the probability distributions for β (gNFW) and r_core (cNFW) obtained by marginalizing over the other parameters, again weighting the samples to incorporate our joint constraint on α_SPS. Results for the individual clusters are listed in Table 1. Every cluster prefers β < 1, i.e., an inner slope shallower than an NFW model. Thick black lines show constraints on the mean: 〈β〉 = 0.50  ±  0.13 and 〈log r_core/kpc〉 = 1.14 ± 0.13; the method for deriving these is outlined in the Appendix. We note that while the typical r_core ≈ 14 kpc is small, the cNFW profile turns over rather slowly at small radii. Thus, while r_core is the radius where the density falls to half of the corresponding NFW profile, significant deviations extend to r ≃ (3–4)r_core.

Zoom In Zoom Out Reset image size

Table 1. Parameters Describing the Inner DM Profile

Cluster	β (gNFW)	log rcore/kpc (cNFW)
MS2137	0.65+0.23−0.30	0.45+0.38−0.48
A963	0.50+0.27−0.30	0.87+0.61−0.71
A383	0.37+0.25−0.23	0.37+0.72−0.64
A611	0.79+0.14−0.19	0.47+0.39−0.50
A2537	0.23+0.18−0.16	1.67+0.24−0.23
A2667	0.42+0.23−0.25	1.29+0.49−0.49
A2390	0.82+0.13−0.18	0.30+0.53−0.34
Ensemble Average
All clusters	0.50 ± 0.13	1.14 ± 0.13
βaniso = +0.2	0.38+0.09−0.07	1.11+0.14−0.10
βaniso = −0.2	0.64+0.05−0.09	0.96+0.24−0.11
Separate αSPS	0.62 ± 0.14	1.09+0.12−0.21

Notes. Median parameters are shown, obtained after weighting samples to incorporate our joint constraint on α_SPS, as described in the text. Error bars encompass the 16–84th percentiles and account for random uncertainties only; see Section 4.3 for a discussion of systematic errors. Results are shown for individual clusters (top) and for the ensemble mean (bottom), including for several alternative assumptions described in Section 4.3.

Download table as:  ASCIITypeset image

We can also ask whether there is evidence for intrinsic variation in the inner DM profiles. This can be quantified by assuming that the parent distributions of β and log r_core are Gaussian, and using the method described in Section 3.1 to infer its dispersion. We find some evidence for intrinsic scatter with σ_β = 0.22^+0.15_−0.11 and $\sigma _{\log r_{{\rm core}}} = 0.57^{+0.33}_{-0.21}$. Its statistical significance can be assessed with the ΔP statistic (Equation (4)): we derive ΔP = 1.5 and 2.6 for β and log r_core, respectively. This indicates a ≃ 2σ preference for the presence of intrinsic scatter in the inner DM profile shape. While we have focused on relaxed clusters, we expect that this variation would increase if a broader sample of clusters that includes recent mergers were considered.

A possible physical origin of this scatter is illustrated in Figure 5. The gray points in the top panel show the total density slope γ_tot. As described in Paper I, these show mild scatter around the mean slope measured in CDM-only simulations (dashed line; Gao et al. 2012b) over the same radial interval (r/r₂₀₀ = 0.003–0.03). Here we see signs of a correlation with the size of the BCG, with more extended BCGs corresponding to shallower total slopes. The effect on the DM slope (colored points) appears stronger: larger BCGs are hosted by clusters with shallower DM slopes β, or equivalently larger core radii r_core (bottom panel). Such a correlation is necessary for the dark and stellar mass to combine to a similar total density profile. The significance can be assessed using the Spearman rank correlation test. We find a probabilities P₀ = 0.18 and 0.07 of obtaining an equally strong correlation between R_e and β or r_core, respectively, in the null hypothesis of uncorrelated data (see caption to Figure 5).

Zoom In Zoom Out Reset image size

Figure 5 suggests that the DM profile in the cluster core is connected to the build-up of stars in the BCG. We return to this point in Section 6 and discuss the physical scenarios that may explain this. Although the correlations with R_e are most convincing, they are not unique: we find correlations between β or r_core and the stellar mass or luminosity with nearly equal statistical significance. There is no sign of a correlation with the virial mass M₂₀₀ (ρ = 0.11 and 0.04 for the gNFW and cNFW models; see caption to Figure 5).⁸

We emphasize that it is preferable to compare directly to the physical density profiles (Figure 3) when possible, rather than only to the marginalized distributions for β. These results do not imply, for example, that a CDM density profile should be modified simply by maintaining the same r_s and changing β = 1 to β = 0.5. Rather, r_s also shifts in our fits such that significant changes in ρ_DM are kept within r ≲ 30 kpc. This degeneracy is simply a result of the gNFW parameterization.

4.3. Systematic Uncertainties

5.1. Stellar Mass-to-light Ratio

5.2. The Total Inner Density Slope

5.3. The Dark Matter Inner Density Slope

By combining strong lensing, weak lensing, and stellar kinematic observations that extend from ≃ 3 kpc to beyond the virial radius, thus spanning the baryon- to DM-dominated regimes, we constrained flexible, physically motivated models of the dark and stellar mass distributions in seven massive, relaxed galaxy clusters. As discussed extensively in Paper I, the density profiles of stars and DM sum to produce a slope close to CDM-only simulations, at least outside the very central ≈5–10 kpc where stars strongly dominate the mass. In this paper, we isolated the dark and stellar density profiles to quantify the behavior of the DM on small scales, finding a mean asymptotic inner power-law slope of 〈β〉 = 0.50 ± 0.13^+0.14_−0.13, or equivalently a mean DM core radius 〈log r_core〉 = 1.14 ± 0.13^+0.14_−0.22. We also presented evidence for possible variation in the inner DM profile from cluster to cluster, which correlates with the size and mass of the BCG (Figure 5). This suggests a connection between the DM profile in cluster cores and the assembly of stars in the BCG.

The conclusion that the inner DM profile is shallower than that of pure CDM halos is fully consistent with our previous claims (Sand et al. 2002, 2004, 2008; Newman et al. 2009, 2011). We have improved upon these earlier works by collecting improved data for a larger sample of clusters and refining our analysis techniques, as discussed in Section 5.3. A particular advance enabled by this enlarged, improved sample was a joint constraint on the stellar mass-to-light ratio ϒ_* of the BCGs in our sample, which we found to be elevated by 〈log α_SPS〉 = 0.27 ± 0.05^+0.10_−0.16 dex relative to fits to SPS models that assume a Chabrier IMF. Our measurements are instead consistent with a Salpeter IMF (or any equivalently "heavy" IMF). As reviewed in Section 5.1, this agrees with recent, independent studies of massive, early-type galaxies based on lensing, dynamics, and detailed spectroscopy (Treu et al. 2010; Auger et al. 2010b; Cappellari et al. 2012a, 2012b; van Dokkum & Conroy 2010, 2012; Conroy & van Dokkum 2012). These rapid developments in our understanding of stellar populations promise significant advances in disentangling the distributions of dark and baryonic mass across a range of systems.

Figure 6 compares our measurements to high-resolution CDM cluster simulations from the Phoenix project (Gao et al. 2012b), clearly demonstrating these DM-only simulations are a better match to the total density profile than that of the DM alone. In assessing the role of baryons on their host halos, much of the focus of the theoretical literature has been on the halo contraction (e.g., Blumenthal et al. 1986; Gnedin et al. 2004, 2011) expected to result from a central dissipative build-up of baryons. The dot-dashed lines in Figure 6 show the effect of applying the modified adiabatic contraction model of Gnedin et al. (2011) to an NFW halo and BCG with parameters typical of our sample and of the Phoenix simulations (see details in caption). As expected, the DM profile steepens (bottom panel), only worsening the disagreement with our observations. It has been argued that this increase in central DM density from adiabatic contraction will boost the gamma-ray flux from DM annihilation in clusters (Ando & Nagai 2012).¹² However, our results suggest that adiabatic contraction is not the main process that sets the density profile and that the net effect on the halo is actually opposite to its predictions. (As discussed in Paper I, this does not necessarily imply that the same theory cannot make valid predictions at the galaxy scale, where the star formation efficiency and assembly history are very different.)

Zoom In Zoom Out Reset image size

A possible formation scenario is that early, dissipative star formation in the main BCG progenitor creates a steep total density slope in the inner ≃ 5–10 kpc, where stars dominate the mass. This size scale is indeed similar to the observed sizes of very massive galaxies at z ≳ 2.5 (e.g., Trujillo et al. 2006; van Dokkum et al. 2008; Newman et al. 2012). The subsequent assembly of the extended stellar envelope of the BCG—thought to be dominated by low mass, dry accretion of satellites (e.g., Naab et al. 2009; Laporte et al. 2012)—then mostly replaces the DM already in place with satellite material, roughly maintaining the density.

Controlled simulations have indeed shown that dynamical friction between infalling satellites and the DM halo can heat the cusp and reduce the central DM density (e.g., El-Zant et al. 2001, 2004; Nipoti et al. 2004; Jardel & Sellwood 2009; Cole et al. 2011). This process is dissipationless, since the orbital energy lost by the satellites is transferred to the halo, and thus contrasts with the AC picture, in which the baryons' energy is radiated away (Lackner & Ostriker 2010). A connection between the assembled stellar mass and the central DM density is naturally expected. Indeed, Nipoti et al. (2004) find an anti-correlation between the amount of stellar mass assembled in the BCG and the inner DM density slope β, similar to our observations (Figure 5; we note that the satellites in their simulations included no DM). Del Popolo (2012) discusses a similar anti-correlation arising in their analytic models for the same physical reason, with higher central baryon fractions corresponding to shallower DM density cusps. The strength of the dynamical friction effect depends on the density of the satellites and their resistance to stripping. Laporte et al. (2012) showed that when a stellar mass–size relation in line with z ≳ 2 observations is imposed in their simulations (offset by 3–5 times in size from the local relation), the central DM cusp is flattened to β ≃ 0.3–0.7, comparable to our observations. It is important to realize that numerical experiments investigating this dynamical effect have generally lacked a fully realistic and consistent treatment of the satellites, so improved simulations are needed. Nonetheless, the current results are promising.

Until the last few years, full hydrodynamical, cosmological cluster simulations that include cooling, star formation, and feedback did not produce shallow DM cusps or cores, which probably reflected overcooling effects. Mead et al. (2010) and Martizzi et al. (2012) showed that the inclusion of active galactic nucleus (AGN) feedback greatly improves this situation (see discussion and references in Paper I, Section 10) and may also play a key role in lowering the central DM density. Understanding the impact that gas cooling, dynamical friction from stellar "clumps," and AGN feedback have on the small-scale DM distribution is an important avenue for future simulations, and the data we have presented provide strong constraints.

In addition to the effect of baryons on the halo, various DM particle scenarios have also been proposed to reduce tension between CDM and observations on small scales, including the "missing satellites" problem and evidence for central DM cores or shallow cusps (for a recent review, see Primack 2009). These include warm sterile neutrinos at the ∼keV scale (e.g., Abazajian et al. 2001; Boyarsky et al. 2009; Macciò et al. 2012; Menci et al. 2012), "fuzzy" CDM composed from an ultralight scalar particle (Hu et al. 2000; Woo & Chiueh 2009), DM produced from early decays (Kaplinghat 2005), and DM that itself decays with a long timescale (Peter et al. 2010), among many other possibilities. The goal is to preserve the large-scale successes of CDM, while allowing for modifications at higher densities where the detailed properties of the DM particle might manifest. A scenario for which halo density profiles has been worked out in detail is a self-interacting DM particle (Spergel & Steinhardt 2000; Yoshida et al. 2000; Davé et al. 2001). Rocha et al. (2012) and Peter et al. (2012) showed that a cross-section σ ∼ 0.1 cm² g⁻¹ can produce ≈20 kpc cores in clusters without violating any current constraints, e.g., from the asphericity of cluster cores or the Bullet Cluster (Randall et al. 2008). Only the dense central regions of the halo are affected, where scattering can occur within a Hubble time.

These ≈20 kpc core sizes are intriguingly similar to our observations. On the other hand, they are also very similar to the scale of the baryons, i.e., the size of the BCG. It is unclear why the total density profile should then match the shape expected of collisionless CDM. In these scenarios, the core size arises from the microphysics of the DM particle and presumably should not "know" about the size of the central galaxy (Figure 5), for example. Thus, observations of clusters alone cannot provide unambiguous support for alternative DM theories. Global comparisons across a wide range of mass scales (for instance, a cross-section that also produces correct core sizes and densities in dwarf galaxies) remain an essential test for attempts to explain low central halo densities in terms of the DM particle.

It is a pleasure to acknowledge helpful conversations with Annika Peter. We thank Liang Gao for providing the Phoenix simulation results. The anonymous referee is thanked for a helpful report. R.S.E. acknowledges financial support from DOE grant DE-SC0001101. Research support by the Packard Foundation is gratefully acknowledged by T.T. The authors recognize and acknowledge the cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

In Section 4.2, we described how posterior probability distributions P(β) and P(log r_core) are derived for each cluster by weighting the samples in the Markov chains derived in Paper I. The weights

$$\begin{equation} w = \frac{1}{\sqrt{2\pi } \sigma } \exp \left[-\frac{1}{2} \left(\frac{\log \alpha _{{\rm SPS}} - \langle \log \alpha _{{\rm SPS}} \rangle }{\sigma }\right)^2\right] \end{equation} \tag{ A1 }$$

effectively convert a flat prior on log α_SPS (Paper I, Section 7) to a Gaussian with mean 〈log α_SPS〉 = 0.27 and a dispersion σ = (σ²_α + σ_SPS²)^1/2. This dispersion accounts for two sources of error: the uncertainty σ_α = 0.05 dex in the global systematic offset 〈log α_SPS〉 from SPS estimates ϒ^SPS_*V, and the random photometric uncertainty σ_SPS = 0.07 dex in ϒ^SPS_*V for each cluster. The first uncertainty is correlated across the entire sample, while the second is not.

Therefore, to obtain constraints on the mean 〈β〉 and 〈log r_core〉, the probability distributions derived for each cluster in this manner cannot simply be multiplied, since they are not independent. Instead, we calculate the posterior probability of 〈β〉 as

$$\begin{equation} P(\langle \beta \rangle) \propto \int P(\langle \beta \rangle | \log \alpha _{{\rm SPS}}) P(\log \alpha _{{\rm SPS}})\,{d}\alpha _{{\rm SPS}}. \end{equation} \tag{ A2 }$$

Here, P(〈β〉|log α_SPS) is the posterior distribution of 〈β〉 at a fixed value of log α_SPS. It is obtained by multiplying the probability densities P(β|log α_SPS) for the seven clusters in our sample, which are each computed with Gaussian weights centered at the fixed value of log α_SPS and a dispersion σ_SPS (i.e., σ = σ_SPS in Equation (A1); we now account for only the random photometric errors in ϒ^SPS_*V since log α_SPS is fixed). P(log α_SPS), which represents our constraint on the common stellar mass scale, is simply a Gaussian with mean 〈log α_SPS〉 = 0.27 and dispersion σ_α = 0.05 dex, as derived in Section 3.1 for isotropic orbits.

We estimate the intrinsic scatter in β (Section 4.2) using the posterior probability densities P(β|log α_SPS = 0.27) for each cluster. That is, we evaluate the cluster-to-cluster scatter in β at a fixed value of log α_SPS. All of the above comments apply equally to our study of the cNFW models, simply replacing β by log r_core.