CLASH: MASS DISTRIBUTION IN AND AROUND MACS J1206.2-0847 FROM A FULL CLUSTER LENSING ANALYSIS^*

We analyze deep BVR_cI_cz' images of MACS1206 observed with the wide-field camera Suprime-Cam (34' × 27'; Miyazaki et al. 2002) at the prime focus of the 8.3 m Subaru telescope. The observations are available in the Subaru archive, SMOKA.³² The seeing FWHM in the co-added mosaic image is 101 in B (2.4 ks), 095 in V (2.2 ks), 078 in R_c (2.9 ks), 071 in I_c (3.6 ks), and 058 in z' (1.6 ks) with 020 pixel⁻¹, covering a field of approximately 36' × 34'. The limiting magnitudes are obtained as B = 26.5, V = 26.5, R_c = 26.2, I_c = 26.0, and z' = 25.0 mag for a 3σ limiting detection within a 2'' diameter aperture. The observation details of MACS1206 are listed in Table 2. Figure 1 shows Subaru BVR_cI_cz' composite color images of the cluster field, produced automatically using the publicly available Trilogy software (Coe et al. 2012).³³

Zoom In Zoom Out Reset image size

Standard reduction steps were performed using the mscred task in IRAF.³⁴ We closely follow the data reduction procedure outlined in Nonino et al. (2009) to create a co-added mosaic of Subaru Suprime-Cam images, incorporating additional reduction steps, such as automated masking of bleeding of bright saturated stars.

To obtain an accurate astrometric solution for Subaru observations, we retrieved processed MegaCam griz images from the Canada–France–Hawaii Telescope (CFHT) archive³⁵ and used MegaCam r data (Filter Number: 9601) as a wide-field reference image. A source catalog was created from the co-added MegaCam r image, using the 2MASS catalog³⁶ as an external reference catalog. The extracted r catalog has been used as a reference for the SCAMP software (Bertin 2006) to derive an astrometric solution for the Suprime-Cam images.

The photometric zero points for the co-added Suprime-Cam images were bootstrapped from a suitable set of reference stars identified in common with the calibrated MegaCam data. These zero points were refined in two independent ways: first by comparing cluster elliptical-type galaxies with the HST/ACS images, and subsequently by fitting SED (spectral energy distribution) templates with the BPZ code (Bayesian photometric redshift estimation; Benítez 2000; Benítez et al. 2004) to Subaru photometry of 1163 galaxies having measured spectroscopic redshifts from VLT/VIMOS (P. Rosati et al. 2012, in preparation). This leads to a final photometric accuracy of ∼0.01 mag in all five passbands (see also Section 3.5). The five-band BVR_cI_cz' photometry catalog was then measured using SExtractor (Bertin & Arnouts 1996) in point-spread function (PSF) matched images created by ColorPro (Coe et al. 2006), where a combination of all five bands was used as a deep detection image. The stellar PSFs were measured from a combination of 100 stars per band and modeled using IRAF routines.

For the weak-lensing shape analysis (Section 3.2), we use the I_c-band data taken in 2009 January, which have the best image quality in our data sets (in terms of the stability and coherence of the PSF anisotropy pattern, taken in fairly good seeing conditions). Two separate co-added I_c-band images, each with a total exposure time of 1.1 ks, were produced based on the imaging obtained at two different camera orientations separated by 90°, in order not to degrade the shape measurement quality.

For shape measurements of background galaxies, we use our weak-lensing analysis pipeline based on the IMCAT package (Kaiser et al. 1995, hereafter KSB), incorporating modifications and improvements outlined in Umetsu et al. (2010). Our KSB+ implementation has been applied extensively to Subaru cluster observations (e.g., Broadhurst et al. 2005a, 2008; Umetsu et al. 2007, 2009, 2010, 2011a, 2011b; Umetsu & Broadhurst 2008; Okabe & Umetsu 2008; Medezinski et al. 2010, 2011; Zitrin et al. 2011c; Coe et al. 2012).

We measure components of the complex image ellipticity, e_α = {Q₁₁ − Q₂₂, Q₁₂}/(Q₁₁ + Q₂₂), from the weighted quadrupole moments of the surface brightness $I(\mbox{\boldmath $\theta $})$ of individual objects,

$$\begin{equation} Q_{\alpha \beta } = \int \!d^2\theta \, W({\theta })\theta _{\alpha }\theta _{\beta } I({\mbox{\boldmath $\theta $}}) \ \ \ (\alpha ,\beta =1,2), \end{equation} \tag{ 8 }$$

where W(θ) is a Gaussian window function matched to the size (r_g) of the object, and the weighted object centroid is chosen as the coordinate origin, which is iteratively refined to accurately measure the object shapes.

Next, we correct observed ellipticities e_α for the PSF anisotropy using a sample of stars in the field as references. We select bright (18 ≲ I_c ≲ 22), unsaturated stellar objects identified in a branch of the object half-light radius (r_h) versus I_c diagram and measure the PSF anisotropy kernel of the KSB algorithm as a function of the object size r_g. Figure 2 shows the distributions of stellar ellipticity components (e*_α) before and after the PSF anisotropy correction. From the rest of the object catalog, we select as a weak-lensing galaxy sample those objects with ν > 10, $r_h > \overline{r_h^*} + 1.5 \sigma (r_h^*)$, and r_g > mode(r*_g), where ν is the KSB detection significance and $\overline{r_h^*}$ and σ(r*_h) are median and rms dispersion values of stellar sizes r*_h. The anisotropy-corrected ellipticities e'_α are then corrected for the isotropic smearing effect as g_α = e'_α/P_g.

Zoom In Zoom Out Reset image size

For each galaxy we assign the statistical weight

$$\begin{equation} w_{(k)} \equiv \frac{1}{\sigma _{g(k)}^2+\alpha _g^2}, \end{equation} \tag{ 9 }$$

where σ²_g(k) is the variance for the reduced shear estimate of the kth galaxy computed from 50 neighbors identified in the r_g–I_c plane and α²_g is the softening constant variance (e.g., Hamana et al. 2003; Umetsu & Broadhurst 2008; Oguri et al. 2009; Okabe et al. 2010). This weighting scheme is essential to downweight faint and small objects that have noisy shape measurements (see Figure 4 of Umetsu et al. 2010). We choose α_g = 0.4, which is a typical value of the mean rms σ_g over the background sample (see Table 3; Umetsu & Broadhurst 2008; Umetsu et al. 2009; Okabe et al. 2010).

We follow the shear calibration strategy of Umetsu et al. (2010) to improve the precision in shear recovery. This is motivated by the general tendency of KSB+ to systematically underestimate the shear signal in the presence of measurement noise (see Umetsu et al. 2010; Okura & Futamase 2012).

First, we select as a sample of shear calibrators those galaxies with ν > ν_c and P_g > 0. Here we take ν_c = 20. Note that the shear calibrator sample is a subset of the target galaxy sample. Second, we divide the calibrator r_g–I_c plane into a grid of 2 × 10 cells, each containing approximately equal numbers of calibrators, and compute a median value of P_g at each cell. Then, each object in the target sample is matched to the nearest point on the (r_g, I_c) calibration grid to obtain a filtered measurement, 〈P_g〉. Finally, we use the calibrated estimator g_α = e'_α/〈P_g〉 for the reduced shear.

We have analyzed the two I_c mosaic images separately to construct a composite galaxy shape catalog, by properly weighting and combining the calibrated distortion measurements (g_α) for galaxies in the overlapping region.

We have tested our analysis pipeline using simulated Subaru Suprime-Cam images (see Section 3.2 of Oguri et al. 2012; Massey et al. 2007). We find that we can recover the weak-lensing signal with good precision, typically, |m| ≲ 5% of the shear calibration bias, where the range of m-values shows a modest dependence of calibration accuracy on seeing conditions and PSF properties, and c ∼ 10⁻³ of the residual shear offset, which is about one order of magnitude smaller than the typical distortion signal in cluster outskirts (|g| ∼ 10⁻²). This level of performance is comparable to other similarly well-tested methods (Heymans et al. 2006; Massey et al. 2007).

A careful background selection is critical for a weak-lensing analysis so that unlensed cluster members and foreground galaxies do not dilute the true lensing signal of the background (Broadhurst et al. 2005a; Medezinski et al. 2007, 2010; Umetsu & Broadhurst 2008). This dilution effect is simply to reduce the strength of the lensing signal when averaged over a local ensemble of galaxies (by a factor of 2–5 at R ≲ 400 kpc h⁻¹; see Figure 1 of Broadhurst et al. 2005a), particularly at small cluster radius where the cluster is relatively dense, in proportion to the fraction of unlensed galaxies whose orientations are randomly distributed.

We use the background selection method of Medezinski et al. (2010) to define undiluted samples of background galaxies, which relies on empirical correlations for galaxies in color–color–magnitude space derived from the deep Subaru photometry, by reference to evolutionary tracks of galaxies (for details, see Medezinski et al. 2010; Umetsu et al. 2010), as well as to the deep photometric-redshift survey in the COSMOS field (Ilbert et al. 2009).

For MACS1206, we have a wide wavelength coverage (BVR_cI_cz') of Subaru Suprime-Cam. We therefore make use of the (B − R_c) versus (R_c − z') color–color (CC) diagram to carefully select two distinct background populations that encompass the red and blue branches of galaxies. We limit the data to z' = 24.6 mag in the reddest band, corresponding approximately to a 5σ limiting magnitude within a 2'' diameter aperture. Beyond this limit incompleteness creeps into the bluer bands, complicating color measurements, in particular of red galaxies.

To do this, we first identify in CC space an overdensity of galaxies with small projected distance <3' (≲ 1 Mpc at z_l = 0.439) from the cluster center. Then, all galaxies within this distinctive region define the green sample (see the green outlined region in Figure 3), comprising mostly the red sequence of the cluster and a blue trail of later type cluster members (Medezinski et al. 2010; Umetsu et al. 2010), showing a number density profile that is steeply rising toward the center (Figure 4, green crosses). The weak-lensing signal for this population is found to be consistent with zero at all radii (Figure 5, green crosses), indicating the reliability of our procedure. For this population of galaxies, we find a mean photometric redshift of 〈z_phot〉 ≈ 0.44 (see Section 3.5), consistent with the cluster redshift. Importantly, the green sample marks the region that contains a majority of unlensed galaxies, relative to which we select our background samples, as summarized below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the background samples, we define conservative color limits, where no evidence of dilution of the weak-lensing signal is visible, to safely avoid contamination by unlensed cluster members and foreground galaxies. The color boundaries for our blue and red background samples are shown in Figure 3. For the blue and red samples, we find a consistent, clearly rising weak-lensing signal all the way to the center of the cluster, as shown in Figure 5.

For validation purposes, we compare in CC space our color samples with a spectroscopic sample of cluster galaxies in MACS1206. Figure 3 shows that the background selection procedure established in our earlier work (Medezinski et al. 2010, 2011; Umetsu et al. 2010) successfully excludes all spectroscopically confirmed cluster members found within the projected cluster virial radius (r_vir ≈ 1.6 Mpc h⁻¹; see Section 6). The cluster members are determined from the ongoing survey with VLT/VIMOS, part of the VLT-CLASH Large Programme 186.A-0798 (P. Rosati et al. 2012, in preparation), using the algorithm of Mamon et al. (2010) in the dynamical analysis that will be presented in a forthcoming paper (A. Biviano et al. 2012, in preparation). We find that about 70% of the cluster members overlap with our CC-selected green galaxies; the rest are cluster members with bluer colors. We note that there is a statistically inevitable fraction of interlopers even in the dynamically selected cluster membership as discussed in Wojtak et al. (2007, their Table 1) and Mamon et al. (2010, their Figure 13).

