A PRECISE CLUSTER MASS PROFILE AVERAGED FROM THE HIGHEST-QUALITY LENSING DATA^*

Clusters of galaxies represent the largest gravitationally bound objects in the universe, which contain a wealth of astrophysical and cosmological information, related to the nature of dark matter, primordial density perturbations, and the emergence of structure over cosmic time. Observational constraints on the properties and evolution of clusters provide independent and fundamental tests of any viable cosmology, structure formation scenario, and possible modifications of the laws of gravity, complementing large-scale cosmic microwave background and galaxy clustering measurements (e.g., Komatsu et al. 2011; Percival et al. 2010).

A key ingredient of cluster-based cosmological tests is the mass and internal mass distribution of clusters. In this context, the current cosmological paradigm of structure formation, the standard Λ cold (i.e., non-relativistic) dark matter (ΛCDM) model, provides observationally testable predictions for CDM-dominated halos over a large dynamical range in density and radius. Unlike galaxies where substantial baryonic cooling is present, massive clusters are not expected to be significantly affected by gas cooling (e.g., Blumenthal et al. 1986; Broadhurst & Barkana 2008). This is because the majority of baryons (∼80%) in massive clusters comprise a hot, X-ray-emitting diffuse intracluster medium (ICM), in which the high temperature and low density prevent efficient cooling and gas contraction, and hence the gas pressure roughly traces the gravitational potential produced by the dominant dark matter (see Kawaharada et al. 2010; Molnar et al. 2010). The ICM represents only a minor fraction of the total mass near the centers of clusters (Lemze et al. 2008; Umetsu et al. 2009). Consequently, for clusters in a state of quasi equilibrium, the form of their total mass profiles reflects closely the distribution of dark matter (Mead et al. 2010).

High-resolution N-body simulations of collisionless CDM exhibit an approximately "universal" form for the spherically averaged density profile of virialized dark-matter halos (Navarro et al. 1997, hereafter Navarro–Frenk–White (NFW)), with some intrinsic variance in the mass assembly histories of individual halos (Jing & Suto 2000; Tasitsiomi et al. 2004; Navarro et al. 2010). The predicted logarithmic gradient γ_3D(r) ≡ −dln ρ/dln r of the NFW form flattens progressively toward the center of mass, with a central cusp slope flatter than a purely isothermal structure (γ_3D = 2) interior to the inner characteristic radius r_s(≲ 300 kpc h⁻¹ for cluster-sized halos), providing a distinctive prediction for the empirical form of CDM halos in gravitational equilibrium. A useful index of the degree of concentration is c_vir = r_vir/r_s, which compares the virial radius r_vir to the characteristic radius r_s of the NFW profile. This empirical NFW profile is characterized by the total mass within the virial radius, M_vir, and the halo concentration c_vir. Theoretical progress has been made in the understanding of the form of this profile in terms of the dynamical structure using both numerical and analytical approaches (Taylor & Navarro 2001; Lapi & Cavaliere 2009; Navarro et al. 2010), though we must currently rely on the quality of N-body simulations when making comparisons with CDM-based predictions for cluster mass profiles.

In the context of standard hierarchical clustering models, the halo concentration should decline with increasing halo mass as dark-matter halos that are more massive collapse later when the mean background density of the universe is correspondingly lower (Bullock et al. 2001; Zhao et al. 2003; Neto et al. 2007). This prediction for the halo mass–concentration relation and its evolution has been established thoroughly with detailed simulations (e.g., Navarro et al. 1997; Bullock et al. 2001; Neto et al. 2007; Duffy et al. 2008; Klypin et al. 2010), although sizable scatter around the mean relation is present due partly to variations in formation epoch of halos (Wechsler et al. 2002; Neto et al. 2007; Zhao et al. 2009). Massive clusters are of particular interest in this context because they are predicted to have a relatively shallow mass profile with a pronounced radial curvature. Gravitational lensing of background galaxies offers a robust way of directly obtaining the mass distribution of galaxy clusters (see Bartelmann & Schneider 2001; Umetsu 2010, and references therein) without requiring any assumptions on the dynamical and physical state of the system (Clowe et al. 2006; Okabe & Umetsu 2008). A detailed examination of this fundamental prediction has been the focus of our preceding work (Broadhurst et al. 2005a, 2008; Medezinski et al. 2007; Umetsu & Broadhurst 2008; Umetsu et al. 2009, 2010, 2011; Lemze et al. 2009).

Systematic cluster lensing surveys are in progress to obtain mass profiles of representative clusters over a wide range of radius by combining high-quality strong- and weak-lensing data. Deep multicolor images of massive cluster cores from Advanced Camera for Surveys (ACS) observations with the Hubble Space Telescope (HST) allow us to identify many sets of multiple images spanning a wide range of redshifts for detailed strong-lens modeling (e.g., Broadhurst et al. 2005b; Zitrin et al. 2009, 2010, 2011a, 2011c, 2011d). The wide-field prime-focus cameras of Subaru and Canada–France–Hawaii Telescope (CFHT) have been routinely producing data of sufficient quality to obtain accurate measurements of the weak-lensing signal, providing model-independent cluster mass profiles out to beyond the virial radius (e.g., Broadhurst et al. 2005a, 2008; Limousin et al. 2007; Umetsu & Broadhurst 2008; Umetsu et al. 2009, 2010, 2011; Coe et al. 2010). Our earlier work has demonstrated that without adequate color information, the weak-lensing signal can be heavily diluted particularly toward the cluster center by the presence of unlensed cluster members, leading to biased cluster mass profile measurements with underestimated concentrations and internal inconsistency, with the weak-lensing-based profile underpredicting the observed Einstein radius (Broadhurst et al. 2005a; Umetsu & Broadhurst 2008; Medezinski et al. 2010).

