CLASH: WEAK-LENSING SHEAR-AND-MAGNIFICATION ANALYSIS OF 20 GALAXY CLUSTERS^*

Since a given flux limit corresponds to different intrinsic luminosities at different source redshifts, count measurements of distinctly different background populations probe different regimes of magnification-bias effects. Deep multi-band photometry allows us to explore the faint end of the intrinsic luminosity function of red galaxies at z ∼ 1 (e.g., Figures 11–13 of Ilbert et al. 2010). For a flux-limited sample of the faint red background population, the effect of magnification bias is dominated by the geometric area distortion because relatively few fainter galaxies can be magnified into the sample, thus resulting in a net count depletion (Taylor et al. 1998; Broadhurst et al. 2005; Umetsu & Broadhurst 2008; Umetsu et al. 2011a, 2012; Coe et al. 2012; Medezinski et al. 2013).

In the present work, we perform count-depletion measurements using flux-limited samples of red background galaxies. If the magnitude shift δm = 2.5log₁₀μ induced by magnification is small compared to that on which the logarithmic slope of the source luminosity function varies, the source number counts can be locally approximated by a power law at the limiting magnitude m. The magnification bias at redshift z is then given by (Broadhurst et al. 1995)

$$\begin{eqnarray} N_\mu ({\boldsymbol {\theta }},z; <m) &=& \overline{N}_\mu \,\mu ^{2.5s-1}({\boldsymbol {\theta }},z) \equiv \overline{N}_\mu b_\mu ({\boldsymbol {\theta }},z; s), \end{eqnarray} \tag{ 8 }$$

where $\overline{N}_\mu {=}\overline{N}_\mu (z;<m)$ is the unlensed mean source counts and s is the logarithmic count slope evaluated at m, $s=d\log _{10} \overline{N}_\mu (z; <m)/dm$. In the regime where 2.5s ≪ 1, the effect of magnification bias is dominated by the geometric expansion of the sky area and hence is not sensitive to the exact form of the intrinsic source luminosity function.

Taking into account the spread of $\overline{N}_\mu (z)$, we express the population-averaged magnification bias as $\langle b_\mu \rangle =\left[\int _0^\infty \!{dz}\,b_\mu (W[z])\overline{N}_{\mu }(z) \right] \left[\int _0^\infty \!dz\,\overline{N}_{\mu }(z) \right]^{-1}$. In this work, we use the following equation to interpret the observed count-depletion measurements (Umetsu 2013):

$$\begin{equation} \langle b_\mu ({\boldsymbol {\theta }})\rangle =N_\mu ({\boldsymbol {\theta }}; <m)/\overline{N}_\mu ({<}m)\approx \langle \mu ^{-1}({\boldsymbol {\theta }})\rangle ^{1-2.5s_{\rm eff}}\protect, \protect \end{equation} \tag{ 9 }$$

with $s_{\rm eff}=d\log _{10}\overline{N}_{\mu }({<}m)/dm$ the effective count slope defined in analogy to Equation (8). Equation (9) is exact for s_eff = 0 and gives a good approximation for depleted background populations with s_eff ≪ 0.4 (Umetsu 2013). Furthermore, the source-averaged inverse magnification is approximated as (Umetsu 2013)

$$\begin{eqnarray} \langle \mu ^{-1}\rangle &=& (1-\langle \kappa \rangle)^2 - |\langle \gamma \rangle |^2 +(f_{W,\mu }-1) (\langle \kappa \rangle ^2 - |\langle \gamma \rangle |^2)\nonumber \\ &\approx & (1-\langle \kappa \rangle)^2 - |\langle \gamma \rangle |^2, \end{eqnarray} \tag{ 10 }$$

where $f_{W,\mu }\equiv \langle W^2\rangle _\mu /\langle W\rangle _\mu ^2$ is of the order unity, 〈κ〉 = 〈W〉_μκ_∞, and 〈γ〉 = 〈W〉_μγ_∞. The error associated with the approximation above is 〈Δμ⁻¹〉 = (f_{W, μ} − 1)(〈κ〉² − |〈γ〉|²) ≡ Δf_{W, μ}(〈κ〉² − |〈γ〉|²), which is much smaller than unity for source populations of our concern (Δf_{W, μ} ∼ O(10⁻²)) in the regime where 〈κ〉 ∼ |〈γ〉| ∼ O(10⁻¹).

In practical observations, the nonvanishing and unresolved angular correlation on small angular scales can lead to a significant increase in the variance of counts in cells (van Waerbeke 2000), which can be much larger than the lensing signal in a given cell. To obtain a clean measure of the lensing signal, such intrinsic variance needs to be downweighted (Umetsu & Broadhurst 2008) and averaged over a sufficiently large sky area.

This local clustering noise can be largely overcome by performing an azimuthal average around the cluster (Umetsu et al. 2011a; Umetsu 2013). Here we calculate the surface number density of background sources n_μ(θ) = dN_μ(θ)/dΩ as a function of cluster radius, by azimuthally averaging the source counts in cells, $N_\mu ({\boldsymbol {\theta }}; <m)$. The source-averaged magnification bias is then expressed as $\langle n_\mu (\theta)\rangle = \overline{n}_\mu \langle \mu ^{-1}(\theta)\rangle ^{1-2.5s_{\rm eff}}$ with $\overline{n}_\mu =d\overline{N}_\mu (<m)/d\Omega$ the unlensed mean surface number density of background sources.

In this work, we adopt the following prescription.

• 1.  
  A positive tail of >νσ cells is excluded in each annulus by iterative σ clipping with ν = 2.5 to reduce the effect of intrinsic angular clustering of source galaxies. We take the systematic change between the mean counts estimated with and without νσ clipping as a systematic error, $\sigma _\mu ^{\rm sys}(\theta)=|n_\mu ^{(\nu)}(\theta)-n_\mu ^{(\infty)}(\theta)|/\nu$, where $n_\mu ^{(\nu)}(\theta)$ and $n_\mu ^{(\infty)}(\theta)$ are the clipped and unclipped mean counts in the annulus θ, respectively. The statistical Poisson uncertainty $\sigma _{\mu }^{\rm stat}(\theta)$ is estimated from the clipped mean counts.
• 2.  
  An additional contribution to the uncertainty from the intrinsic clustering of source galaxies, $\sigma _\mu ^{\rm int}(\theta)$, is estimated empirically from the variance in each annulus due to variations of the counts $N_\mu ({\boldsymbol {\theta }})$ along the azimuthal direction.
• 3.  
  Each grid cell is weighted by the fraction of its area lying within the respective annular bins. We use Monte Carlo integration to calculate the area fractions for individual cells.
• 4.  
  Masking of background galaxies due to cluster galaxies, foreground objects, and saturated pixels is properly taken into account and corrected for, by calculating at each annulus f_mask(θ) = ΔΩ_mask(θ)/ΔΩ_tot(θ), with ΔΩ_mask(θ) the area of masked regions in the annulus and ΔΩ_tot(θ) the total area of the annulus. In our analysis, we use Method B of Appendix A developed by Umetsu et al. (2011a).

The errors are combined in quadrature as $\sigma _\mu ^2 = (\sigma _\mu ^{\rm stat})^2 + (\sigma _\mu ^{\rm int})^2 + (\sigma _\mu ^{\rm sys})^2$. We note that the $\sigma _\mu ^{\rm sys}$ contribution may account for (1) strong contamination by background clusters projected near the line of sight and (2) spurious excess counts of red galaxies due perhaps to spatial variation of the photometric zero point and/or to residual flat-field errors.

Finally, we apply the correction to the number counts for the masking effects by n_μ(θ) → n_μ(θ)/[1 − f_mask(θ)] and σ_μ(θ) → σ_μ(θ)/[1 − f_mask(θ)]. Similarly, this correction is applied to the mean background counts $\overline{n}_\mu$ and its total uncertainty. The typical level of this correction is about 8% in our weak-lensing observations (see Section 5.2).

In the Bayesian framework of Umetsu et al. (2011a), the signal is described by a vector ${\boldsymbol {s}}$ of parameters containing the binned convergence profile $\lbrace \kappa _{\infty,i}\rbrace _{i=1}^N$ given by N binned κ values and the average convergence enclosed by the innermost aperture radius θ_min, κ_{∞, min} ≡ κ_∞(< θ_min), so that

