The Projected Dark and Baryonic Ellipsoidal Structure of 20 CLASH Galaxy Clusters^*

We use the 2D reduced shear field averaged on a grid as the primary constraint from our wide-field weak-lensing observations. From shape measurements of background galaxies, the source-averaged reduced shear g_n = g(${\boldsymbol{\theta }}$_n) is measured on a regular Cartesian grid of pixels (n = 1,2, ..., ${N}_{\mathrm{pix}}$) as

$$\begin{eqnarray}&&{g}_{n}=\left[\displaystyle \sum _{k}S({{\boldsymbol{\theta }}}_{(k)},{{\boldsymbol{\theta }}}_{n}){w}_{(k)}{g}_{(k)}\right]{\left[\displaystyle \sum _{k}S({{\boldsymbol{\theta }}}_{(k)},{{\boldsymbol{\theta }}}_{n}){w}_{(k)}\right]}^{-1},\end{eqnarray} \tag{ 9 }$$

where S(${\boldsymbol{\theta }}$, ${\boldsymbol{\theta }}$') is a spatial window function, g_(k) is an estimate of g(${\boldsymbol{\theta }}$) for the kth galaxy at ${\boldsymbol{\theta }}$_(k), and w_(k) is its statistical weight, ${w}_{(k)}=1/({\sigma }_{g(k)}^{2}+{\alpha }_{g}^{2})$, with ${\sigma }_{g(k)}^{2}$ the error variance of g_(k). Following Umetsu et al. (2014), we choose α_g = 0.4, a typical value of the mean rms ${(\overline{{\sigma }_{g}^{2}})}^{1/2}$ found in Subaru observations (Umetsu et al. 2009; Oguri et al. 2010; Okabe & Smith 2016).

The theoretical expectation (denoted by a hat symbol) for the estimator (9) is approximated by (Seitz & Schneider 1997; Umetsu et al. 2015)

$$\begin{eqnarray}&&\widehat{g}({{\boldsymbol{\theta }}}_{n})\approx \displaystyle \frac{\langle W{\rangle }_{g}{\gamma }_{\infty }({{\boldsymbol{\theta }}}_{n})}{1-{f}_{W,g}\langle W{\rangle }_{g}{\kappa }_{\infty }({{\boldsymbol{\theta }}}_{n})},\end{eqnarray} \tag{ 10 }$$

where $\langle W{\rangle }_{g}$ is the source-averaged relative lensing strength (Section 2.2), and ${f}_{W,g}=\langle {W}^{2}{\rangle }_{g}/\langle W{\rangle }_{g}^{2}$ is a dimensionless correction factor of order unity. The error variance σ_g,n² for g_n is expressed as

$$\begin{eqnarray}&&{\sigma }_{g,n}^{2}=\displaystyle \frac{{\displaystyle \sum }_{k}{S}^{2}({{\boldsymbol{\theta }}}_{(k)},{{\boldsymbol{\theta }}}_{n}){w}_{(k)}^{2}{\sigma }_{g(k)}^{2}}{{\left[{\displaystyle \sum }_{k}S({{\boldsymbol{\theta }}}_{(k)},{{\boldsymbol{\theta }}}_{n}){w}_{(k)}\right]}^{2}}.\end{eqnarray} \tag{ 11 }$$

We adopt the top-hat window of radius θ_f (Merten et al. 2009; Umetsu et al. 2015), $S({\boldsymbol{\theta }},{\boldsymbol{\theta }}^{\prime} )=H({\theta }_{{\rm{f}}}-| {\boldsymbol{\theta }}-{\boldsymbol{\theta }}^{\prime} | )$, with H(x) the Heaviside function defined such that H(x) = 1 if x ≥ 0 and H(x) = 0 otherwise. The shape-noise covariance matrix for g_α,n = g_α(${\boldsymbol{\theta }}$_n) is then given as (Oguri et al. 2010)³¹

$$\begin{eqnarray}&&{({C}_{g}^{\mathrm{shape}})}_{\alpha \beta ,{mn}}=\displaystyle \frac{1}{2}{\delta }_{\alpha \beta }{\sigma }_{g,m}{\sigma }_{g,n}{\xi }_{H}(| {{\boldsymbol{\theta }}}_{m}-{{\boldsymbol{\theta }}}_{n}| ),\end{eqnarray} \tag{ 12 }$$

where the indices α and β run over the two components of the reduced shear (α, β = 1, 2), δ_αβ denotes the Kronecker delta, and ξ_H(x) is the autocorrelation of a pillbox of radius θ_f (White et al. 1999; Park et al. 2003; Umetsu et al. 2015),

$$\begin{eqnarray}&&{\xi }_{H}(x)=\displaystyle \frac{2}{\pi }\left[{\cos }^{-1}\left(\displaystyle \frac{x}{2{\theta }_{{\rm{f}}}}\right)-\left(\displaystyle \frac{x}{2{\theta }_{{\rm{f}}}}\right)\sqrt{1-{\left(\displaystyle \frac{x}{2{\theta }_{{\rm{f}}}}\right)}^{2}}\right]\end{eqnarray} \tag{ 13 }$$

for $| x| \leqslant 2{\theta }_{{\rm{f}}}$ and ξ_H(x) = 0 for $| x| \gt 2{\theta }_{{\rm{f}}}$.

In this study, we directly measure 2D reduced shear maps in 20 CLASH cluster fields using the wide-field shear catalogs obtained by the CLASH collaboration (Umetsu et al. 2014). The shear measurement pipeline of Umetsu et al. (2014) is based on the IMCAT package (Kaiser et al. 1995, hereafter KSB) with modifications incorporating the improvements developed by Umetsu et al. (2010).

Briefly summarizing, the key feature in our analysis pipeline is that only those galaxies detected with sufficiently high significance, ν_g > 30, are used to model the isotropic PSF correction as a function of object size and magnitude. Here, ν_g is the peak detection significance given by the IMCAT peak-finding algorithm hfindpeaks. A very similar procedure was employed by the LoCuSS collaboration in their weak-lensing study of 50 clusters based on Subaru/Suprime-Cam data (Okabe & Smith 2016). Another key feature is that we select those galaxies isolated in projection for the shape measurement, reducing the impact of crowding and blending (for details, see Umetsu et al. 2014). After the close-pair rejection, objects with low detection significance ν_g < 10 were excluded from our analysis. All galaxies with usable shape measurements are matched with those in our CC-selected samples of background galaxies (Section 3.2), ensuring that each galaxy is detected in both the reddest CC-selection band and the shape measurement band. Applying conservative selection criteria (Section 3.2), Umetsu et al. (2014) find a typical surface number density of n_g ≃ 12 galaxies arcmin⁻² for their weak-lensing-matched background catalogs (Umetsu et al. 2014, their Table 3).

In Appendix A, we show the results of our shape measurement test based on simulated Subaru/Suprime-Cam images provided by M. Oguri. From this test, we find that the reduced shear signal can be recovered up to $| {g}_{\alpha }| \simeq 0.3$ with m_α ≃ −0.05 of the multiplicative calibration bias and $| {c}_{\alpha }| \lt {10}^{-3}$ of the residual shear offset (for details, see Appendix A), where the observed and true values of the reduced shear are related by (Heymans et al. 2006; Massey et al. 2007)

$$\begin{eqnarray}&&{g}_{\alpha }^{\mathrm{obs}}=(1+{m}_{\alpha }){g}_{\alpha }^{\mathrm{true}}+{c}_{\alpha }.\end{eqnarray} \tag{ 14 }$$

Umetsu et al. (2014) included for each galaxy a shear calibration factor with m₁ = m₂ = −0.05 to account for residual calibration. The degree of multiplicative bias m_α depends on the seeing conditions and the PSF properties (Figure 10), so that the variation with the PSF properties limits the shear calibration accuracy to δm_α ∼ 0.05 (Umetsu et al. 2012, Section 3.3). We note that the same simulation data set was used by the LoCuSS collaboration (Okabe & Smith 2016) to test their shape measurement pipeline, and Okabe & Smith (2016) found a similar level of shear calibration bias (m_α ∼ −0.03) from their mock observations.

For all cluster fields in our sample, the estimated values for $\langle \beta {\rangle }_{g}$ and f_W,g are summarized in Table 3 of Umetsu et al. (2014). The calibration uncertainty in $\langle \beta {\rangle }_{g}$ is marginalized over in our joint analysis of shear and magnification data (Section 3.5).

Lensing magnification can influence the observed surface number density of background sources, amplifying their apparent fluxes and expanding the area of sky (Hildebrandt et al. 2011; Umetsu et al. 2011b; Morrison et al. 2012; Garcia-Fernandez et al. 2016; Tudorica et al. 2017). The former effect increases the number of sources detectable above the limiting flux, whereas the latter reduces the effective observing area in the source plane, reducing the number of sources per solid angle. The net effect is known as magnification bias (Broadhurst et al. 1995) and depends on the slope of the intrinsic source luminosity function. Since a given flux limit corresponds to different luminosities at different source redshifts, number counts of distinctly different source populations probe different regimes of magnification bias (Umetsu 2013).

Deep multiband photometry can be used to sample the faint end of the luminosity function of red quiescent galaxies lying at z ∼ 1 (e.g., Ilbert et al. 2010). For such a source population, the effect of magnification bias is dominated by the geometric area distortion, because there are relatively few fainter galaxies that can be magnified into the flux-limited sample. This effect results in a net depletion of source counts, a phenomenon known as negative magnification bias or weak-lensing depletion (e.g., Broadhurst 1995; Taylor et al. 1998; Broadhurst et al. 2005b; Umetsu & Broadhurst 2008; Umetsu et al. 2011b, 2012, 2014, 2015; Coe et al. 2012; Ford et al. 2012; Medezinski et al. 2013; Radovich et al. 2015; Ziparo et al. 2016; Wong et al. 2017). A practical advantage of this technique, at the expense of very deep multicolor imaging, is that the effect is not sensitive to the exact form of the source luminosity function (Umetsu et al. 2014).

In the weak-lensing regime, the shift in magnitude δm = 2.5 log₁₀ μ due to magnification is small compared to the range in which the slope of the luminosity function varies. The number counts can then be approximated by a power law at the limiting magnitude m_lim. The expectation value (denoted by a hat symbol) for the lensed counts at source redshift z is then expressed as (Broadhurst et al. 1995)

$$\begin{eqnarray}&&{\widehat{N}}_{\mu }({\boldsymbol{\theta }},z| \lt {m}_{\mathrm{lim}})={\overline{N}}_{\mu }(z| \lt {m}_{\mathrm{lim}})\,{{\rm{\Delta }}}_{\mu }{({\boldsymbol{\theta }},z)}^{1-2.5s},\end{eqnarray} \tag{ 15 }$$

with ${\overline{N}}_{\mu }(z| \lt {m}_{\mathrm{lim}})$ the unlensed mean counts per cell and s the logarithmic count slope evaluated at m = m_lim,³²

$$\begin{eqnarray}&&s=\displaystyle \frac{d{\mathrm{log}}_{10}\overline{N}(z| \lt m)}{{dm}}{| }_{{m}_{\mathrm{lim}}}.\end{eqnarray} \tag{ 16 }$$

A count depletion (enhancement) results when s < 0.4 (>0.4).

Following Umetsu et al. (2014), we interpret the observed source-averaged magnification bias as

$$\begin{eqnarray}{\widehat{b}}_{\mu }({\boldsymbol{\theta }}) & \equiv & \displaystyle \frac{{\widehat{N}}_{\mu }({\boldsymbol{\theta }}| \lt {m}_{\mathrm{lim}})}{{\overline{N}}_{\mu }(\lt {m}_{\mathrm{lim}})}\approx {{\rm{\Delta }}}_{\mu }{({\boldsymbol{\theta }})}^{1-2.5{s}_{\mathrm{eff}}},\\ {{\rm{\Delta }}}_{\mu }({\boldsymbol{\theta }}) & = & \displaystyle \frac{{\displaystyle \int }_{0}^{\infty }{dz}\,{\overline{N}}_{\mu }(z| \lt {m}_{\mathrm{lim}}){{\rm{\Delta }}}_{\mu }({\boldsymbol{\theta }},z)}{{\displaystyle \int }_{0}^{\infty }{dz}\,{\overline{N}}_{\mu }(z| \lt {m}_{\mathrm{lim}})}\\ & \approx & {[1-\langle W{\rangle }_{\mu }{\kappa }_{\infty }({\boldsymbol{\theta }})]}^{2}-\langle W{\rangle }_{\mu }^{2}| {\gamma }_{\infty }({\boldsymbol{\theta }}){| }^{2},\end{eqnarray} \tag{ 17 }$$

where ${\overline{N}}_{\mu }(\lt {m}_{\mathrm{lim}})={\int }_{0}^{\infty }{dz}\,{\overline{N}}_{\mu }(z| \lt {m}_{\mathrm{lim}})$, ${s}_{\mathrm{eff}}=d{\mathrm{log}}_{10}{\overline{N}}_{\mu }(\lt m)/{dm}{| }_{{m}_{\mathrm{lim}}}$, and $\langle W{\rangle }_{\mu }$ is the source-averaged relative lensing strength (Section 2.2). Equation (17) gives a good approximation for depleted populations with s_eff ≪ 0.4 (for details, see Appendix A.2 of Umetsu 2013). For simplicity, we write ${N}_{\mu }({\boldsymbol{\theta }})\,={N}_{\mu }({\boldsymbol{\theta }}| \lt {m}_{\mathrm{lim}})$ and ${\overline{N}}_{\mu }={\overline{N}}_{\mu }(\lt {m}_{\mathrm{lim}})$. In the weak-lensing limit, Equation (17) reads ${\widehat{b}}_{\mu }-1\approx (5{s}_{\mathrm{eff}}-2)\langle W{\rangle }_{\mu }{\kappa }_{\infty }$.

We azimuthally average the observed counts N_μ(${\boldsymbol{\theta }}$) in clustercentric, circular annuli and calculate the surface number density $\{{n}_{\mu ,i}\}{}_{i=1}^{{N}_{\mathrm{bin}}}$ of background galaxies as (Umetsu et al. 2015, 2016)

$$\begin{eqnarray}&&{n}_{\mu ,i}=\displaystyle \frac{1}{(1-{f}_{\mathrm{mask},i}){{\rm{\Omega }}}_{\mathrm{cell}}}\displaystyle \sum _{m}{{ \mathcal P }}_{{im}}{N}_{\mu }({{\boldsymbol{\theta }}}_{m}),\end{eqnarray} \tag{ 18 }$$

with Ω_cell the solid angle per cell and ${{ \mathcal P }}_{{im}}={({\sum }_{m}{A}_{{mi}})}^{-1}{A}_{{mi}}$ the projection matrix normalized in each annulus by ${\sum }_{m}{{ \mathcal P }}_{{im}}=1$. Here, A_mi is the area fraction of the mth cell lying within the ith radial bin (0 ≤ A_mi ≤ 1), and f_mask,i is the mask correction factor for the ith radial bin due to saturated objects, foreground and cluster galaxies, and bad pixels (for details, see Section 3.4 of Umetsu et al. 2016).

The theoretical expectation for the estimator (18) is

$$\begin{eqnarray}&&{\widehat{n}}_{\mu ,i}={\overline{n}}_{\mu }\displaystyle \sum _{m}{{ \mathcal P }}_{{im}}{{\rm{\Delta }}}_{\mu }{({{\boldsymbol{\theta }}}_{m})}^{1-2.5{s}_{\mathrm{eff}}}\end{eqnarray} \tag{ 19 }$$

with ${\overline{n}}_{\mu }={\overline{N}}_{\mu }/{{\rm{\Omega }}}_{\mathrm{cell}}$.

The choice of annular constraints rather than pixelated ones is mainly because magnification constraints are by far noisier than shear measurements, especially due to the local clustering noise (see Umetsu & Broadhurst 2008). An optimal choice of the resolution is to have per-element signal-to-noise ratio of order unity or above. This is satisfied by azimuthally averaging noisy count measurements, while it allows us to estimate the variance due to angular clustering at large clustercentric distances.

We use the CLASH weak-lensing magnification measurements obtained using flux-limited samples of red background galaxies as published in Umetsu et al. (2014). They measured the magnification effects in ${N}_{\mathrm{bin}}$ = 10 log-spaced circular annuli centered on the cluster. The radial bins range from 09 to 16' for all clusters, except 09 to 14' for RX J2248.7−4431 observed with ESO/WFI. Our magnification analysis begins at θ_min = 09, which is sufficiently large compared to the range of effective Einstein radii for our sample (Zitrin et al. 2015; Umetsu et al. 2016). The magnification profiles ${\{{n}_{\mu ,i}\}}_{i=1}^{{N}_{\mathrm{bin}}}$ used in the present work are presented in Figure 2 of Umetsu et al. (2014).

Here we briefly describe the magnification analysis performed in Umetsu et al. (2014). Their analysis was limited to the 24' × 24' region centered on the cluster. They accounted for the Poisson, intrinsic clustering, and additional systematic contributions to the total uncertainty σ_μ. The clustering noise term ${\sigma }_{\mu ,i}^{\mathrm{int}}$ was estimated in each circular annulus from the variance due to variations of the counts along the azimuthal direction. Besides, a positive tail of >νσ cells with ν = 2.5 was removed in each circular annulus by iterative σ clipping to reduce the bias due to intrinsic angular clustering of red galaxies. We checked that the clipping threshold chosen is sufficiently high compared to the maximum variations of the magnification signal due to halo ellipticity (typically, $| \delta \kappa (\theta )| /\langle \kappa (\theta )\rangle \lesssim 0.5$ for a projected halo axis ratio of ≥0.6). The Poisson noise term ${\sigma }_{\mu ,i}^{\mathrm{stat}}$ was estimated from the clipped mean counts in each annulus. The difference between the mean counts estimated with and without σ clipping was taken as a systematic error, ${\sigma }_{\mu ,i}^{\mathrm{sys}}=| {n}_{\mu ,i}^{(\nu )}-{n}_{\mu ,i}^{(\infty )}| /\nu$, where ${n}_{\nu ,i}^{(\nu )}$ and ${n}_{\mu ,i}^{(\infty )}$ denote the clipped and unclipped mean counts in the ith annulus, respectively. Finally, these errors were combined in quadrature as

$$\begin{eqnarray}&&{\sigma }_{\mu ,i}^{2}={({\sigma }_{\mu ,i}^{\mathrm{int}})}^{2}+{({\sigma }_{\mu ,i}^{\mathrm{stat}})}^{2}+{({\sigma }_{\mu ,i}^{\mathrm{sys}})}^{2}.\end{eqnarray} \tag{ 20 }$$

Our magnification bias measurements are stable and insensitive to the particular choice of ν because of the inclusion of the ${\sigma }_{\mu }^{\mathrm{sys}}$ term in the error analysis.

Masking of observed sky was accounted and corrected for using the method of (Umetsu et al. 2011b, Method B of Appendix A), which can be fully automated once the configuration parameters of SExtractor (Bertin & Arnouts 1996) are optimally tuned (Umetsu et al. 2011b, 2014). Chiu et al. (2016) adopted this method to estimate the masked area fraction in their magnification analysis and found that the SExtractor configuration of Umetsu et al. (2011b) is optimal for their data taken with Megacam on the Magellan Clay telescope.

The count normalization and slope parameters $({\overline{n}}_{\mu },{s}_{\mathrm{eff}})$ were estimated in the cluster outskirts (Umetsu et al. 2014).³³ The mask-corrected magnification bias profile ${b}_{\mu ,i}={n}_{\mu ,i}/{\overline{n}}_{\mu }$ is proportional to (1 − f_mask,back)/(1 − f_mask,i) ≡ 1 + Δf_mask,i, with f_mask,back estimated in the background region (see Umetsu et al. 2014, 2016). The effect of the mask correction is thus sensitive to the difference of the f_mask values, which is insensitive to the particular choice of the SExtractor configuration parameters. The typical variation of Δf_mask,i across the full radial range is ∼5% (see also Chiu et al. 2016), much smaller than the typical magnification signal $\delta {n}_{\mu }/{\overline{n}}_{\mu }\sim -0.3$ in the innermost bin. Accordingly, the systematic uncertainty on the mask correction is likely negligible.

The estimated values and errors for $\langle \beta {\rangle }_{\mu }$, ${\overline{n}}_{\mu }$, and s_eff are summarized in Table 4 of Umetsu et al. (2014). The values of s_eff span the range [0.11, 0.20] with a mean of $\langle {s}_{\mathrm{eff}}\rangle =0.153$ and a typical fractional uncertainty of 33% per cluster field. We marginalize over the calibration parameters ($\langle \beta {\rangle }_{\mu },{\overline{n}}_{\mu },{s}_{\mathrm{eff}}$) for each cluster in our joint likelihood analysis of shear and magnification (Section 3.5).

The log-likelihood function ${l}_{g}\equiv -\mathrm{ln}{{ \mathcal L }}_{g}$ for 2D shear data can be written (ignoring constant terms) as (Oguri et al. 2010)

$$\begin{eqnarray}{l}_{g}({\boldsymbol{\lambda }}) & = & \displaystyle \frac{1}{2}\displaystyle \sum _{m,n=1}^{{N}_{\mathrm{pix}}}\displaystyle \sum _{\alpha ,\beta =1}^{2}[{g}_{\alpha ,m}-{\widehat{g}}_{\alpha ,m}({\boldsymbol{\lambda }})]{({{ \mathcal W }}_{g})}_{\alpha \beta ,{mn}}\\ & & \times [{g}_{\beta ,n}-{\widehat{g}}_{\beta ,n}({\boldsymbol{\lambda }})],\end{eqnarray} \tag{ 21 }$$

where ${\widehat{g}}_{\alpha ,m}({\boldsymbol{\lambda }})$ is the theoretical expectation for g_α,m = g_α(${\boldsymbol{\theta }}$_m), and ${({{ \mathcal W }}_{g})}_{\alpha \beta ,{mn}}$ is the shear weight matrix,

$$\begin{eqnarray}&&{({{ \mathcal W }}_{g})}_{\alpha \beta ,{mn}}={M}_{m}{M}_{n}{({C}_{g}^{-1})}_{\alpha \beta ,{mn}}.\end{eqnarray} \tag{ 22 }$$

Here, M_m is a mask weight, defined such that M_m = 0 if the mth cell is masked out and M_m = 1 otherwise, and C_g is the shear covariance matrix. We account for contributions from the shape covariance ${C}_{g}^{\mathrm{shape}}$ and the cosmic covariance ${C}_{g}^{\mathrm{lss}}$ due to uncorrelated LSS projected along the line of sight as

$$\begin{eqnarray}&&{({C}_{g})}_{\alpha \beta ,{mn}}={({C}_{g}^{\mathrm{shape}})}_{\alpha \beta ,{mn}}+{({C}_{g}^{\mathrm{lss}})}_{\alpha \beta ,{mn}},\end{eqnarray} \tag{ 23 }$$

where ${({C}_{g}^{\mathrm{lss}})}_{\alpha \beta ,{mn}}={\xi }_{\alpha \beta }^{\mathrm{lss}}(| {{\boldsymbol{\theta }}}_{m}-{{\boldsymbol{\theta }}}_{n}| )$ with ${\xi }_{\alpha \beta }^{\mathrm{lss}}={\xi }_{\beta \alpha }^{\mathrm{lss}}$ (α, β = 1, 2) the cosmic shear correlation function (Hu & White 2001; Oguri et al. 2010). We compute the elements of the ${C}_{g}^{\mathrm{lss}}$ matrix for a given source population, 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). For each cluster, we use the effective mean source redshift (${\overline{z}}_{\mathrm{eff}};$ Table 3 of Umetsu et al. 2014) estimated with our multiband photometric redshifts.