Careful lensing work on individual clusters has shown that full mass profiles constructed from combined strong- and weak-lensing measurements show a continuous steepening radial trend consistent with the predicted form for the family of collisionless CDM halos (Gavazzi et al. 2003; Broadhurst et al. 2005a, 2008; Umetsu & Broadhurst 2008; Umetsu et al. 2010, 2011). Intriguingly these initial results from combined strong- and weak-lensing measurements reveal a relatively high degree of halo concentration in lensing clusters (e.g., Gavazzi et al. 2003; Kneib et al. 2003; Broadhurst et al. 2008; Oguri et al. 2009; Zitrin et al. 2011d), lying well above the mass–concentration relation for cluster-sized halos predicted by the ΛCDM model, despite careful attempts to correct for potential projection and selection biases inherent to lensing (Hennawi et al. 2007; Meneghetti et al. 2010b). This apparent overconcentration of lensing clusters is also indicated by the generally large Einstein radii determined from strong-lensing data (Broadhurst & Barkana 2008; Meneghetti et al. 2010a; Zitrin et al. 2011a).

In this paper, we explore in greater depth the utility of high-quality lensing data for obtaining highest-precision cluster mass profiles by combining all possible lensing information available in the cluster regime. This extends our recent weak-lensing work by Umetsu et al. (2011), where a Bayesian method was developed for a direct reconstruction of the projected cluster mass profile from complementary weak-lensing distortion and magnification effects (Umetsu & Broadhurst 2008), the combination of which can be used to unambiguously determine the absolute mass normalization. For a massive cluster acting as a super-critical lens, the strong- and weak-lensing regimes contribute equal logarithmic coverage of the radial profile (Umetsu et al. 2011), so that here we concentrate on those clusters for which we have high-quality data in both these regimes. The high quality of our data is such that we have now become significantly limited by the cosmic noise from large-scale structure behind the cluster center, where magnified sources lie at greater distances. This noise is correlated between radial bins, and so can be overcome by stacking clusters, along independent sight lines. Stacking also helps average over the effects of cluster asphericity and substructure (Mandelbaum et al. 2006; Johnston et al. 2007; Okabe et al. 2010; Umetsu et al. 2011), allowing a tighter comparison of the averaged profile with theoretical models. Our aim here is to develop and apply comprehensive methods to a sample of four similarly high-mass lensing clusters (A1689, A1703, A370, and Cl0024+17), for which we have previously identified multiply lensed images and measured weak magnification and distortion effects from deep HST and Subaru observations (Broadhurst et al. 2005b; Umetsu & Broadhurst 2008; Umetsu et al. 2010, 2011; Zitrin et al. 2010; Medezinski et al. 2010, 2011).

The paper is organized as follows. In Section 2, we briefly summarize the basis of cluster gravitational lensing. In Section 3, we outline our comprehensive methods for obtaining projected cluster mass profiles from weak-lensing distortion, magnification, and strong-lensing measurements. In Section 4, we apply our methodology to deep HST and Subaru observations of four massive clusters to derive a mean radial mass profile over the entire radial range, demonstrating how stacking the weak- and strong-lensing signals improves upon the statistical precision of the mass profile determination; we then examine the radial dependence of the stacked cluster mass profile. Finally, we discuss our results and conclusions in Section 5. Throughout this paper, we adopt a concordance ΛCDM cosmology with Ω_m = 0.3, Ω_Λ = 0.7, and h ≡ H₀/(100 km s⁻¹ Mpc⁻¹) = 0.7, unless otherwise noted.

The gravitational deflection of light rays by a cluster can be described by the thin lens equation, which relates the angular position of a lensed image $\mbox{\boldmath $\theta $}$ to the angular position of the intrinsic source $\mbox{\boldmath $\beta $}$ as

$$\begin{equation} \mbox{\boldmath $\beta $}=\mbox{\boldmath $\theta $}-\mbox{\boldmath $\nabla $}\psi, \end{equation} \tag{ 1 }$$

where $\mbox{\boldmath $\alpha $}\equiv \mbox{\boldmath $\nabla $}\psi (\mbox{\boldmath $\theta $})$ is the deflection field and $\psi (\mbox{\boldmath $\theta $})$ is the effective lensing potential, which is defined by the two-dimensional Poisson equation as $\triangle \psi (\mbox{\boldmath $\theta $})=2\kappa (\mbox{\boldmath $\theta $})$ with the lensing convergence κ given as a source term. This equation can be readily inverted to yield: $\psi (\mbox{\boldmath $\theta $})=2\int \!d^2\theta ^{\prime }\,\triangle ^{-1}(\mbox{\boldmath $\theta $},\mbox{\boldmath $\theta $}^{\prime })\kappa (\mbox{\boldmath $\theta $}^{\prime })=(1/\pi)\int \!d^2\theta ^{\prime }\,\ln |\mbox{\boldmath $\theta $}-\mbox{\boldmath $\theta $}^{\prime }|\kappa (\mbox{\boldmath $\theta $}^{\prime })$, so that the deflection field is expressed in terms of κ as

$$\begin{equation} \mbox{\boldmath $\alpha $}(\mbox{\boldmath $\theta $})=\frac{1}{\pi }\int \!d^2\theta ^{\prime }\,\frac{\mbox{\boldmath $\theta $}-\mbox{\boldmath $\theta $}^{\prime }}{|\mbox{\boldmath $\theta $}-\mbox{\boldmath $\theta $}^{\prime }|^2} \kappa (\mbox{\boldmath $\theta $}^{\prime }). \end{equation} \tag{ 2 }$$

For gravitational lensing in the cluster regime (e.g., Umetsu 2010), κ is expressed as $\kappa (\mbox{\boldmath $\theta $})=\Sigma _{\rm crit}^{-1}\Sigma (\mbox{\boldmath $\theta $})$, namely, the projected mass density $\Sigma (\mbox{\boldmath $\theta $})$ in units of the critical surface mass density for gravitational lensing, defined as

$$\begin{equation} \Sigma _{\rm crit} = \frac{c^2}{4\pi G D_l} \beta ^{-1}; \quad \beta (z_s) \equiv {\rm max}\left[ 0,\frac{D_{ls}(z_s)}{D_s(z_s)}\right], \end{equation} \tag{ 3 }$$

where D_s, D_l, and D_ls are the proper angular-diameter distances from the observer to the source, from the observer to the lens, and from the lens to the source, respectively, and β = D_ls/D_s is the angular-diameter distance ratio associated with the population of background sources.