$$\begin{equation} {\boldsymbol {s}}=\lbrace \kappa _{\infty, {\rm min}}, \kappa _{\infty,i}\rbrace _{i=1}^N \equiv \Sigma _{{\rm c},\infty }^{-1}{\boldsymbol {\Sigma }}, \end{equation} \tag{ 11 }$$

specified by (N + 1) parameters.

The joint likelihood function ${\cal L}({\boldsymbol {s}})$ for multi-probe lensing observations is given as a product of their separate likelihoods, ${\cal L}={\cal L}_g{\cal L}_\mu$, with ${\cal L}_{g}$ and ${\cal L}_{\mu }$ the likelihood functions for $\lbrace \langle g_{+,i}\rangle \rbrace _{i=1}^N$ and $\lbrace \langle n_{\mu,i}\rangle \rbrace _{i=1}^N$, defined as

$$\begin{eqnarray} \ln {\cal L}_{g}({\boldsymbol {s}})&=& {-} \frac{1}{2}\sum _{i=1}^{N}\frac{[\langle g_{+,i}\rangle -\hat{g}_{+,i}({\boldsymbol {s}},{\boldsymbol {c}})]^2}{\sigma ^2_{+,i}}, \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} \ln {\cal L}_\mu ({\boldsymbol {s}})&=& {-} \frac{1}{2}\sum _{i=1}^{N}\frac{[\langle n_{\mu,i}\rangle -\hat{n}_{\mu,i}({\boldsymbol {s}},{\boldsymbol {c}})]^2}{\sigma ^2_{\mu,i}}, \end{eqnarray} \tag{ 13 }$$

with ($\hat{g}_+,\hat{n}_\mu$) the theoretical predictions for the corresponding observations and ${\boldsymbol {c}}$ the calibration nuisance parameters to marginalize over,

$$\begin{equation} {\boldsymbol {c}}=\lbrace \langle W\rangle _{g}, f_{W,g}, \langle W\rangle _\mu, \overline{n}_\mu, s_{\rm eff}\rbrace. \end{equation} \tag{ 14 }$$

For each parameter of the model ${\boldsymbol {s}}$, we consider a flat uninformative prior with a lower bound of ${\boldsymbol {s}}=0$. Additionally, we account for the calibration uncertainty in the observational parameters ${\boldsymbol {c}}$.

We implement our method using a Markov Chain Monte Carlo (MCMC) approach with Metropolis–Hastings sampling, by following the prescription outlined in Umetsu et al. (2011a). The shear+magnification method has been tested (Umetsu 2013) against synthetic weak-lensing catalogs from simulations of analytical NFW lenses performed using the public package glafic (Oguri 2010). The results suggest that, when the mass-sheet degeneracy is broken, both maximum-likelihood (ML) and Bayesian marginal maximum a posteriori probability (MMAP) solutions produce reliable reconstructions with unbiased profile measurements, so that this method is not sensitive to the choice and form of priors. In the presence of a systematic bias in the background-density constraint ($\overline{n}_\mu$), the global ML estimator is less sensitive to systematic effects than MMAP and provides robust reconstructions (Section 7.4.2).

On the basis of our simulations, we thus use the global ML estimator for determination of the mass profile. In our error analysis we take into account statistical, systematic, and cosmic-noise contributions to the total covariance matrix C_ij ≡ Cov(s_i, s_j) as

$$\begin{equation} C= C^{\rm stat} + C^{\rm sys} + C^{\rm lss}, \end{equation} \tag{ 15 }$$

where C^stat is estimated from the posterior MCMC samples, C^sys accounts for systematic errors due primarily to the mass-sheet uncertainty,

$$\begin{equation} C^{\rm sys}_{ij}=({\boldsymbol {s}}_{\rm ML}-{\boldsymbol {s}}_{\rm MMAP})_i^2\delta _{ij}, \end{equation} \tag{ 16 }$$

with ${\boldsymbol {s}}_{\rm ML}$ and ${\boldsymbol {s}}_{\rm MMAP}$ the ML and MMAP solutions, respectively, and C^lss is the cosmic-noise covariance due to uncorrelated LSS projected along the line of sight (Hoekstra 2003; Hoekstra et al. 2011; Umetsu et al. 2011b). For a given depth of weak-lensing observations, the impact of cosmic noise is most important when the cluster signal itself is small (Hoekstra 2003), that is, when nearby clusters are studied, or when data at large cluster radii are examined.

The C^lss matrix is computed for a given source population as outlined in Umetsu et al. (2011b),²³ using the nonlinear matter power spectrum of Smith et al. (2003) for the Wilkinson Microwave Anisotropy Probe (WMAP) seven-year cosmology (Komatsu et al. 2011).

First, we derive an expression for the averaged ΔΣ₊ profile from stacking of tangential distortion signals around the cluster centers, following the general procedure of Umetsu et al. (2011b). For a given cluster sample, we center the shear catalogs on the respective cluster centers and construct their individual distortion profiles $\langle {\boldsymbol {g}}_+\rangle$ in physical proper length units across the range R = D_lθ = [R_min, R_max]. As we shall see below, our choice of stacking in physical length units is to reduce systematic errors in determining the ensemble-averaged cluster mass profile from stacked-lensing measurements (Okabe et al. 2013), although this choice is not optimized for the signal-to-noise ratio (S/N) of the stacked signal.

For each cluster, we express the covariance matrix (C₊)_ij of $\langle {\boldsymbol {g}}_+\rangle$ as a sum of the contributions from the shape noise ($C_+^{\rm stat}$) and the cosmic shear ($C_+^{\rm lss}$) due to uncorrelated LSS projected along the line of sight (Hoekstra 2003), $C_+ = C_+^{\rm stat}+C_+^{\rm lss}$, where $(C_+^{\rm stat})_{ij}=\sigma _{+,i}^2\delta _{ij}$ is estimated from bootstrap resampling of the background source catalog and $C_+^{\rm lss}$ is computed for a given source population (see Section 3.3.2). This cosmic noise is correlated between radial bins but can be overcome by stacking an ensemble of clusters along independent lines of sight (Umetsu et al. 2011b).

Since the noise in different clusters is uncorrelated, the tangential distortion profiles of individual clusters can be co-added according to (Umetsu et al. 2011b)

$$\begin{equation} \langle \!\langle \widehat{{\boldsymbol {\Delta \Sigma }}_+}\rangle \!\rangle = \left(\displaystyle \sum _n {{\cal W}_{+,n}} \right)^{-1} \, \left(\displaystyle \sum _n{ {{\cal W}_{+,n}} \widehat{\boldsymbol {\Delta \Sigma }}_{+,n} } \right), \end{equation} \tag{ 17 }$$

where 〈〈...〉〉 denotes the sensitivity-weighted average over all clusters in the sample, $\widehat{{\boldsymbol {\Delta \Sigma }}_{+,n}} = \Sigma _{{\rm c},n} \, \langle {\boldsymbol {g}}_{+,n}\rangle$, and ${{\cal W}_{+,n}}$ is the shear-sensitivity matrix of the nth cluster

$$\begin{equation} ({{\cal W}_{+,n}})_{ij} \equiv \Sigma _{{\rm c},n}^{-2}\,({C_{+,n}}^{-1})_{ij}\protect, \protect \end{equation} \tag{ 18 }$$

with C_{+, n} the covariance matrix of the nth $\langle {\boldsymbol {g}}_+\rangle$ profile.

The statistical covariance matrix ${\cal C}_+^{\rm stat}$ for $\langle \!\langle \widehat{\boldsymbol {\Delta \Sigma }}_+\rangle \!\rangle$ is estimated from bootstrap resampling of the cluster sample in Equation (17), which accounts for the statistical total variation of the mean mass profile averaged over the sample. Additionally, we include in our error analysis the photo-z uncertainties on the mean depth calibration, ${\cal C}^{\rm sys}$, and the net cosmic-noise covariance,

$$\begin{equation} {\cal C}_+^{\rm lss} =\left[ \displaystyle \sum _n \Sigma _{{\rm c},n}^{-2} \left({C^{\rm lss}_{+,n}} \right)^{-1} \right]^{-1}. \end{equation} \tag{ 19 }$$