As a further consistency check, we also plot in Figure 4 the galaxy surface number density as a function of radius, n(θ), for the blue and red samples. As can be seen, no clustering is observed toward the center for the background samples, which demonstrate that there is no significant contamination by cluster members in the samples. The red sample reveals a systematic decrease in their projected number density toward the cluster center, caused by the lensing magnification effect (Section 2). A more quantitative magnification analysis is given in Section 5.2.2.

To summarize, our CC-selection criteria yielded a total of N = 13, 252, 1638, and 4570 galaxies, for the red, green, and blue photometry samples, respectively (Table 3). For our weak-lensing distortion analysis, we have a subset of 8969 and 4154 galaxies in the red and blue samples (with usable I_c shape measurements), respectively (Table 4).

The lensing signal depends on the source redshift z_s through the distance ratio β(z_s) = D_ls/D_s. We thus need to estimate and correct for the respective depths 〈β〉 of the different galaxy samples, when converting the observed lensing signal into physical mass units.

For this we used BPZ (Section 3.1) to measure photometric redshifts (photo-zs) z_phot for our deep Subaru BVR_cI_cz' photometry (Section 3.1). BPZ employs a Bayesian inference where the redshift likelihood is weighted by a prior probability, which yields the probability density P(z, T|m) of a galaxy with apparent magnitude m of having certain redshift z and spectral type T. In this work we used a new library (N. Benitez 2012, in preparation) composed of 10 SED templates originally from PEGASE (Fioc & Rocca-Volmerange 1997) but recalibrated using the FIREWORKS photometry and spectroscopic redshifts from Wuyts et al. (2008) to optimize its performance. This library includes five templates for elliptical galaxies, two for spiral galaxies, and three for starburst galaxies. In our depth estimation we utilize BPZ's ODDS parameter, which measures the amount of probability enclosed within a certain interval Δz centered on the primary peak of the redshift probability density function (pdf), serving as a useful measure to quantify the reliability of photo-z estimates (Benítez 2000).³⁷ We used our VLT/VIMOS sample of 1163 galaxies with spectroscopic redshifts z_spec(≲ 1.5) to assess the performance of our photo-z estimation. From the whole sample, we find an rms scatter of σ(δ_z) ≈ 0.027 in the fractional error δ_z ≡ (z_phot − z_spec)/(1 + z_spec), with a small mean offset μ(δ_z) = −0.0021 and a 5σ outlier fraction of ≈5.5%. Using a subsample of ∼510 galaxies with 0.3 < z_spec < 0.5, we find σ(δ_z) ≈ 0.031 with ≈1.5% of outliers.

For a consistency check, we also make use of the COSMOS catalog (Ilbert et al. 2009) with robust photometry and photo-z measurements for the majority of galaxies with i' < 25 mag. For each sample, we apply the same CC selection to the COSMOS photometry and obtain the redshift distribution N(z) of field galaxies.

For each background population, we calculate weighted moments of the distance ratio β as

$$\begin{equation} \langle \beta ^n \rangle = \frac{\int \!dz\,w(z) N(z)\beta ^n(z)}{\int \!dz\,w(z)N(z)}, \end{equation} \tag{ 10 }$$

where w(z) is a weight factor, w is taken to be the Bayesian ODDS parameter for the BPZ method, and w = 1 otherwise. The sample mean redshift 〈z_s〉 is defined similarly to Equation (10). The first moment 〈β〉 represents the mean lensing depth.³⁸ It is useful to define the effective single-plane source redshift, z_{s, eff}, such that (Umetsu & Broadhurst 2008; Umetsu et al. 2009, 2010)

$$\begin{equation} \beta (z_{s,{\rm eff}})= \langle \beta \rangle. \end{equation} \tag{ 11 }$$

In Table 4 we summarize the mean depths 〈β〉 and the effective source redshifts z_{s, eff} for our background samples. For each background sample, we obtained consistent mean-depth estimates 〈β〉 (within 2%) using the BPZ- and COSMOS-based methods. In the present work, we adopt a conservative uncertainty of 5% in the mean depth for the combined blue and red sample of background galaxies, 〈β(back)〉 = 0.54 ± 0.03, which corresponds to z_{s, eff} = 1.15 ± 0.1. We marginalize over this uncertainty when fitting parameterized mass models to our weak-lensing data.

Here we briefly summarize our well-tested approach to strong-lens modeling, developed by Broadhurst et al. (2005b) and optimized further by Zitrin et al. (2009), which has previously uncovered large numbers of multiply lensed galaxies in HST images of many clusters (e.g., Broadhurst et al. 2005b; Zitrin et al. 2009, 2010, 2011a, 2011b, 2011c). In the present work, we use a new Markov Chain Monte Carlo (MCMC) implementation of the Zitrin et al. (2009) method, where also the BCG mass is allowed to vary.⁴⁰

Our flexible mass model consists of four components, namely, the BCG, cluster galaxies, a smooth DM halo, and the overall matter ellipticity (corresponding to a coherent external shear; for details, see Zitrin et al. 2009), described by seven free parameters in total.⁴¹ The basic assumption adopted is that cluster galaxy light approximately traces the DM; the latter is modeled as a smoothed version of the former (see, for details, Zitrin et al. 2009). 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.

Zitrin et al. (2012) identified 47 new multiple images of 12 distant sources (including three candidate systems; Systems 9–11 therein), in addition to the known giant arc system at z_s = 1.03 (Ebeling et al. 2009), bringing the total known for this cluster to 50 multiply lensed images of 13 sources, spanning a wide redshift range of 1 ≲ z_s ≲ 5.5, spread fairly evenly over the central region, 3'' ≲ θ ≲ 1'. Zitrin et al. (2012) used the position and redshift of 32 secure multiple images of nine systems to constrain the mass model. Following Zitrin et al. (2012), we adopt an image positional error of 2'' (≈14 in each dimension), which is a typical value in the presence of uncorrelated large-scale structure (LSS) along the line of sight (for details, see Zitrin et al. 2012; Host 2012; Jullo et al. 2010). Including the BCG mass as an additional free parameter, we find here an acceptable fit with the minimized χ² value (χ²_min) of 22.8 for 39 dof, with an image-plane reproduction error of 176. The new MCMC results are in good agreement with the results of Zitrin et al. (2012), as shown here in Figure 6, with only some minor differences at the innermost radii ≲ 2'' (∼8 kpc h⁻¹) dominated by the BCG (see Newman et al. 2009). The detailed central mass map reveals a fairly elliptical outer critical curve (see Figure 1 of Zitrin et al. 2012). For a source at z_s = 2.54, the outer critical curve encloses an area with an effective Einstein radius of θ_Ein = 28'' ± 3''; for the lower-redshift system with z_s = 1.03, the effective Einstein radius of the critical area is θ_Ein = 17'' ± 3'' (Table 1).

Zoom In Zoom Out Reset image size

To obtain meaningful radial profiles, one must carefully define the cluster center. It is often assumed that the cluster mass centroid coincides with the BCG position, whereas BCGs can be offset from the mass centroids of the corresponding DM halos (Johnston et al. 2007; Oguri & Takada 2011; Umetsu et al. 2011a, 2011b).

Here we utilize our detailed mass model of Zitrin et al. (2012), which allows us to locate the peak position of the smooth DM component, providing an independent mass centroid determination (e.g., Umetsu et al. 2010, 2011b). In this method, 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 to represent the DM distribution. The success of this simple model in describing the projected mass distributions of lensing clusters, as well as identifying many sets of multiply lensed images, assures us that the effective DM center can be determined using multiple images and the distribution of cluster member galaxies. In this context, the DM peak location is primarily sensitive to the degree of smoothing (S) and the index of the power law (q) of Zitrin et al. (2009).

We find only a small offset of ∼1'', or a projected offset distance of d_off = 4 kpc h⁻¹ at z_l = 0.439, between the BCG and the DM peak of mass, well within the uncertainties. The BCG position also coincides well with the peak of X-ray emission within 2'' in projection (Table 1). This level of cluster centering offset is fairly small as compared to those found in other high-mass clusters, say, d_off ≈ 20 kpc h⁻¹ in RX J1347-11 (Umetsu et al. 2011b), often implied by other massive bright galaxies in the vicinity of the BCG. In the present work, we thus adopt the BCG position as the cluster center and limit our analysis to radii greater than 4'' (≈16 kpc h⁻¹), which is approximately the location of the innermost strong-lensing constraint (see Section 4.1) and sufficiently large to avoid the BCG contribution. This inner radial limit corresponds roughly to 4d_off(> 2d_off), beyond which smoothing from the cluster miscentering effects on the Σ profile is sufficiently negligible (Johnston et al. 2007; Umetsu et al. 2011a; Sereno & Zitrin 2012).

The MCMC approach allows for a full parameter-space extraction of the underlying lensing signal. We construct from MCMC samples a central mass profile κ_i and its covariance matrix ${\cal C}_{ij}$ in linearly spaced radial bins, spanning from θ = 1'' to the limit of our ACS data, θ ∼ 100''. Note that multiple-image constraints are available out to a radius of ≈1' (Section 4.1), so that the mass model beyond this radius is constrained essentially by the light distribution of cluster member galaxies, and hence the constraints there are driven by the prior. We find that the mass profile is positively correlated from bin to bin, especially at radii beyond θ_Ein ≈ 28'' (z_s = 2.5). Accordingly, the ${\cal C}$ matrix is nearly singular, with very small eigenvalues associated with large-scale modes where the constraints are weaker, leading to underestimated diagonal errors at θ ≳ θ_Ein ≈ 28'' (z_s = 2.5).

Here, we use a regularization technique with a single dof to calibrate the ${\cal C}$ matrix and obtain conservative errors for strong lensing, accounting for possible systematic errors introduced by the prior assumptions in the modeling. We first perform an eigenvalue decomposition as ${\cal C}=U\Lambda U^t$, where Λ is a diagonal matrix of eigenvalues and U is a unitary matrix of eigenvectors. Then, we determine our regularization constant, the minimum eigenvalue Λ_min, by conservatively requiring that the outermost κ value, κ_min = κ(100'') ≈ 0.22, is consistent with a null detection, i.e., Λ_min = κ²_min = (0.22)². Replacing those less than Λ_min by Λ_min and restoring the ${\cal C}$ matrix with the regularized Λ yields the desired, self-calibrated ${\cal C}$ matrix. All points at ≳ 1' are then excluded from our analysis. We find that a weaker regularization with Λ_min = (0.1)² only affects the halo parameters (M_vir, c_vir) by less than 4%.

In Figure 6 we show our strong-lensing constraints on the central κ profile using the self-calibrated ${\cal C}$ matrix, where the outer radial boundary is conservatively set to θ = 53'' (≈2θ_Ein at z_s = 2; see Zitrin et al. 2012). This calibration scheme produces conservative error estimates. Overall, the level of correction applied to the ${\cal C}$ matrix increases with increasing radius. We introduce here an estimator for the total signal-to-noise ratio (S/N) for detection, integrated over the radial range considered, and quantify the significance of the reconstruction, by the following equation (Umetsu & Broadhurst 2008):