Similarly, the log-likelihood function for magnification bias data ${l}_{\mu }\equiv -\mathrm{ln}{{ \mathcal L }}_{\mu }$ is written as (Umetsu et al. 2015)

$$\begin{eqnarray}&&{l}_{\mu }({\boldsymbol{\lambda }})=\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{{N}_{\mathrm{bin}}}[{n}_{\mu ,i}-{\widehat{n}}_{\mu ,i}({\boldsymbol{\lambda }})]{({{ \mathcal W }}_{\mu })}_{{ij}}[{n}_{\mu ,j}-{\widehat{n}}_{\mu ,j}({\boldsymbol{\lambda }})],\end{eqnarray} \tag{ 24 }$$

where ${\widehat{n}}_{\mu ,i}({\boldsymbol{\lambda }})$ is the theoretical expectation for the observed counts n_μ,i, and ${({{ \mathcal W }}_{\mu })}_{{ij}}$ is the magnification weight matrix, ${{ \mathcal W }}_{\mu }={C}_{\mu }^{-1}$, with C_μ the corresponding covariance matrix,

$$\begin{eqnarray}&&{({C}_{\mu })}_{{ij}}={\sigma }_{\mu }^{2}{\delta }_{{ij}}+{({C}_{\mu }^{\mathrm{lss}})}_{{ij}}.\end{eqnarray} \tag{ 25 }$$