Finally, the full covariance matrix for $\langle \!\langle \widehat{\boldsymbol {\Delta \Sigma }}_+\rangle \!\rangle$ is obtained as

$$\begin{equation} {\cal C}_+ = {\cal C}_+^{\rm stat} + {\cal C}_+^{\rm sys} + {\cal C}_+^{\rm lss}. \end{equation} \tag{ 20 }$$

The relation between the observable lens distortion and the lensing fields is nonlinear (Equation (3)), so that the stacked $\langle \!\langle \widehat{\Delta \Sigma _+}\rangle \!\rangle$ profile is nonlinearly related to the averaged lensing fields. Expanding the right-hand side of Equation (17) and taking the ensemble average, we obtain the next-to-leading-order correction as

$$\begin{eqnarray} \langle \!\langle \widehat{{\boldsymbol {\Delta \Sigma }}_{+}}\rangle \!\rangle &=& {\rm E}[{\boldsymbol {\Delta \Sigma }}_{+}] + \langle \!\langle \Sigma _{\rm c}^{-1} \rangle \!\rangle {\rm E}[{\boldsymbol {\Sigma }}{\boldsymbol {\Delta \Sigma }}_+]\nonumber \\ &\approx & {\rm E}[{\boldsymbol {\Delta \Sigma }}_{+}] /(1-\langle \!\langle \Sigma _{\rm c}^{-1} \rangle \!\rangle {\rm E}[{\boldsymbol {\Sigma }}]), \end{eqnarray} \tag{ 21 }$$

where E[...] denotes the ensemble average, and we have used the trace approximation $({{\cal W}_{+,n}})_{ij} \propto \delta _{ij}{\rm tr}({{\cal W}_{+,n}})$ in the first line of Equation (21); $\langle \!\langle \Sigma _{\rm c}^{-1}\rangle \!\rangle$ is defined as

$$\begin{equation} \langle \!\langle \Sigma _{\rm c}^{-1}\rangle \!\rangle = \frac{\sum _n {\rm tr}({{\cal W}_{+,n}}) \Sigma _{{\rm c},n}^{-1}}{\sum _n {\rm tr}({{\cal W}_{+,n}})}, \end{equation} \tag{ 22 }$$

where the statistical weight ${\rm tr}({{\cal W}_{+,n}})$ is proportional to $1/\Sigma _{{\rm c},n}^2$ (Equation (18)). We note that ${\rm tr}({{\cal W}_{+,n}})$ is independent of the cluster mass when stacking in physical length units (Okabe et al. 2010, 2013; Umetsu et al. 2011b; Oguri et al. 2012). On the other hand, as discussed by Okabe et al. (2013), stacking in length units scaled to r_Δ weights the contribution of each cluster to each radial bin in a nonlinear and model-dependent manner, such that ${\rm tr}{\cal W}_+ \propto r_\Delta ^2 \propto M_\Delta ^{2/3}$ when C₊ is dominated by the shape noise contribution ($C_+^{\rm stat}$).

In this work, we shall use Equations (21) and (22) to obtain a best-fit model for a given $\langle \!\langle \widehat{\boldsymbol {\Delta \Sigma }}_+\rangle \!\rangle$. The $\langle \!\langle {\boldsymbol {\Delta \Sigma }}_+\rangle \!\rangle$ is then obtained as

$$\begin{equation} \langle \!\langle {\boldsymbol {\Delta \Sigma }}_+\rangle \!\rangle \approx \langle \!\langle \widehat{{\boldsymbol {\Delta \Sigma }}_+}\rangle \!\rangle - \langle \!\langle \Sigma _{\rm c}^{-1}\rangle \!\rangle ({\boldsymbol {\Sigma }}{\boldsymbol {\Delta \Sigma }}_+), \end{equation} \tag{ 23 }$$

where the nonlinear correction $({\boldsymbol {\Sigma }}{\boldsymbol {\Delta \Sigma }}_+)$ is calculated using the best-fit solution to $\langle \!\langle \widehat{\boldsymbol {\Delta \Sigma }}_+\rangle \!\rangle$ in Equation (21).

Having obtained the mass density profiles ${\boldsymbol {\Sigma }}$ of individual clusters from combined weak-lensing shear and magnification measurements (Section 3), we can stack the clusters to produce an averaged radial mass profile.

Following Umetsu et al. (2011b), we re-evaluate the mass profiles of individual clusters in proper length units across the range R = [R_min, R_max] and construct ${\boldsymbol {\Sigma }}=\lbrace \Sigma ({<}R_{\rm min}),\Sigma (R_i)\rbrace _{i=1}^N$ on the same radial grid for all clusters. Stacking an ensemble of clusters is expressed as

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

where ${\cal W}_n$ is the nth sensitivity matrix defined as $({\cal W}_n)_{ij} \equiv \Sigma _{({\rm c},\infty)n}^{-2} \, ({C}^{-1}_n)_{ij}$ with C_n the total covariance matrix (Equation (15)) of the nth cluster.²⁴

Finally, the full covariance matrix ${\cal C}$ for the stacked $\langle \!\langle {\boldsymbol {\Sigma }}\rangle \!\rangle$ profile can be obtained in a similar manner as for $\langle \!\langle {\boldsymbol {\Delta \Sigma }}_+\rangle \!\rangle$ (Equation (20)), accounting for the profile variations in individual clusters, observational uncertainties, photo-z uncertainties on the mean-depth calibration, and the net cosmic-noise contribution.

A robust determination of the concentration parameter requires sensitive lensing measurements over a wide range of cluster radii because the concentration is sensitive to the radial curvature in the mass profile. In practice, such a measurement can be achieved by combining strong- and weak-lensing data (Merten et al. 2009), performing wide-field weak-lensing observations of nearby clusters (e.g., A2142 in Umetsu et al. 2009), or stacking the lensing signal of a statistical sample of clusters (Johnston et al. 2007; Umetsu et al. 2011b; Oguri et al. 2012).

In this work we have not attempted to determine the concentration for each cluster because the weak-lensing profiles for individual clusters are highly degenerate in M and c: the observed lensing profile can be explained by halo density profiles with larger M and smaller c than the best-fit values and vice versa. This c–M degeneracy is more significant for low-mass, high-redshift systems, which are unresolved by weak-lensing observations and for which the scale radius r_s is unconstrained by the data. In such a case, the constraints on c are essentially imposed by prior information. This inherent covariance can potentially bias the slope of the c(M) relation determined from weak lensing (Hoekstra et al. 2011).

Unlike the shear, which is sensitive to the mean interior density, the majority of the κ signal (with respect to the noise) comes from relatively inner regions and the outer profile exhibits a high degree of positive correlation (∼50% in the last few bins). The relative contribution of projected uncorrelated LSS noise increases with increasing radius (Figure 3), so that the effect of including C^lss is to further downweight the lensing signal in the outer regions especially beyond θ ∼ 10' (Hoekstra 2003). On average, we find that cosmic noise contributes ∼25% to the total error budget $\sqrt{C_{ii}}$ for the reconstructed κ profile.

Without restricting the radial range for NFW fitting, we find ∼2% lower virial masses (M_vir) relative to our baseline results obtained with a maximum fitting radius of R = 2 Mpc h⁻¹. This effect is less significant at higher overdensities, Δ_c ⩾ 200. We thus conclude that our cluster mass estimates $M_\Delta$ are statistically insensitive to the choice of the outer fitting radius.

Lensing mass measurements are sensitive to the halo triaxiality (Oguri et al. 2005). In the context of ΛCDM, prolate halo shapes are expected to develop along filaments at early stages of halo assembly, so that dynamically young cluster halos tend to have a prolate morphology (Shaw et al. 2006). Accordingly, a large fraction of clusters are expected to be elongated in the plane of the sky. On average, this will lead to an underestimation of the cluster mass when spherical symmetry is assumed (Rasia et al. 2012). Numerical simulations show that, for a sample of randomly oriented prolate clusters, their mass estimates are biased low by ≲ 5% on average when the masses are recovered from the projected mass distributions Σ under the assumption of spherical symmetry (M14).

