SCALING RELATIONS AND OVERABUNDANCE OF MASSIVE CLUSTERS AT z ≳ 1 FROM WEAK-LENSING STUDIES WITH THE HUBBLE SPACE TELESCOPE^*

Our reduction started with the FLT images, which are the products of the standard STScI CALACS pipeline (Hack et al. 2003). We used the "apsis" pipeline (Blakeslee et al. 2003) to determine shift and rotation, correct geometric distortion, remove cosmic rays, and create final mosaic images. The determination of the shift and rotation between different pointings is the most critical step among these for a weak-lensing analysis because any potential misalignment induces a coherent distortion of object shapes, mimicking gravitational lensing signals. The "apsis" pipeline uses high signal-to-noise (S/N) astronomical objects iteratively to measure and refine the alignment through the "match" program developed by Richmond (2002). We find that the average uncertainties for shift and rotation are ∼0.02 pixels and ∼02, respectively, and thus the systematics, if any, due to the alignment errors is negligible. The Lanczos3 (windowed sinc function) kernel was used for drizzling (Fruchter & Hook 2002) with the native 005 pixel scale of ACS/WFC.

A detection image was created for each cluster by weight averaging the two passband images via apsis. Then, SExtractor (Bertin & Arnouts 1996) was run in dual image mode so that objects were identified from the detection image while photometry was measured on an individual filter. This provides common isophotal areas to both filters, which enables a robust object color measurement. We filtered the detection image with a Gaussian kernel that approximately matches the PSF of the instrument before looking for five or more contiguous pixels above 1.5 times the sky rms. Although our criteria were experimentally determined to minimize the return of false detections by the software, inevitably manual identification of spurious detections (e.g., saturated stars, diffraction spikes, uncleaned cosmic rays near field boundaries, H ii regions inside nearby galaxies, etc.) is always required. We divided the 22 clusters into several groups and a few authors were assigned to each group to carry out this manual procedure. A possible concern is whether or not different groups might have non-negligible biases in this false object identification, which may propagate into cluster mass determination. To examine the possibility of bias in the manual procedure, several clusters were randomly selected by one author and compared to the two false object catalogs created by this author and others. We find that about ∼60% of false objects are shared by the two catalogs, and the disagreement is mostly attributed to very faint objects, whose fluxes are just above the sky background. Because these extremely faint objects were discarded anyway by our magnitude and shape measurement error criteria, the resulting difference in mass estimation and two-dimensional mass reconstruction is far below statistical errors.

We define source galaxies as objects whose i₇₇₅–z₈₅₀ colors are bluer than the color of the red sequence in each cluster while their total z₈₅₀ magnitudes (MAG_AUTO) are fainter than z₈₅₀ = 24. We also require the source galaxies to have ellipticity measurement error less than 0.25. The total exposure times for the 22 clusters vary substantially and thus, so does the number density of source galaxies as shown in Table 1.

Our shape measurement method is detailed in Jee et al. (2009). Briefly, we fit a PSF-convolved elliptical Gaussian function to source galaxies to determine their semi-major and -minor axes. Convolution with an elliptical Gaussian is required to measure galaxies' ellipticity before photons reach the instrument. Because an elliptical Gaussian profile is not an unbiased representation of true galaxies, the measurement inevitably introduces a bias. On average, the bias is toward underestimation of ellipticity because a Gaussian profile is steeper than that of real objects; more extended objects are subject to larger underestimates of ellipticity. In future analyses, it may be worth experimenting with more generalized functional forms and examining the effect based on simulated images. Nevertheless, we stress that the amount of bias induced by the current shape measurement is small and to first order was corrected here using the simulation results of Jee et al. (2007a), which shows a ∼6% bias for γ ≳ 0.5.

We model the PSF using our library constructed from stellar field observations (Jee et al. 2007a). In order to find the matching PSF template, we use ∼10 high S/N stars present in each cluster field. This matching is done for each visit, and the final PSF model is obtained after the model for each visit is shifted, rotated, and stacked. The mean residual ellipticity is less than ∼0.01, and this level of accuracy is sufficient for the current cluster lensing analysis. The detailed results for each cluster can be found in Appendix A.

Together with the PSF, CTI may be a potential source of systematics in weak lensing. For bright stars, the trailing from deferred charges is visually apparent. We quantify the effect using cosmic rays and warm pixels, and find that the effect on the ellipticity of the galaxies that we use to extract the signal is in fact much less than the statistical noise set by the finite number of galaxies. Thus, we confirm our previous argument (Jee et al. 2009) that the ACS CTI effect is negligible for cluster lensing analysis although it is a critical source of systematics in cosmic shear studies. We present the details in Appendix B.

In the weak-lensing regime, the characteristic length of the lensing signal variation is larger than the object size. Thus, the resulting shape distortion can be linearized as follows:

$$\begin{equation} \textbf {A}(\textbf {x}) = \delta _{ij} - \frac{\partial ^2 \Psi (\textbf {x})}{\partial x_i \partial x_j} = \left(\begin{array}{@{}c c@{}}1 - \kappa - \gamma _1 & -\gamma _2 \\ \ms {-}\gamma _2 & 1- \kappa + \gamma _1 \end{array} \right), \end{equation} \tag{ 1 }$$