The deformation of the image for a background source can be described by the Jacobian matrix $\cal {A}_{\alpha \beta }\equiv (\partial \mbox{\boldmath $\beta $}/\partial \mbox{\boldmath $\theta $})_{\alpha \beta }=\delta _{\alpha \beta }-\psi _{,\alpha \beta }$ (α, β = 1, 2) of the lens mapping, where δ_αβ is Kronecker's delta.⁷ The real, symmetric Jacobian ${\cal A}_{\alpha \beta }$ can be decomposed as ${\cal A}_{\alpha \beta } = (1-\kappa)\delta _{\alpha \beta } -\Gamma _{\alpha \beta }$, where $\Gamma _{\alpha \beta }(\mbox{\boldmath $\theta $}) \equiv (\partial _\alpha \partial _\beta -\delta _{\alpha \beta }\nabla ^2/2)\psi (\mbox{\boldmath $\theta $})$ is the trace-free, symmetric shear matrix,

$$\begin{equation} \Gamma _{\alpha \beta }= \left(\begin{array}{@{}cc@{}}+{\gamma }_1 & {\gamma }_2 \\ \ms {\gamma }_2 & -{\gamma }_1 \end{array} \right)\!, \end{equation} \tag{ 4 }$$

with γ_α being the components of spin-2 complex gravitational shear γ := γ₁ + iγ₂. In the strict weak-lensing limit where κ, |γ| ≪ 1, Γ_αβ induces a quadrupole anisotropy of the background image, which can be observed from ellipticities of background galaxy images. Given an arbitrary circular loop of radius ϑ on the sky, the average tangential shear γ₊(ϑ) around the loop satisfies the following identity (e.g., Kaiser 1995):

$$\begin{equation} \gamma _+ (\vartheta) = \bar{\kappa }({<}\vartheta)-\kappa (\vartheta), \end{equation} \tag{ 5 }$$

where κ(ϑ) is the azimuthal average of $\kappa (\mbox{\boldmath $\theta $})$ around the loop, and $\bar{\kappa }(<\vartheta)$ is the average convergence within the loop.

The local area distortion due to gravitational lensing, or magnification, is given by the inverse Jacobian determinant

$$\begin{equation} \mu = \frac{1}{|{\rm det}{\cal A}|} = \frac{1}{|(1-\kappa)^2-|\gamma |^2|}, \end{equation} \tag{ 6 }$$

which can influence the observed surface density of background sources, expanding the area of sky, and enhancing the observed flux of background sources (Broadhurst et al. 1995). The former effect reduces the effective observing area in the source plane, decreasing the number of background sources per solid angle; on the other hand, the latter effect amplifies the flux of background sources, increasing the number of sources above the limiting flux. The net effect is known as magnification bias and depends on the intrinsic slope of the luminosity function of background sources.

In general, the observable quantity for quadrupole weak lensing is not the gravitational shear γ but the complex reduced shear

$$\begin{equation} g(\mbox{\boldmath $\theta $})=\frac{\gamma (\mbox{\boldmath $\theta $})}{1-\kappa (\mbox{\boldmath $\theta $})} \end{equation} \tag{ 7 }$$

in the subcritical regime where ${\rm det}\,{\cal A}>0$ (or 1/g* in the negative parity region with ${\rm det}{\cal A}<0$). The spin-2 reduced shear g is invariant under the following global linear transformation:

$$\begin{equation} \kappa (\mbox{\boldmath $\theta $}) \rightarrow \lambda \kappa (\mbox{\boldmath $\theta $}) + 1-\lambda, \ \ \ \gamma (\mbox{\boldmath $\theta $}) \rightarrow \lambda \gamma (\mbox{\boldmath $\theta $}) \end{equation} \tag{ 8 }$$

with an arbitrary scalar constant λ ≠ 0 (Schneider & Seitz 1995). This transformation is equivalent to scaling the Jacobian matrix ${\cal A}(\mbox{\boldmath $\theta $})$ with λ, $\cal {A}(\mbox{\boldmath $\theta $}) \rightarrow \lambda {\cal A}(\mbox{\boldmath $\theta $})$, and hence leaves the critical curves ${\rm det}{\cal A}(\mbox{\boldmath $\theta $})=0$ invariant. Furthermore, the curve $\kappa (\mbox{\boldmath $\theta $})=1$, on which the gravitational distortions disappear, is left invariant under the transformation (Equation (8)).

This mass-sheet degeneracy can be unambiguously broken by measuring the magnification effects, because the magnification μ transforms under the invariance transformation (Equation (8)) as

$$\begin{equation} \mu (\mbox{\boldmath $\theta $}) \rightarrow \lambda ^2 \mu (\mbox{\boldmath $\theta $}). \end{equation} \tag{ 9 }$$

In practice, the lens magnification μ can be measured from characteristic variations in the number density of background galaxies due to magnification bias (Broadhurst et al. 1995; Umetsu et al. 2011) as

$$\begin{equation} n_\mu (\mbox{\boldmath $\theta $})=n_0\mu (\mbox{\boldmath $\theta $})^{2.5s-1}, \end{equation} \tag{ 10 }$$

where n₀ = dN₀(< m_cut)/dΩ is the unlensed number density of background sources for a given magnitude cutoff m_cut, approximated locally as a power-law cut with slope s = dlog₁₀N₀(< m)/dm (s > 0). In the strict weak-lensing limit, the magnification bias is δn_μ/n₀ ≈ (5s − 2)κ. For red background galaxies, the intrinsic count slope s at faint magnitudes is relatively flat, s ∼ 0.1, so that a net count depletion results (Broadhurst et al. 2005a; Umetsu & Broadhurst 2008; Umetsu et al. 2010, 2011). On the other hand, the faint blue background population tends to have a steeper intrinsic count slope close to the lensing invariant slope (s = 0.4). Alternatively, the constant λ can be determined such that the mean κ averaged over the outermost cluster region vanishes, if a sufficiently wide sky coverage is available.⁸

In this section, we outline our methods for obtaining cluster mass profiles in a continuous radial coverage from the central region to beyond the virial radius, by combining independent weak-lensing distortion, magnification, and strong-lensing measurements.

3.1. Cluster Weak Lensing

The relation between observable distortion (g) and underlying convergence (κ) is nonlocal, and the convergence derived from distortion data alone suffers from a mass-sheet degeneracy (Section 2). However, by combining the complementary distortion (g) and magnification (μ) measurements the convergence can be obtained unambiguously with the correct mass normalization.