Measuring the shear and magnification effects provides a consistency check of weak-lensing measurements, thereby allowing us to assess the robustness of our cluster mass estimates. Here we compare our mass estimates based on the ${\boldsymbol {\Sigma }}$ profiles recovered from the joint shear+magnification analysis with those obtained using the standard shear-only approach. Since background samples defined in different color regions (Section 4.4) suffer different degrees of (if any) contamination by cluster members, this comparison is sensitive to the presence of residual contamination of the background and to the shear calibration uncertainty. To do this, we adopt the Bayesian inference approach described in Section 5.3 and fit the NFW model to the tangential reduced-shear profiles $\langle {\boldsymbol {g}}_+\rangle$ in the range R ⩽ 2 Mpc h⁻¹.

In Figure 4, we show for our sample the ratio of cluster masses obtained using these two weak-lensing methods as a function of spherical radius. At each cluster radius, we compute the unweighted geometric mean and median of the shear-only to shear+magnification mass ratios. Note that we use geometric averaging, which satisfies 〈Y/X〉 = 1/〈X/Y〉 (see also Donahue et al. 2014). We see that the averaged mass ratio is consistent with unity within the errors but increases monotonically with cluster radius from 0.95 ± 0.08 at r ≃ 200 kpc h⁻¹ to 1.08 ± 0.09 at r ≃ 2 Mpc h⁻¹. This systematic radial trend can be explained by the c–M degeneracy discussed in Section 5.4.1. On the basis of this comparison, we estimate the systematic uncertainty in the overall mass calibration to be of the order of ±8%.

Zoom In Zoom Out Reset image size

The azimuthally averaged tangential distortion is a measure of the radially modulated surface mass density and is insensitive to sheet-like mass overdensities, which resemble the projected two-halo term within a couple of virial radii (Oguri & Hamana 2011). Hence, the stacked tangential-distortion signal around a large sample of clusters is a sensitive probe of the cluster-only (one-halo term) mass distribution.

In Figure 5 we show the stacked tangential-shear profile $\langle \!\langle {\boldsymbol {\Delta \Sigma }}_+\rangle \!\rangle$ derived for our sample where individual clusters and background galaxies are weighted by the shear-sensitivity kernel (${\rm tr}{\cal W}_+$ in Section 3.4.1). The individual profiles are co-added in physical length units across the range R = [R_min, R_max] = [200, 3500] kpc h⁻¹, in 11 log-spaced bins. Here, the radial limits [R_min, R_max] of our stacking analysis represent approximately the respective median values of the radial boundaries [θ_min, θ_max] covered by the data for our clusters at 0.19 ≲ z_l ≲ 0.69. For individual clusters, we impose their respective radial cuts [θ_min, θ_max] on the background samples, to be consistent with our individual cluster analysis. For our sample, we find a sensitivity-weighted average redshift of 〈〈z_l〉〉 = 0.345, in close agreement with the median redshift of $\overline{z}_{\rm l}=0.352$. The effective lensing sensitivity $\langle \!\langle \Sigma _{\rm c}^{-1}\rangle \!\rangle$ (Equation (22)) is $1/\langle \!\langle \Sigma _{\rm c}^{-1}\rangle \!\rangle \simeq 3.88\times 10^{15}\, h\ M_\odot$ Mpc⁻².

Zoom In Zoom Out Reset image size

We detect the stacked-lensing signal at a total S/N of ≃ 25 using the full covariance matrix ${\cal C}_+$ (Equation (20)) to take into account intrinsic ellipticity and projected uncorrelated LSS noise, photo-z uncertainties in the mean-depth calibration, and profile variations in individual clusters. The 45° rotated × component $\langle \!\langle {\boldsymbol {\Delta \Sigma }}_\times \rangle \!\rangle$ is consistent with a null signal within 2σ at all radii, with a total S/N of ≃ 2.8, indicating that residual systematic errors are at least an order of magnitude smaller than the measured lensing signal.

Here we quantify and characterize the ensemble-averaged mass distribution of our cluster sample using the stacked tangential-distortion signal. We examine the following five models for the halo mass density profile, ρ(r) = dM_3D(< r)/dr/(4πr²), each described by N_p parameters.

• 1.  
  Singular isothermal sphere (SIS) model with N_p = 1:

  $$\begin{equation} \rho _{\rm SIS}(r)=\frac{\sigma _v^2}{2\pi Gr^2}\protect, \protect \end{equation} \tag{ 28 }$$

  with σ_v the one-dimensional isothermal velocity dispersion.
• 2.  
  Isothermal-β model with N_p = 2 (Hattori et al. 1999):

  $$\begin{equation} \rho _{\rm iso}(r)=\frac{M_{\rm c}}{2\pi r_{\rm c}^3}\frac{3+(r/r_{\rm c})^2}{[1+(r/r_{\rm c})^2]^2}\protect, \protect \end{equation} \tag{ 29 }$$

  with M_c = M_3D(< r_c) the total mass enclosed within the core radius r_c. Note ρ_iso(r)∝r⁻² at r ≫ r_c.
• 3.  
  NFW model with N_p = 2: ρ_NFW(r) by Equation (25).
• 4.  
  Baltz–Marshall–Oguri truncated-NFW model with N_p = 2 (Baltz et al. 2009, BMO):

  $$\begin{equation} \rho _{\rm BMO}(r)=\frac{\rho _{\rm s}}{(r/r_{\rm s})(1+r/r_{\rm s})^2}\left[\frac{1}{1+(r/r_{\rm t})^2}\right]^2\protect, \protect \end{equation} \tag{ 30 }$$

  with r_t = 2.6r_vir(∼ 3r_200c) the truncation radius, where the ratio of the truncation to virial radius, τ_v ≡ r_t/r_vir, is fixed to the typical value for cluster-sized halos in the ΛCDM cosmology (Oguri & Hamana 2011).
• 5.  
  Einasto model with N_p = 3 (Einasto 1965):

  $$\begin{equation} \rho _{\rm E}(r)=\rho _{\rm s}\exp \left[-\frac{2}{\alpha _{\rm E}} \left(\frac{r}{r_{\rm s}}\right)^{\alpha _{\rm E}}\right]\protect, \protect \end{equation} \tag{ 31 }$$

  with α_E the shape parameter describing the degree of curvature. An Einasto profile with $\alpha _{\rm E}\hspace{0.5pt}\approx\hspace{0.5pt}0.18$ closely resembles the NFW profile over roughly two decades in radius (Ludlow et al. 2013). The logarithmic density gradient equals −2 at r = r_s.

The NFW, BMO, and Einasto density profiles represent a family of phenomenological models for DM halos motivated by simulations and observations. The SIS profile with ρ(r)∝r⁻² is often adopted as a lens model for its simplicity. The isothermal-β model describes the total gravitating mass profile for isothermal gas with a β density profile (Cavaliere & Fusco-Femiano 1978).³³

For the NFW, BMO, and Einasto models, we define the halo concentration parameter by c_200c = r_200c/r_s. We specify the NFW and BMO models with (M_200c, c_200c), the Einasto model with (M₂₀₀, c₂₀₀, α_E), the isothermal-β model with (M_200c, r_c), and the SIS model with M_200c. For all cases, we can ignore the two-halo contribution to ΔΣ₊, which is estimated to be γ₊ ≲ 10⁻³ within our radial range R ≲ 2r_vir and is an order of magnitude smaller than the observed lensing signal (see Figure 5).

We constrain the model parameters with our $\langle \!\langle \widehat{{\boldsymbol {\Delta \Sigma }}_{+}}\rangle \!\rangle$ profile and its full covariance matrix ${\cal C}_+$ (Section 3.4.1). We use Equation (21) to make predictions for$\langle \!\langle \widehat{{\boldsymbol {\Delta \Sigma }}_{+}}\rangle \!\rangle$. The χ² minimization is performed using the minuit minimization package from the CERN program libraries. The best-fit parameters are reported in Table 7, along with the reduced χ² and corresponding probability-to-exceed (PTE) values. The NFW, BMO, and Einasto models provide excellent fits to the data with PTE values of 0.66, 0.58, and 0.51, respectively. The isothermal β model yields a poorer but acceptable fit with a PTE of 0.33, while the SIS model gives an unacceptable fit of PTE = 4.7 × 10⁻³. The SIS model is disfavored at 4.3σ significance from a likelihood-ratio test based on the differenced χ² values $\Delta \chi ^2 \equiv \chi ^2_{\rm SIS, min}-\chi ^2_{\rm NFW,min} \simeq 18.6$ relative to the NFW model. If the truncation radius of the BMO model is allowed to vary, we obtain τ_v = (2 ± 106) × 10² (PTE = 0.56), which, however, makes the model essentially equivalent to the NFW profile.