where $\textbf {A}(\textbf {x})$ is the transformation matrix $\textbf {x}^\prime = \textbf {A} \textbf {x}$, which relates a position $\textbf {x}$ in source plane to a position $\textbf {x}^\prime$ in image plane, and Ψ is the two-dimensional lensing potential. The convergence κ is the surface mass density in units of critical surface mass density

$$\begin{equation} \Sigma _c= \frac{c^2}{4\pi G} \frac{D(z_s)}{D(z_l) D(z_l,z_s)}, \end{equation} \tag{ 2 }$$

where D(z_s), D(z_l), and D(z_l, z_s) are the angular diameter distance from the observer to the source, from the observer to the lens, and from the lens to the source, respectively. The convergence κ and the shears γ₁₍₂₎ are related to the lensing potential Ψ via

$$\begin{eqnarray} &&\kappa =\frac{1}{2} (\psi _{11}+\psi _{22}), \quad \gamma _1=\frac{1}{2} (\psi _{11} -\psi _{22}),\nonumber\\ &&\qquad \mbox{and} \quad \gamma _2=\psi _{12}=\psi _{21}, \end{eqnarray} \tag{ 3 }$$

where the subscripts on ψ_i(j) denote partial differentiation with respect to x_i(j).

A useful quantity to estimate masses is the reduced tangential shear defined by

$$\begin{equation} g_T = \left< - g_1 \cos 2\phi - g_2 \sin 2\phi \right>, \end{equation} \tag{ 4 }$$

where ϕ is the position angle of the object with respect to the cluster center and g is the reduced shear γ/(1 − κ) (valid only in the weak-lensing regime). This reduced tangential shear profile informs us of how the mass density of the lens changes as a function of radius. Of course, without any lensing signal, the resulting shear profile should vanish with fluctuations consistent with shot noise.

Hoekstra (2001, 2003) demonstrates that a cosmic shear (lensing by large-scale structure, LSS) is an important limiting factor in the accuracy of cluster masses derived from the cluster tangential shears. Following the formalism of Hoekstra (2003), we estimate that the cosmic shear contribution on average is γ ∼ 0.01 for our sample. Although this is small compared to the shot noise set by the finite number of source galaxies, we include the effect in our final mass uncertainty.

A potential ambiguity is the choice of the cluster center when there is disagreement between centroids determined from the distribution of cluster galaxies, lensing mass peaks, and X-ray emission. If we adopt a location that maximizes the amplitude of the tangential shear or a peak in the convergence map as the cluster center, the derived mass will be always biased high. Fortunately, in previous high S/N weak-lensing studies, where the number density of background galaxies is high (>120 per square arcmin), the weak-lensing mass peaks are in good spatial agreement with cluster galaxies. However, in many cases the X-ray peaks are conspicuously offset, indicative of merging activity or the collisional properties of the intracluster medium (ICM). In this paper, we adopt the centroid defined by the cluster galaxies as the center for the construction of the tangential shear profile. For the clusters in our sample where the statistical significance of the lensing signal is high, these centroids also agree well with the mass peaks. In several cases, where the lensing signal is weak (the cluster is not massive and/or the image is not sufficiently deep), we occasionally observe rather large (>20'') offsets between mass and galaxies. The choice of centroid will lead to an underestimate of the cluster mass if the offsets reflect real features in the system. We do not use the inner most data points in our mass estimation because the tangential shears at very small radii are highly sensitive to the centroid choice. Furthermore, both observationally and theoretically, it is not clear how the cluster mass profile behaves near the center.

There are two popular ways of deriving the lensing mass from tangential shear: aperture-mass densitometry (Fahlman et al. 1994) and parameterized model fitting. Aperture-mass densitometry is useful when one attempts to estimate the total projected mass within some aperture radius without requiring an assumption on the behavior of the cluster mass profile. However, this approach is not practical for the current data set, which in most cases provides only ∼3' × 3' areas smaller than the virial radii of the clusters. Therefore, we use the second method to derive the cluster masses. This method determines the parameters of analytic halo models by comparing the observed tangential shears with the expected values. Of course, a scatter is introduced because the real cluster profile is different from the model. Nevertheless, many studies show that this parametric approach gives consistent results with the values obtained by the former (e.g., Jee et al. 2005a; Hoekstra 2007).