$$\begin{equation} {\rm (S/N)^2} = \displaystyle \sum _{i,j} \kappa _i {\cal C}^{-1}_{ij} \kappa _j = \mbox{\boldmath $\kappa $}^{t} {\cal C}^{-1} \mbox{\boldmath $\kappa $}. \end{equation} \tag{ 12 }$$

With the calibrated ${\cal C}$ matrix, we find a total S/N of ≈18 for our strong-lensing κ profile in the radial range θ ⩽ 53''. We check that our results are insensitive to the choice of radial binning scheme when the self-calibration technique is applied.

We have performed complementary semi-independent strong-lensing analyses (Lenstool, LensPerfect, Pixelens, SaWLens), using as input the sets (or subsets) of multiple images identified by Zitrin et al. (2012) and the same spectroscopic and photometric redshift information.

In our Lenstool analysis, we parameterize the lens mass distribution $\Sigma (\mbox{\boldmath $\theta $})$ as a multi-component model consisting of an elliptical NFW potential and truncated elliptical halos (Kassiola & Kovner 1993) for the 86 brightest cluster members. All nine of the secure image systems are included as observational constraints. Our best solution reproduces all arc systems included and the critical lines at z_s = 2.54 and 1.03 derived in Zitrin et al. (2012), with an image-plane rms of 19, very similar to the value of ∼18 obtained by Zitrin et al. (2012) and typical to parametric mass models for clusters with many multiple images (Broadhurst et al. 2005b; Halkola et al. 2006; Limousin et al. 2007; Zitrin et al. 2009).

In the Pixelens analysis we model the lens mass distribution on a circular grid of 52'' radius divided into 18 pixels. We consider 200 models with decreasing projected mass profiles (i.e., Σ(R)∝R^−α with α > 0). We use as constraints the spectroscopically confirmed Systems 1–4 (Zitrin et al. 2012) of 14 multiple images, spanning the range 35–46'' in radius. We check that adding other multiple-image systems identified in Zitrin et al. (2012) does not significantly affect the Pixelens mass reconstruction.

In the LensPerfect analysis, we assume a prior that the projected mass is densest near the center of the BCG and decreases outward. Other priors include overall smoothness and approximate azimuthal symmetry (for details, see Coe et al. 2010). All secure image systems are used in this modeling, where including the three candidate systems (9–11) does not change the results significantly.

The SaWLens method combines central strong-lensing constraints from multiple-image systems with weak-lensing distortion constraints in a non-parametric manner to reconstruct the underlying lensing potential on an adaptively refined mesh. For this cluster we use two levels of refinement, providing a 6'' pixel resolution in the strong-lensing regime covered by CLASH imaging and a ≈22'' resolution in the Subaru weak-lensing field where the background source galaxies are sparsely sampled. The field size for the reconstruction is 25' on a side. All image systems except 10 and 11 are included as strong-lensing constraints. The lens distortion measurements for the blue+red sample are used as weak-lensing constraints. The reconstruction errors are derived from 1000 bootstrap realizations of the weak-lensing background catalog and 1000 samples of the redshift uncertainties in the catalog of strong-lensing features. The number of realizations is limited by runtime constraints.

Figure 7 shows and compares the resulting projected integrated mass profiles M_2D(< θ) derived from our comprehensive strong-lensing analyses, along with our primary strong-lensing results and model-independent Einstein-radius constraints based on Zitrin et al. (2012). All these models are broadly consistent with the Einstein-radius constraints. The calibrated error bars of Zitrin et al. (2012) are roughly consistent with the spread of the semi-independent mass profiles derived here. This comparison shows clear consistency among a wide variety of analysis methods with different assumptions and systematics, which firmly supports the reliability of our strong-lensing analyses and calibration.

Zoom In Zoom Out Reset image size

Weak-lensing distortion measurements (g) can be used to reconstruct the underlying projected mass density field $\Sigma (\mbox{\boldmath $\theta $})$ (see Equation (3)). Here we use the linear map-making method outlined in Section 4.4 of Umetsu et al. (2009) to derive the projected mass distribution from the Subaru distortion data presented in Section 3.

In the left panel of Figure 8, we show the $\Sigma (\mbox{\boldmath $\theta $})$ field in the central 24' × 24' region, reconstructed from the blue+red sample (Section 3.4), where for visualization purposes the mass map is smoothed with a Gaussian with 15 FWHM. A prominent mass peak is visible in the cluster center. This first maximum in the mass map is detected at a significance level of 9.5σ and coincides well with the optical/X-ray cluster center within the statistical uncertainty: ΔR.A. = 70  ±  72, Δdecl. = −14  ±  76, where ΔR.A. and Δdecl. are right ascension and declination offsets, respectively, from the BCG center.

Zoom In Zoom Out Reset image size

Also compared in Figure 8 are member galaxy distributions in the MACS1206 field, Gaussian smoothed to the same resolution of θ_FWHM = 15. The middle and right panels display the number and (K-corrected) R_c-band luminosity density fields, respectively, of green cluster galaxies (see Table 3).

Overall, mass and light are similarly distributed in the cluster. The cluster is fairly centrally concentrated in projection and associated with elongated LSS running northwest–southeast (NW–SE), both in the projected mass and galaxy distributions. A more quantitative characterization of the 2D matter distribution around the cluster will be given in Section 6.

Now we derive azimuthally averaged lens distortion and magnification profiles from the Subaru data. We calculate the weak-lensing profiles in N discrete radial bins from the cluster center (Section 4.2), spanning the range [θ_min, θ_max] with a constant logarithmic radial spacing Δln θ = ln (θ_max/θ_min)/N, where the inner radial boundary θ_min is taken to be θ_min = 08 (>θ_Ein). The outer radial boundary θ_max is chosen to be θ_max = 16' (R_max ≈ 3.8 Mpc h⁻¹), sufficiently larger than the typical virial radius r_vir of high-mass clusters (r_vir ≈ 1.6 Mpc h⁻¹ for MACS1206; see Section 6), but sufficiently small with respect to the size of the Suprime-Cam's field of view so as to ensure accurate PSF anisotropy correction. The number of radial bins is set to N = 8, chosen such that the detection S/N (defined as in Equation (12)) is of the order of unity per pixel.

5.2.1. Lens Distortion

For each galaxy, we define the tangential distortion g₊ and the 45° rotated component, with respect to the cluster center, from linear combinations of the distortion coefficients (g₁, g₂) as g₊ = −(g₁cos 2ϕ + g₂sin 2ϕ) and g_× = −(g₂cos 2ϕ − g₁sin 2ϕ), with ϕ being the position angle of an object with respect to the cluster center. In the absence of higher-order effects, weak lensing only induces curl-free tangential distortions, while the azimuthal averaged × component is expected to vanish. In practice, the presence of × modes can be used to check for systematic errors.

For each galaxy sample, we calculate the weighted average of g₊ in a set of radial bins (i = 1, 2, ..., N) as

$$\begin{equation} g_{+,i} \equiv g_+(\theta _i)= \left[ \displaystyle \sum _{k\in i} w_{(k)}\, g_{+(k)} \right] \left[ \displaystyle \sum _{k\in i} w_{(k)}\right]^{-1}, \end{equation} \tag{ 13 }$$

where the index k runs over all objects located within the ith annulus, θ_i is the weighted center of the ith radial bin, and the weight factor w_(k) is defined by Equation (9). We use the continuous limit of the area-weighted center for θ_i (see Appendix A of Umetsu & Broadhurst 2008). We perform a bootstrap error analysis to assess the uncertainty σ_{+, i} in the tangential distortion profile g_{+, i} (Umetsu et al. 2010).

In Figure 5, we compare azimuthally averaged radial profiles of g₊ and g_× as measured from our red, blue, green, and blue+red galaxy samples (Section 3.4). For all samples, the × component is consistent with a null detection well within 2σ at all radii, indicating the reliability of our distortion analysis. The red and blue populations show a very similar form of the radial g₊ profile that declines smoothly from the cluster center. The observed tangential distortion signal is significant with a total detection S/N of 8.1 and 5.1 for the red and the blue sample, respectively, both remaining positive to the limit of our data, θ_max = 16'. The detection significance is improved to 9.3σ using a full composite sample of Subaru blue+red background galaxies (see the top panel of Figure 9).

Zoom In Zoom Out Reset image size

In Figure 5 we also compare the Subaru data with the results obtained from CFHT/MegaCam data (Section 3.1) using the same analysis pipeline as described in Section 3. For this we identified 15,875 background galaxies (n_g ≈ 4.5 galaxies arcmin⁻²) with MegaCam grz photometry using our CC background selection method (Section 3.4) and estimated a mean depth of z_{s, eff} ≈ 1.09, comparable to that of the Subaru full background sample (z_{s, eff} = 1.15 ± 0.1, n_g ≈ 13 galaxies arcmin⁻²; Section 3.5). This comparison shows excellent agreement where the data overlap, demonstrating the robustness of our analysis.

5.2.2. Magnification Bias

The relation between observable distortion (g) and underlying convergence (κ) is non-local. Hence, the mass distribution derived from distortion data alone suffers from a mass-sheet degeneracy (Section 2).

Here we construct a radial mass profile from complementary lens distortion and magnification measurements, {g_{+, i}}^N_{i = 1} and {n_{μ, i}}^N_{i = 1}, following the Bayesian prescription given by Umetsu et al. (2011b), effectively breaking the mass-sheet degeneracy. A brief summary of this Bayesian method is provided in Appendix A.1. The model is described by a vector $\mbox{\boldmath $s$}$ of parameters containing the discrete convergence profile {κ_i}^N_{i = 1} in the subcritical regime (θ_i > θ_Ein) and the average convergence within the inner radial boundary θ_min of the weak-lensing data, $\overline{\kappa }_{\rm min}\equiv \overline{\kappa }({<}\theta _{\rm min})$, so that $\mbox{\boldmath $s$}=\lbrace \overline{\kappa }_{\rm min},\kappa _i\rbrace _{i=1}^{N}$, being specified by (N + 1) parameters.

We find a consistent mass profile solution $\mbox{\boldmath $s$}$ based on a joint Bayesian fit to the observed distortion and magnification measurements, as shown in Figure 9. The detection significance has been improved from 9.3σ to 11.4σ by adding the magnification measurement, corresponding to an improvement by ∼23%, compared to the lensing distortion signal (Umetsu et al. 2011b; Coe et al. 2012).

The resulting mass profile $\mbox{\boldmath $s$}$ is shown in Figure 6, along with our primary strong-lensing model (Sections 4.1–4.3). Our independent strong- and weak-lensing mass profiles are in good agreement where they overlap, and together they form a well-defined mass profile. The outer mass profile derived from weak lensing exhibits a fairly shallow radial trend with a nearly isothermal logarithmic density slope in projection, γ_2D ≡ −dln Σ/dln R ∼ 1. Note that this flat behavior is not clearly evident in the tangential distortion profile, which is insensitive to sheet-like mass overdensities (Section 2). To constrain the cluster properties from the composite halo+LSS mass profile, this LSS contribution needs to be taken into account and corrected for. We will come back to this point in Sections 6.2 and 6.4.