Table 7. Best-fit Models for the Stacked Distortion Profile of the X-Ray-selected Subsample

Model	M200c	c200c	Structural Parameter	$\chi ^2_{\rm min}/{\rm dof}$a	PTEb
	($10^{14}\ M_{\odot }\,h_{70}^{-1}$)
SIS	$10.4^{+0.6}_{-0.6}$	⋅⋅⋅	⋅⋅⋅	25.4/10	0.00
Isothermal β	$14.2^{+1.2}_{-1.2}$	⋅⋅⋅	$r_{200{\rm c}}/r_c = 16.5^{+2.4}_{-1.6}$	10.3/9	0.33
NFW	$13.4^{+1.0}_{-0.9}$	$4.01^{+0.35}_{-0.32}$	⋅⋅⋅	6.8/9	0.66
NFW+pm	$13.1^{+1.2}_{-1.1}$	$3.66^{+0.68}_{-0.50}$	Mp = (11 ± 33) × 1012 M☉	6.7/8	0.57
BMO	$13.1^{+0.9}_{-0.9}$	$3.73^{+0.33}_{-0.31}$	⋅⋅⋅	7.6/9	0.58
Einasto	$13.2^{+1.0}_{-1.0}$	$3.73^{+0.43}_{-0.52}$	$\alpha _E= 0.191^{+0.071}_{-0.068}$	7.3/8	0.51

Notes. ^aMinimum χ² per degrees of freedom (dof). ^bProbability to exceed (PTE) the given $\chi ^2_{\rm min}/{\rm dof}$ based on the standard χ² probability distribution function.

Download table as:  ASCIITypeset image

In what follows, we will focus on our best models among those studied here, namely, the NFW, Einasto, and BMO density profiles. These models constrain c_200c in the range (1σ) 3.2 ≲ c_200c ≲ 4.4 ($c_{200{\rm c}}=4.01^{+0.35}_{-0.32}$ for NFW), for our 16 X-ray-selected clusters with $M_{200{\rm c}}=(1.3\pm 0.1)\times 10^{15}\ M_\odot \,h_{70}^{-1}$. For the NFW model, we find r_s ≃ 360 kpc h⁻¹, so that our data cover the range 0.6 ≲ R/r_s ≲ 10.

To examine the impact of the uncertainty in profile shape, we compare in Figure 6 the spherical mass profiles $M_{\Delta _{\rm c}}$ for the best-fit NFW, BMO, and Einasto models at several overdensities. These profiles are nearly identical at Δ_c ⩾ 200. At Δ_c = 100, the BMO-to-NFW and Einasto-to-NFW mass ratios are 0.95  ±  0.10 and 1.03 ± 0.14, respectively, both consistent with unity.

Zoom In Zoom Out Reset image size

Also shown in Figure 6 is the composite profile 〈〈M_Δ〉〉 for our sample constructed from a weighted average of the individual mass estimates (Section 5.3) constrained by the combined shear+magnification measurements (Figures 2 and 3). Specifically, 〈〈M_Δ〉〉 is defined by

$$\begin{equation} \langle \!\langle M_{\Delta }\rangle \!\rangle = \frac{\sum _{n=1}^{N_{\rm cl}} {\rm tr}({\cal W}_{+,n})\, M_{\Delta,n}}{\sum _{n=1}^{N_{\rm cl}} {\rm tr}({\cal W}_{+,n})} \end{equation} \tag{ 32 }$$

(N_cl = 16), using the shear-sensitivity kernel ${\rm tr}({\cal W}_{+})$. At Δ_c = 200, $\langle \!\langle M_{200{\rm c}}\rangle \!\rangle\hspace{0.5pt}=\hspace{0.5pt}(1.32\pm 0.08)\times 10^{15}\hspace{0.5pt} M_\odot \,h_{70}^{-1}$, in excellent agreement with the best-fit NFW halo mass, $M_{200{\rm c}}=1.34^{+0.10}_{-0.09}\times 10^{15}\ M_\odot \,h_{70}^{-1}$, from the stacked-shear-only analysis (Table 7). These comparisons show that our results are robust against different ensemble-averaging techniques and data combinations, and that the uncertainty in profile shape does not significantly affect our cluster mass measurement within the virial radius. Similarly, the composite NFW prediction for $\langle \!\langle {\boldsymbol {\Delta \Sigma }}_+\rangle \!\rangle$ is shown in Figure 5, demonstrating consistency between the shear and magnification measurements.

On the other hand, the effects of baryonic physics can in principle affect the mass density profile in the central high-density region of the cluster. Here we turn to assess the possible impact of adiabatic contraction on the total measured mass profile (Gnedin et al. 2004; Okabe et al. 2013) by introducing a central point mass into our modeling. The results for the combined NFW and point-mass (NFW+pm) model are listed in Table 7. For the central point-mass component, we obtain $M_{\rm p}=(11\pm 33)\times 10^{12}\ M_\odot \,h_{70}^{-1}$, or a 1σ upper limit of $M_{\rm p} < 44\times 10^{12}\ M_\odot \,h_{70}^{-1}$ within R_min = 200 kpc h⁻¹. Including this additional degree of freedom to account for the central unresolved mass component does not significantly change the best-fit mass and concentration parameters for our data (Table 7); however, it reduces the lower bound on c₂₀₀, allowing for somewhat lower concentrations (c_200c ≳ 3.2 at 1σ).

Now we compare our estimates of M_vir obtained by our uniform analysis with those from our previous work on individual CLASH clusters, namely, A2261 (Coe et al. 2012), MACS J1206.2−0847 (Umetsu et al. 2012), and MACS J0717.5+3745 (Medezinski et al. 2013). The major differences between this and previous weak-lensing analyses are summarized as follows.

• 1.  
  We have applied a shear calibration correction factor of 1/0.95 ≃ 1.05 (Section 4.3), which was not included in our previous studies of A2261 and MACS J1206.2−0847.
• 2.  
  The nonlinear correction given by Equation (7) was not included in our previous work.
• 3.  
  All mean depth estimates have been uniformly performed here using a self-consistent photo-z method as described in Section 4.4.
• 4.  
  A Bayesian inference approach has been used here to measure cluster masses at several overdensities, by limiting the fitting range to R ⩽ 2 Mpc h⁻¹. The χ² minimization was performed in our previous work to derive the best-fit virial mass and concentration parameters for the full radial range.

Besides, in earlier papers, we assumed a slightly different cosmology with $(\Omega _{\rm m},\Omega _\Lambda)=(0.3,0.7)$. As a result, these changes lead to a typical increase of ∼10% in our M_vir estimates relative to our previous work. This increase is primarily due to the inclusion of the shear calibration correction, which translates into a typical increase of ∼8% in M_vir. On the other hand, the present analysis takes into account systematic uncertainties C^sys (Equation (16)) due primarily to the residual mass-sheet degeneracy, providing more conservative error estimates for clusters with poorer magnification constraints.

The following is a detailed comparison of each cluster.

A2261. Coe et al. (2012) obtained $M_{\rm vir}=2.21^{+0.25}_{-0.23}\times 10^{15}\ M_\odot \,h_{70}^{-1}$ from an NFW fit to the mass profile from their shear+magnification data (Constraints (8) of their Table 4), with an estimated correction of $\Delta M_{\rm vir}\simeq 0.15\times 10^{15}\ M_\odot \,h_{70}^{-1}$ for projection effects due to LSS as specifically observed behind A2261. This is compared to $M_{\rm vir}=(2.58\pm 0.54)\times 10^{15}\ M_\odot \,h_{70}^{-1}$ in this work, corresponding to an increase of 17% in the best-fit M_vir, in which the revised estimates of mean depth result in a ∼5% increase in mass. The large increase in the uncertainty of M_vir is primarily caused by the inclusion of C^sys (Equation (16)) due to the residual mass-sheet uncertainty, which dominates the error budget for the four outermost Σ bins of A2261. A large contribution from C^sys is generally expected for low-z, high-mass clusters owing to their large angular extent on the sky.