The diagonal term in Equation (25) is responsible for the observational errors, and the bin-to-bin covariance matrix ${C}_{\mu }^{\mathrm{lss}}$ accounts for the cosmic noise contribution due to projected uncorrelated LSS, ${({C}_{\mu }^{\mathrm{lss}})}_{{ij}}={[(5{s}_{\mathrm{eff}}-2){\overline{n}}_{\mu }]}^{2}{({C}_{\kappa }^{\mathrm{lss}})}_{{ij}}$, where ${C}_{\kappa }^{\mathrm{lss}}$ is the cosmic convergence matrix (Umetsu et al. 2011a). We compute the elements of the ${C}_{\kappa }^{\mathrm{lss}}$ matrix for a given source redshift (${\overline{z}}_{\mathrm{eff}};$ Table 4 of Umetsu et al. 2014) in a similar manner to those of the ${C}_{g}^{\mathrm{lss}}$ matrix (Section 3.5.1). We evaluate the ${C}_{\mu }^{\mathrm{lss}}$ matrix by fixing the values of ${\overline{n}}_{\mu }$ and s_eff to the observed ones (Section 3.4.2).

The l_μ function sets azimuthally integrated constraints on the projected mass distribution and provides the otherwise unconstrained normalization of Σ(${\boldsymbol{R}}$) over a set of concentric annuli where magnification measurements are available. In this algorithm, no assumption is made about the azimuthal symmetry or isotropy of Σ(${\boldsymbol{R}}$). We use Monte Carlo integration to compute the projection matrix ${{ \mathcal P }}_{{im}}$ (Equation (18)) of size ${N}_{\mathrm{bin}}$ × ${N}_{\mathrm{pix}}$, which is needed to predict ${\{{\widehat{n}}_{\mu ,i}({\boldsymbol{\lambda }})\}}_{i\,=\,1}^{{N}_{\mathrm{bin}}}$ for a given model ${\boldsymbol{\lambda }}$ = (${\boldsymbol{s}}$, ${\boldsymbol{c}}$).