Also shown in Figures 6 and 7 is a purely shear-based reconstruction using the one-dimensional (1D) method of Umetsu & Broadhurst (2008; see also Umetsu et al. 2010), based on the nonlinear extension of aperture mass densitometry (Clowe et al. 2000). Here we have adopted a zero-density boundary condition in the outermost radial bin, 16' ⩽ θ ⩽ 18'. The total S/N in the recovered mass profile is ≈9.2, which agrees well with ≈9.3 in the g₊ profile (Section 5.2.1). Our results with different combinations of lensing measurements and boundary conditions, having different systematics, are in agreement with each other. This consistency demonstrates that our results are robust and insensitive to systematic errors.

The projected cumulative mass profile M_2D(< θ) is given by integrating the density profile $\mbox{\boldmath $s$}=\lbrace \overline{\kappa }_{\rm min},\kappa _i\rbrace _{i=1}^{N}$ (see Appendices A and B of Umetsu et al. 2011b) as

$$\begin{equation} M_{\rm 2D}({<}\theta _i)=\pi (D_l\theta)^2 \Sigma _{\rm crit}\overline{\kappa }_{\rm min} +2\pi D_l^2 \Sigma _{\rm crit} \int _{\theta _{\rm min}}^{\theta _i}\! d\ln \theta \,\theta ^2\kappa (\theta). \end{equation} \tag{ 14 }$$

We compare in Figure 7 the resulting M_2D profiles derived here from a wide variety of strong- (Section 4) and weak-lensing analyses, along with the model-independent Einstein-radius constraints of M_2D(< 17'') = 5.8^+1.3_{− 1.4} × 10¹³ M_☉ h⁻¹ at θ_Ein = 17'' ± 2'' (z_s = 1.03) and M_2D(< 28'') = 1.1^+0.2_{− 0.3} × 10¹⁴ M_☉ h⁻¹ at θ_Ein = 28'' ± 3'' (z_s = 2.54).⁴² Again, we find good agreement in the regions of overlap among the results obtained from a variety of lensing analyses, ensuring consistency of our lensing analysis and methods.

Unlike the non-local distortion effect, the magnification falls off sharply with increasing distance from the cluster center. For MACS1206, we find κ ≲ 1% at radii ≳ 10', where the expected level of the depletion signal is n_μ/n₀ − 1 ≈ −2κ for a maximally depleted sample with s = 0, indicating a depletion signal of ≲ 2% in the outer region where we have estimated the unlensed background counts, n₀. This level of signal is smaller than the fractional uncertainties in estimated unlensed counts n₀ of 3% (Section 5.2.2), thus consistent with the assumption. Note that the calibration uncertainties in our observational parameters (n₀, s, ω) have been marginalized over in our Bayesian analysis (the Appendix).

In the presence of magnification, one probes the number counts at an effectively fainter limiting magnitude: m_cut + 2.5log₁₀μ(θ). The level of magnification is on average small in the weak-lensing regime but reaches μ ≈ 1.6 (at z_{s, eff} ≈ 1.1) for the innermost bin in this cluster. Hence, we have implicitly assumed in our analysis that the power-law behavior (Equation (7)) persists down to ∼0.5 mag fainter than m_cut where the count slope may be shallower. For a given level of count depletion, an underestimation of the effective count slope could lead to an underestimation of μ, thus biasing the resulting mass profile. However, the count slope for our data flattens only slowly with depth varying from s ∼ 0.13 to ∼0.05 from a limit of z' = 24.6–25.1 mag, so that this introduces a small correction of only ∼10% for the most magnified bins (μ ∼ 2). In fact, we have found a good consistency between the results with and without the magnification data.

First, we constrain the NFW model parameters $\mbox{\boldmath $p$}\equiv (M_{\rm vir},c_{\rm vir})$ by combining model-independent weak-lensing distortion, magnification, and strong-lensing Einstein-radius measurements, whose systematic errors are well understood from numerical simulations (e.g., Meneghetti et al. 2011; Rasia et al. 2012). The χ² function for the combined Einstein-radius and weak-lensing constraints is expressed as

$$\begin{equation} \chi ^2 = \chi ^2_{\rm Ein} + \chi ^2_{\rm WL}, \end{equation} \tag{ 16 }$$

where the χ²_Ein for the Einstein-radius constraints is defined by (see Umetsu & Broadhurst 2008; Umetsu et al. 2010)

$$\begin{eqnarray} &&\chi ^2_{\rm Ein} = \displaystyle \sum _{i=1}^{N_{\rm Ein}} \frac{[1-\hat{g}_{+,i}(\mbox{\boldmath $p$},z_{s,i})]^2}{\sigma _{+,i}^2}, \end{eqnarray} \tag{ 17 }$$

with NEin being the number of independent Einstein-radius constraints {θEin, i}NEini = 1 from sources with different redshifts {zs, i}NEini = 1 and $\hat{g}_{+,i}(\mbox{\boldmath $p$},z_{s,i})=\hat{g}(\theta _{{\rm Ein},i}|\mbox{\boldmath $p$},z_{s,i})$ being the NFW model prediction for the reduced tangential shear at θ = θEin, i, evaluated at the source redshift zs = zs, i. Note that the Einstein radius marks the point of maximum distortion, $g_+=(\overline{\kappa }-\kappa)/(1-\kappa)=1$, i.e., $\overline{\kappa }=1$ within θEin. The χ2 function for our full weak-lensing analysis (Section 5.3) is described by

$$\begin{equation} \chi ^2_{\rm WL} = \displaystyle \sum _{i,j} [s_{i}-\hat{s}_{i}(\mbox{\boldmath $p$},z_{s,{\rm eff}})] \left({\cal C}_{\rm WL}\right)_{ij}^{-1} [s_{j}-\hat{s}_{j}(\mbox{\boldmath $p$},z_{s,{\rm eff}})], \end{equation} \tag{ 18 }$$

where $\mbox{\boldmath $s$}=\lbrace \overline{\kappa }_{\rm min},\kappa _i\rbrace _{i=1}^{N}$ is the mass profile reconstructed from the combined lens distortion and magnification measurements, $\hat{\mbox{\boldmath $s$}}(\mbox{\boldmath $p$},z_{s,{\rm eff}})$ is the NFW model prediction for $\mbox{\boldmath $s$}$, and ${\cal C}_{\rm WL}$ is the full covariance matrix of $\mbox{\boldmath $s$}$ defined as

$$\begin{equation} {\cal C}_{\rm WL} = {\cal C} + {\cal C}^{\rm lss}, \end{equation} \tag{ 19 }$$

with ${\cal C}$ being responsible for statistical measurement errors (Appendix A.1) and ${\cal C}^{\rm lss}$ being the cosmic covariance matrix responsible for the effect of uncorrelated LSS along the line of sight (Hoekstra 2003; Hoekstra et al. 2011; Umetsu et al. 2011a; Oguri & Takada 2011).⁴³ In all modeling below, the effective source redshift z_{s, eff} = 1.15 ± 0.1 of our full background sample is treated as a nuisance parameter, and its uncertainty is marginalized over. In order to evaluate ${\cal C}^{\rm lss}$, we assume the concordance ΛCDM cosmological model of Komatsu et al. (2011) and use the fitting formula of Peacock & Dodds (1996) to compute the nonlinear matter power spectrum. We project the matter spectrum out to an effective source redshift of z_{s, eff} = 1.15 to calculate ${\cal C}^{\rm lss}$ for weak-lensing observations. For details, see Umetsu et al. (2011a). For Einstein-radius measurements, we conservatively assume an rms displacement of 2'' due to uncorrelated LSS, as predicted by recent theoretical work (∼2'' for a distant source at z_s ∼ 2.5; see Host 2012; Jullo et al. 2010). This is combined in quadrature with the measurement error in θ_Ein (Table 1) to estimate a total uncertainty σ_{+, i}.⁴⁴

For strong lensing, we use double Einstein-radius constraints (N_Ein = 2) from the multiple-image systems at z_s = 1.03 and z_s = 2.54 (Table 1). For weak lensing, the cluster mass profile $\mbox{\boldmath $s$}$ is measured in N + 1 = 9 bins. Hence, we have a total of 11 constraints.

The resulting constraints on the NFW model parameters are summarized in Table 5.

6.1.1. Weak-lensing Constraints

6.1.2. Combining Einstein-radius Constraints with Weak Lensing

The presence of surrounding LSS in MACS1206 has a non-negligible impact on the determination of cluster mass profile especially at large radii (Sections 5.3 and 6.1). It is therefore necessary to assess and correct for their effects on the projected mass profile. Here we use two different methods to quantify the ellipticity and orientation of the projected mass distribution in and around the cluster.

First, following the prescription given by Umetsu et al. (2009), we introduce mass-weighted quadrupole shape moments around the cluster center, in analogy to Equation (8), defined as

$$\begin{eqnarray} &&Q_{\alpha \beta } = \int _{\Delta \theta \le \theta _{\rm max}} \!d^2\theta \, \Delta \theta _\alpha \Delta \theta _\beta \, \Sigma (\mbox{\boldmath $\theta $}) \ \ \ (\alpha ,\beta =1,2), \end{eqnarray} \tag{ 20 }$$

where θ_max is the circular aperture radius and Δθ_α is the angular displacement vector from the cluster center. We construct with {Q_αβ} a spin-2 ellipticity measure $e_\Sigma =|e_\Sigma |e^{2i\phi _e}$, where the ellipticity is defined such that, for an ellipse with major and minor axes a and b, it reduces to $|e_\Sigma |=1-b/a$ and ϕ_e is the position angle of the major axis (Bertin & Arnouts 1996), measured north of west here. Similarly, the spin-2 ellipticity for the cluster galaxies is defined using the surface number and R_c-band luminosity density fields of CC-selected cluster galaxies (Section 5.1), $\Sigma _n(\mbox{\boldmath $\theta $})$ and $\Sigma _l(\mbox{\boldmath $\theta $})$. We calculate weighted moments using only those pixels above the 2σ threshold with respect to the background level (estimated with the biweight scale and location; see Beers et al. 1990). Practical shape measurements are done using pixelized maps shown in Figure 8.

Next, we constrain the ellipticity and orientation of the projected mass distribution by directly fitting a 2D shear map with a single elliptical lens model. Here, we closely follow the prescription given by Oguri et al. (2010) to construct an elliptical NFW (eNFW, hereafter) model (see also Oguri et al. 2012), by introducing the mass ellipticity $|e_\Sigma |=1-b/a$ in the isodensity contours of the projected NFW profile Σ(R) as $R^2\rightarrow X^2 (1-|e_\Sigma |) + Y^2 /(1-|e_\Sigma |)$ (for details, see Oguri et al. 2010).⁴⁶ The model shear field is computed by solving the 2D Poisson equation (Keeton 2001). We then construct from Subaru data a lens distortion map $(g_1(\mbox{\boldmath $\theta $}),g_2(\mbox{\boldmath $\theta $}))$ and its covariance matrix ${\cal C}_g$ (Equation (A5)) on a 2D Cartesian grid with 1' spacing, centered at the BCG. We exclude from our analysis the five innermost cells lying in the central region, θ < 1', to avoid systematic errors (see Appendix A.2). The halo centroid is fixed to the BCG position. Accordingly, the eNFW model is specified by four model parameters, $\mbox{\boldmath $p$}=(M_{\rm vir},c_{\rm vir},|e_\Sigma |,\phi _e)$. The constraints on individual parameters are obtained by projecting the 2D shear likelihood function (Equation (A6) in Appendix A.2) to the parameter space (or, minimizing χ²).