MACS J1206.2-0847. Umetsu et al. (2012) found $M_{\rm vir}=1.64^{+0.49}_{-0.40}\times 10^{15}\ M_\odot \,h_{70}^{-1}$ based on their shear+magnification measurements (Method (2) of their Table 5), compared to our estimate of $M_{\rm vir}=(1.87\pm 0.46)\times 10^{15}\ M_\odot \,h_{70}^{-1}$, which represents an increase of 14%. Their primary NFW model from full-lensing constraints (their Method (7)), including a correction for the surrounding LSS, yields $M_{\rm 500c}=(1.01 \pm 0.15)\times 10^{15}\ M_\odot \,h_{70}^{-1}$ (see also Rozo et al. 2014), which agrees well with our estimate of $M_{\rm 500c}=(1.06 \pm 0.21)\times 10^{15}\ M_\odot \,h_{70}^{-1}$.

MACS J0717.5+3745. The cluster was recently studied by Medezinski et al. (2013), who derived $M_{\rm vir}=3.19^{+0.63}_{-0.54}\times 10^{15}\ M_\odot \,h_{70}^{-1}$ based on their two-dimensional shear and azimuthally averaged magnification-bias measurements (their Table 5),³⁵ assuming a composite model of an NFW halo and a constant component, where the latter accounts for the surrounding LSS. Their estimated M_vir is 11% lower than our estimate of $M_{\rm vir}=(3.53\pm 0.60)\times 10^{15}\ M_\odot \,h_{70}^{-1}$, as the correlated LSS contributions are partly absorbed into their mass-sheet parameter.

Finally, we compare our mass estimates of RX J2248.7−4431 with the results from a weak-shear analysis of Gruen et al. (2013) based on ESO/WFI data. For this cluster we use the same co-added images created by Gruen et al. (2013) but adopt substantially different approaches in our analysis (Section 4.2). They derived mass estimates of this cluster at several overdensities, $M_{\rm 200m}\hspace{0.5pt}=\hspace{0.5pt}33.1^{+9.6}_{-6.8}\times 10^{15}\hspace{0.5pt}M_\odot \,h_{70}^{-1}$, $M_{\rm 101c}=32.2^{+9.3}_{-6.6}\times 10^{15}\ M_\odot \,h_{70}^{-1} (\approx M_{\rm 100c})$, $M_{\rm 200c}=22.8^{+6.6}_{-4.7}\times 10^{15}\ M_\odot \,h_{70}^{-1}$, and $M_{\rm 500c}=12.7^{+3.7}_{-2.6}\times 10^{15}\ M_\odot \,h_{70}^{-1}$, which are all consistent within errors with ours (Table 6).

The majority of the CLASH clusters were targeted as part of the WtG project. Recently, the WtG collaboration published results of their weak-lensing shear mass measurements of 51 X-ray-luminous clusters at 0.15 ≲ z ≲ 0.7 based on deep multi-color Subaru/Suprime-Cam and CFHT/MegaPrime optical imaging (von der Linden et al. 2014; Kelly et al. 2014; Applegate et al. 2014).

Figure 8 shows a comparison of Subaru shear-based mass estimates (M_WtG) from Applegate et al. (2014) to our Subaru weak-lensing results. There are 17 clusters in common between the two studies, including 14 X-ray-selected and 3 high-magnification CLASH clusters. To make this comparison, we measure the mass (M_CLASH) within a fixed physical radius of r = 1.5 Mpc $h_{70}^{-1}$ assuming the spherical NFW model and using the same cosmology (h, Ω_m, Ω_Λ) = (0.7, 0.3, 0.7) adopted by Applegate et al. (2014). Applegate et al. (2014) derived cluster masses from NFW fits to their tangential reduced-shear data over the radial range $R\hspace{0.5pt}=\hspace{0.5pt}0.75$–3 Mpc $h_{70}^{-1}$, whereas our mass measurements are based on the lensing convergence (κ) data over the range R ⩽ 2 Mpc h⁻¹ ≃ 2.9 Mpc $h_{70}^{-1}$, constrained by the combined shear+magnification measurements. For this comparison, we use their results based on the photo-z posterior probability distributions of individual galaxies where available; otherwise, we use those from their color-cut method. For these 17 clusters, we find a median M_WtG/M_CLASH ratio of 1.02 and an unweighted geometric mean (Section 5.4.3) of 〈M_WtG/M_CLASH〉 = 1.10 ± 0.09. For the error estimation, we have neglected the correlation between the two axes due to overlap of source galaxies used to measure the shear because our analysis includes independent constraints from the magnification effects on red galaxy counts. Considering the systematic uncertainty of ±8% in our overall mass calibration (Section 5.4.4), we find no significant offset between our and WtG mass measurements.

Zoom In Zoom Out Reset image size

Rozo et al. (2014) performed a self-consistent calibration of cluster scaling relations between X-ray, SZE, and optical observables. They used their calibrated mass–X-ray luminosity relation to predict masses (M_R14) within the overdensity radius r_500c for the CLASH clusters at z ⩽ 0.4 and MACS J1206.2−0847 at z = 0.44. For each cluster we calculate the mass (M_CLASH) within the r_500c of Rozo et al. (2014) assuming the spherical NFW model.

For 13 clusters in common with Rozo et al. (2014), we obtain a median mass ratio of 1.12 and an error-weighted geometric mean of 〈M_R14/M_CLASH〉 = 1.13 ± 0.10, corresponding to a mass offset of Δln M ≡ 〈ln M_R14〉 − 〈ln M_CLASH〉 = 0.12 ± 0.09. This offset is not significant compared to the overall calibration uncertainty of ±8% (Section 5.4.4). Excluding two obvious outliers with M_R14/M_CLASH ∼ 2 (RX J 2129.7+0005 and RX J 1532.9+3021), we find a median mass ratio of 1.11 and a weighted geometric mean ratio of 1.03  ±  0.10, which corresponds to a mass offset of Δln M = 0.03  ±  0.09.

We compare our results with ΛCDM predictions from numerical simulations. Specifically, we examine the c–M relations obtained by Duffy et al. (2008), Bhattacharya et al. (2013), De Boni et al. (2013), and M14 using the standard NFW model. Hence, we use our NFW results for a baseline comparison with the numerical simulations.

To make a quantitative comparison between theory and observations, we calculate, for a given c(M,z) relation, the sensitivity-weighted average 〈〈c_200c〉〉 of concentrations for the observed mass and redshift distribution of our cluster sample, $\lbrace M_{200{\rm c},n}\rbrace _{n=1}^{N_{\rm cl}}$ (Table 6) and $\lbrace z_{{\rm l},n}\rbrace _{n=1}^{N_{\rm cl}}$ (Table 1), by

$$\begin{equation} \langle \!\langle c_{200{\rm c}}\rangle \!\rangle = \frac{\sum _{n=1}^{N_{\rm cl}} {\rm tr}({\cal W}_{+,n})\, c_{200{\rm c}}(M_{200{\rm c},n},z_{{\rm l},n})}{\sum _{n=1}^{N_{\rm cl}} {\rm tr}({\cal W}_{+,n})}\protect, \protect \end{equation} \tag{ 33 }$$

which is marked at the average cluster mass 〈〈M_200c〉〉 (Section 6.1.2) in Figure 9. Note that the masses $\lbrace M_{200{\rm c},n}\rbrace _{n=1}^{N_{\rm cl}}$ are constrained by the combined shear+magnification data. We use Monte Carlo simulations to estimate the total uncertainty in 〈〈c_200c〉〉, taking into account both the observational uncertainties in the cluster mass measurements (Section 5.3) and the intrinsic scatter in the concentration. Here we assume a Gaussian intrinsic scatter of dispersion σ_c = 0.33c_200c as found by Bhattacharya et al. (2013).