We construct a discrete convergence profile in the weak-lensing regime from observable lens distortion and magnification profiles, g₊(θ) = γ₊(θ)/[1 − κ(θ)] and n_μ(θ) = n₀μ(θ)^{2.5s − 1} (see Section 3 and Appendix B of Umetsu et al. 2011, for details of weak-lensing profile measurements), following the Bayesian prescription given by Umetsu et al. (2011). The Bayesian approach allows for a full parameter-space extraction of model and calibration parameters. A proper Bayesian statistical analysis is of particular importance to explore the entire parameter space and investigate the parameter degeneracies, arising in part from the mass-sheet degeneracy.

In the Bayesian framework, we sample from the posterior probability density function (PDF) of the underlying signal $\mbox{\boldmath $s$}$ given the data $\mbox{\boldmath $d$}$, $P(\mbox{\boldmath $s$}|\mbox{\boldmath $d$})$. Expectation values of any statistic of the signal $\mbox{\boldmath $s$}$ shall converge to the expectation values of the a posteriori marginalized PDF, $P(\mbox{\boldmath $s$}|\mbox{\boldmath $d$})$. The covariance matrix C of $\mbox{\boldmath $s$}$ is obtained from the resulting posterior sample. In our problem, the signal $\mbox{\boldmath $s$}$ is a vector containing the discrete convergence profile, κ_i ≡ κ(θ_i) with i = 1, 2, ..., N^wl in the weak-lensing regime (θ_i > θ_Ein), and the average convergence within the inner radial boundary θ^wl_min of the weak-lensing data, $\overline{\kappa }_{\rm min}\equiv \overline{\kappa }(<\theta _{\rm min}^{\rm wl})$, so that $\mbox{\boldmath $s$}=\lbrace \overline{\kappa }_{\rm min},\kappa _i\rbrace _{i=1}^{N^{\rm wl}}$, being specified by (N^wl + 1) parameters. The Bayes' theorem states that

$$\begin{equation} P(\mbox{\boldmath $s$}|\mbox{\boldmath $d$}) \propto P(\mbox{\boldmath $s$}) P(\mbox{\boldmath $d$}|\mbox{\boldmath $s$}), \end{equation} \tag{ 11 }$$

where ${\cal L}(\mbox{\boldmath $s$})\equiv P(\mbox{\boldmath $d$}|\mbox{\boldmath $s$})$ is the likelihood of the data given the model ($\mbox{\boldmath $s$}$), and $P(\mbox{\boldmath $s$})$ is the prior probability distribution for the model parameters. The ${\cal L}(\mbox{\boldmath $s$})$ function for combined weak-lensing observations is given as a product of the two separate likelihoods, ${\cal L}_{\rm wl}={\cal L}_g{\cal L}_\mu$, where ${\cal L}_g$ and ${\cal L}_\mu$ are the likelihood functions for distortion and magnification, respectively, as given in Umetsu et al. (2011). The log-likelihood for combined weak-lensing distortion and magnification observations, {g_{+, i}}^Nwl_{i = 1} and {n_{μ, i}}^Nwl_{i = 1}, is given as

$$\begin{equation} {-}2\ln {\cal L_{\rm wl}}=\displaystyle \sum _{i=1}^{N^{\rm wl}} \frac{[g_{+,i}-\hat{g}_{+,i}(\mbox{\boldmath $s$})]^2}{\sigma _{+,i}^2}+ \displaystyle \sum _{i=1}^{N^{\rm wl}} \frac{[n_{\mu,i}-\hat{n}_{\mu,i}(\mbox{\boldmath $s$})]^2}{\sigma _{\mu,i}^2}, \end{equation} \tag{ 12 }$$

where $(\hat{g}_{+,i},\hat{n}_{\mu,i})$ are the theoretical predictions for the corresponding observations; the errors σ_{+, i} for g_{+, i} (i = 1, 2, ..., N^wl), due primarily to the variance of the intrinsic source ellipticity distribution, can be conservatively estimated from the data using bootstrap techniques; the errors σ_{μ, i} for n_{μ, i} (i = 1, 2, ..., N^wl) include both contributions from Poisson errors in the counts and contamination due to intrinsic clustering of red background galaxies (Umetsu et al. 2011).

For each parameter of the model $\mbox{\boldmath $s$}$, we consider a simple flat prior with a lower bound of $\mbox{\boldmath $s$}=0$, that is, $\overline{\kappa }_{\rm min}>0$ and κ_i > 0. Additionally, we account for the calibration uncertainty in the observational parameters, such as the normalization and slope parameters (n₀, s) of the background counts and the relative lensing depth due to population-to-population variations between the background samples used for the magnification and distortion measurements (see Umetsu et al. 2011).

3.2. Cluster Strong Lensing

We apply our well-tested approach to strong-lens modeling, which has previously uncovered large numbers of multiply lensed galaxies in ACS images of many clusters, such as A1689 at z = 0.183 (Broadhurst et al. 2005b), Cl0024+17 at z = 0.395 (Zitrin et al. 2009), 12 high-z MACS clusters (Zitrin et al. 2011a), MS 1358+62 at z = 0.33 (Zitrin et al. 2011c), and A383 at z = 0.188 (Zitrin et al. 2011d). Briefly, the basic assumption adopted is that mass approximately traces light, so that the photometry of the red cluster member galaxies is used as the starting point for our model. Cluster member galaxies are identified as lying close to the cluster sequence by HST multiband photometry.

In the strong-lensing regime, we approximate the large-scale distribution of cluster mass by assigning a power-law mass profile to each cluster galaxy, the sum of which is then smoothed. The degree of smoothing (S) and the index of the power law (q) are the most fundamental parameters determining the cluster mass profile dominated by dark matter. A worthwhile improvement in fitting the location of the lensed images is generally found by expanding to first order the gravitational potential of this smooth component, equivalent to a coherent external shear Γ^ex_αβ (α, β = 1, 2) describing the overall matter ellipticity. The direction ϕ_ex of the spin-2 external shear Γ^ex_αβ and its amplitude |γ_ex| are free parameters, allowing for some flexibility in the relation between the distribution of dark matter and the distribution of galaxies, which cannot be expected to trace each other in detail.