In Table 6, we summarize our cluster ellipticity and orientation measurements. In this analysis, we are mainly interested in the orientation of the ellipticity, in order to correct for the effects of LSS along the axis of elongation. An overall agreement is found between the shapes of mass, light, and galaxy distributions in MACS1206, especially in terms of orientation (Figure 8), within large uncertainties (Table 6). The mass distribution in and around the cluster is aligned well with the luminous galaxies in the green sample, composed mostly of cluster member galaxies (Section 3.4). For all cases, the position angle ϕ_e of the major axis is found to be fairly constant with radius θ_max and lies in the range 15° ≲ ϕ_e ≲ 30°.

Table 6. Ellipticity and Position Angle Measurements

Method	θmaxa	Ellipticityb	PAc
	(')		($\deg$)
BCG	10''	0.53+0.03− 0.03	15.0+2.3− 2.3
Chandra X-ray	15	0.30+0.03− 0.03	21.9+1.7− 1.7
Galaxy density	8'	0.53+0.04− 0.04	15.7+1.3− 5.9
Galaxy light	8'	0.41+0.06− 0.06	19.0+5.9− 5.4
WL mass map	8'	0.37+0.13− 0.13	19.4+8.5− 17.7
WL 2D shear fit	8'	0.68+0.18− 0.28	28.6+5.8− 7.9

Notes. ^aCircular aperture radius. ^bEllipticity modulus defined such that, for an ellipse with major and minor axes a and b, it reduces to 1 − b/a. ^cPosition angle of the major axis measured north of west.

Download table as:  ASCIITypeset image

In the central region, we find a projected mass ellipticity of |e_Σ| ∼ 0.3 and a position angle of ϕ_e ∼ 14° from the Pixelens analysis; we obtain consistent values for both $|e_\Sigma |$ and ϕ_e from a different strong-lensing analysis (C. Grillo et al. 2012, in preparation) using only System 7 of Zitrin et al. (2012). A similar value is found for the projected mass ellipticity of $|e_\Sigma |=0.26\pm 0.16$ (ϕ_e ∼ 19°) at θ_max = 4' using the weak-lensing Σ map. From an elliptical King model fit to Chandra X-ray data (Figure 10; for details, see Section 7.4), we find an ellipticity of 0.30 ± 0.03 (a/b ≈ 1.5) and ϕ_e = 219 ± 17 at θ_max = 15.

Zoom In Zoom Out Reset image size

On the other hand, we obtained higher values of ellipticity on large angular scales beyond the cluster virial radius, θ_vir ≡ r_vir/D_l ∼ 7'. We find $|e_\Sigma |\sim 0.4$–0.5 at θ_max = 8' using the pixelized cluster mass, galaxy, and light distributions. From the 2D shear fitting to a single eNFW model, the projected mass ellipticity is constrained in the range $|e_\Sigma | =0.68^{+0.18}_{-0.28}$ ($|e_\Sigma |\gtrsim 0.4$ or a/b ≳ 1.7 at 1σ) at θ_max = 8'. This apparent increase in ellipticity with radius could be partly explained by the additional contribution from the surrounding LSS that is extended along the cluster major axis. Note that the observed tendency for the shear-based method to yield higher ellipticity estimates, compared to the mass-map-based method, could be due to the non-local nature of the shear field, in conjunction with our single-component assumption in the 2D shear fitting analysis. Overall, this level of ellipticity is consistent within large errors with the mean cluster ellipticity $\langle |e_\Sigma |\rangle =0.46{\,\pm\,} 0.04$ obtained by Oguri et al. (2010) from a 2D weak-lensing analysis of 25 X-ray-luminous clusters.

In what follows, we fix the position angle of the NW–SE cluster-LSS major axis to a reference value of ϕ_e = 20°, which is close to the values derived from the Chandra X-ray data, Σ_l and Σ maps. We note that, in principle, the X-ray structure in a triaxial system is expected to be tilted with respect to the total matter in projection, even in the absence of intrinsic misalignments (see Romanowsky & Kochanek 1998). In the present work, we define the NW and SE excess regions, respectively, as NW and SE outer cone regions with θ > 4' centered on the cluster center, with opening angle 90° and position angle ϕ_e = 20°, as defined by the pair of gray solid lines in each panel of Figure 8.

We have also obtained CLASH constraints on the mean BCG ellipticity and position angle derived from the ACS F814W image. For this we performed a detailed structural analysis on the BCG using the snuc task in the XVista software package.⁴⁷ In Figure 11 we show the ACS F814W image, best-fit model, and image residuals after subtraction of the model. No systematic deviations are seen in the residuals between the data and the model, suggesting that the BCG has not undergone any major merger recently. The radial profiles of ellipticity and position angle were measured in several independent radial bins (02 ≲ θ ≲ 10''), and their respective (sensitivity weighted) mean values were obtained as $\langle e_\Sigma \rangle =0.53\pm 0.03$ and $\phi _e=(15.0\pm 2.3)\deg$ (Table 6). Consistent results were found in several other HST bands (ACS F475W to F814W and WFC3 F105W to F160W). The mean BCG ellipticity is found to lie in the range 0.46–0.53 with a small scatter of 0.02 across the ACS and WFC3 bands. The BCG position angle is constrained to be $\phi _e=(15.2\pm 0.4)\deg$, which is in excellent agreement especially with that derived independently from the large-scale distribution Σ_n of galaxies.

Zoom In Zoom Out Reset image size

In this subsection, we look into the azimuthal dependence of the radial projected mass distribution, Σ(R, ϕ), to assess and correct for the effect of surrounding LSS on the cluster mass profile measurement. Because of the non-local nature and inherent insensitivity to sheet-like overdensities of the shear field, it is essential to use the combination of lens magnification and distortion to reconstruct the projected cluster mass distribution embedded in LSS. For this purpose, we extend the 1D Bayesian method of Umetsu et al. (2011b) into a 2D mass distribution by combining the 2D shear pattern $g(\mbox{\boldmath $\theta $})$ with the azimuthally averaged magnification measurements n_μ(θ). In the 2D analysis, our model $\mbox{\boldmath $s$}$ is a vector of parameters containing a set of discrete mass elements on a grid of N_cell independent cells, $\mbox{\boldmath $s$}=\lbrace \kappa _m\rbrace _{m=1}^{N_{\rm cell}}$. A brief summary of this 2D method is given in Appendix A.2. The details of the method will be presented in our forthcoming paper (K. Umetsu et al. 2012, in preparation).

By combining Subaru distortion and magnification data, we construct here a mass map over a 30 × 30 grid with 08 spacing, covering a 24' × 24' field around the cluster (N_cell = 900). We have 2 × 896 distortion constraints $\lbrace g_1(\mbox{\boldmath $\theta $}_m)\rbrace _{m=1}^{N_{\rm cell}}$ and $\lbrace g_2(\mbox{\boldmath $\theta $}_m)\rbrace _{m=1}^{N_{\rm cell}}$ over the mass grid, excluding the four innermost cells lying in the cluster central region (see Appendix A.2), and N = 8 radial magnification constraints {n_μ(θ_i)}^N_{i = 1}. Hence, we have a total of 1800 constraints (900 dof). Additionally, we marginalize over the calibration uncertainties in the observational parameters (n₀, s, ω; Section 3.5). The best solution $\mbox{\boldmath $s$}$ has been obtained with χ²_min/dof = 1058/900. We then follow Umetsu & Broadhurst (2008) to calculate the radial mass distribution 〈Σ(R)〉 and its covariance matrix from a weighted projection of the κ map, where we conservatively limit our 2D analysis to radii smaller than θ = 12' (R ≈ 2.9 Mpc h⁻¹). We check that the azimuthally averaged radial mass profile constructed from the κ map reproduces our corresponding 1D results (Section 5.3).

We show in Figure 6 the radial mass distribution obtained excluding the NW and SE excess regions (defined in Section 6.4; see also Figure 8). This weak-lensing mass profile, corrected for the effect of surrounding LSS, exhibits a steeper radial trend than that averaged over all azimuthal angles. We note that a slight remaining excess is seen at θ ≳ 5' (R ≳ 1.2 Mpc h⁻¹). By fitting the "LSS-corrected" mass profile with an NFW profile, we find a higher concentration c_vir = 7.5^+2.5_{− 1.8} with M_vir = 1.15^+0.25_{− 0.20} × 10¹⁵ M_☉ h⁻¹ (χ²_min/dof = 10.6/6). This model predicts an Einstein radius of θ_Ein ≈ 32'' for z_s = 2.5, comparable to the observed value, θ_Ein = 28'' ± 3''.

As shown in Figures 6 and 7, our weak- and strong-lensing data agree well in their region of overlap. Here we further improve the statistical constraints on the halo parameters $\mbox{\boldmath $p$}=(M_{\rm vir},c_{\rm vir},\alpha)$ by combining the joint weak-lensing distortion and magnification constraints $\chi ^2_{\rm WL}(\mbox{\boldmath $p$},z_{s,{\rm eff}})$ (Section 6.1) with the inner mass profile κ_i based on the detailed strong-lensing analysis of Zitrin et al. (2012).

We write the combined χ² function of our full-lensing constraints as

$$\begin{equation} \chi ^2=\chi ^2_{\rm WL} + \chi ^2_{\rm SL} \end{equation} \tag{ 21 }$$

with χ²_SL for strong lensing being defined as

$$\begin{eqnarray} &&\chi ^2_{\rm SL} = \displaystyle \sum _{i,j} [\kappa _{i}-\hat{\kappa }_{i}(\mbox{\boldmath $p$})] \left({\cal C}_{\rm SL}\right)^{-1}_{ij} [\kappa _{j}-\hat{\kappa }_{j}(\mbox{\boldmath $p$})], \end{eqnarray} \tag{ 22 }$$

where κ_i is defined in 25 discrete bins over the radial range [4'', 53''] (see Section 4) and scaled to a fiducial depth z_s = 2.54 of the strong-lensing observations, matched to the spectroscopically confirmed five-image system (System 4 of Zitrin et al. 2012); $\hat{\kappa }_i$ is the theoretical prediction for κ_i; and ${\cal C}_{\rm SL}={\cal C}+{\cal C}_{\rm lss}$ is the bin-to-bin covariance matrix for the discrete κ profile, with ${\cal C}$ being the self-calibrated covariance matrix derived in Section 4.3 and ${\cal C}_{\rm lss}$ being the cosmic noise contribution. We use a consistent single source plane at z_s = 2.54 to evaluate ${\cal C}^{\rm lss}$.

The resulting NFW and gNFW fits are summarized in Tables 5 and 7, respectively. For both models, we show the respective fits derived with and without the LSS correction for the outer weak-lensing profile (R ≳ 1 Mpc h⁻¹). We find that, when the detailed strong-lensing information is combined with weak lensing, the LSS correction does not significantly affect the fitting results with the adopted NFW/gNFW form. Moreover, all these models properly reproduce the observed location of the Einstein radius, θ_Ein ≈ 28''.

Here we summarize our primary results obtained with the LSS correction. The confidence contours on the NFW parameters (M_vir, c_vir) are shown in Figure 12. The constraints are strongly degenerate when only the inner or outer mass profile is included in this fit. Combining complementary weak- and strong-lensing information significantly narrows down the statistical uncertainties on the NFW model parameters, placing tighter constraints on the entire mass profile (Model 7 of Table 5): M_vir = 1.07^+0.20_{− 0.16} × 10¹⁵ M_☉ h⁻¹ and c_vir = 6.9^+1.0_{− 0.9} with χ²_min/dof = 18.0/31, corresponding to a Q-value goodness of fit of Q = 0.970. Next, when α is allowed to vary (Table 7), we find M_vir = 1.06^+0.23_{− 0.18} × 10¹⁵ M_☉ h⁻¹, c₋₂ = 7.0^+1.5_{− 1.4}, and α = 0.97^+0.28_{− 0.23} with χ²_min/dof = 18.0/30 and Q = 0.960, being consistent with the simple NFW form with α = 1. Thus, the addition of the α parameter has little effect on the fit, as shown by the quoted χ² and Q values. The two-dimensional marginalized constraints on (M_vir, α) and (c₋₂, α) are shown in Figure 13.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this subsection, we address the impact of the choice of strong-lensing models on the determination of the halo mass and concentration parameters in a joint weak- and strong-lensing analysis. As an alternative choice to the Zitrin et al. (2012) model, we consider here Pixelens (non-parametric) and Lenstool (parametric) models, in combination with our LSS-corrected weak-lensing mass model (Section 6.4). For each case, we define the χ² function for strong lensing as in Equation (22) and minimize the total χ² function (Equation (21)) to estimate the NFW parameters (M_vir, c_vir).

The resulting model constraints are tabulated in Table 8. We find that both parameters based on different strong-lensing profiles are consistent with each other within the statistical errors. This also indicates consistency between these strong-lensing models and our weak-lensing measurements, as shown in Figure 7. We find a tendency for Pixelens to yield somewhat higher mass estimates compared to other strong-lens modeling methods, as discussed by Grillo et al. (2010, their Appendix).

When the NFW (gNFW) form is assumed, the Zitrin et al. (2012) model predicts a somewhat higher concentration and a lower mass than other models as implied by its correspondingly higher central density at ≲ 02 (see Figure 7). When the inner fitting radius is increased from 4'' to 12'' (∼50 kpc h⁻¹), we find a fractional increase of ∼9% in M_vir and a fractional decrease of ∼17% in c_vir (6.9^+1.0_{− 0.9} → 5.7^+1.4_{− 1.1}). Including these variations as systematic uncertainties in our mass-concentration determination, the spherical NFW model for MACS1206 is constrained as M_vir = (1.07^+0.20_{− 0.16} ± 0.10) × 10¹⁵ M_☉ h⁻¹ and c_vir = 6.9^+1.0_{− 0.9} ± 1.2 (statistical followed by systematic uncertainty). Similarly, when the central 50 kpc h⁻¹ region is excluded from the fit, we have M_vir = (1.17^+0.25_{− 0.20} ± 0.10) × 10¹⁵ M_☉ h⁻¹ and c_vir = 5.7^+1.4_{− 1.1} ± 1.2.

Motivated by the apparently shallow projected density profile in the outer regions (cf. XMMU J2235.3−2557 at z = 1.4; Jee et al. 2009), we consider here a softened power-law sphere (SPLS) model (Grogin & Narayan 1996) as an alternative to the NFW profile and perform profile fitting analyses on our full-range mass profile data (derived from Methods 6 and 7 in Table 5; see Sections 6.4 and 6.5).

The SPLS model has the same number of free parameters as gNFW, namely, three. The SPLS density profile is given by ρ(r) = ρ₀(1 + r²/r²_c)^{(η − 3)/2}, where ρ₀ = ρ(0) is the central density, r_c is the core radius, and the power-law index η is restricted to lie in the range 0 ⩽ η ⩽ 2 (Grogin & Narayan 1996). At r ≫ r_c, M(< r)∝r^η. This reduces to a non-singular isothermal sphere (NIS) model when η = 1. The fitting results with and without the outer LSS correction (Methods 7 and 6, respectively) are summarized in Table 9.

First, when η is fixed to unity (NIS), the NIS model provides acceptable fits, but with larger residuals (χ²) compared to the corresponding NFW fits with the same degrees of freedom (31): Δχ² = χ²_{min, NIS} − χ²_{min, NFW} = 2.3 (Method 6) and 5.9 (Method 7) between the best-fit NIS and NFW models. Note that because of the asymptotic M(< r)∝r behavior, the assumed NIS form leads to substantially higher masses at large radius (r ≫ r_c) than what the NFW model predicts (∼35% higher than the NFW values at r = 1.6 Mpc h⁻¹).

Next, when the outer slope is allowed to vary, the fit is noticeably improved for the results with the outer LSS correction (Method 7), corresponding to a difference of Δχ² = χ²_{min, NIS} − χ²_{min, SPLS} = 4.4 between NIS and SPLS for 1 additional dof. For this, the best-fitting slope parameter is obtained as η = 0.77^+0.13_{− 0.17} (χ²_min = 19.5 for 30 dof), corresponding to 2.1 ⩽ γ_3D(r ≫ r_c) ⩽ 2.4. This SPLS model yields a virial mass of M_vir = (1.26 ± 0.37) × 10¹⁵ M_☉ h⁻¹ (r_vir ≈ 1.73 Mpc h⁻¹).

Gravitational lensing probes the total mass projected onto the sky along the line of sight, so that the lensing-based cluster mass measurements are sensitive to projection effects arising from (1) additional mass overdensities (underdensities) along the line of sight (Meneghetti et al. 2010b; Rasia et al. 2012) and (2) halo triaxiality (Hennawi et al. 2007; Oguri & Blandford 2009; Meneghetti et al. 2010b; Rasia et al. 2012).

7.1.1. Projection of Additional Mass Structures

7.1.2. Halo Triaxiality

Complementary multiwavelength observations serve as a useful guide to the likely degree of lensing bias. Here we retrieved and analyzed archival Chandra and XMM-Newton data of MACS1206 to obtain an independent cluster mass estimate, as well as to constrain the physical properties of the X-ray gas.

We perform a simultaneous fit to Chandra and XMM data sets under the assumption that the intracluster gas is in hydrostatic equilibrium (HSE) with the overall cluster potential of the NFW form. The tool used for this analysis is Joint Analysis of Cluster Observations (JACO; Mahdavi et al. 2007b); we refer the reader to this paper for the details of the X-ray analysis procedure, which we briefly summarize below.

We use Chandra ObsID 3277 and XMM-Newton observation 0502430401. We screen periods of flaring background according to standard procedure, resulting in usable exposure times of 23 ks and 26 ks, respectively. Appropriate co-added blank-sky fields allow us to subtract particle background spectra for both telescopes, and the residual (positive or negative) astrophysical background is included and marginalized over in the global cluster gas model. Spectra are extracted over seven annular bins for both Chandra and XMM-Newton. The extracted spectra extended out to a distance of 37 (1.26 Mpc) and contain an average of 1500 counts each.

The model for the gas density distribution is a single β-model multiplied by a power law of slope γ:

$$\begin{equation} \rho _g(r) = \rho _0 \left(\frac{r_c}{r}\right)^\gamma \left(1 + \frac{r^2}{r_c^2} \right)^{-3 \beta /2}. \end{equation} \tag{ 23 }$$

The power-law component is required to capture the steep increase of the density toward the center of the cluster; all parameters of the gas distribution are fit to the data. The metallicity is allowed to vary with radius as well, as are the parameters of the NFW mass profile. Model spectra are generated self-consistently in concentric spherical shells and forward projected onto the annular sky regions matching the extracted annuli. The resulting spectra are mixed using in-orbit energy- and position-dependent PSFs for both Chandra and XMM-Newton. Systematic calibration uncertainties between Chandra and XMM-Newton spectra are taken into account by adding a 4% error (a typical correction used in Mahdavi et al. 2008) in quadrature to each spectral bin used for the joint fits. This brings the joint χ² into the acceptable range (χ² = 1603 for 1541 dof). An MCMC procedure is used to estimate errors on the best-fit quantities. After marginalizing over all other parameters, we measure a total mass M₂₅₀₀ = (4.45 ± 0.28) × 10¹⁴ M_☉, a gas mass M_{gas, 2500} = (0.54 ± 0.02) × 10¹⁴ M_☉, an NFW concentration parameter of c₂₀₀ = 3.5 ± 0.5, an inner gas density profile slope of 0.7 ± 0.03, and a central cooling time of 2.1 ± 0.1 Gyr. In what follows, the examination of the X-ray results is conservatively limited to r < 1 Mpc.

In Figure 14 we plot the resulting X-ray-based total mass profile, M(< r), shown along with our NFW model from the full-lensing analysis. The results of the NFW fit are also reported in Table 10. This X-ray model yields a total mass of M_X = (4.6 ± 0.2) × 10¹⁴ M_☉ at the lensing-derived overdensity radius of r₂₅₀₀ ≈ 0.60 Mpc. This is in excellent agreement with the lensing mass at the same radius, M_lens = (4.9 ± 0.9) × 10¹⁴ M_☉, which corresponds to the X-ray-to-lensing mass ratio, a₂₅₀₀ = M_X(< r₂₅₀₀)/M_lens(< r₂₅₀₀) = 0.95^+0.23_{− 0.25}. The a₂₅₀₀ value obtained here is in good agreement with results from mock observations of 20 ΛCDM clusters by Rasia et al. (2012): a₂₅₀₀ = 0.94 ± 0.02. At this overdensity, no significant bias was observed in detailed observational studies by Zhang et al. (2008) and Mahdavi et al. (2008), who performed a systematic comparison of weak-lensing and X-ray mass measurements for sizable cluster samples. In the bottom panel of Figure 14, we show the X-ray-to-lensing mass ratio $a_\Delta$ as a function of cluster radius, in the radial range where X-ray observations are sufficiently sensitive. Overall, the mass ratio is consistent with unity especially at r ∼ r₂₅₀₀.

Zoom In Zoom Out Reset image size

Ebeling et al. (2009) obtained a hydrostatic mass estimate of M_X = (1.7 ± 0.1) × 10¹⁵ M_☉ at r = 2.3 Mpc (their estimate for r₂₀₀) assuming an isothermal β-model with β = 0.57 ± 0.02 and their estimated temperature k_BT = 11.6 ± 0.7 keV in the radial range [70, 1000] kpc (M_X∝β^1/2T), which is high but consistent within the errors with M_lens(< 2.3 Mpc) = (1.4 ± 0.3) × 10¹⁵ M_☉ obtained with our best NFW model based on the full-lensing analysis.

Our full-lensing results, when combined with X-ray gas mass measurements (M_gas), yield a direct estimate for the cumulative gas mass fraction, f_gas(< r) ≡ M_gas(< r)/M(< r), free from the HSE assumption. For this we use reduced Chandra X-ray data presented in the Archive of Chandra Cluster Entropy Profile Tables (ACCEPT; Cavagnolo et al. 2009). In Figure 15, we plot our f_gas measurements as a function of cluster radius. We find a gas mass fraction of f_gas(< r) = 13.7^+4.5_{− 3.0}% at a radius of r = 1 Mpc ≈ 1.7 r₂₅₀₀(≈ 0.8 r₅₀₀), a typical value observed for high-mass clusters (Umetsu et al. 2009; Zhang et al. 2009). When compared to the cosmic baryon fraction f_b = Ω_b/Ω_m = 0.1675 ± 0.006 constrained from the Wilkinson Microwave Anisotropy Probe (WMAP) seven-year data (Jarosik et al. 2011), this indicates f_gas/f_b = 0.82^+0.27_{− 0.18} at this radius. At the innermost measurement radius r ≈ 40 kpc where the lensing and X-ray data overlap, we have f_gas(< r) = 3.4^+1.2_{− 0.8}%. Thus, the hot gas represents only a minor fraction of the total lensing mass near the cluster center, as found for other high-mass clusters (Lemze et al. 2008; Umetsu et al. 2009).

Zoom In Zoom Out Reset image size

Additionally, we derive a mass profile from simulated annealing fits of the ACCEPT pressure profile (Cavagnolo et al. 2009), adopting the Arnaud et al. (2010, A10) "universal profile" (M. Donahue et al. 2012, in preparation). This Chandra-only mass profile is shown to be in good agreement with the lensing as well as joint Chandra+XMM results (Figure 14). The Chandra-only gas mass fraction profile is also shown in Figure 15.

We conclude, on the basis of these results and comparison with detailed statistical studies, that the level of orientation bias in this cluster is not significant given the large uncertainties in our lensing/X-ray observations, as well as the possible contribution from non-thermal pressure in the cluster core (e.g., Kawaharada et al. 2010).

We have also compared our lensing-derived results to mass estimates determined from the SZE data. Using Bolocam at the Caltech Submillimeter Observatory, we observed MACS1206 for approximately 11 hr in 2011 April. These data were collected with Bolocam configured at an SZE-emission-weighted band center of 140 GHz. Further details of the Bolocam instrument are given in Haig et al. (2004). We detect the cluster with an S/N value of 21.1 and a white noise rms of 24.9 μK_CMB arcmin. We reduced these data according to the procedure described in detail in Sayers et al. (2011), but with the updated calibration model reported in Sayers et al. (2012) and some other minor modifications.

The key steps involved in our Bolocam data reduction and cluster modeling are summarized as follows. We first remove sky noise from the time streams by subtracting a template of the correlated signal over the field of view followed by a high-pass filter (Sayers et al. 2011). This process results in a filtered image of the true SZE signal (see the left panel in Figure 16), where the filtering is weakly dependent on the cluster shape due to the correlated template removal. To characterize this filtering, we process a beam-smoothed, initial best-fit cluster profile by reverse mapping it using our pointing information. These data are then processed iteratively with a new best-fit profile using our full reduction pipeline, until the procedure converges. For this analysis we use the A10 "universal pressure profile," which adopts a form of the Nagai et al. (2007) pressure profile with its slopes fixed to the values given in A10, allowing the overall normalization and scale radius to vary.

Zoom In Zoom Out Reset image size

We have derived cluster mass estimates from our SZE data alone using the method outlined in Mroczkowski (2011). The key innovation of this method is that, in addition to assuming HSE, the virial theorem is used, which is no stronger an assumption than HSE and can be derived from HSE and thermodynamics. This method determines the underlying total mass profile from an SZE-determined pressure profile, with the added assumption of a constant gas mass fraction f_gas. Cluster mass estimates derived with this method have been shown to be consistent with X-ray-derived results using data from the SZA (Mroczkowski 2011) and SZA follow-up of blind SZE detections using the Atacama Cosmology Telescope (Reese et al. 2012). The SZE-only mass estimates for MACS1206 are given in Table 11, which presents $M_\Delta$ and $r_\Delta$ values at overdensities Δ = 2500 and 500 derived from our Bolocam data alone, under the assumptions made. This table also contains Bolocam-derived mass estimates at the lensing-derived values of r₂₅₀₀ and r₅₀₀. We note that the values in Table 11 include an estimate of our systematic errors on the SZE-derived masses, which we describe in detail below.

Table 11. Bolocam SZE-derived Cluster Mass Estimates

Overdensitya	Bolocam-derived $r_\Delta$b		Lensing-derived $r_\Delta$c
Δ	$r_\Delta$	$M({<}r_\Delta)$	$r_\Delta$	$M({<}r_\Delta)$
	(Mpc)	(1014 M☉)	(Mpc)	(1014 M☉)
2500	0.63+0.01 + 0.06− 0.02 − 0.05	5.8+0.4 + 1.7− 0.4 − 1.4	0.60	5.3+0.2 + 0.8− 0.2 − 0.7
500	1.67+0.09 + 0.12− 0.08 − 0.12	21.2+3.7 + 5.1− 3.0 − 4.3	1.31	15.7+1.2 + 2.3− 1.1 − 2.1

Notes. For each value the first error estimate represents our measurement uncertainty and the second error estimate represents our uncertainty due to systematics in our fitting method, flux calibration, and choice of parameterization. See for details Section 7.3. All quantities here are given in physical units assuming the concordance ΛCDM cosmology ($h=0.7, \Omega _m=0.3, \Omega _\Lambda =0.7$). ^aMean interior overdensity with respect to the critical density of the universe at the cluster redshift z = 0.439. ^bBolocam cluster mass estimates at the Bolocam-SZE-derived values of overdensity radius $r_\Delta$. ^cBolocam cluster mass estimates at the lensing-derived values of overdensity radius $r_\Delta$.

Download table as:  ASCIITypeset image

The dominant source of uncertainty in our mass estimates, as discussed in Mroczkowski (2011), stems from the uncertainty in the assumed value of a radially constant f_gas(r). Masses derived under this assumption scale as ∝f^−1/2_gas. We adopt the value f_gas = 0.13 and marginalize over uncertainties for a range f_gas = [0.1, 0.17], consistent with our X-ray-determined gas fraction measurements at radii near r₂₅₀₀ (see Figure 15). An additional source of systematic uncertainty is the absolute calibration of the Bolocam maps, which is about 5% and results in a ≲ 5% uncertainty in our derived masses. Finally, we include a ±1.5% systematic at r₂₅₀₀ and ±5% systematic at r₅₀₀, due to our particular choice of parameterization for the pressure profile. These values are roughly consistent with those shown in Mroczkowski (2011) for different parameterization of the exponents in the pressure profiles.

By comparison to the lensing mass estimates, we find an SZE-to-lensing mass ratio of a₂₅₀₀ = 1.08 ± 0.29 ± 0.22 (statistical followed by systematic at 68% confidence) at the lensing-derived overdensity radius r₂₅₀₀ of 0.60 Mpc (Table 11). Hence, our lensing mass estimate is in agreement with both the X-ray and SZE mass estimates at r₂₅₀₀. At a lower overdensity of Δ = 500, we find an SZE-to-lensing mass ratio of a₅₀₀ = 1.55 ± 0.30 ± 0.26 at the lensing-derived radius r₅₀₀ of 1.3 Mpc, roughly consistent with unity within large errors.

MACS1206 is an X-ray-luminous cluster at a redshift of z = 0.439, or a cosmic time of t ∼ 9 Gyr. The cluster appears relatively relaxed in projection in both optical and X-ray images, with a pronounced X-ray peak at the BCG position (Ebeling et al. 2009; Postman et al. 2012). This cluster was classified to be relaxed by Gilmour et al. (2009) on the basis of a visual examination of its X-ray morphology. Our detailed morphology analysis shows no sign of significant recent merging activity around the BCG, which is also supported by our strong-lensing analysis, finding no significant offset between the DM center of mass, BCG, and X-ray peak (Section 4.2). A good agreement between the lensing and X-ray mass estimates (Section 7.2) indicates that the hot gas is not far from a state of HSE in the cluster potential well.

However, some evidence of merger activity along the line of sight was suggested by the high velocity dispersion of σ_v ≈ 1580 km s⁻¹ based on 38 redshift measurements (Ebeling et al. 2009). Recently, a much larger spectroscopic sample of cluster members has been obtained for this cluster using VLT/VIMOS (P. Rosati et al. 2012, in preparation). Defining membership is crucial for a dynamical analysis since interlopers by projection effects can largely bias the derived projected velocity dispersion, especially at large radii where the number density of cluster members is low (Wojtak et al. 2007). Using a secure sample of >400 cluster members identified in the projected phase space (e.g., Biviano & Salucci 2006; Lemze et al. 2009), we find that the velocity dispersion profile decreases outward fairly rapidly from ∼1500 km s⁻¹ in the central region to ∼800 km s⁻¹ at a projected distance of R ∼ 2 Mpc. Accordingly, the dynamical mass estimate is in agreement with the lensing estimate (A. Biviano et al. 2012, in preparation). This may argue against a strong deviation from dynamical relaxation.

The present Chandra analysis yields a gas temperature of 10.8 ± 0.7 keV averaged in the radial range [70, 700] kpc. Assuming that the galaxies and the gas are confined in the same gravitational potential well, this is consistent with a line-of-sight velocity dispersion of ∼1300 km s⁻¹, which is again in agreement with the observed value. This may also suggest that the cluster is not far from equilibrium. The cluster appears fairly round in both Chandra and XMM images at large distances from the cluster center, as demonstrated in Figure 10. X-ray emission is concentrated around and peaked on the BCG but shows some elongation within θ ≲ 1' at a position angle around 120° east of north, aligned with the orientation of the projected mass distribution. The surface brightness profile is fairly smooth, but there might be some tiny hints of discontinuities (see the ACCEPT catalog).⁴⁸ However, a much deeper observation is required to confirm them.

Finally, morphological analysis of Bolocam data has been performed in an identical way to the procedure used in Sayers et al. (2011). We find an ellipticity of (10 ± 7)% with a position angle of 55° ± 27° north of west from elliptical A10 model fits to our Bolocam SZE data. The fits include all data within a 14 × 14 arcmin square, corresponding to a fairly large aperture of θ_max ≈ 9' > θ_vir ∼ 7'. Of the approximately 50 clusters observed with Bolocam and fit with an elliptical A10 model, MACS1206 is one of the more circularly symmetric model fits.

In Figure 17, we summarize our full-lensing constraints on the mass and concentration parameters of MACS1206, along with recent ΛCDM predictions for relaxed cluster-sized halos based on N-body simulations (Duffy et al. 2008; Klypin et al. 2011; Bhattacharya et al. 2011; Prada et al. 2011). Our range of allowed concentration values (4.6 ⩽ c_vir ⩽ 7.9 at 1σ; see Section 6.6) span the high end and average expectations (4 ≲ 〈c_vir〉 ≲ 7) from ΛCDM simulations (Duffy et al. 2008; Zhao et al. 2009; Klypin et al. 2011; Bhattacharya et al. 2011; Prada et al. 2011). Average concentrations for relaxed clusters are found to be ∼10% higher and have lower scatter than those for the full population of halos (Duffy et al. 2008; Bhattacharya et al. 2011). A relatively high concentration of MACS1206 may also be indicated by the large Einstein radius θ_Ein ≈ 28'' (17'') at z_s = 2.5 (1.0) (Ebeling et al. 2009; Zitrin et al. 2012).

Zoom In Zoom Out Reset image size

Care must be taken when comparing these predictions for spherically averaged halo structure with our lensing results, which are obtained from an NFW fit to the projected lensing measurements assuming a spherical halo. In the previous subsection (Section 7.1), we have shown that our lensing results are in good agreement with the X-ray-derived mass profiles (see Figures 14 and 15) in the region of overlap (≲ 1 Mpc), as well as with the Bolocam SZE mass estimates (Section 7.3), suggesting that the level of orientation bias (see Section 7.1.2) is not significant in this cluster.

Additionally, the effects of baryonic physics can impact the inner halo profile (at r ≲ 0.0 5r_vir; Duffy et al. 2010) and thus modify the gravity-only c–M relation, especially for less massive halos (M_vir ≲ 4 × 10¹⁴ M_☉ h⁻¹; see Bhattacharya et al. 2011). Using cosmological hydrodynamical simulations including the back-reaction of baryons on DM, Duffy et al. (2010) found a <20% increase in the halo concentration for cluster-sized halos (M_vir < 6 × 10¹⁴ M_☉ h⁻¹ at z = 0). When excluding the central 50 kpc h⁻¹ (≈0.03 r_vir) region from our primary strong-lensing mass model (Zitrin et al. 2012), we find a ≈17% decrease in the best-fit concentration parameter derived from our full-lensing analysis (Section 6.6), as demonstrated in Figure 17. We note that the CLASH clusters are massive (5 × 10¹⁴ < M_vir/M_☉ < 3 × 10¹⁵; see Postman et al. 2012) and hence expected to be less affected by baryonic effects.

For this cluster, the lensing-derived total mass distribution is consistent with the NFW form (α = γ_3D(r → 0) = 0.96^+0.31_{− 0.49}), as found for several relaxed clusters: A611 (Newman et al. 2009); A383 (Zitrin et al. 2011c); A1703 (α ≈ 0.9 Richard et al. 2009; Oguri et al. 2009); a stacked full-lensing analysis of A1689, A1703, A370, and Cl 0024+17 (α = 0.89^+0.27_{− 0.39}; Umetsu et al. 2011a). Multiwavelength observations can be used to measure gas and stellar density profiles for subtraction from lensing-derived total mass profiles to yield DM-only mass profiles (Lemze et al. 2008; Newman et al. 2009), allowing for a more direct comparison with CDM predictions from gravity-only simulations. We defer this analysis to a forthcoming paper.

We first summarize the Bayesian method of Umetsu et al. (2011b) for a direct reconstruction of the cluster mass profile from combined radial distortion and magnification profiles.

In the Bayesian framework, we sample from the posterior 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$})$. For a mass profile analysis, $\mbox{\boldmath $s$}$ is a vector containing the discrete convergence profile, κ_i ≡ κ(θ_i) with i = 1, 2, ..., N in the subcritical regime (θ_i > θ_Ein), and the average convergence within the inner radial boundary θ_min of the weak-lensing data, $\overline{\kappa }_{\rm min}\equiv \overline{\kappa }({<}\theta _{\rm min})$, so that $\mbox{\boldmath $s$}=\lbrace \overline{\kappa }_{\rm min},\kappa _i\rbrace _{i=1}^{N}$, being specified by (N + 1) parameters.

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{ A1 }$$

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}={\cal L}_{g_+}{\cal L}_\mu$, where ${\cal L}_{g_+}$ and ${\cal L}_\mu$ are the likelihood functions for tangential distortion and magnification bias, respectively. The log-likelihood functions for the weak-lensing observations {g_{+, i}}^N_{i = 1} and {n_{μ, i}}^N_{i = 1} are given, respectively (ignoring constant terms), as