The resulting predictions are summarized in Table 8 and shown in Figure 9. Overall, the predicted concentrations 〈〈c_200c〉〉 for relaxed halos overlap well with our range of allowed c values, 3.2 ≲ c_200c ≲ 4.4 at 1σ (NFW and NFW+pm). Our NFW results are in excellent agreement with the predicted concentration from nonradiative hydrodynamical N-body simulations in a WMAP seven-year cosmology (M14, σ₈ = 0.82). Our measurements are slightly higher than, but consistent with, the DM-only predictions of Bhattacharya et al. (2013, σ₈ = 0.8) and Duffy et al. (2008, σ₈ = 0.796). We find a discrepancy of about 1.8σ between our NFW results and the DM-only prediction of De Boni et al. (2013, σ₈ = 0.776). These differences can be partly explained by different cosmological settings, such as the adopted values of σ₈, as discussed in detail by M14. Alternatively, this discrepancy can be reconciled if the NFW+pm model (Figure 9) is adopted to account for the possible impact of unresolved adiabatic contraction (Section 6.1.2).

Using the c–M relations derived for the full population of halos (including both relaxed and unrelaxed halos), we find that the predicted concentrations are 4%–9% lower than those for the relaxed halos (Table 8). The full-sample predictions of the Bhattacharya et al. (2013) and M14 models are consistent with our results within ∼1σ.

Finally, we examine the M14 predictions for a sample of CLASH-like X-ray-regular halos, including the effects of the CLASH X-ray selection function and possible biases due to projection effects. M14 characterized a sample of halos that reproduces the observed distribution of X-ray regularity characteristics in the CLASH X-ray-selected subsample. They found that a large fraction of this sample (∼70%) is composed of relaxed halos, but it also contains a non-negligible fraction of unrelaxed halos. Matching their simulations to the individual CLASH clusters based on the X-ray morphology, M14 predict that the NFW concentrations recovered from the lensing analysis of the CLASH X-ray-selected clusters are in the range 3 ≲ c_200c ≲ 6. We find this model provides an excellent match to the observed concentration (Table 8).

For a massive cluster acting as a supercritical lens, the strong- and weak-lensing regimes contribute similar logarithmic radial coverage (Umetsu et al. 2011b). Hence, adding strong-lensing information to weak-lensing allows us to provide tighter constraints on the inner density profile. Merten et al. (2014) conducted a strong- and weak-lensing (SaWLenS hereafter; Merten et al. 2009) analysis of 19 CLASH X-ray-selected clusters, including all clusters in our X-ray-selected subsample, by combining data from 16 band HST imaging with the wide-field weak-shear data analyzed in the present study. Their work is thus complementary to our wide-field weak-lensing analysis, which combines the shear and magnification effects. Merten et al. (2014) derived a c–M relation for their clusters, finding a moderately significant trend of decreasing c_200c with increasing halo mass, which is in good agreement with ΛCDM predictions.

The square in Figure 9 represents the average concentration predicted for our X-ray-selected subsample using their M_200c masses and best-fit c(M, z) relation, demonstrating consistency between the results obtained with different lensing methods.

Our findings are in agreement with the results obtained by Okabe et al. (2013), who carried out a stacked weak-lensing analysis of 50 X-ray clusters (0.15 < z_l < 0.3) from R = 0.1 to 2.8 Mpc h⁻¹, using two-band imaging from Subaru/Suprime-Cam observations. For their full sample, Okabe et al. (2013) found a mean concentration of $c_{\rm 200c}=4.22^{+0.40}_{-0.36}$ at their effective halo mass of $M_{200{\rm c}} =(8.5\,{\pm}\, 0.6)\times 10^{14}\ M_\odot \,h_{70}^{-1}$, which exceeds the predicted concentration from numerical simulations (Duffy et al. 2008; Zhao et al. 2009; Bhattacharya et al. 2013; De Boni et al. 2013).

More recently, Covone et al. (2014) used the CFHT Lensing Survey (CFHTLenS; Heymans et al. 2012) shear catalog to measure the stacked shear signal of 1176 optically selected clusters ($0.1 \hspace{0.5pt}<\hspace{0.5pt} z_{\rm l} \hspace{0.5pt}<\hspace{0.5pt} 0.6$) from R = 0.1 to 20 Mpc h⁻¹, in six bins of optical richness corresponding to the mass range 0.7 ≲ M_200c/(10¹⁴ M_☉ $h_{70}^{-1}) \lesssim 5$ at a median redshift of $\overline{z}_{\rm l}=0.36$. The normalization of their c(M) relation is higher than but marginally consistent with the prediction by Duffy et al. (2008) and is in closer agreement with recent simulations by Bhattacharya et al. (2013), which is qualitatively consistent with our results.

Umetsu et al. (2011b) obtained an averaged total mass profile $\langle \!\langle {\boldsymbol {\Sigma }}\rangle \!\rangle$ of four similar-mass, strong-lensing clusters at 〈〈z_l〉〉 ≃ 0.32 using combined strong-lensing, weak-lensing shear and magnification measurements from high-quality HST and Subaru observations. These clusters display prominent strong-lensing features, characterized by large Einstein radii of θ_Ein ≳ 30'' for a fiducial source redshift of z_s = 2. Umetsu et al. (2011b) found that their stacked total mass profile is well described by a single NFW profile over a wide radial range, R = 40–2800 kpc h⁻¹. The mean concentration is $c_{\rm vir}=7.68^{+0.42}_{-0.40}$ at $M_{\rm vir}=2.20^{+0.16}_{-0.14}\times 10^{15} \ M_\odot \,h_{70}^{-1}$, corresponding to an Einstein radius of θ_Ein ≃ 36'' (z_s = 2), which is compared to our CLASH X-ray-selected subsample with c_vir = 5.0  ±  0.4, $M_{\rm vir}=(1.58 \pm 0.12)\times 10^{15}\ M_\odot \,h_{70}^{-1}$, and θ_Ein = 15''  ±  4'' (z_s = 2). Intriguingly, the results from these two stacking studies are in good agreement with the ΛCDM prediction for the c_vir–θ_Ein relation based on semi-analytic calculations of Oguri et al. (2012, see their Figure 10), in which clusters with larger Einstein radii are statistically more concentrated in projection.

Here we examine the robustness of our results against the inclusion of less relaxed clusters. To do this, we perform a stacked-shear-only analysis by excluding those clusters likely with a lesser degree of dynamical relaxation, as indicated by their relatively higher degree of substructure (Postman et al. 2012; Lemze et al. 2013), namely, A209, A2261, RX J2248.7−4431, MACS J0329.7−0211, RX J1347.5−1145, and MACS J0744.9+3927. For the remaining subset of 10 X-ray regular clusters, the best-fit NFW parameters are obtained as c_200c = 3.9 ± 0.4 and $M_{\rm 200c}=(1.30 \pm 0.12)\times 10^{15}\ M_\odot \,h_{70}^{-1}$, at 〈〈z_l〉〉 = 0.334. We thus find an only negligible change in the best-fit NFW parameters compared to the total uncertainties.

We note that this subset exhibits a smaller level of cluster-to-cluster variations in the tangential distortion signal, and the total uncertainties in the stacked-lensing signal are dominated by statistical noise. Accordingly, we find that the uncertainties in the derived parameters here are comparable to those found for our total X-ray-selected subsample of 16 clusters.

Our CLASH X-ray selection (T_X > 5 keV) is designed to minimize the impact of baryonic effects on the cluster mass distribution. The effects of baryonic physics can in principle impact the inner halo profile (r ≲ 0.1r_vir; Duffy et al. 2010) and thus modify the gravity-only c–M relation, especially for less massive halos (Duffy et al. 2010; Bhattacharya et al. 2013).

In the present stacked-shear-only analysis, we examined the possible impact of adiabatic contraction by introducing a central point mass into our modeling. From this we obtained a 1σ upper limit on the unresolved central mass of M_p < 44 × 10¹² M_☉ $h_{70}^{-1}$ within R_min = 200 kpc h⁻¹ (NFW+pm in Table 7). We find that this does not significantly impact the best-fit M_200c and c_200c values for our data, but allows for somewhat lower concentrations, c_200c ≳ 3.2 at 1σ. Our findings are consistent with the conclusions of Okabe et al. (2013), who obtained a tighter limit of $M_{\rm p}<17\times 10^{12}\ M_\odot \,h_{70}^{-1}$ within R_min = 80 kpc h⁻¹ from their stacked weak-shear analysis of 50 X-ray-luminous clusters at 0.15 < z < 0.3.