We account for the calibration uncertainty in the observational nuisance parameters,

$$\begin{eqnarray}&&{\boldsymbol{c}}=(\langle W{\rangle }_{g},{f}_{W,g},\langle W{\rangle }_{\mu },{\overline{n}}_{\mu },{s}_{\mathrm{eff}}).\end{eqnarray} \tag{ 26 }$$

To this end, we include in our joint-likelihood analysis Gaussian priors on ${\boldsymbol{c}}$ with mean values and errors estimated from data.

The log-posterior function $F({\boldsymbol{\lambda }})=-\mathrm{ln}P({\boldsymbol{\lambda }}| {\boldsymbol{d}})$ is expressed as a linear sum of the log-likelihood and prior terms. For each cluster, we find the global maximum of the joint posterior probability distribution over ${\boldsymbol{\lambda }}$ by minimizing F(${\boldsymbol{\lambda }}$) with respect to ${\boldsymbol{\lambda }}$. In the clumi-2D implementation of Umetsu et al. (2015), we use the conjugate-gradient method (Press et al. 1992) to find the global solution $\widehat{{\boldsymbol{\lambda }}}$. We employ an analytic expression for the gradient function ${\boldsymbol{\nabla }}$F(${\boldsymbol{\lambda }}$) obtained in the nonlinear, subcritical regime (Appendix B).

To quantify the reconstruction errors, we evaluate the Fisher matrix at ${\boldsymbol{\lambda }}=\widehat{{\boldsymbol{\lambda }}}$ as

$$\begin{eqnarray}&&{{ \mathcal F }}_{{pp}^{\prime} }=\left\langle \displaystyle \frac{{\partial }^{2}F({\boldsymbol{\lambda }})}{\partial {\lambda }_{p}\partial {\lambda }_{p^{\prime} }}\right\rangle {| }_{\widehat{{\boldsymbol{\lambda }}}},\end{eqnarray} \tag{ 27 }$$