We employ two kinds of halo models: singular isothermal sphere (SIS) and Navarro–Frenk–White (NFW; Navarro et al. (1997) profiles. Although the SIS profile is considered inappropriate at large radii (many lines of evidence suggest the mass density drops faster than 1/r²), the resulting parameter can be conveniently translated to the Einstein radius θ_E or the velocity dispersion σ_v of the system. We use this predicted velocity dispersion to compare with the dynamical value. When we fit NFW profiles, it is assumed that the cluster virial mass M₂₀₀ (the total mass at the radius inside of which the mean density is 200 times the critical density of the universe at the cluster redshift) is tied to the concentration c via the following relation from Duffy et al. (2008):

$$\begin{equation} c=5.71 \left(\frac{M_{200}}{2 \times 10^{12}\, h^{-1}\, M_{\odot }} \right) ^{-0.084} (1+z)^{-0.47}. \end{equation} \tag{ 5 }$$

Therefore, effectively, our NFW model is described by a single parameter. The relation between projected mass density and observed shear is simple for SIS with an Einstein radius r_E: κ = 0.5r_E/r and g = κ/(1 − κ) (in the weak-lensing regime). For NFW, the relation is rather complicated, and we refer readers to Bartelmann (1996). We present tangential shears and fitting results for individual clusters in Figure 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In principle, the two-dimensional projected mass distribution κ can be obtained by convolving the shear γ as follows (Kaiser & Squires 1993):

$$\begin{equation} \kappa (\textbf {x}) = \frac{1}{\pi } \int D^*(\textbf {x}-\textbf {x}^\prime) \gamma (\textbf {x}^\prime) d^2 \textbf {x}, \end{equation} \tag{ 6 }$$

where $D(\textbf {x}) = - 1/ (x_1 - i x_2)^2$ is the convolution kernel. However, the direct application of Equation (6) gives rises to some practical problems. First, the ellipticity of individual galaxies should be smoothed, and we do not know the optimal smoothing scale beforehand. Second, the smoothed galaxy ellipticity gives only an estimate for the reduced shear g = γ/(1 − κ). Third, the relation g = γ/(1 − κ) is only valid in the |γ/(1 − κ)| < 1 regime. Fourth, it is not obvious how to weight shape measurement noise properly via Equation (6).

Certainly, a more robust two-dimensional mass reconstruction algorithm is required to account for the aforementioned subtleties and minimize numerous artifacts. Many suggestions are present in the literature (Bridle et al. 1998; Seitz et al. 1998), and the consensus is that a reliable mass reconstruction should be achieved through iterations. In addition, it is desirable to use the ellipticity of individual galaxies rather than the smoothed values so as to reveal small-scale, but significant, features. In the present paper, we reconstruct two-dimensional mass maps using the entropy-regularized, maximum likelihood code of Jee et al. (2007a). This method of mass reconstruction accounts for the aforementioned subtleties and is effective in revealing low-contrast structures. We present our mass reconstructions in Figures 2–23. Results on individual clusters are discussed in Section 3.3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Because Σ_c scales as ∝D(z_s)/[D(z_l)D(z_l, z_s)] (Equation (2)), the mass density unit rises sharply as the redshift of the lens approaches that of the source. Thus, in the weak-lensing analysis of high-redshift clusters, small errors in our estimate of the source redshift distribution result in large errors in the translation of the signal into the physical unit. For example, a 10% systematic uncertainty in the effective source redshift leads to 1%–2% error in mass for a z = 0.2 object (e.g., A1689). The same amount of uncertainty would give ∼30% error for the mass of a z = 1.4 lens (e.g., XMM2235–2557).¹⁷

We estimate the redshift distribution of the source population using the publicly available Ultra Deep Field (UDF) photometric redshift catalog (Coe et al. 2006). The ultra deep images in six filters from F435W to F160W provide unparalleled high-quality photometric redshift information well beyond the limiting magnitude of the cluster images that we analyze here. The catalog has been extensively used in our previous studies (e.g., Jee et al. 2009), and thus we only briefly summarize the procedure.

First, we bin our source galaxies in 0.5 mag intervals. Next, we randomly draw the galaxies that match the selection criteria of each bin from the UDF photo-z catalog. Finally, we combine these simulated catalogs for all bins to create the redshift catalog for the entire source population. Because the UDF is much deeper than the cluster images and also because there is sample variance, it is important to take into account the difference in number density in this step. Consequently, we are utilizing the magnitude–color–redshift relation measured in the UDF data, rather than simplistically imposing the UDF redshift distribution on our images, which would cause much greater systematics.

As noted in Jee et al. (2009), the total error in z_eff due to the sample variance, the resampling error, and the difference among the photo-z estimation codes is δz_eff ≃ 0.06. This causes an ∼11% (∼3%) uncertainty in mass for a z ∼ 1.4 (z ∼ 0.9) cluster. We include these errors in our final uncertainty. The assumption that there exists a single source plane while in reality each source galaxy is at a different redshift also creates a bias in the interpretation of the lensing signal. We correct this bias to first order using the equation derived by Seitz & Schneider (1997).

We here comment on specific features in the measurement of each cluster. All of the key numbers are summarized in Table 2.

3.3.1. XMMXCS J2215−1738 (z = 1.46)

3.3.2. XMMU J2205−0159 (z = 1.12)

3.3.3. XMMU J1229+0151 (z = 0.98)

3.3.4. WARPS J1415+3612 (z = 1.03)

3.3.5. ISCS J1432+3332 (z = 1.11)

3.3.6. ISCS J1429+3437 (z = 1.26)

3.3.7. ISCS J1434+3427 (z = 1.24)

3.3.8. ISCS J1432+3436 (z = 1.35)

3.3.9. ISCS J1434+3519 (z = 1.37)

3.3.10. ISCS J1438+3414 (z = 1.41)

3.3.11. RCS 0220−0333 (z = 1.03)

3.3.12. RCS 0221−0321 (z = 1.02)

3.3.13. RCS 0337−2844 (z = 1.10)

3.3.14. RCS 0439−2904 (z = 0.95)

3.3.15. RCS 2156−0448 (z = 1.07)

3.3.16. RCS 1511+0903 (z = 0.97)

3.3.17. RCS 2345−3632 (z = 1.04)

3.3.18. RCS 2319+0038 (z = 0.91)

3.3.19. XLSS 0223−0436 (z = 1.22)

3.3.20. RDCS J0849+4452 (z = 1.26)

3.3.21. RDCS J0910+5422 (z = 1.11)

3.3.22. RDCS J1252−2927 (z = 1.24)

Table 2 lists the dynamical velocity dispersions and the predicted values from the lensing analysis for the combined sample. There are a total of 23 clusters whose velocity dispersions are available either in the literature or through our collaborations (Meyers et al. 2011; D. G. Gilbank et al. 2011, in preparation).

Figure 24 shows the comparison of the lensing prediction from the SIS fit with the dynamical measurement for these 23 clusters. The filled and open circles represent the clusters analyzed in the present paper and in our previous studies, respectively. The dotted line is a fit to the data, whose slope α = 1.12 ± 0.31 and intercept b = 28 ± 260 km s⁻¹ are consistent with the line of equality (solid); the shadow represents the 1σ range of the slope. Individual data points have rather large scatters around these lines mostly due to their large statistical (e.g., small number of known spectroscopic members and source galaxies) and systematic (e.g., cluster mass profile being different from SIS) uncertainties.

Zoom In Zoom Out Reset image size

Scaling relations between X-ray properties and lensing masses not only provide invaluable insight into the physical mechanism of cluster formation, but also help us to calibrate mass estimates based on X-ray data. For a virialized cluster, we expect the mass to scale with the gas temperature via the power-law relation M∝T^3/2 (Kaiser 1986). In this paper, we use the following specific form to fit the data:

$$\begin{equation} E(z)M_{\Delta } = M_5 \left(\frac{T}{5 \, \mbox{keV}} \right)^{\alpha }, \end{equation} \tag{ 7 }$$

where E(z) is the redshift-dependent Hubble parameter

$$\begin{equation} E(z)=\frac{H(z)}{H_0} = \sqrt{ \Omega _M (1+z)^3 + \Omega _{\Lambda }}. \end{equation} \tag{ 8 }$$

In Equation (7), the mass M_Δ is usually defined as the total mass within the radius, at which the mean density becomes Δ times the critical density. We quote the results for M₂₅₀₀ for an easy comparison with previous studies. In our combined sample, 14 clusters have published X-ray temperatures (Table 2) spanning the 1.7 keV < T < 10.4 keV range with reasonable constraints on their uncertainties. For these clusters, we present the mass–temperature relation in Figure 25. The thick red line shows the best-fit power-law relation with α = 1.54 ± 0.23 and M₅ = (9.13 ± 0.85) × 10¹³ h⁻¹ M_☉. The cluster RCS 0439−2904 (marked with red circle) is a significant outlier from this relation and the fit shown here is obtained without including this cluster. The exclusion of this cluster in our estimation of the M–T relation is justified because there is a strong indication that RCS 0439−2904 might be a line-of-sight superposition of multiple clusters (Section 3.3). However, we stress that adding RCS 0439−2904 does not change our results (reducing the slope to 1.48 ± 0.23).

Zoom In Zoom Out Reset image size

Allen et al. (2001) studied six relatively relaxed clusters at z ≲ 0.46 observed with Chandra and obtained the mass–temperature relation α = 1.51 ± 0.27 and M₅ = (1.32 ± 0.13) × 10¹⁴ h⁻¹ M_☉. Arnaud et al. (2005) used six nearby (z ≲ 0.15) T > 3.5 keV clusters observed with XMM-Newton and measured the slope α = 1.51 ± 0.11 and the normalization M₅ = (1.25 ± 0.04) × 10¹⁴ h⁻¹ M_☉, which are consistent with the results of Allen et al. (2001). A larger sample from Chandra has been studied by Vikhlinin et al. (2006), who analyzed 13 low-redshift clusters in the temperature range 0.7–9 keV. They also agreed that the observed scaling relation is in good accordance with the theoretical prediction, albeit with slightly higher slope 1.64 ± 0.06 than previous X-ray results. Their normalization M₅ = (1.25 ± 0.05) × 10¹⁴h⁻¹ M_☉ is in good agreement with previous X-ray results.

Because X-ray studies derive cluster masses from temperature with hydrostatic equilibrium assumption, these results showing that the correlation is tight and the power-law slope of the M–T_X relation is close to the theoretical prediction 3/2 may not be a direct proof for the cluster self-similarity. It is critical to investigate if the relation still holds when the masses are given by an independent estimator.

The weak-lensing mass versus X-ray temperature relation is investigated by Hoekstra (2007), who derived α = 1.34^+0.30_−0.28 and M₅ = (1.4 ± 0.2) × 10¹⁴ h⁻¹ M_☉ from the 17 clusters at 0.17 ≲ z ≲ 0.54, which are in good agreements with the values from X-ray studies. Mahdavi et al. (2008) updated the redshift distribution used by Hoekstra (2007) and found that the new n(z) from the photometric redshift catalog of Canada–France–Hawaii Telescope Legacy Survey derived by Ilbert et al. (2006) gives on average ∼10% smaller masses. We show this revised result in Figure 25. We further note that Hoekstra (2007) used the ASCA temperatures that do not correct for cool cores. The absence of this cool–core correction can bias the temperatures low.

We verified that the above M–T_X relation is not sensitive to a mass–concentration relation. When we recompute our M₂₅₀₀ using the mass–concentration relation of Bullock et al. (2001), the power-law slope α virtually remains the same and the normalization (M₅) decreases ∼3% only. The X-ray and Hoekstra (2007) results were obtained without any assumption on mass–concentration relation. Nevertheless, we tested if the result of Hoekstra (2007) changes when the masses are derived from NFW profile fitting with the Duffy et al. (2008) mass–concentration relation as is done in this paper. Again, we found that neither the normalization nor the slope of the M–T_X relation of Hoekstra (2007) changes because of this mass determination method.

It is easy to see in Figure 25 that there is a 20%–30% discrepancy in normalization between our result and previous results. The discrepancy may imply that the normalization decreases with redshift or our weak-lensing masses are biased low. One cause for possible weak-lensing bias is the redshift estimation bias, similar to the case in Hoekstra (2007). We estimate that the sample variance of the UDF may be responsible for the 5%–12% shift in mass depending on the cluster redshift. However, although we cannot exclude this possibility, this would aggravate the already problematic existence of the most massive clusters in the high-z universe (see Section 5 for details). Therefore, it seems more plausible that there is some evolution in the normalization of the mass–temperature relation to the highest redshifts.

Finally, we compare the masses derived from X-rays and lensing analyses. The cluster masses at r₂₀₀ obtained from the X-ray results are listed in Table 2. We use an isothermal β assumption to derive these values in order to take advantage of the published parameters in the literature. For the clusters that do not have published measurements of their X-ray surface brightness profile (because the S/N of the existing data do not constrain the shape of the X-ray profile well) we assume β = 0.7. Figure 26 shows the result for 14 clusters in our combined sample. The agreement between the lensing and X-ray estimates is good except for RCS0439−2904, which, as mentioned already, might be a line-of-sight superposition of multiple components (Cain et al. 2008). It is important to remember that an SIS profile is assumed for the derivation of M₂₀₀ from X-ray data. Therefore, it is possible that this excellent agreement might be a coincidence to some extent and a systematic difference might appear if masses are evaluated at different radii.

Zoom In Zoom Out Reset image size

The mass function is often given in the following form:

$$\begin{equation} \frac{dn}{d \ln M} = \frac{\rho _{m,0}}{M} \left| \frac{d \ln \sigma (M,z)}{d \ln M} \right| f, \end{equation} \tag{ 9 }$$

where σ is the rms variation of the density field when smoothed on scale M. ρ_{m, 0} is the present matter density. The function f determines the shape of the mass function.

Jenkins et al. (2001) proposed a redshift-independent form f = f(σ), which makes the above mass function agree with their numerical simulation results with high precision. Tinker et al. (2008) improve the agreement by considering the redshift-dependence of f = f(σ, z), which is expressed as

$$\begin{equation} f(\sigma,z)=A \left[ \left(\frac{a}{b} \right) ^{-a} + 1 \right] e^{-c/\sigma ^2}, \end{equation} \tag{ 10 }$$

where A = 0.186(1 + z)^−0.14, a = 1.47(1 + z)^−0.06, b = 2.57(1 + z)^−0.011, and c = 1.19 when the mass is defined as the total mass at the radius, inside of which the mean density is 200 times the mean density of the universe at the cluster redshift. We refer to it as M_200 m when necessary to distinguish it from M₂₀₀, which is more commonly used to refer to the mass at the radius, where the mean density becomes 200 times the critical density. The deviation of the universality of the mass function by Jenkins et al. (2001) from the numerical simulation results is 20%–50% (Tinker et al. 2008); Bhattacharya et al. (2011) argue that the difference between different fitting formulae is probably larger than the scatter seen in the largest simulations.

The cluster survey area that one chooses to adopt to evaluate the cosmological significance for discovery of unusually massive clusters is somewhat controversial. For example, Mortonson et al. (2011) argue that ∼300 deg² should be used for the search area of XMM2235−2557 at z = 1.4 instead of the pilot survey area of ∼11 deg² (Jee et al. 2009). The argument for this ∼300 deg² area is based on the study of Hoyle et al. (2011), who assumed ∼168 deg² for XCS, ∼64 deg² for XMM-LSS, ∼11 deg² for XDCP, ∼17 deg² for WARPS, and ∼33 deg² for RDCS. However, the different XMM-Newton-based surveys overlap significantly. In addition, because the flux limit that enables a secure cluster detection at z > 1 is ∼10⁻¹⁴ erg  cm⁻² s⁻¹, the effective search area must be reduced significantly. For example, the XDCS survey has continued since the end of the 11 deg² pilot survey and now reaches ∼50 deg² (R. Fassbender 2011, in preparation) for a flux limit of ∼10⁻¹⁴ erg  cm⁻² s⁻¹. The RDCS survey covers a geometrical area of 50 deg² only at relatively high fluxes. However, the area suitable for z > 1 cluster detection is ∼5 deg², corresponding to fluxes below ∼3 × 10⁻¹⁴ erg  cm⁻¹ s⁻¹ (Rosati et al. 1998). Considering the argument listed above, the combined survey area for the X-ray-selected high-z clusters should be ∼100 deg².

In the full sky survey, the number of clusters with mass and redshift greater than M_min and z_min, respectively, is given by

$$\begin{equation} N(M,z) = \int _{z_{{\rm min}}}^{z_{{\rm max}}} \frac{dV(z)}{dz} dz \int _{M_{{\rm min}}}^{M_{{\rm max}}} \frac{dn}{{\it dM}} {\it dM}, \end{equation} \tag{ 11 }$$

where dV/dz is the volume element and dn/dM is the mass function. For massive high-redshift clusters, the result is sensitive to the mass function dn/dM near z_min and M_min.

We compute the probability of finding each cluster by marginalizing over the measurement errors, assuming a log-normal distribution for P(M). This is different from the approach of Jee et al. (2009), where we simply used a threshold mass M_thr to evaluate the integral (Equation (11)). Nevertheless, we emphasize that the new method yields values very close to those obtained from the previous method when the mass uncertainty is relatively small (≲ 20%). For example, marginalizing over the measurement error gives an abundance of ∼2 × 10⁻³ for XMM J2235−2557 in the original 11 deg² XDCP survey, which is in good agreement with the previous value of ∼5 × 10⁻³. As we use a fixed value of z_max, the resulting abundance is slightly conservative because in principle the maximum redshift should decrease at the low end of P(M); remember that at the high end of P(M), the steep mass function cancels the effect of increased z_max. For the conversion from expected number to probability we assume Poissonian distribution. That is, if the expected number is N within a particular survey, the discovery probability is simply P = 1 − e^−N.

Mortonson et al. (2011) noted that an Eddington bias (Eddington 1913) is an important factor in the discussion of the cosmological significance using most massive clusters. Eddington bias is the combined effect of a steep mass function and measurement uncertainty. For example, consider the weak-lensing mass estimate of XMM2235 of M₂₀₀ = (7.3 ± 1.3) × 10¹⁴ M_☉. Because the slope of the mass function near the cluster mass is steep, there should be more ≲ 7.3 × 10¹⁴ M_☉ objects scattering up than ≳ 7.3 × 10¹⁴ M_☉ objects scattering down. As a matter of course, this does not mean that XMM2235's mass should be quoted with a lower mass (for individual clusters the scatter is symmetric). However, when we discuss the existence of the cluster in probabilistic terms, this additional chance of bias should be considered from a Bayesian point of view.

We adopt the convenient Mortonson et al. (2011) prescription, which is to replace the observed mass with the following reduced mass M' in the estimation of the abundance: M' = exp (1/2γσ_ln M)M, where γ is the local power-law slope (dn/dln M ∼ M^γ) and σ_ln M is the 1σ uncertainty of ln M (log-normal distribution is assumed for mass errors).

In Figure 27 and Table 3, we show the discovery probability of all clusters (i.e., similarly massive clusters) studied in this paper. We compute the probabilities after marginalizing over the Wilkinson Microwave Anisotropy Probe 7 year (WMAP7) Monte Carlo Markov Chain data.¹⁸ The probabilities estimated within the parent survey may lead to a somewhat aggressive interpretation for some clusters. For example, most of the RDCS clusters have significantly low probabilities mainly due to the small (∼5 deg²) survey area. However, most non-X-ray-selected clusters (ISCS and RCS) have discovery probability equal or close to unity even within their parent surveys.

Zoom In Zoom Out Reset image size

As mentioned in Section 5.2, we believe that the survey area of 100 deg² is still a conservative choice for most of the clusters shown here (the exceptions are MS1054−0321 and CL J0152−1357, which were discovered in the EMSS ∼735 deg² survey). With this survey area, there remain four clusters whose discoveries are still unusual given our current understanding of cosmological parameters (the names of these clusters are shown in Figure 27). The cluster with the lowest probability (∼1%) is the z ∼ 0.9 cluster CLJ1226+3332, whose weak-lensing mass and X-ray temperature are M_200 m ∼ 1.6 × 10¹⁵ M_☉ and T_X ∼ 10 keV, respectively. Because the cluster was discovered in a 72 deg² survey (WARPS), using 100 deg² does not increase the probability significantly. In addition, XMM-Newton archival search programs look for regions where no known X-ray luminous clusters exist, and thus the search volume in the XMM-Newton archival survey excludes the z ≲ 1 volume by observer's selection. This makes our probability estimation for CLJ1226+3332 even more conservative. The second lowest probability is found with XMM2235−2557 at z = 1.4. The probability of ∼3% is much higher (after the area and Eddington bias correction) than our previous value of ∼0.5% within the 11 deg² pilot XDCP survey, but still gives ∼2σ tension with the predicted number of clusters from ΛCDM. RDCS1252−2927 at z = 1.24 is ranked as third in rarity, giving P ∼ 6%, which is followed by XLSS0223−0436 at z = 1.22 (P ∼ 18%).

Combining discovery probabilities of all our clusters is somewhat tricky. Because of Poissonian fluctuation, the discovery probability is always lower than the estimated abundance (i.e., P = 1 − e^−N). For example, if the expected number of a given cluster is N = 1 within 100 deg² survey, the probability is ∼0.63. Obviously, simply multiplying probabilities of all clusters potentially leads to an unreasonably low final value if there are many slightly less-than-one-probability clusters. Nevertheless, if we include only rare clusters (P ≪ 1), multiplying individual probabilities provides a useful measure on how much the existence of these multiple clusters give rise to tension with the current cosmological parameters.

One set of cosmological parameters that increase the predicted abundance of one massive cluster can boost another similarly massive cluster's abundance. Therefore, multiplication of probabilities must be done separately for each chain, and the final probability should be obtained after properly taking into account each chain's weight. We find that the final combined probability obtained in this way is ∼0.03%.

Mortonson et al. (2011) presented fitting formulae, which provide an exclusion mass as a function of redshift for given sample and parameter variance confidence limits (CLs). Even a single cluster equal to or above the exclusion mass rules out both ΛCDM and quintessence.

In Figure 28, we compare our clusters to the exclusion curves of Mortonson et al. (2011). Although we believe that the proper choice for the sky area is 100 deg² or less, we also plot the result for the 300 deg² assumption, which was adopted in Mortonson et al. (2011). For simplicity, we only consider the case where the sample variance CL is equal to the parameter variance CL, which is referred to as joint CL in Mortonson et al. (2011).

Zoom In Zoom Out Reset image size

The central mass of the cluster CL J1226+3332 after the Eddington bias correction lies above the joint CL 95% exclusion curve when we adopt 100 deg² for the survey area. All the masses of the four clusters shown here are approximately on or above the joint CL 80% exclusion curves. When we consider 300 deg² for the search area instead, we note that the two clusters CLJ1226+3332 and XMM J2235−2557 are above the 80% exclusion curves. We clarify here that Mortonson et al. (2011) used the mass of XMM J2235−2557 reported by Stott et al. (2011), whose value is similar to our weak-lensing mass, but has significantly larger uncertainty (∼40%); Stott et al. (2011) derived the cluster mass from their X-ray scaling relation. This large uncertainty leads to ∼40% Eddington bias, and Mortonson et al. (2011) concluded that the resulting mass of the cluster provides only a negligible tension with ΛCDM when they assume 300 deg² for the sky area. On the other hand, our mass uncertainty is ∼18%, and the corresponding Eddington bias is ∼12%, which makes the corrected mass still close to the ∼85% exclusion curve for Mortonson et al. (2011)'s 300 deg² assumed area.

The above exclusion curve test provides weaker constraints when one is using more than a single cluster to examine their cosmological significance because the tool adopts the lowest redshift for the sample in the estimation of the cluster abundance (Mortonson et al. 2011). Therefore, we do not attempt to estimate the tension by combining the four most massive clusters.

Assuming the current well-accepted ΛCDM cosmology, we have demonstrated above that the discovery of the most massive clusters in our sample are rare events. Among the most popular interpretation would be that the primordial density fluctuation is non-Gaussian. For example, Jimenez & Verde (2009) claim that the local f_NL should be in the range of 150–260 to accommodate XMMU2235−2557 within the pilot 11 deg² survey. We suspect that their estimate of f_NL would decrease somewhat when both the Eddington bias and revised area of 50–100 deg² are taken into account. A tighter constraint is possible if one utilizes the weak-lensing masses of the four most massive clusters studied here. A similar study is performed by Hoyle et al. (2011) from the investigation of the 14 high-redshift cluster masses reported in the literature, giving a significantly large f_NL > 375 (after marginalization over WMAP5 parameters).

One crucial test that however should accompany these non-Gaussianity studies is the investigation of the mass function at the high mass end. The existing limited-volume simulations do not constrain the number density of extremely massive clusters accurately, and the commonly used fitting functions are simply extrapolated results in this regime. Another potentially important contribution from numerical simulation is the predicted mass function specific for a given survey, which takes into account the various aspects of selection limits and projection effects.

The projection effect is always a concern in the cosmological interpretation of extremely massive clusters. Certainly, a superposition of two moderately massive clusters or a long filament viewed along the line of sight can be identified with an extreme object. The former case is easily detected by a redshift histogram if the field is spectroscopically surveyed. Initially, RDCS1252−2927 at z = 1.24 was such a candidate because the redshift catalog showed three possible foreground groups at z = 0.47, 0.68, and 0.74 (Lombardi et al. 2005). However, the spatial distributions of these groups are not compact, and they appear to be just loose galaxy groups commonly found in any random field. The latter case (a filament viewed along the line of sight) is very difficult to prove (or disprove) directly. Nevertheless, there are a few circumstantial lines of evidence, which help us to infer if a system possesses extreme triaxiality. Comparing X-ray temperature with weak-lensing mass is useful because a severe triaxiality will give rise to much higher weak-lensing mass than what the temperature predicts. In addition, a line-of-sight elongation inflates velocity dispersion, and the result will be significantly larger than what the lensing predicts. XLSS0223−0436 is potentially such an object because our lensing mass of M₂₀₀ = 7.4^+2.5_−1.8 × 10¹⁴ M_☉ is much larger than what the X-ray temperature of ∼3.8 keV (Bremer et al. 2006) suggests. However, the X-ray photon statistics are poor, and Bremer et al. (2006) could not determine an upper limit of the X-ray temperature. Deeper Chandra observations would be necessary to robustly examine this large (a factor of 2–3) discrepancy. We note that the existing data for CLJ1226+3332 and XMM J2235−2557, the two clusters with the largest tension with ΛCDM in our sample, do not suggest the possibility that weak-lensing masses are significantly boosted by any projection effect.

We characterize the CTI with the ellipticity bias of the SSFs. Then, a nontrivial question is how to quantitatively relate this ellipticity bias to galaxy ellipticity change. Obviously, the ellipticity bias in SSF is a significant exaggeration of what occurs to galaxy shapes; a small CTI trail induces much larger ellipticity to a δ-function-like feature than to, for example, a 05 galaxy. The exact translation requires our knowledge of the shape of the CTI-induced trail and the surface brightness profile of galaxies. Consequently, the recovery of the galaxy shape prior to the CTI effect involves the same delicacies in PSF-effect correction, which is a more familiar problem.

Indeed, the CTI trailing can be treated as convolution with a one-dimensional kernel, which further smears the post-seeing objects, but only along the readout direction,

$$\begin{equation} o(x,y)=i(x,y)\otimes p(x,y) \otimes c(x,y) = i(x,y) \otimes r(x,y), \end{equation} \tag{ B1 }$$

where i(x, y) and o(x, y) are the intrinsic and observed images, respectively. p(x, y) and c(x, y) are the convolution kernel representing the PSF and the CTI trailing, respectively. Because the associative law holds for convolution, the smearing of the object can be viewed as a single convolution by the kernel r(x, y), which includes both the PSF and CTI trailing. One caveat here is that the CTI kernel c(x, y) in practice depends on flux, and, because a galaxy image consists of multiple pixels of varying flux, it is essential to approximate the effect with a single kernel, representing the collective effect of CTI on each galaxy.

One can determine the shape of the kernel c(x, y) by stacking many CTI trails. We find that on average at short distances, CTI trails are well approximated by an exponential tail ∝e^−y/τ, where the time constant τ (expressed in units of pixel) determines how quickly the trail drops. Upon close examination, however, we observe that the measured CTI trails differ in that (1) at large distances (> a few pixels) the measured CTI-trail slope is shallower and (2) occasionally a couple of small bumps appear, which suggests that different species of charge traps might be present. Nevertheless, we chose the simple functional form of a single exponential because object shapes are most sensitive to the information in the central few pixels. Based on the laboratory experiments of Dawson et al. (2008), Rhodes et al. (2010) also concluded that object shapes are most sensitive to the first few pixels and a single exponential term is adequate to describe the trails.

Now the remaining question is how to determine the relation between the ellipticity bias δe₊ in the SSFs and the CTI time constant τ. For this investigation, we rely on image simulations, where we create many artificial SSFs features and trail them with the ∼e^−y/τ kernel. Because we select only the r_h ∼ 1 SSFs to perform the study above, no exhaustive efforts are required to match the artificial SSFs to the real ones. Figure 32 shows the τ versus SSF ellipticity bias calibration that we derived from this image simulation. The thick solid line represents the relation when the size of the simulated SSFs matches the observed ones. The other lines illustrate how much the slope of the CTI trail depends on object sizes. The flowchart in Figure 33 summarizes the procedure to correct the PSF and CTI effects.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although the CTI correction described here is complicated and the result of time-consuming effort, we emphasize that the effect is on average small (also perhaps negligible in many cases) for cluster lensing analysis. Figure 34 shows that very few galaxies are subject to large elongation (≳ 0.05). A majority of galaxies, which fall to the CTI-mitigation regime are elongated by δe₊ ≪ 0.01.

Zoom In Zoom Out Reset image size