Since the tangential shear is a measure of the radially modulated surface mass density, which is insensitive to sheet-like structures, the stacked-shear-only analysis can provide powerful constraints on the halo structure. We find that the shape of the stacked shear profile exhibits a steepening radial trend across the radial range 200–3500 kpc h⁻¹, which is well described by the NFW, BMO (truncated-NFW), and Einasto models (Section 6.1.2), whereas the two-halo contribution to ΔΣ₊ is negligible across the radial range of our observations.

The Einasto shape parameter is constrained as $\alpha _{\rm E}=0.191^{+0.071}_{-0.068}$, which agrees fairly well with numerical simulations: α_E = 0.175 ± 0.046 (Gao et al. 2012; the best-fit for the averaged profile of their nine cluster-sized halos), α_E = 0.172 ± 0.032 (Navarro et al. 2004, the average and dispersion of their 19 dwarf- to cluster-sized halos), α_E ∼ 0.17 (Merritt et al. 2006; the median of their six cluster-sized halos), α_E = 0.21 ± 0.07 (M14; α_E = 0.24 ± 0.09 when fitted to the projected total mass density profiles). The fitting formula given by Gao et al. (2008) yields α_E ∼ 0.29 for our high-mass clusters at 〈〈z_l〉〉 ≃ 0.35, which is consistent with our results within 1.3σ. Our results are also consistent with the recent stacked weak-shear analysis of 50 X-ray-luminous clusters by Okabe et al. (2013), who obtained $\alpha _{\rm E}=0.188^{+0.062}_{-0.058}$ for their sample with $M_{200{\rm c}} =(8.5 \pm 0.6)\times 10^{14}\ M_\odot \,h_{70}^{-1}$.

Misidentification of cluster centers is a potential source of systematic errors for stacked weak-lensing measurements on small scales. George et al. (2012) examined the impact of the choice for the cluster center on the stacked weak-shear signal based on 129 X-ray-selected galaxy groups at $0\hspace{0.5pt}<\hspace{0.5pt}z\hspace{0.5pt}<\hspace{0.5pt}1$ detected in the COSMOS field. They show that the brightest or most massive galaxies near the X-ray centroids appear to best trace the center of mass of halos. Zitrin et al. (2012) analyzed the strong-lensing signature of 10,000 clusters from the Gaussian Mixture Brightest Cluster Galaxy (Hao et al. 2010) catalog, finding a small mean offset of ≃ 13 kpc h⁻¹ between the BCG and the smoothed optical light that is assumed to trace the DM in their analysis.

Johnston et al. (2007) demonstrated that the smoothing effects of miscentering on ΔΣ₊ are much larger than on Σ and produce a noticeable effect on ΔΣ₊ out to 10 times the typical positional offset from the cluster mass centroid (Johnston et al. 2007). This is not surprising because ΔΣ₊ is insensitive to flat sheet-like structures. Here our CLASH X-ray-selection criteria ensure well-defined cluster centers, reducing the smoothing effects of cluster miscentering. Assuming that the BCG–X-ray positional offset is a good proxy for the offset from the mass centroid, the smoothing effects on ΔΣ₊ vanish at R ≳ 10σ_off ∼ 110 kpc h⁻¹ (Section 5.1), which is sufficiently smaller than the innermost measurement radius, R_min = 200 kpc h⁻¹, for our stacked shear analysis.

We compare in Figure 10 the ensemble-averaged projected mass density profile $\langle \!\langle {\boldsymbol {\Sigma }}\rangle \!\rangle$ of our X-ray-selected subsample with the results of Umetsu et al. (2011b), who analyzed combined strong-lensing, weak-lensing shear and magnification measurements of four strong-lensing-selected clusters (Section 7.2.3), characterized by $c_{\rm vir}=7.68^{+0.42}_{-0.40}$ and $M_{\rm vir}=2.20^{+0.16}_{-0.14}\times 10^{15} \ M_\odot \,h_{70}^{-1}$ at 〈〈z_l〉〉 ≃ 0.32. This is translated into the halo bias factor, b_h ≃ 10.9, which is only ∼20% higher than that estimated for our X-ray-selected subsample, b_h ≃ 9.0 (Section 6.2). In fact, the two cluster samples have similar "peak heights" in the linear (primordial) density field (Tinker et al. 2010): ≃ 3.8σ for our X-ray-selected subsample and ≃ 4.1σ for the Umetsu et al. (2011b) sample.

Zoom In Zoom Out Reset image size

Figure 10 shows that the two samples have very similar outer mass profiles at R ≳ 1 Mpc h⁻¹ ∼ 0.5r_vir, which are sensitive to the underlying mass accretion rate or halo peak height (Diemer & Kravtsov 2014). At R ≲ 400 kpc h⁻¹(∼ r_s), on the other hand, we start to see systematic deviations between the two profiles, reflecting different degrees of projected mass concentration.

As shown by high-resolution cluster simulations of Gao et al. (2012), the asphericity of clusters can lead to large variations of up to a factor of three in the projected density Σ(R) at a given radius R, depending on projection, especially within R ∼ 500 kpc h⁻¹ (see their Figure 9). Such projection effects due to cluster asphericity could explain the high apparent concentration and high central surface mass density of these four strong-lensing clusters (see also Oguri et al. 2012).

In the presence of magnification, one probes the number counts at an effectively fainter limiting magnitude of m_lim + 2.5log₁₀μ. The level of magnification is on average small in the weak regime but for our innermost bin θ = [09, 12] reaches a median factor of μ ∼ 1.7, corresponding to a magnitude shift of δm ∼ 0.6. Hence, we have implicitly assumed that the power-law behavior in Equation (8) persists down to ∼0.6 mag fainter than m_lim, where the count slope may be shallower. For a given level of count depletion, an overestimation of the count slope could lead to an overestimation of the magnification, thus biasing the resulting mass profile. However, the count slope s_eff for our data flattens slowly with depth varying from s_eff ∼ 0.15 to s_eff ∼ 0.1 from a typical magnitude limit of m_lim = 25.4 to m_lim + δm (see also Umetsu et al. 2011a), so that this introduces a small correction of only 7% for the innermost bin, much smaller compared to our noisy depletion measurements with a ∼54% median uncertainty, corresponding to $54\%/\sqrt{N_{\rm cl}}\sim 14\%$ when all clusters are combined. Therefore, we conclude that the effect of this correction is subdominant with respect to the total uncertainty.

The background density parameter $\overline{n}_\mu$ for the count-depletion analysis has been estimated from the red galaxy counts in the outer region θ = [10', θ_max] (Table 5). We find from the stacked mass profile that the mean convergence at θ ⩾ 10', where we have estimated $\overline{n}_\mu$, is κ = (8 ± 4) × 10⁻³ at 〈〈z_l〉〉 ≃ 0.35. This corresponds to a depletion signal of $\delta n_\mu /\overline{n}_\mu \approx (5\langle s_{\rm eff}\rangle -2)\kappa \sim -0.01$ with the mean count slope 〈s_eff〉 ∼ 0.15 of our sample, indicating that the estimated $\overline{n}_\mu$ is biased low by 1%. This level of the signal offset, however, is smaller than the calibration uncertainties of $\sigma (\overline{n}_\mu)/\overline{n}_\mu = (2\hbox{--}8)\%$ for an individual cluster. Hence, for all clusters in our sample, the offset signal lies within the prior range considered, $[\overline{n}_\mu -\sigma (\overline{n}_\mu),\overline{n}_\mu +\sigma (\overline{n}_\mu)]$ (Section 3.3.1). In fact, we find that the ML (best-fit) estimates of $\overline{n}_\mu$, as constrained by the combined shear+magnification data, are on average (1.0 ± 0.6)% larger than the values estimated from the counts at [10', θ_max], so that the lensing signal is consistently recovered.

This analysis demonstrates that the background level determination is not critically sensitive to our calibration prior on the background density parameter $\overline{n}_\mu$, but driven by the combined shear+magnification data.³⁶