The total deflection field $\mbox{\boldmath $\alpha $}(\mbox{\boldmath $\theta $})=\sum _j\mbox{\boldmath $\alpha $}_j(\mbox{\boldmath $\theta $})=(\Sigma _{\rm crit}^{-1}/\pi)\int \!d^2\theta ^{\prime } (\mbox{\boldmath $\theta $}-\mbox{\boldmath $\theta $}^{\prime })/|\mbox{\boldmath $\theta $}-\mbox{\boldmath $\theta $}^{\prime }|^2\sum _j\Sigma _j(\mbox{\boldmath $\theta $}^{\prime })$ consists of the galaxy component $\mbox{\boldmath $\alpha $}_{\rm gal}(\mbox{\boldmath $\theta $})$, scaled by a factor K, the smooth cluster dark-matter component $\mbox{\boldmath $\alpha $}_{\rm DM}(\mbox{\boldmath $\theta $})$, scaled by (1 − K), and the external-shear component $\mbox{\boldmath $\alpha $}_{\rm ex}(\mbox{\boldmath $\theta $})$

$$\begin{equation} \mbox{\boldmath $\alpha $}(\mbox{\boldmath $\theta $})=K\mbox{\boldmath $\alpha $}_{\rm gal}(\mbox{\boldmath $\theta $})+(1-K)\mbox{\boldmath $\alpha $}_{\rm DM}(\mbox{\boldmath $\theta $}) +\mbox{\boldmath $\alpha $}_{\rm ex}(\mbox{\boldmath $\theta $}), \end{equation} \tag{ 13 }$$

where $\alpha _{{\rm ex},\alpha }(\mbox{\boldmath $\theta $})=\left(\Gamma ^{\rm ex}\right)_{\alpha \beta }\Delta \theta _\beta$ with $\Delta \mbox{\boldmath $\theta $}$ being the displacement vector of the angular position $\mbox{\boldmath $\theta $}$ with respect to a fiducial reference position. The overall normalization ${\cal N}$ of the model and the relative scaling K of the smooth dark-matter component versus the galaxy contribution bring the total number of free parameters in the model to 6. This approach to strong lensing is sufficient to accurately predict the locations and internal structure of multiple images, since in practice the number of multiple images uncovered readily exceeds the number of free parameters, so that the fit is fully constrained.

We use this preliminary model to delens the more obvious lensed galaxies back to the source plane by subtracting the derived deflection field. We then relens the source plane in order to predict the detailed appearance and location of additional counter images, which may then be identified in the data by morphology, internal structure, and color. The best-fit strong-lensing model is assessed by minimizing the χ² value in the image plane:

$$\begin{equation} \chi ^2_{\rm sl}=\displaystyle \sum _i{ \frac{[\mbox{\boldmath $\theta $}_i-\hat{\mbox{\boldmath $\theta $}}_i(q,S,{\cal N},K,\Gamma ^{\rm ex})]^2}{\sigma _i^2}}, \end{equation} \tag{ 14 }$$

where i runs over all lensed images, $\hat{\mbox{\boldmath $\theta $}}_i(q,S,{\cal N},\Gamma ^{\rm ex})$ is the position given by the model, $\mbox{\boldmath $\theta $}_i$ is the observed image position, and σ_i is the positional measurement error. For each model parameter, we estimate the 1σ uncertainty by Δχ² ≡ χ² − χ²_min = 1 in the six-parameter space. The uncertainties for the $\Sigma (\mbox{\boldmath $\theta $})$ field and the Σ(θ) profile are estimated by propagating the errors on the strong-lens model parameters, $(q,S,{\cal N},K,\Gamma ^{\rm ex})$.

3.3. Combining Weak and Strong Lensing

We derive a full-radial mass profile on an individual cluster basis by combining independent weak- and strong-lensing data, which can be compared for consistency in the region of overlap. In order to obtain meaningful radial profiles, one must carefully define the center of the cluster. It is often assumed that the cluster mass centroid coincides with the position of the brightest cluster galaxy (BCG), whereas the BCGs can be offset from the mass centroids of the corresponding dark-matter halos (Johnston et al. 2007; Oguri et al. 2010; Oguri & Takada 2011). Umetsu et al. (2011) adopted the location of the BCG as the cluster center in their one-dimensional profile analysis of five massive clusters. A small offset of typically ≲ 20 kpc h⁻¹ ≡ d_off is found by Umetsu et al. (2011) between the BCG and the dark-matter center of mass recovered from strong-lens modeling (Section 3.2). In the following, we will adopt the BCG position as the cluster center, and limit our analysis to radii greater than R_min ≡ 2d_off = 40 kpc h⁻¹, beyond which the cluster miscentering effects on the Σ profile are negligible (see Section 10 of Johnston et al. 2007).

Having defined the cluster center, we can construct a joint discrete mass profile $\mbox{\boldmath $\Sigma $}=\lbrace \Sigma (R_i)\rbrace _{i=1}^{N}$ as a function of the projected radius R = D_lθ by combining the weak- and strong-lensing κ profiles: Σ(R_i) = w⁻¹_iκ(θ_i) (i = 1, 2, ..., N), where w_i is the lensing efficiency function, or the inverse critical surface mass density, w_i ≡ (Σ_{crit, i})⁻¹ = (4πG/c²)D_lβ_i, Note, the i dependence arises because strong- and weak-lensing profiles with different depths are combined together. To simplify the analysis, we exclude the strong-lensing data points in the region of overlap (typically, θ_Ein ≲ θ ≲ 2θ_Ein) as well as the central weak-lensing bin $\overline{\kappa }_{\rm min}$, when defining the joint Σ profile.

The formulation thus far allows us to derive covariance matrices C^stat_ij of statistical measurement errors for individual cluster $\mbox{\boldmath $\kappa $}$ profiles. Here, we take into account the effect of uncorrelated large-scale structure projected along the line of sight on the error covariance matrix C^lss_ij as C = C^stat + C^lss, where C^lss is given as (Schneider et al. 1998; Hoekstra 2003; Dalal et al. 2005; Hoekstra et al. 2011; Oguri & Takada 2011)

$$\begin{equation} C^{\rm lss}_{ij} = \int \!\frac{l\,dl}{2\pi }\, C^{\kappa \kappa }(l)\, \hat{J}_0(l\theta _i)\hat{J}_0(l\theta _j). \end{equation} \tag{ 15 }$$