$$\begin{eqnarray} l_{g_+}(\mbox{\boldmath $s$}) &\equiv & {-}\ln {{\cal L}_g}= \frac{1}{2} \displaystyle \sum _{i=1}^{N} \frac{[g_{+,i}-\hat{g}_{+,i}(\mbox{\boldmath $s$})]^2}{\sigma _{+,i}^2}, \end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray} l_\mu (\mbox{\boldmath $s$}) &\equiv & {-}\ln {{\cal L}_\mu }= \frac{1}{2} \displaystyle \sum _{i=1}^{N} \frac{[n_{\mu ,i}-\hat{n}_{\mu ,i}(\mbox{\boldmath $s$})]^2}{\sigma _{\mu ,i}^2}, \end{eqnarray} \tag{ A3 }$$

where $\lbrace \hat{g}_{+,i}\rbrace _{i=1}^{N}$ and $\lbrace \hat{n}_{\mu ,i}\rbrace _{i=1}^{N}$ are the theoretical predictions for the corresponding observations. The total likelihood $l_{\rm 1D}(\mbox{\boldmath $s$})\equiv -\ln {\cal L}$ of the combined observations is obtained as

$$\begin{equation} l_{\rm 1D} = l_{g_+} +l_\mu. \end{equation} \tag{ A4 }$$