where the angular brackets denote an ensemble average, and the indices (p, p') run over all model parameters ${\boldsymbol{\lambda }}$ = (${\boldsymbol{s}}$, ${\boldsymbol{c}}$). The error covariance matrix of the reconstructed parameters is obtained by $C={{ \mathcal F }}^{-1}$. We note that the reconstructed mass pixels are correlated primarily because the relation between the shear and convergence is nonlocal (Equation (2)). Additionally, the effects of spatial averaging (Equation (13)) and cosmic noise (Sections 3.5.1 and 3.5.2) produce a covariance between different pixels. In our analysis, the effects of correlated errors are modeled analytically (i.e., ${\xi }_{H},{\xi }^{\mathrm{lss}},{C}_{\kappa }^{\mathrm{lss}}$).

We model the radial mass distribution in galaxy clusters with a Navarro−Frenk−White (NFW) density profile, motivated by cosmological N-body simulations (Navarro et al. 1996, 1997) as well as by direct lensing observations (e.g., Newman et al. 2013; Okabe et al. 2013; Umetsu et al. 2014, 2016; Niikura et al. 2015; Okabe & Smith 2016; Umetsu & Diemer 2017). The radial dependence of the spherical NFW density profile is given by (Navarro et al. 1996)

$$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{{\rho }_{{\rm{s}}}}{(r/{r}_{{\rm{s}}}){(1+r/{r}_{{\rm{s}}})}^{2}},\end{eqnarray} \tag{ 28 }$$

with ρ_s the characteristic density parameter and r_s the characteristic scale radius at which the logarithmic slope of the density profile equals −2. We specify the spherical NFW model with the halo mass, M_200c, and the concentration parameter, c_200c ≡ r_200c/r_s.

The surface mass density Σ(R) as a function of projected clustercentric radius R is given by projecting ρ(r) along the line of sight. We employ an analytic expression given by Wright & Brainerd (2000) for the radial dependence of the projected NFW profile, ${\rm{\Sigma }}(R| {M}_{200{\rm{c}}},{c}_{200{\rm{c}}})$, which provides a good approximation for the projected halo model within a couple of virial radii (Oguri & Hamana 2011) and an excellent description of the projected mass distribution in clusters at R ≲ r_200m (Umetsu et al. 2016; Umetsu & Diemer 2017).

We follow the prescription given by Oguri et al. (2010, 2012) to construct an elliptical NFW (eNFW hereafter) model, which can be used to characterize the morphology of projected triaxial ellipsoids. To this end, we introduce the mass ellipticity, , in isodensity contours of the projected NFW profile ${\rm{\Sigma }}(R| {M}_{200{\rm{c}}},{c}_{200{\rm{c}}})$ as (Evans & Bridle 2009; Oguri et al. 2010; Umetsu et al. 2012; Medezinski et al. 2016)

$$\begin{eqnarray}&&{R}^{2}=X{{\prime} }^{2}(1-\epsilon )+Y{{\prime} }^{2}/(1-\epsilon ),\end{eqnarray} \tag{ 29 }$$

where our definition of the ellipticity is  = 1 − q_⊥, with q_⊥ ≤ 1 the projected minor-to-major axis ratio of isodensity contours, and we have chosen the coordinate system (X', Y') centered on the cluster halo, such that the X' axis is aligned with the major axis of the projected ellipse. Note that the isodensity area is πR², so that R represents the geometric mean radius of the isodensity ellipse. Accordingly, the M_200c and c_200c parameters in the eNFW model can be interpreted as respective spherical equivalent quantities.³⁵ In this work, we adopt the observer's coordinate system in which the X- and Y-axes are aligned with the west and north, respectively. With this coordinate system, the position angle (PA) of the projected major axis is measured east of north. An alternative definition for the projected ellipticity is $e=(1-{q}_{\perp }^{2})/(1+{q}_{\perp }^{2})$ (e.g., Evans & Bridle 2009).

We use a Bayesian Markov Chain Monte Carlo (MCMC) method to obtain an accurate inference of the eNFW parameters from our 2D weak-lensing data (Section 4.1). In this study and subsequent companion papers  (Chiu et al. 2018; Sereno et al. 2018), we perform model fitting to the 2D surface mass density data, rather than fitting directly to the combined shear and magnification data sets.³⁶ This allows consistency checks with existing codes used in previous work (e.g., Oguri et al. 2005; Umetsu & Broadhurst 2008; Morandi et al. 2011; Sereno & Umetsu 2011; Sereno et al. 2013; Umetsu et al. 2015; Sereno et al. 2017b) that also used weak-lensing Σ map data for 2D and 3D mass modeling. The mass maps obtained in this work are also useful for further studies of substructures in the context of the multiwavelength CLUMP-3D program.

The projected eNFW model is specified by four parameters, ${\boldsymbol{p}}$ = (M_200c, c_200c, q_⊥, PA). We use uniform prior distributions for the projected axis ratio and position angle in the range 0.1 ≤ q_⊥ ≤ 1 and −90° ≤ PA < 90° (Oguri et al. 2010). Following Umetsu et al. (2014, 2016), we assume log-uniform priors for M_200c and c_200c in the range $0.1\leqslant {M}_{200{\rm{c}}}/({10}^{15}\,{M}_{\odot }\,{h}^{-1})\leqslant 10$ and 0.1 ≤ c_200c ≤ 10. As found by Sereno et al. (2015) and Umetsu et al. (2016), the mass and concentration estimates for the CLASH sample are not sensitive to the choice of the priors, thanks to the deep high-quality weak-lensing observations. The χ² function for our observations is

$$\begin{eqnarray}&&{\chi }^{2}({\boldsymbol{p}})=\displaystyle \sum _{m,n=1}^{{N}_{\mathrm{pix}}}[{{\rm{\Sigma }}}_{m}-{\widehat{{\rm{\Sigma }}}}_{m}({\boldsymbol{p}})]{({C}^{-1})}_{{mn}}[{{\rm{\Sigma }}}_{n}-{\widehat{{\rm{\Sigma }}}}_{n}({\boldsymbol{p}})],\end{eqnarray} \tag{ 30 }$$

where ${\widehat{{\rm{\Sigma }}}}_{m}({\boldsymbol{p}})$ denotes the surface mass density at the grid position (X_m, Y_m) predicted by the model ${\boldsymbol{p}}$. For all clusters, we restrict the fitting to a square region of side 4 $\mathrm{Mpc}\,{h}^{-1}$ centered on the cluster, where half the side length corresponds to the typical r_200m radius of CLASH clusters. This is to minimize the impact of the 2-halo term and local substructures that are abundant in cluster outskirts, which otherwise can lead to bias in cluster mass estimates (Meneghetti et al. 2010; Becker & Kravtsov 2011; Rasia et al. 2012). Similarly, our previous CLASH studies performed fitting to azimuthally averaged lensing profiles by restricting the fitting range to R ≤ 2 Mpc h⁻¹ (Umetsu et al. 2014; Merten et al. 2015; Umetsu et al. 2016).

We also fit the data with a spherical NFW halo (q_⊥ = 1) using the same log-uniform priors on M_200c and c_200c in order to examine the consistency of our results with those of Umetsu et al. (2014).

In Table 1, we list the marginalized posterior constraints on each of the model parameters ${\boldsymbol{p}}$ = (M_200c, c_200c, q_⊥, PA) for 20 individual clusters of our sample. In this work, we employ the robust biweight estimators of Beers et al. (1990) for the center location (C_BI) and scale (S_BI) of the marginalized 1D posterior distributions (e.g., Stanford et al. 1998; Sereno & Umetsu 2011; Biviano et al. 2013). The biweight estimator is insensitive to and stable against outliers, as it assigns higher weight to points that are closer to the center of the distribution (Beers et al. 1990). In the table, we also report, for each cluster, the minimum χ² value per degree of freedom (dof) as an indicator of goodness of fit. Here, the number of dof is defined as the difference between the number of mass pixels within the fitting region and the number of free parameters. We find that the minimum χ²/dof values for our sample range from 0.41 (MACS J0744.9+3927) to 1.32 (RX J1347.5−1145), with a median of 0.80. This indicates that complex morphologies in the projected cluster mass distribution (e.g., substructures and deviations from elliptical isodensity contours) are not statistically significant in individual clusters, and the eNFW model provides an adequate description of our 2D weak-lensing data.

In Appendix C, we show, for each of the clusters, the 2D marginalized posterior distributions of the eNFW parameters, with the contours enclosing 68% and 95% of the posterior probability (Figure 11). For each parameter, we also present the 1D marginalized distribution, in which the C_BI location is marked with a vertical line. We see that the mass and concentration parameters are well constrained by the data, except for MACS J1931.8−2635 and MACS J0429.6−0253, for which the posterior distribution on c_200c is largely informed by the prior, in the sense that the likelihood extends outside of the prior range. Overall, the 2D weak-lensing constraints on the halo shape parameters (q_⊥, PA) are not strongly degenerate with halo mass and concentration, as found by Oguri et al. (2010).

We examine here the ensemble distribution of projected axis ratios using the results from the Bayesian inference of the eNFW parameters. Figure 3 shows the resulting distribution function for our sample of 20 clusters (blue squares), constructed from posterior point estimates (C_BI) of the q_⊥ parameter (Table 1). The median axis ratio for our sample, measured within a radial scale of 2 $\mathrm{Mpc}\,{h}^{-1}$, is $\langle {q}_{\perp }\rangle \,=0.67\pm 0.07$, or $\langle \epsilon \rangle =0.33\pm 0.07$ and $\langle e\rangle =0.38\pm 0.08$ in terms of the projected halo ellipticity (Table 2), where the errors were estimated by bootstrap resampling the cluster sample. To check at which radius the constraint on the halo ellipticity effectively comes from, we calculate the mass-weighted, projected clustercentric radius R_eff for our sample, averaged within the central 2 $\mathrm{Mpc}\,{h}^{-1}$ region. Using the best-fit NFW model based on the stacked lensing analysis by Umetsu et al. (2016), we find R_eff ≃ 0.89 $\mathrm{Mpc}\,{h}^{-1}$.

Zoom In Zoom Out Reset image size

Table 2.  Weak-lensing Halo Ellipticity Measurements

Sample	N	$\langle \epsilon \rangle$	$\langle e\rangle$
Full sample	20	0.33 ± 0.07	0.38 ± 0.08
X-ray selected	16	0.28 ± 0.07	0.33 ± 0.06
High-magnification	4	0.45 ± 0.11	0.53 ± 0.14

Note. Median values and 1σ errors of weak-lensing cluster shape measurements derived for the full sample, the X-ray-selected subsample, and the high-magnification subsample. We adopt the halo ellipticity defined in two ways:  = 1 − q_⊥ and $e=(1-{q}_{\perp }^{2})/(1+{q}_{\perp }^{2})$ with ${q}_{\perp }\leqslant 1$ the projected halo axis ratio.

Download table as:  ASCIITypeset image

For a consistency check, we also create a composite probability distribution function (PDF) of q_⊥ by stacking the marginal posterior distributions (Figure 11) for all clusters in our sample. This yields a median axis ratio of $\langle {q}_{\perp }\rangle \,=0.65\pm 0.05$, or $\langle \epsilon \rangle =0.35\pm 0.05$ and $\langle e\rangle =0.41\pm 0.07$, where the errors are based on 50,000 random samples that are drawn from the posterior distributions for individual clusters. These results are in good agreement with those from the posterior point estimates.

It is important to note that the shape of the clusters in the CLASH sample is expected to be rounder on average due to a high fraction of relaxed clusters (Section 3.1; Meneghetti et al. 2014). Accordingly, there could be a bias toward higher values of the projected axis ratio in our sample. For our X-ray-selected subsample (excluding the four high-magnification clusters with complex morphologies), we find the median axis ratio to be $\langle {q}_{\perp }\rangle =0.72\pm 0.07$ from the posterior point estimates and $\langle {q}_{\perp }\rangle =0.68\pm 0.06$ from the stacked PDF. Although the amplitude of the shift is not statistically significant, the direction of the shift is consistent with the effect of the CLASH X-ray-selection function (Section 3.1).

Donahue et al. (2016) found the typical axis ratio for all 25 CLASH clusters in the Chandra X-ray brightness distribution to be 0.88 ± 0.06 within an aperture radius of 350 $\mathrm{kpc}\,{h}^{-1}$ (∼1.1r_2500c), consistent with the value of 0.90 ± 0.06 inferred from the SZE maps observed with Bolocam operating at 140 GHz (Sayers et al. 2013; Czakon et al. 2015). This comparison for the CLASH sample, albeit at different radial scales, indicates that the shape of the projected mass distribution as measured from weak lensing is more elongated than the gas distribution. This is qualitatively consistent with the theoretical expectation that the intracluster gas in hydrostatic equilibrium is rounder than the underlying matter distribution (Lee & Suto 2003). However, we note that the values of cluster morphological parameters for the X-ray and SZE maps could be different if measured over larger radial scales (see the discussion in Section 5.2.2). On the other hand, our ground-based weak-lensing data do not sufficiently resolve morphological structures within a radial scale of 350 $\mathrm{kpc}\,{h}^{-1}$ (∼16 at the median sample redshift of z = 0.377) as they are limited by the small number density of background galaxies (Umetsu et al. 2014).

These results can be compared with predictions from ΛCDM N-body simulations, in which DM halos are often modeled as triaxial ellipsoids (e.g., Jing & Suto 2000; Bett et al. 2007; Despali et al. 2014; Bonamigo et al. 2015; Despali et al. 2016; Vega-Ferrero et al. 2017).³⁷ Here we restrict the comparison to predictions for the shape of the DM distribution measured within the virial radius, which is close to the maximum fitting radius of our analysis (2 $\mathrm{Mpc}\,{h}^{-1}$ ∼ r_200m ∼ 1.1r_vir). In particular, the triaxial model of Bonamigo et al. (2015) is of particular interest because they extended the original work of Jing & Suto (2000) to a wider mass range with higher precision. Bonamigo et al. (2015) selected and analyzed those halos for which the offset between the center of mass and geometrical center of the ellipsoid is less than 5% of their virial radius. The fraction of selected clusters is ∼50% at z = 0 (see Figure 2 of Bonamigo et al. 2015). On the other hand, applying stringent relaxation selection criteria typically results in a much smaller fraction of selected halos (e.g., ∼15% at z = 0.25 as found by Meneghetti et al. 2014).

In Figure 3, we show the theoretical expectation P(q_⊥) for our composite clusters (red shaded histogram) obtained using the fitting formula of Bonamigo et al. (2015), which describes the intrinsic distribution of triaxial axis ratios for their DM halos. Here, the predicted distribution has been constructed as follows: first, we compute for each cluster the intrinsic axis ratio distribution as a function of halo virial mass ${M}_{\mathrm{vir}}$ and redshift (Bonamigo et al. 2015) by accounting for the uncertainty in the mass determination. Next, we construct the distribution of projected axis ratios by projecting triaxial halos onto the sky plane assuming random orientations. The q_⊥ distribution predicted for each cluster is then convolved with a Gaussian of width equal to the observed uncertainty σ(q_⊥). The composite PDF is finally obtained by adding the renormalized distributions of all 20 clusters. In Figure 3, we also show the theoretical distribution P(q_⊥) without the Gaussian convolution (red open histogram). Because of the projection effect, the overall shape of the distribution of the projected axis ratio is broader and shifted to values (rounder) higher than those of the intrinsic minor-to-major axis ratio (Suto et al. 2016). The expected median value is $\langle {q}_{\perp }\rangle =0.59$ in both cases with and without the Gaussian convolution, in agreement with our measurement within the uncertainty.We reiterate that our sample is expected to contain a high fraction of relaxed clusters (∼60% in our full composite sample, compared to ∼70% in the X-ray-selected subsample; see Section 3.1), and hence the average projected axis ratio is likely biased high to some degree. In particular, if we restrict our comparison to the X-ray-selected subsample, the observed median projected axis ratio (Table 2) is higher at the 1.8σ level than predicted from the triaxial model of Bonamigo et al. (2015).

We plot in Figure 4 the projected axis ratio q_⊥ as a function of the virial mass ${M}_{\mathrm{vir}}$ for our sample of 20 CLASH clusters (black squares). The blue triangle in Figure 4 represents the (unweighted) median values for the sample, $\langle {q}_{\perp }\rangle =0.67\,\pm 0.07$ at $\langle {M}_{\mathrm{vir}}\rangle =(15.2\pm 2.8)\times {10}^{14}{M}_{\odot }\,{h}_{70}^{-1}$. Again, our CLASH weak-lensing measurements are in good agreement with ΛCDM expectations at the median sample redshift of z = 0.377. All 20 CLASH clusters lie within the 2σ distribution predicted with the triaxial model of Bonamigo et al. (2015) assuming random orientations of the clusters.

Zoom In Zoom Out Reset image size

Suto et al. (2016) studied the mass and radial dependence and the redshift evolution of the non-sphericity of cluster-size halos using DM-only simulations. They found that the average 3D minor-to-major axis ratio of simulated halos has a strong radial dependence as a function of enclosed mass, indicating that the internal structure of halos deviates from a self-similar geometry of concentric ellipsoidal surfaces. Suto et al. (2016) thus constructed the PDF of the projected axis ratio, P(q_⊥), directly from the projected density distributions of simulated halos, without involving triaxial modeling.³⁸ Using their fitting formula for the PDF at z = 0.4, the mean, standard deviation, and median are found to be 0.57, 0.17, and 0.58, respectively. The predicted median agrees with our full-sample result at the 1.3σ level. On the other hand, when halos in multicomponent systems are excluded from their analysis, the mean, standard deviation, and median of the PDF are 0.59, 0.16, and 0.60, respectively. The observed median axis ratio for the CLASH X-ray-selected subsample is 1.6σ above this predicted median.

More recently, it was shown by Suto et al. (2017) that baryonic physics operating in the cluster central region, such as cooling and feedback effects, can have a substantial impact on the non-sphericity of cluster halos up to half the virial radius, even though these baryonic effects have little impact on the spherically averaged DM density profile. They found that the DM distribution becomes more spherical, depending on the distance from the cluster center, when the effects of baryons are included. Since our sample comprises highly relaxed and highly disturbed clusters (Section 3.1), a more quantitative comparison with theoretical expectations would require a detailed modeling of baryonic physics by accounting for the selection function. In fact, Figure 4 shows a slight tendency for lower-mass clusters to have projected axis ratios that are higher than the predicted distribution. This tendency is qualitatively consistent with the combination of the CLASH selection function and the baryonic effects. Nevertheless, our analysis is currently limited by the small number of clusters.

Our CLASH lensing results can be compared to the weak-lensing measurements obtained by the LoCuSS collaboration (Oguri et al. 2010), who performed a 2D shear-fitting analysis on a sample of 25 X-ray-luminous clusters at $\langle z\rangle \sim 0.2$. Oguri et al. (2010) modeled the projected mass distribution in each individual cluster assuming a single eNFW halo, as done in our work. They fitted an eNFW profile to the grid-averaged 2D reduced shear field. For a subset of 18 clusters that were specifically chosen to give good ellipticity constraints, they found a ∼7σ detection of the mean cluster ellipticity, $\langle \epsilon \rangle =0.46\pm 0.04$, for their clusters with $\langle {M}_{\mathrm{vir}}\rangle \sim 10\,\times {10}^{14}{M}_{\odot }\,{h}_{70}^{-1}$. Their mean ellipticity is higher than, but consistent within the errors with, our full-sample measurement (Table 2). This difference is not statistically significant, but could be due in part to the CLASH X-ray selection function (Section 3.1).

Similarly, Oguri et al. (2012) found from their 2D shear analysis of 25 strong-lensing-selected clusters that the stacked cluster ellipticity is nearly constant, $\langle \epsilon \rangle \sim 0.45$, with cluster radius within their errors.

On group/cluster scales, several authors have constrained the halo ellipticity from stacked weak-lensing measurements (Evans & Bridle 2009; Clampitt & Jain 2016; van Uitert et al. 2017; Shin et al. 2018) using quadrupole shear estimators and their variants (e.g., Natarajan & Refregier 2000; Adhikari et al. 2015; Clampitt & Jain 2016). Their stacked weak-lensing measurements span a range of mean q_⊥ values from ${0.48}_{-0.09}^{+0.14}$ (Evans & Bridle 2009) to ∼0.78 (Clampitt & Jain 2016), broadly consistent with ΛCDM predictions. We note that, in their approach, one probes halo quadrupoles with respect to the major axis of the light distribution (e.g., BCGs) chosen as a reference orientation. On the other hand, we have directly measured the shape and orientation of individual cluster halos from deep high-quality weak-lensing data (Umetsu et al. 2014). We present in Section 5.3.2 our stacked quadrupole shear measurement for our full sample of 20 CLASH clusters using the shape of the Chandra X-ray brightness distribution (see Section 5.3.1) as a reference orientation.

Donahue et al. (2015, 2016) presented a detailed study of morphologies and alignments of BCGs, intracluster gas, and mass at small cluster radii for the CLASH sample. Donahue et al. (2015) measured PAs for all 20 X-ray-selected CLASH clusters from the rest-frame ultraviolet (UV) and near-infrared (NIR) light distributions of the BCGs using CLASH HST data (see their Table 3). Similarly, Donahue et al. (2016) examined the morphological properties of all 25 CLASH clusters with Chandra X-ray data (Donahue et al. 2014) and compared them with those inferred from Bolocam SZE maps of the intracluster gas (Sayers et al. 2013; Czakon et al. 2015) and those from mass maps based on HST strong and weak lensing (hereafter HST-GL; Zitrin et al. 2015). The lensing maps used are based on parametric modeling assuming that light traces mass for the galaxy-scale mass components (see Meneghetti et al. 2017). They found a strong correlation between PAs as measured from the X-ray, SZE, and HST-GL maps inside a consistent metric aperture of radius 350 $\mathrm{kpc}\,{h}^{-1}$ (∼0.2r_vir for these clusters). They also found a strong alignment of the cluster shapes at this scale, as measured from the X-ray, SZE, and HST-GL maps, with that of the NIR BCG light at 10 kpc scales. In particular, they obtained a median misalignment angle of $| {\rm{\Delta }}\mathrm{PA}| \simeq 11^\circ$ between the BCG and X-ray orientations for the 20 X-ray-selected CLASH clusters.

Table 3.  Weak-lensing Halo Misalignment Statistics

Data Sets	N	$\langle {\rm{\Delta }}\mathrm{PA}\rangle$	$\langle | {\rm{\Delta }}\mathrm{PA}| \rangle$	p-value	pBin
		(degrees)	(degrees)
BCG/WL	16	−25 ± 14	34 ± 13	1.5 × 10−1	1.8 × 10−1
X-ray/WL	20	−3 ± 9	21 ± 7	4.1 × 10−3	4.6 × 10−3
SZE/WL	18	7 ± 17	42 ± 10	3.8 × 10−1	1.9 × 10−1
HST-GL/WL	20	−12 ± 6	16 ± 4	3.1 × 10−4	1.1 × 10−3

Note. The basic statistics of misalignment angles are listed. The results here are based on the point estimates of weak-lensing PAs from individual cluster posterior distributions (Figure 11). Column 1: combination of data sets. Column 2: number of clusters in the overlapping sample. Column 3: median of the distribution of misalignment angles, ΔPA  ∈ [−90°, 90°). Column 4: median of the distribution of absolute misalignment angles, $| {\rm{\Delta }}\mathrm{PA}| \in [0^\circ ,90^\circ ]$. Column 5: probability of finding the value $\langle | {\rm{\Delta }}\mathrm{PA}| \rangle$ smaller than or equal to the observed value when the null hypothesis of a uniform distribution in [0°, 90°] is true. Column 6: probability, obtained with the binomial test, that the observed distribution in $| {\rm{\Delta }}\mathrm{PA}|$ has random orientations (West et al. 2017).

Download table as:  ASCIITypeset image

Now we compare the PAs determined from our 2D weak-lensing analysis (WL) to those from three baryonic tracers, namely the BCGs, X-ray maps, and SZE maps, as well as from the HST-GL maps. To this end, we use the BCG PAs derived from the HST NIR images (Table 3 of Donahue et al. 2015), and the X-ray, SZE, and HST-GL PAs measured within an aperture radius of 350 $\mathrm{kpc}\,{h}^{-1}$ (Tables 3, 6, and 5 of Donahue et al. 2016). To be conservative, we estimate errors on the BCG PAs from differences between the HST UV and NIR measurements of Donahue et al. (2015). There are 16, 20, 18, and 20 clusters available for BCG/WL, X-ray/WL, SZE/WL, and HST-GL/WL comparisons, respectively. The typical uncertainty in the weak-lensing PA measurements is ∼20° per cluster, comparable to that in the SZE measurements, while those in the BCG, X-ray, and HST-GL measurements are much smaller, ∼10°, ∼4°, and ∼4° per cluster, respectively.

In Figure 5, the PAs of the clusters determined from our weak-lensing analysis are plotted against those from the BCGs, X-ray maps, SZE maps, and HST-GL maps. To quantify the degree of correlation between PAs of different tracers, we calculate the Spearman rank correlation coefficient r_P and the corresponding probability (p) for the null hypothesis of random orientations. The test indicates a similarly good correlation in all four of these comparisons, where a low probability indicates high significance of correlation: r_S = 0.69 (p = 2.9 × 10⁻³) for BCG/WL, r_S = 0.66 (p = 1.7 × 10⁻³) for X-ray/WL, r_S = 0.58 (p = 1.1 × 10⁻²) for SZE/WL PAs, and r_S = 0.74 (p = 1.7 × 10⁻⁴) for HST-GL/WL PAs.

Zoom In Zoom Out Reset image size

Next, we quantify the alignment of the different components in the clusters by constructing the probability distributions of the absolute misalignment angles $| {\rm{\Delta }}\mathrm{PA}|$ (e.g., Faltenbacher et al. 2009; Zhang et al. 2013; West et al. 2017) with respect to the weak-lensing halo shape. In the absence of any alignment, $| {\rm{\Delta }}\mathrm{PA}|$ follows a uniform distribution between 0° and 90°, with a mean (median) of $\langle | {\rm{\Delta }}\mathrm{PA}| \rangle =45^\circ$. Values of $\langle | {\rm{\Delta }}\mathrm{PA}| \rangle \lt 45^\circ$ indicate that the PAs of the two distributions are, on average, aligned parallel with each other.

We show in Figure 6 the histogram distributions of $| {\rm{\Delta }}\mathrm{PA}|$ between the BCG/WL, X-ray/WL, SZE/WL, and HST-GL/WL orientations derived for our sample. Table 3 lists, for the respective comparisons, the median values of $| {\rm{\Delta }}\mathrm{PA}|$ and the corresponding significance probabilities based on the posterior point estimates of weak-lensing PAs. As evident from Table 3 and Figure 6, among the three baryonic tracers studied here, the X-ray morphology is best aligned with the weak-lensing mass distribution. For this comparison, we find a median misalignment angle of $\langle | {\rm{\Delta }}\mathrm{PA}| \rangle =21^\circ \pm 7^\circ$ (Table 3), corresponding to 3.6σ significance with respect to the null hypothesis of random orientations. However, these median values are biased high relative to their intrinsic values because they are estimated from the noisy distributions. If the intrinsic PA difference ΔPA =PA(X-ray) − PA(WL) follows a (truncated) Gaussian distribution with a zero mean and a dispersion σ_int, the median of the intrinsic distribution of $| {\rm{\Delta }}\mathrm{PA}|$ is ≃σ_int/1.483. In the presence of noise, the apparent dispersion from observations is increased. We correct for the effect of this noise bias assuming that the PA errors follow a Gaussian distribution. Adopting the typical uncertainties in the PA measurements, we find the bias-corrected median misalignment angle to be $\langle | {\rm{\Delta }}\mathrm{PA}| \rangle \,=16^\circ \pm 7^\circ$.

Zoom In Zoom Out Reset image size

For the BCG/WL and SZE/WL comparisons, we find weaker alignment signals (<45°), which are consistent with a null detection within the errors. A weak constraint on the SZE/WL alignment is in line with expectations accounting for the large errors in both measurements. On the other hand, the weak signal in the BCG/WL alignment appears to be contradictory to the high degree of correlation found between the BCG/WL PAs (Figure 5).³⁹ This is largely because the distribution of ΔPA = PA(BCG) − PA(WL) is not symmetric about ΔPA =0° (Figure 5), with a median offset of $\langle {\rm{\Delta }}\mathrm{PA}\rangle =-25^\circ \pm 14^\circ$. For the X-ray/WL, SZE/WL, and HST-GL/WL comparisons, the median offsets are found to be $\langle {\rm{\Delta }}\mathrm{PA}\rangle =-3^\circ \pm 9^\circ$, 7° ± 17°, and −12° ± 6°, respectively (Table 3).

On the other hand, a strong alignment is found between the WL and HST-GL shapes at different radial scales, with a median misalignment angle of $\langle | {\rm{\Delta }}\mathrm{PA}| \rangle =16^\circ \pm 4^\circ$. This corresponds to a significance level of 7.7σ with respect to the null hypothesis. Applying the noise-bias correction yields $\langle | {\rm{\Delta }}\mathrm{PA}| \rangle =8^\circ \pm 4^\circ$. This is consistent with the results of Despali et al. (2016), who found from N-body simulations that, for cluster-scale halos, the innermost and outermost mass ellipsoids are aligned with each other within 10°. We emphasize that this strong alignment signal has been found despite using two independent data sets (HST versus ground-based observations) and substantially different modeling methods.

Here we present a complementary quadrupole shear analysis to test the consistency of our cluster ellipticity measurements for our sample of 20 CLASH clusters. To this end, we employ the Cartesian estimators of Clampitt & Jain (2016) that null the purely tangential, monopole shear contribution. Specifically, we measure the stacked quadrupole shear signal with respect to a coordinate system with the x-axis aligned with the X-ray major axis of each cluster. We adopt the same sign convention for the Cartesian g₁ and g₂ components as defined in Clampitt & Jain (2016; see their Figure 1) and use θ to denote the azimuthal angle relative to the x-axis. Following Clampitt & Jain (2016), we group together the first and second shear components of background galaxies in the regions where cos 4θ and sin 4θ have the same sign (Figure 7), respectively, and define the following estimator:

$$\begin{eqnarray}&&{\rm{\Delta }}{{\rm{\Sigma }}}_{\alpha }^{(s)}(R)={{\rm{\Sigma }}}_{{\rm{c}}}\left[\displaystyle \sum _{k}{w}_{(k)}\,{g}_{\alpha ,k}\right]{\left[\displaystyle \sum _{k}{w}_{(k)}\right]}^{-1},\end{eqnarray} \tag{ 31 }$$

where we have introduced the notation in analogy to the tangential shear, which probes the differential surface mass density ΔΣ (e.g., Umetsu et al. 2014). Here, ${{\rm{\Sigma }}}_{{\rm{c}}}\propto \langle \beta {\rangle }_{g}^{-1}{D}_{{\rm{l}}}^{-1}$ is the source-averaged critical surface density for a given cluster (Section 2.1), w_(k) is the statistical weight for the kth background galaxy (Section 3.3.1), and k runs over all background galaxies that fall in the specified bin, different for each shear component α and sign s (Clampitt & Jain 2016): α = 1, s = −, −π/8 ≤ θ_j < π/8; α = 1, s = +, π/8 ≤ θ_j< 3π/8; α = 2, s = −, 0 ≤ θ_j < π/4; α = 2, s = +, π/4 ≤ θ_j < π/2. For each case, the summation in Equation (31) also includes background galaxies lying in symmetrical regions shifted by π/2, π, and 3π/2, as illustrated in Figure 7.

Zoom In Zoom Out Reset image size

In this work, we first measure for each cluster the quadrupole shear profiles ${\rm{\Delta }}{{\rm{\Sigma }}}_{\alpha }^{(s)}(R)$ from the background-selected shear catalog according to Equation (31), and then stack all clusters together by

$$\begin{eqnarray} \langle\langle {\rm{\Delta }}{{\rm{\Sigma }}}_{\alpha }^{(s)}(R) \rangle\rangle & = & \left[\displaystyle \sum _{n}{W}_{n}\,{\rm{\Delta }}{{\rm{\Sigma }}}_{\alpha ,n}^{(s)}(R)\right]{\left[\displaystyle \sum _{n}{W}_{n}\right]}^{-1},\\ {W}_{n} & = & 1/{[{\sigma }_{\alpha ,n}^{(s)}(R)]}^{2},\end{eqnarray} \tag{ 32 }$$

where $\langle\langle ... \rangle\rangle$ denotes the sensitivity-weighted average over the cluster sample (Umetsu et al. 2014, 2016), n runs over all 20 clusters in our sample, and ${\sigma }_{\alpha ,n}^{(s)}(R)$ is the uncertainty of ${\rm{\Delta }}{{\rm{\Sigma }}}_{\alpha ,n}^{(s)}(R)$ estimated from bootstrap resampling of the background galaxies. We estimate the errors of the stacked $\langle\langle {\rm{\Delta }}{{\rm{\Sigma }}}_{\alpha }^{(s)}(R) \rangle\rangle$ profiles by bootstrap resampling the cluster sample. The resulting stacked quadrupole profiles for our cluster sample are shown in Figure 8.

Zoom In Zoom Out Reset image size

Clampitt & Jain (2016) modeled the quadrupole shear signal using a multipole expansion of the surface mass density of elliptical halos (Adhikari et al. 2015). However, this method can only be applied to the case with a small halo ellipticity, so that the higher-order terms can be safely ignored.

In order to accurately model the observed signal and to make a direct comparison with our 2D cluster ellipticity measurements (Section 5.2), we forward-model the stacked quadrupole shear profiles by assuming an eNFW halo with the major axis aligned with the X-ray major axis, that is, ΔPA = PA(X-ray) − PA(WL) = 0°. Therefore, any misalignment $| {\rm{\Delta }}{\rm{P}}{\rm{A}}| \gt {0}^{\circ }$ will lead to dilution of the quadrupole signal and hence underestimation of the halo ellipticity. We use a Bayesian MCMC approach (see Section 4.2.2) to simultaneously fit an eNFW model to the four stacked quadrupole profiles, namely $\langle\langle {\rm{\Delta }}{{\rm{\Sigma }}}_{1}^{(-)} \rangle\rangle$, $\langle\langle {\rm{\Delta }}{{\rm{\Sigma }}}_{1}^{(+)} \rangle\rangle$, $\langle\langle {\rm{\Delta }}{{\rm{\Sigma }}}_{2}^{(-)} \rangle\rangle$, and $\langle\langle {\rm{\Delta }}{{\rm{\Sigma }}}_{2}^{(+)} \rangle\rangle$, each measured in four radial bins (Figure 8). We use a uniform prior distribution for the projected axis ratio in the range 0.1 ≤ q_⊥ ≤ 1 (Section 5.2). We marginalize over the mass and concentration parameters using Gaussian priors of ${M}_{200{\rm{c}}}=(10.9\pm 0.7)\times {10}^{14}\,{M}_{\odot }\,{h}^{-1}$ and c_200c = 3.3 ± 0.2 based on the joint strong-lensing, weak-lensing shear and magnification analysis of Umetsu et al. (2016). Marginalized posterior constraints (C_BI ± S_BI) on the projected axis ratio are obtained as q_⊥ = 0.67 ± 0.10 (Figure 8), or  = 0.33 ± 0.10 and e = 0.38 ± 0.12 in terms of the halo ellipticity. These are in excellent agreement with the results from the 2D weak-lensing analysis of 20 individual clusters (Table 2), supporting the robustness of our results and indicating that the effect of dilution due to X-ray/WL misalignment is not significant for our cluster sample.

Okumura et al. (2009) measured intrinsic ellipticity correlation functions for a large sample of luminous red galaxies (LRGs) in the redshift range 0.16–0.47 selected from Data Release 6 (DR6) of the Sloan Digital Sky Survey (SDSS), finding a clear signal up to a scale of 30 $\mathrm{Mpc}\,{h}^{-1}$. To model the observed ellipticity correlation, they populated galaxies into DM halos in cosmological N-body simulations using a halo occupation distribution approach. In this context, the fraction of synthetic central LRGs is 93.7%, and they are hosted by cluster-scale DM halos with $M\gt 8\times {10}^{13}\,{M}_{\odot }$. The ellipticity correlation is predicted to have an amplitude that is about four times higher than their measurement when assuming a perfect alignment between the central LRGs and their host DM halos. Assuming a Gaussian misalignment between the major axes of central LRGs and host halos, they found a misalignment dispersion of $\sigma ({\rm{\Delta }}\mathrm{PA})={35.4}_{-3.3}^{+4.0}$ degrees, or an absolute median of $\langle | {\rm{\Delta }}\mathrm{PA}| \rangle ={24}_{-2}^{+3}$ degrees. A similar constraint was derived by Okumura & Jing (2009) from the gravitational shear–intrinsic ellipticity correlation function measured using the LRG sample. Their constraints on the misalignment angle are in agreement with our BCG/WL results within the errors (Table 3).

Huang et al. (2016) studied central galaxy alignments with respect to the spatial distribution of member galaxies using the red-sequence Matched-filter Probabilistic Percolation (redMaPPer) cluster catalog based on the SDSS DR8, finding an absolute mean misalignment angle of ∼35°, corresponding to an absolute median angle of $\langle | {\rm{\Delta }}\mathrm{PA}| \rangle \sim 30^\circ$. This is in agreement with our estimate for the median misalignment angle between the BCG/WL major axes (Table 3) if assuming that cluster member galaxies are an unbiased probe of the underlying mass distribution.

van Uitert et al. (2017) presented a stacked quadrupole shear analysis of ∼2600 galaxy groups selected from the Galaxy And Mass Assembly survey using the galaxy shear catalog from the Kilo Degree Survey. On small scales (<250 kpc), they found the major axis of the BCG to be the best proxy of the orientation of the underlying mass distribution. On larger scales, a much weaker correlation was found between the orientations of the BCG and the mass distribution, while the distribution of member galaxies provides a better proxy for the orientation of the overall mass distribution in groups.

More recently, Shin et al. (2018) studied the stacked halo ellipticity for a sample of ∼10⁴ SDSS redMaPPer clusters. Stacking the quadrupole shear signal along the major axis of the cluster member distribution, they found a mean axis ratio of 0.558 ± 0.086 ± 0.026 (statistical followed by systematic uncertainty). This agrees well with the mean axis ratio of the member distribution, 0.573 ± 0.002 ± 0.039, indicating that cluster galaxies trace the shape of the cluster mass distribution within their errors. On the other hand, they found an rms misalignment angle of 30° ± 10° between the central galaxy and the cluster mass distribution.