Here, C^κκ(l) is the weak-lensing power spectrum as a function of angular multipole l evaluated for a given source population and a cosmology, and $\hat{J}_0(l\theta _i)$ is the Bessel function of the first kind and order zero averaged over the ith annulus between θ_{i, 1} and θ_{i, 2}(> θ_{i, 1}), given as

$$\begin{equation} \hat{J}_0(l\theta _i)=\frac{2}{(l\theta _{i,2})^2-(l\theta _{i,1})^2} [ l\theta _{i,2}J_1(l\theta _{i,2})-l\theta _{i,1}J_1(l\theta _{i,1}) ]. \end{equation} \tag{ 16 }$$

We will assume the concordance ΛCDM cosmological model of Komatsu et al. (2011) and use the fitting formula of Peacock & Dodds (1996) to compute the nonlinear mass power spectrum that enters in Equation (15).

3.4. Stacked Lensing Analysis

The utility of high-quality data is ultimately limited by the cosmic noise from large-scale structure along the line of sight, producing covariance between radial bins, particularly behind the cluster center, where magnified sources lie at greater distances. This noise is correlated between radial bins, but can be overcome by stacking an ensemble of clusters along independent lines of sight. Stacking also helps average over the effects of cluster asphericity and substructure inherent in projected lensing measurements. The statistical precision can be greatly improved by stacking together a number of clusters, especially on small angular scales (see Okabe et al. 2010), allowing a tighter comparison of the averaged profile with theoretical models.

With the full mass profiles of individual clusters from combined weak and strong lensing (Section 3.3), we can stack the clusters to produce an averaged radial mass profile. Here, we re-evaluate the mass profiles of the individual clusters in M logarithmically spaced radial bins in the range R = [R_min, R_max] following the prescription given in Umetsu et al. (2011). Since the noise in different clusters is uncorrelated, the mass profiles of individual clusters can be co-added according to (Umetsu et al. 2011)

$$\begin{equation} \langle \mbox{\boldmath $\Sigma $}\rangle = \left(\displaystyle \sum _n {\cal W}_n \right)^{-1} \, \left(\displaystyle \sum _n{ {\cal W}_n \mbox{\boldmath $\Sigma $}_n} \right), \end{equation} \tag{ 17 }$$

where the index n runs over all clusters, $\mbox{\boldmath $\Sigma $}_n$ is a vector containing the discrete surface mass density profile for the nth cluster, and ${\cal W}_n$ is the window matrix defined as

$$\begin{equation} ({\cal W}_n)_{ij} \equiv \big(C_n^{-1}\big)_{ij} (w_n)_i (w_n)_j \end{equation} \tag{ 18 }$$

with (C_n)_ij and (w_n)_i (i = 1, 2, ..., M) being the full covariance matrix and the lensing efficiency function for the nth cluster, respectively. The error covariance matrix for the stacked mass profile $\langle \mbox{\boldmath $\Sigma $}\rangle$ is obtained as

$$\begin{equation} {\cal C} =\left(\displaystyle \sum _n {\cal W}_n\right)^{-1}. \end{equation} \tag{ 19 }$$

4.1. Cluster Sample and Lensing Data

4.2. Results

Our weak- and strong-lensing data together cover a wide range of radius ranging typically from R ∼ 10 kpc h⁻¹ to 2000–3500 kpc h⁻¹ (Umetsu et al. 2011), depending on the cluster redshift as limited by the field of view of Subaru/Suprime-Cam (34' × 27'). Table 1 lists for each cluster the radial ranges R = [R^sl_min, R_max^sl] and [R^wl_min, R_max^wl] of strong- and weak-lensing measurements, respectively, used to define a joint discrete mass profile $\mbox{\boldmath $\Sigma $}=\lbrace \Sigma (R_i)\rbrace _{i=1}^{N}$, given in a total of N radial bins spanning from R_min = R^sl_min to R_max = R^wl_max. In Table 1, we also quote values of the total signal-to-noise ratio (S/N) in our joint cluster mass profiles ($\mbox{\boldmath $\Sigma $}$) obtained using the full covariance matrix C. We find that ignoring the cosmic noise contribution (Equation (15)) will underestimate the errors by ∼30%–40%. To evaluate C^lss for strong-lensing observations, we projected the matter power spectrum out to a fiducial depth of z_s = 2, which is a typical source redshift of strongly lensed arcs in clusters at intermediate redshifts. We used the estimated mean source redshifts given in Table 3 of Umetsu et al. (2011) for weak lensing.

We show in the top panel of Figure 1 the resulting averaged radial mass profile 〈Σ(R)〉 in M = 15 logarithmically spaced bins with its statistical 1σ uncertainty (given as the square root of the diagonal part of the full covariance matrix ${\cal C}$), obtained by stacking the four clusters using Equations (17) and (19). Note, no scaling has been applied to match the mass normalizations between the four clusters, which span a relatively narrow range in mass, 1.3 ≲ M_vir/(10¹⁵ M_☉ h⁻¹) ≲ 2.3 (see Table 6 of Umetsu et al. 2011). For our sample, we find a sensitivity-weighted average cluster redshift of 〈z_l〉 ≃ 0.32, which is fairly close to the simple average of $\overline{z_l}=0.31$ due to the narrow redshift coverage of our cluster sample. The stacked mass profile exhibits a smooth radial trend with a clear radial curvature over a wide range of radius from R = 40 kpc h⁻¹ to 2800 kpc h⁻¹ ≈ 1.4r_vir, and is detected at a high significance level of 58σ, with the contribution from cosmic covariance included. Here, the maximum radius for the stacking analysis represents approximately the average maximum radius 〈R_max〉 covered by our data. Also shown in Figure 1 is the cosmic noise contribution, which increases toward the cluster center. A noticeable increase of the stacked cosmic noise is seen at R ≲ 150 kpc h⁻¹, within which the averaged profile is dominated by strong-lensing measurements with greater depth. In the bottom panel of Figure 1, we plot the logarithmic density slope γ_2D(R) ≡ −dln 〈Σ〉/dln R of the stacked mass profile. The logarithmic gradient of the average profile shows a clear continuous steepening with increasing radius in projection.

Zoom In Zoom Out Reset image size