Here we consider a simple flat prior with a lower bound of $\mbox{\boldmath $s$}=0$. Additionally, we account for the uncertainty in the calibration parameters, $\mbox{\boldmath $c$}=(n_0,s,\omega)$, namely, the normalization and slope parameters (n₀, s) of the background counts and the relative lensing depth ω ≡ 〈β(red)〉/〈β(back)〉 between the background samples used for the magnification and distortion measurements.

We use the MCMC technique with Metropolis–Hastings sampling to constrain our mass model $\mbox{\boldmath $s$}$. The covariance matrix ${\cal C}$ of $\mbox{\boldmath $s$}$ is obtained from MCMC samples.

Here we extend the 1D method of Umetsu et al. (2011b) to a 2D mass distribution $\kappa (\mbox{\boldmath $\theta $})$, by combining 2D distortion data with the azimuthally averaged magnification information. For this analysis, the signal $\mbox{\boldmath $s$}$ is a vector of parameters containing discrete mass elements on a 2D Cartesian grid of independent cells: $\mbox{\boldmath $s$}=\lbrace \kappa _m\rbrace _{m=1}^{N_{\rm cell}}$. The $\gamma (\mbox{\boldmath $\theta $})$ field can be written as a linear combination of the parameters $\mbox{\boldmath $s$}$ (Equation (3)). Then, the distortion $g(\mbox{\boldmath $\theta $})$ and magnification $\mu (\mbox{\boldmath $\theta $})$ fields can be uniquely specified in the subcritical regime (Section 2).