To quantify and characterize the averaged cluster mass distribution, we compare the 〈Σ〉 profile with the physically and observationally motivated NFW model. Here, we consider a generalized parameterization of the NFW model (hereafter gNFW) of the form (Zhao 1996; Jing & Suto 2000)

$$\begin{equation} \rho (r)=\frac{\rho _s}{(r/r_s)^\alpha (1+r/r_s)^{3-\alpha }}, \end{equation} \tag{ 20 }$$

where ρ_s is the characteristic density, r_s is the characteristic scale radius, and α represents the inner slope of the density profile. This reduces to the NFW model for α = 1. We introduce the radius r₋₂ at which the logarithmic slope of the density is isothermal, i.e., γ_3D = 2. For the gNFW profile, r₋₂ = (2 − α)r_s, and thus the corresponding concentration parameter reduces to c₋₂ ≡ r_vir/r₋₂ = c_vir/(2 − α). We specify the gNFW model with the central cusp slope, α, the halo virial mass, M_vir, and the concentration, c₋₂ = c_vir/(2 − α). We employ the radial dependence of the gNFW lensing profiles given by Keeton (2001).

First, when the central cusp slope is fixed to α = 1 (NFW), the best-fit model for the averaged 〈Σ〉 profile is obtained as M_vir = 1.54^+0.11_{− 0.10} × 10¹⁵ M_☉ h⁻¹ and c₋₂ = c_vir = 7.68^+0.42_{− 0.40} with the minimized χ² value (χ²_min) of 5.8 for 13 degrees of freedom (dof), corresponding to a Q-value goodness-of-fit of Q = 0.952. This model yields an Einstein radius of θ_Ein = 399^+4.4_{− 4.1} for a fiducial source at z_s = 3. The resulting best-fit NFW parameters from the stacked analysis are consistent with the respective sample-weighted means of the individual NFW model fits obtained by Umetsu et al. (2011, Table 6): 〈M_vir〉 = 1.44 ± 0.11 × 10¹⁵ M_☉ h⁻¹ and 〈c_vir〉 = 7.76 ± 0.79. Next, when α is allowed to vary, a gNFW fit to $\langle \mbox{\boldmath $\Sigma $}\rangle$ gives M_vir = 1.50^+0.14_{− 0.13} × 10¹⁵ M_☉ h⁻¹, c₋₂ = 7.91^+0.72_{− 0.75}, and α = 0.89^+0.27_{− 0.39} with χ²_min/dof = 5.7/12 and Q = 0.931 (θ_Ein = 384^+12.2_{− 10.2} at z_s = 3), being consistent with a simple NFW model with α = 1. Thus, the addition of the α parameter does not improve the fit substantially, as shown by the quoted χ² and Q-values (see also Zitrin et al. 2011d). The two-dimensional marginalized constraints (68.3%, 95.4%, and 99.7% confidence levels) on (M_vir, α) and (c₋₂, α) are shown in Figure 2. Finally, a force fit to the singular isothermal sphere (SIS) model (ρ∝r⁻²) yields a poor fit with χ²_min/dof = 78.5/14, so that the SIS model is strongly disfavored at 62σ significance from a likelihood-ratio test, based on the difference between χ² values of the best-fit NFW and SIS models: Δχ² ≡ χ²_{SIS, min} − χ²_{NFW, min} = 72.6 for 1 dof.

Zoom In Zoom Out Reset image size

We have developed a method for improving the statistical precision of cluster mass profiles, combining independent weak-lensing distortion, magnification, and strong-lensing measurements. This extends recent weak-lensing work by Umetsu et al. (2011) to include the central strong-lensing information in a stacking analysis, for full-radial coverage. Our methods take into account the cosmic covariance from uncorrelated large-scale structure projected along the line of sight (Hoekstra 2003; Hoekstra et al. 2011), as well as the effect of different cluster redshifts, so that error propagation in terms of lensing efficiency of individual clusters can be properly averaged.

We have applied our method to a sample of four similarly high-mass lensing clusters (A1689, A1703, A370, and Cl0024+17), for which we have previously identified multiply lensed images and measured weak magnification and distortion effects from deep HST and Subaru observations (Broadhurst et al. 2005b; Umetsu & Broadhurst 2008; Umetsu et al. 2010, 2011; Zitrin et al. 2010; Medezinski et al. 2010, 2011). For our sample of massive clusters, the strong- and weak-lensing regimes contribute equal logarithmic coverage of the radial profile and can be compared for consistency in the region of overlap. We have formed an averaged radial mass profile 〈Σ(R)〉 from stacking the clusters (Figure 1), which shows a progressive steepening with increasing radius from R = 40 kpc h⁻¹ to 2800 kpc h⁻¹. The inner radial boundary is chosen to be sufficiently large to avoid smoothing from cluster miscentering effects (Johnston et al. 2007), where the typical offset between the BCG and the dark-matter center is estimated as d ≲ 20 kpc h⁻¹ for our sample from our detailed strong-lens modeling (see Section 3.3). The stacked full mass profile is detected at a high significance level of 58σ over the entire radial range. It is found here that ignoring the cosmic noise contribution will underestimate the errors by ∼30%–40%. This is due to the correlation of this noise between radial bins and can only be reduced by averaging over independent lines of sight, with uncorrelated line-of-sight structures, i.e., by averaging over well-separated clusters.

Our stacked projected mass profile with a continuously steepening radial trend is very accurately described by the NFW form predicted for the family of CDM-dominated halos, whereas it strongly disfavors the SIS model at 62σ significance. In the context of an assumed gNFW profile, the central cusp slope is constrained as α = 0.89^+0.27_{− 0.39} (at r ≳ 0.02r_vir; see Figures 1 and 2), being consistent with, but slightly shallower than, the simple NFW form with α = 1. Our results are in agreement with recent high-resolution simulations, which find asymptotic inner slopes somewhat shallower than unity, γ_3D(r → 0) ≲ 0.9, for galaxy- and cluster-sized ΛCDM halos (e.g., Merritt et al. 2006; Graham et al. 2006; Navarro et al. 2010). Note NFW define this profile for halos which they identify as in virial equilibrium, in terms of the simulated CDM particles (see Section 2.2.2 of Navarro et al. 1997). The clusters we have selected for our stacked analysis are, in terms of their lensing properties, very well behaved with at most only ∼10% perturbations in mass visible locally in the two-dimensional mass distribution, and otherwise very symmetric over most of the radius (Broadhurst et al. 2005a, 2005b, 2008; Umetsu et al. 2010). Detailed hydrodynamical simulations show that equilibrium is relatively rapidly achieved in only a few sound crossing times after a major merger, though some dynamical and gas disruption may continue for over a Gyr. This is not important in terms of the central relaxation time of the dark matter (Ricker & Sarazin 2001; Umetsu et al. 2010).

An accurate measurement of the cluster mass profile enables us to constrain dark-matter models. Recently, Woo & Chiueh (2009) examined in detail an extremely light bosonic dark matter (ELBDM) model (m ∼ 10²² eV) as an alternative to CDM in the context of nonlinear cosmic structure formation. ELBDM with a de-Broglie wavelength of astronomical length scales, if it exists, may well be in a ground-state Bose–Einstein condensate and hence well described by a coherent wave function, which may naturally account for the perceived lack of small galaxies relative to the ΛCDM model (Klypin et al. 1999; Peebles & Nusser 2010). Woo & Chiueh (2009) showed that, irrespective of whether halos form through accretion or merger, ELBDM halos can form steepening density profiles of the form similar to the standard CDM, but with perhaps a steeper central cusp slope of γ_3D ≃ 1.4 and a shallower outer slope of γ_3D ≃ 2.5. During a merger between condensates interesting large-scale interference occurs which will differ markedly from standard collisionless CDM, and it will be important to explore this class of dark matter further via more extensive and detailed simulations for testing against accurate lensing profiles of both relaxed and merging clusters.

The mean concentration for the four massive lensing clusters considered here is found to be c_vir = 7.68^+0.42_{− 0.40} (at a mean virial mass M_vir = 1.54^+0.11_{− 0.10} × 10¹⁵ M_☉ h⁻¹), which is apparently higher than the standard ΛCDM predictions evaluated at the mean redshift 〈z_l〉 = 0.32 of our sample: c_vir = 4.5^+1.3_{− 1.0} (the errors quoted represent a 1σ lognormal scatter of σ[log₁₀c_vir] = 0.11) for relaxed clusters derived by Duffy et al. (2008) from N-body simulations based on the Wilkinson Microwave Anisotropy Probe (WMAP) five-year data and c_vir ≈ 4.4 by Klypin et al. (2010) from the recent Bolshoi ΛCDM N-body simulation. More recent results with greater mass resolution based on four large N-body simulations (Bolshoi, MultiDark, Millennium-I and II) exhibit a complex mass and redshift dependence of the median concentration, namely, a flattening and upturn of concentration at very high mass and redshift (Prada et al. 2011). Accordingly, their concentrations derived for cluster-sized halos (i.e., rare objects corresponding to high-σ peaks in the primordial density field) are substantially higher than previous results based on smaller simulations. Interestingly, they find a concentration of c_vir ∼ 7 for their most-massive relaxed halos with M_vir ≈ 10¹⁵ M_☉ h⁻¹ at z = 0 (Figure 15 of Prada et al. 2011). A comparison between our results and the ΛCDM predictions (Duffy et al. 2008; Klypin et al. 2010; Prada et al. 2011) is given in Figure 3.

Zoom In Zoom Out Reset image size

An accurate characterization of the observed sample is crucial for any cluster-based cosmological tests. In the extreme case, those clusters identified by the presence of a giant arc represent the most lensing-biased population. Calculations of the enhancement of the projected mass and hence boosted Einstein radii (say, θ_Ein > 20'') find a statistical bias of ∼34% derived from N-body simulations of the ΛCDM model (Hennawi et al. 2007). Semianalytical simulations incorporating idealized triaxial halos yield a ∼50% bias correction (Oguri & Blandford 2009). Applying a conservative 50% bias correction, we find a discrepancy of about 1.8σ with respect to the ΛCDM predictions by the Duffy et al. (2008) model for relaxed clusters (see Figure 3). If this large bias (∼50%) is coupled to a sizable intrinsic scatter in concentration, estimated for the full halo population to be σ[log₁₀c_vir] = 0.11–0.15, then our measurements can come into line with standard ΛCDM.

The results presented here are very favorable in terms of the standard explanation for dark matter, as collisionless and non-relativistic, interacting only via gravity, with a very precise match between our composite mean mass profile, and that of the general form of the mass profile advocated for massive halos in virial equilibrium. The relatively high concentration we obtain for the averaged profile is consistent with previous lensing work which similarly detected a concentration excess in the lensing-based measurements for many individual relaxed strong-lensing clusters (e.g., Gavazzi et al. 2003; Kneib et al. 2003; Comerford & Natarajan 2007; Broadhurst et al. 2008; Oguri et al. 2009). This possibly interesting tension between cluster lensing observations and ΛCDM models can be more definitively addressed with full-lensing data for new cluster surveys, such as CLASH,¹⁰ LoCuSS, Subaru Hyper Suprime-Cam, and XMM-XXL (Pierre et al. 2011), to meaningfully examine the c_vir–M_vir relation over a wider mass and redshift range when applied to sizable samples of relaxed clusters. It is highly desirable to cover the full profile by combining accurate weak- and strong-lensing measurements, requiring several sets of multiple images over a wide range of source redshift, to obtain a meaningful model-independent inner profile and to add weak lensing with sufficient color information to exclude the otherwise sizable dilution effect on the weak-lensing signal from foreground and cluster members. The CLASH survey is in particular designed to generate such useful data free of systematics in both the weak and strong regime, with first results for the substantial smaller mass cluster A383 with M_vir = 5.37^+0.70_{− 0.63} × 10¹⁴ M_☉ h⁻¹ (Zitrin et al. 2011d) showing similar behavior (c_vir = 8.77^+0.44_{− 0.42}).

We thank the anonymous referee for a careful reading of the manuscript and for providing useful comments. We are very grateful for discussions with Nobuhiro Okabe, Sandor Molnar, Jack Sayers, and Alister Graham, whose comments were very helpful. We thank Nick Kaiser for making the IMCAT package publicly available. The work is partially supported by the National Science Council of Taiwan under the grant NSC97-2112-M-001-020-MY3. K.U. acknowledges support from the Academia Sinica Career Development Award.