In analogy to Equation (13), we calculate the weighted average $g_{\alpha ,m}\equiv g_{\alpha }(\mbox{\boldmath $\theta $}_m)$ (α = 1, 2) of individual distortion estimates and its covariance matrix,

$$\begin{equation} {\rm Cov}[g_{\alpha ,m}, g_{\beta ,n}]\equiv (C_g)_{\alpha \beta ,mn} = \frac{1}{2} \sigma ^2_{g}(\mbox{\boldmath $\theta $}_m) \delta _{mn} \delta _{\alpha \beta }, \end{equation} \tag{ A5 }$$

where $\sigma ^2_g(\mbox{\boldmath $\theta $}_m)$ is the standard error of the weighted mean distortion, $g(\mbox{\boldmath $\theta $}_m)$. Accordingly, the 2D shear log-likelihood function $l_g(\mbox{\boldmath $s$})\equiv -\ln {{\cal L}_g}$ is written as

$$\begin{equation} l_g(\mbox{\boldmath $s$}) = \frac{1}{2} \sum _{m,n=1}^{N_{\rm cell}} \sum _{\alpha ,\beta =1}^{2} [g_{\alpha ,m}-\hat{g}_{\alpha ,m}(\mbox{\boldmath $s$})] ({\cal W}_g)_{\alpha \beta ,mn}[g_{\beta ,n}-\hat{g}_{\beta ,n}(\mbox{\boldmath $s$})], \end{equation} \tag{ A6 }$$

where $\hat{g}_{\alpha ,m}(\mbox{\boldmath $s$})$ is the theoretical prediction for g_{α, m} and $({\cal W}_g)_{\alpha \beta ,mn}$ is the shear weight matrix,

$$\begin{equation} ({\cal W}_g)_{\alpha \beta ,mn} = M_m M_n \left({\cal C}_g^{-1}\right)_{\alpha \beta ,mn}, \end{equation} \tag{ A7 }$$

with M_m being a mask weight, defined such that M_m = 0 if the mth cell is masked out and M_m = 1 otherwise. In practice, we exclude from our analysis innermost cells that lie in the cluster central region, where the surface-mass density can be close to or greater than the critical value (i.e., κ ≳ 1). Furthermore, this is crucial to minimize contamination by unlensed cluster member galaxies (see Section 3.4).

Now we combine 2D distortion data with magnification information to obtain the total log-likelihood $l_{\rm 2D}(\mbox{\boldmath $s$})$ as

$$\begin{equation} l_{\rm 2D} = l_{g} + l_\mu , \end{equation} \tag{ A8 }$$

where l_μ, given by Equation (A3), imposes a set of azimuthally integrated constraints on the underlying κ field. Since the degree of magnification is locally related to κ, this will essentially provide the (otherwise unconstrained) normalization of $\kappa (\mbox{\boldmath $\theta $})$ over a set of concentric rings where count measurements n_{μ, i} are available. Note that no assumption is made of azimuthal symmetry or isotropy of the cluster mass distribution.

This 2D inversion problem involves estimation of a large number of parameters $\mbox{\boldmath $s$}$; typically, N_cell ≳ 1000 when distortion data are binned into subarcminute pixels. We use in our implementation the conjugate-gradient method (Press et al. 1992) to find the best solution. We include Gaussian priors on the calibration nuisance parameters $\mbox{\boldmath $c$}=(s,n_0,\omega)$, given by means of quadratic penalty terms with mean values and variances directly estimated from data. The log posterior pdf, $F = -\ln {P(\mbox{\boldmath $s$}|\mbox{\boldmath $d$})}$, is expressed as a linear sum of $l_{\rm 2D}(\mbox{\boldmath $s$})$ and the prior terms on $\mbox{\boldmath $c$}$. The best-fit parameters are determined with a maximum-likelihood estimation, by minimizing the function F with respect to $\mbox{\boldmath $p$}\equiv (\mbox{\boldmath $s$},\mbox{\boldmath $c$})$, a vector containing the mass and calibration parameters. Here we employ an analytic expression for the gradient function $\mbox{\boldmath $\nabla $}F(\mbox{\boldmath $p$})$ obtained in the nonlinear subcritical regime. To quantify the errors on the mass reconstruction, we evaluate the Fisher matrix at the maximum-likelihood estimate $\mbox{\boldmath $p$}=\hat{\mbox{\boldmath $p$}}$ as

$$\begin{equation} {\cal F}_{mn} = \left\langle \frac{\partial ^2 F(\mbox{\boldmath $p$})}{\partial p_m \partial p_n} \right\rangle \Big |_{\mbox{\boldmath $p$}=\hat{\mbox{\boldmath $p$}}}, \end{equation} \tag{ A9 }$$

where the angular brackets represent an ensemble average and the indices (m, n) run over all model parameters. We estimate the covariance matrix ${\cal C}$ of $\mbox{\boldmath $s$}$ as

$$\begin{equation} {\cal C}_{mn} = ({\cal F}^{-1})_{mn}. \end{equation} \tag{ A10 }$$