The Atacama Cosmology Telescope: Weighing Distant Clusters with the Most Ancient Light

Most of the mass in the universe is thought to consist of dark matter that does not interact with electromagnetic radiation other than through the gravitational force. Distortions in background light sources due to the gravitational influence of massive clusters of galaxies can be used to produce maps of the total matter distribution. This technique has served not only as evidence for the existence of dark matter (Trimble 1987; Massey et al. 2010), but also as a method for inferring the total mass of galaxy clusters themselves (e.g., Hoekstra et al. 2013). Such mass measurements are critical for the program of using the abundance of clusters to infer cosmological parameters like the dark energy equation of state or the mass scale of neutrinos (e.g., Allen et al. 2011; Madhavacheril et al. 2017). One key aspect of this program is our ability to constrain the abundance of clusters to high redshifts, directly probing the growth of structure through cosmic time (Voit 2005).

Sensitive large-area optical/near-infrared surveys are now allowing measurements of the mean mass of clusters up to z ∼ 1 (Murata et al. 2019; Chiu et al. 2020) through the lensing effects induced on background galaxies. However, as the distance of the clusters increases, the number of background galaxies that are useful for weak lensing measurements decreases rapidly. Measurements of this galaxy weak lensing effect at large distances are therefore only possible at present through deep targeted observations with the Hubble Space Telescope (Jee et al. 2011; Schrabback et al. 2018), with which observations currently exist for only some dozens of rich, massive clusters at redshifts of z > 0.8. A valuable complementary probe is emerging as cosmic microwave background (CMB) measurements are becoming sensitive enough to allow measurements of weak lensing by galaxy clusters in maps of the temperature and polarization of the CMB (e.g., Baxter et al. 2015; Madhavacheril et al. 2015; Planck Collaboration & Ade 2016; Geach & Peacock 2017; Raghunathan et al. 2019; Zubeldia & Challinor 2019). The high source redshift of the CMB allows weak lensing measurements to higher redshifts than galaxy lensing. Consequently, measurements of CMB lensing by clusters are anticipated to provide more stringent constraints on the masses of high-redshift clusters than enabled by future optical surveys (e.g., Madhavacheril et al. 2017).

We provide a mass estimate using gravitational lensing of the CMB for a blindly selected sample of galaxy clusters whose average mass has not been previously determined using galaxy lensing, which is additionally the highest redshift, $\langle z\rangle =1.08$, where a detection of gravitational lensing by galaxy clusters has been reported for a blindly selected sample to date. As opposed to targeted measurements of the most massive clusters, our work allows for inference of the average mass of a representative cluster sample. We make the code used in this analysis available at https://github.com/ACTCollaboration/madcows_lensing.

We use a combination of CMB data from the ground-based Atacama Cosmology Telescope (ACT) and the Planck satellite at the location of galaxy clusters selected from the Massive and Distant Clusters of WISE Survey (MaDCoWS, Gonzalez et al. 2019). The MaDCoWS clusters were identified as galaxy overdensities in near-infrared imaging (at 3.4 and 4.6 μm) from the WISE all-sky survey (Wright et al. 2010), and a large number of them were followed up with the Spitzer Space Telescope (at 3.6 and 4.5 μm). At declinations >−30°, the addition of optical data (grizy bands) from the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS, Chambers et al. 2016) allowed reliable photometric redshift estimation and subsequent cluster richness measurements. As a first attempt at a mass proxy, Gonzalez et al. (2019) define the richness, $\lambda$, to be the overdensity of red-sequence galaxies (identified using both PanSTARRS and Spitzer data) brighter than 15 μJy in the Spitzer 4.5 μm band and within 1 Mpc of the brightest cluster galaxy. Photometric redshifts were estimated from the Spitzer 3.6 and 4.5 μm bands, aided with the PanSTARRS i-band to remove low-redshift galaxies, with an estimated scatter ${\sigma }_{z}/(1+z)=0.04$ and no significant bias (but with an outlier fraction potentially of order 5%).⁴⁰ We discuss the impact of photometric redshift uncertainties below. We use this subset of the MaDCoWS "WISE-PanSTARRS" sample with available photometric redshifts, imposing an additional $\lambda \gt 20$ cut to reduce contamination by false detections in the MaDCoWS catalog, for a total of 677 clusters after masking point sources and other artifacts in the ACT maps (which is also limited to declinations  ≤20°). We show the richness and redshift distributions of these clusters in Figure 1. The mean redshift of the sample is $\langle z\rangle =1.08$ with a 16th and 84th percentile range of 0.93–1.26.

Zoom In Zoom Out Reset image size

To reconstruct the lensing signal, we use coadded maps of ACT and Planck CMB temperature data prepared separately at 98 and 150 GHz and described in Naess et al. (2020). The coadded maps include night-time data collected during the years 2008–2018 using the MBAC (Swetz et al. 2011), ACTPol (Thornton et al. 2016), and AdvACT (Henderson et al. 2016) receivers. The Planck maps used in the coadded maps are the PR2 (2015) CMB temperature maps at 100 and 143 GHz. We also use the Planck 2018 SMICA tSZ-deprojected maps (Planck Collaboration et al. 2020) as an additional input to the lensing reconstruction (see Appendix A). This is done in order to remove the bias from the thermal Sunyaev–Zeldovich (tSZ) effect due to inverse-Compton scattering of CMB photons off ionized electrons in hot gas in massive clusters (Madhavacheril & Hill 2018). Deprojection of the tSZ is possible because the frequency dependence of the spectral distortion due to the tSZ effect is well understood.

We reconstruct the CMB lensing convergence κ in regions centered on the locations of MaDCoWS clusters. The lensing convergence is related to the line-of-sight integral ϕ of the underlying gravitational potential sourced by a cluster, and to the lensing deflection angle ${\boldsymbol{\alpha }}$ through ${{\rm{\nabla }}}^{2}\phi ={\boldsymbol{\nabla }}\cdot {\boldsymbol{\alpha }}=-2\kappa$. The convergence is related to the surface mass density Σ through κ = Σ/Σ_cr, where Σ_cr is a characteristic mass density for the formation of multiple images that depends on distances to the source and lens (see Appendix B). Since the convergence map is directly proportional to the surface mass density, it can be thought of as a mass map. Gravitational lensing of the CMB by clusters leads to a remapping of the temperature anisotropies $T({\boldsymbol{x}})={T}_{\mathrm{unlensed}}({\boldsymbol{x}}+{\boldsymbol{\alpha }})$. In the 2D Fourier space of the temperature anisotropy image, this remapping corresponds to coupling between previously independent Fourier modes (say with wavenumbers ${\boldsymbol{\ell }}$ and ${\boldsymbol{\ell }}^{\prime}$) that is proportional to the lensing convergence: $\langle T({\boldsymbol{\ell }})T({\boldsymbol{\ell }}^{\prime} )\rangle \propto \kappa ({\boldsymbol{\ell }}+{\boldsymbol{\ell }}^{\prime} )$. This allows us to reconstruct the underlying convergence mode by mode by using a "quadratic estimator," i.e., a weighted sum over products of pairs of image modes (Hu et al. 2007), which in practice can be written as the divergence of the product of the large-scale CMB gradient and the small-scale CMB fluctuations (see Appendix A). The quadratic estimator reconstruction provides an unbiased estimate of the Fourier modes of the cluster mass map within a range of scales set by the band limit of the CMB maps. Thus, the output reconstruction image $\hat{\kappa }({\boldsymbol{x}})$ is effectively filtered, and we forward model this filtering when fitting to theoretical expectations. The signal-to-noise ratio for any given cluster is well below unity, and so our final measurements are a weighted average of the mass maps and of the radial profile (azimuthal average) of the lensing convergence, where the weights are inversely proportional to the noise variance in the lensing convergence reconstruction. We estimate a weighted covariance matrix for the radial profile from the scatter among the profiles. With each cluster profile weighted by w_i, the weighted covariance of the weighted mean of the N clusters in our sample is

$$\begin{eqnarray}&&{\boldsymbol{C}}=\displaystyle \frac{1}{N}\displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}{w}_{i}{\left({{\boldsymbol{p}}}_{i}-{\boldsymbol{d}}\right)}^{T}\left({{\boldsymbol{p}}}_{i}-{\boldsymbol{d}}\right)}{{V}_{1}-({V}_{2}/{V}_{1})}.\end{eqnarray} \tag{ 1 }$$

where ${{\boldsymbol{p}}}_{i}$ are the individual radial profiles of each cluster mass map, ${\boldsymbol{d}}=\tfrac{1}{{V}_{1}}{\sum }_{i=1}^{N}{w}_{i}{{\boldsymbol{p}}}_{i}.$ is the weighted mean profile, ${V}_{1}\,={\sum }_{i=1}^{N}{w}_{i}$ and ${V}_{2}={\sum }_{i=1}^{N}{w}_{i}^{2}$.

Details of the reconstruction process and mitigation of astrophysical foregrounds are provided in Appendix A. We show the resulting stacked CMB lensing convergence reconstruction in Figure 2.

Zoom In Zoom Out Reset image size

After averaging lensing convergence maps across all the clusters in our sample, we detect an excess relative to null within 8' at 4.9σ confidence. We estimate the mean mass of the sample by fitting the binned radial profile of the average mass map assuming the Navarro–Frenk–White (NFW; Navarro et al. 1996) density profile model for the distribution of matter within galaxy clusters with a fixed relation between the profile shape (or concentration parameter) and the mass. We account for the distribution of clusters of different masses across redshift and their clustering (i.e., the two-halo term), and apply an approximate correction for selection effects in order to construct an ensemble average of the lensing signal. This halo model approach, detailed in Appendix B, allows us to estimate the mean mass of the underlying sample from which the MaDCoWS sample is drawn.

We show the model predictions for the filtered lensing signal in Figure 3, color-coded by the excess χ² with respect to the minimum χ² obtained from our fit, corresponding to ${\chi }_{\min }^{2}=1.4$. Since our model has only one free parameter, the 1σ, 2σ, and 3σ credible ranges are given by Δχ² = 1, 4, and 9, respectively. We calculate the probability-to-exceed ${\chi }_{\min }^{2}$ (PTE) by drawing one million random samples with mean zero and using the covariance of our measurements; we find a PTE of 0.93, showing that the model is an adequate description of the data. The best-fit mean mass is $\langle {M}_{500c}\rangle =(1.7\pm 0.4)\,\times {10}^{14}\,{M}_{\odot }$, with a preference for our best fit over zero mass at the 4.2σ level.

Zoom In Zoom Out Reset image size

The clusters in our sample have a mean redshift of $\langle z\rangle =1.08$, the highest-redshift, blindly selected sample for which a gravitational lensing measurement has been obtained to date. Our results highlight the potential of both CMB lensing in general and ACT in particular to constrain the scaling relation between mass and other observables for high-redshift galaxy clusters.

In Appendix B, we explore systematic uncertainties from cluster redshift errors, centering errors, and the concentration parameter of the NFW profile, which are all smaller than 30% of our reported uncertainty. Since our measurement probes the overall lensing amplitude, the reported mass is also robust to assumptions about both the selection function and the scaling relation between mass and richness, which conversely we are not able to constrain.

Because we stack the lensing maps around all MaDCoWS clusters within the ACT footprint with a richness larger than 20, our measurement is representative of the full MaDCoWS sample above this cut. Targeted observations of the tSZ effect (which scales steeply with cluster mass as $\sim {M}^{\tfrac{5}{3}}$) of MaDCoWS subsamples by Gonzalez et al. (2019), Di Mascolo et al. (2020), and Dicker et al. (2020) have naturally picked preferentially the most luminous objects, and thus the mean richness of our sample, $\langle \lambda \rangle =32$, is lower than that probed in those works. The mean mass we estimate is consistent with the scaling relations reported in those works, which however have large uncertainties and are susceptible to additional systematic effects associated with an uncertain selection function.

While looking for galaxy overdensities is an efficient means of finding galaxy clusters, characterizing the selection effects that these techniques are subject to is a challenging endeavor. No estimate of the selection function is available for MaDCoWS, which precludes the use of this sample for cosmological parameter inference. Future CMB lensing measurements of well-defined cluster samples at high redshift will yield tight constraints on the parameters that govern the growth of structure.

We thank Peter Melchior for helpful discussions. Software used for this analysis includes healpy (Zonca et al. 2019), HEALPix (Górski et al. 2005), and Astropy (Astropy Collaboration et al. 2013, 2018). Research at the Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Industry Canada and by the Province of Ontario through the Ministry of Colleges and Universities. C.S. acknowledges support from the Agencia Nacional de Investigación y Desarrollo (ANID) through FONDECYT Iniciación grant No. 11191125. N.B. acknowledges support from NSF grant AST-1910021. E.C. acknowledges support from the STFC Ernest Rutherford Fellowship ST/M004856/2 and STFC Consolidated grant ST/S00033X/1, and from the Horizon 2020 ERC Starting Grant (grant agreement No. 849169). J.D. is supported through NSF grant AST-1814971. R.D. thanks CONICYT for grant BASAL CATA AFB-170002. D.H., A.M., and N.S. acknowledge support from NSF grant numbers AST-1513618 and AST-1907657. M.H.i. acknowledges support from the National Research Foundation. J.P.H. acknowledges funding for SZ cluster studies from NSF grant No. AST-1615657. K.M. acknowledges support from the National Research Foundation of South Africa. This work was supported by the U.S. National Science Foundation through awards AST-1440226, AST0965625, and AST-0408698 for the ACT project, as well as awards PHY-1214379 and PHY-0855887. Funding was also provided by Princeton University, the University of Pennsylvania, and a Canada Foundation for Innovation (CFI) award to UBC. ACT operates in the Parque Astronómico Atacama in northern Chile under the auspices of the Comisión Nacional de Investigación Científica y Tecnológica de Chile (CONICYT). Computations were performed on the GPC and Niagara supercomputers at the SciNet HPC Consortium. SciNet is funded by the CFI under the auspices of Compute Canada, the Government of Ontario, the Ontario Research Fund—Research Excellence, and the University of Toronto. The development of multichroic detectors and lenses was supported by NASA grants NNX13AE56G and NNX14AB58G. Colleagues at AstroNorte and RadioSky provide logistical support and keep operations in Chile running smoothly. We also thank the Mishrahi Fund and the Wilkinson Fund for their generous support of the project. This document was prepared by Atacama Cosmology Telescope using the resources of the Fermi National Accelerator Laboratory (Fermilab), a U.S. Department of Energy, Office of Science, HEP User Facility. Fermilab is managed by Fermi Research Alliance, LLC (FRA), acting under Contract No. DE-AC02-07CH11359.

The quadratic estimator for CMB lensing reconstruction can be recast as the divergence of the product of the small-scale CMB and the large-scale CMB gradient:

$$\begin{eqnarray}&&\hat{\kappa }({\boldsymbol{\theta }})=-{{ \mathcal F }}^{-1}\left\{{A}^{{TT}}({\boldsymbol{L}}){ \mathcal F }\left\{\mathrm{Re}\left[{\rm{\nabla }}\cdot \left[{\boldsymbol{\nabla }}{T}_{g}({\boldsymbol{\theta }}){T}_{h}({\boldsymbol{\theta }})\right]\right]\right\}\right\},\end{eqnarray} \tag{ A1 }$$

where ${ \mathcal F }$ and ${{ \mathcal F }}^{-1}$ denote 2D Fourier and inverse-Fourier transforms, respectively. Our analysis applies this estimator locally to cutouts of CMB data centered on each cluster location that are 128' wide with pixels of width 0.5 arcminutes. The large width relative to the typical arcminute size of clusters allows us to estimate the large-scale gradient.

Following Hu et al. (2007), we use a low-pass (top-hat) filtered temperature anisotropy gradient ${\boldsymbol{\nabla }}{T}_{g}$ with a maximum CMB multipole of ${{\ell }}_{G}=2000$ to mitigate bias in massive clusters as well as to reduce contamination from foregrounds. Furthermore, following Madhavacheril & Hill (2018), for the gradient map ${\boldsymbol{\nabla }}{T}_{g}$, we use a CMB map from which tSZ has been explicitly deprojected using multifrequency information (the Planck PR3 SMICA tSZ-deprojected), so as to null a large tSZ-induced bias. For the high-resolution map T_h, we use an internal linear combination (ILC) of the postage stamp cutouts from the 98 and 150 GHz 2018 coadds of Planck and ACT from Naess et al. (2020), where the high-resolution ACT data dominates the information content. While the tSZ cleaning for ${\boldsymbol{\nabla }}{T}_{g}$ nulls the tSZ bias, a large amount of variance can be induced by the presence of large tSZ decrements in T_h. To reduce this, following Patil et al. (2020), prior to the ILC we subtract best estimates of the tSZ decrements from each of the ∼4000 SZ clusters detected in the coadd (Hilton et al. 2020) since some of these clusters either appear in our sample or may appear in the postage stamps that we perform our reconstruction on. This has the effect of reducing the uncertainty in the first three bins of our measurement by ≈10%. The weights for the ILC are designed to minimize the power spectrum of the combination and are determined after fitting the total 1D power spectrum of each single-frequency stamp to a simple two parameter model that captures both an atmospheric component as well as white noise. We impose a maximum multipole cut of ℓ = 6000 on the high-resolution map. Both maps are inverse-variance filtered as in Hu et al. (2007). For both the gradient and high-resolution map, we do not use scales below ℓ = 200 due to the size of our cutouts. The reconstruction is normalized with an analytic expression ${A}^{{TT}}({\boldsymbol{L}})$ (Hu et al. 2007) that depends on the filters we apply. The final reconstruction is filtered to ensure it only contains modes 200 < L < 5000, and this filter is subsequently propagated in our theoretical model when fitting the profile.

Prior to reconstruction, the gradient and high-resolution maps are multiplied by a tapering cosine window (of approximate width 18') to enforce periodicity required by Fourier transforms. The presence of such a window (as well as other sources of anisotropy such as inhomogeneity of instrument noise) induces a spurious lensing signal referred to as the mean-field. We estimate the mean-field by repeating the above stacking procedure on a large number of random locations (∼300 times the number of clusters in our sample) and subtract this from our main cluster stack. The mean-field profile is roughly constant as a function of distance from the center of the stack and is 20% of the measurement in the first bin, 50% in the second bin, and larger than the measurement in subsequent bins.

We perform a number of systematics tests to ensure the robustness of our detection. This includes a curl null test, where the divergence in Equation (A1) is replaced by the curl (and normalized appropriately), shown in Figure 3. With this test, we obtain a measurement that is consistent with null (PTE of 0.05, corresponding to consistency with null at the $1.9\sigma$ level). As shown in Madhavacheril & Hill (2018), one needs correlated contaminants in the gradient and high-resolution maps of the quadratic estimator in order to be biased by these contaminants. Therefore, to test the possibility that dust emission might be biasing our measurement, we stack the Planck PR3 SMICA tSZ-deprojected map (used for the gradient) at the location of all 1676 MaDCoWS clusters with photometric redshifts (not just the 677 that fall in the ACT footprint used in this analysis). We detect no residual in this stack and thus conclude that dust contamination in the lensing reconstruction itself is unlikely. The presence of contaminants in the high-resolution map will contribute noise, which is captured in our empirically determined covariance matrix.

In order to infer the mean mass of the sample, we model the CMB lensing signal by assuming that all mass is contained within spherical halos, and that these halos cluster together (i.e., using the halo model; Peacock & Smith 2000; Seljak 2000). We adopt the best-fit flat ΛCDM cosmology from Planck Collaboration et al. (2016), with present-day matter density parameter ${{\rm{\Omega }}}_{{\rm{m}}}=0.307$ and a Hubble constant ${H}_{0}=67.7\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$, also adopted by Gonzalez et al. (2019).

Within the halo model formalism, the lensing signal can be decomposed as

$$\begin{eqnarray}&&\kappa (\theta )={\kappa }_{1{\rm{h}}}(\theta )+{\kappa }_{2{\rm{h}}}(\theta ),\end{eqnarray} \tag{ B1 }$$

where the subscripts 1h and 2h represent the one-halo (due to each halo) and two-halo (due to halo clustering) contributions. We model each component as follows.

For the one-halo term, we use an NFW density profile (Navarro et al. 1996). We use a mass definition, M_500c, corresponding to the mass within an overdensity of 500 times the critical density, ${\rho }_{{\rm{c}}}(z)=3{H}^{2}(z)/8\pi G$. We model the concentration, ${c}_{500c}\equiv {r}_{500c}/{r}_{{\rm{s}}}$ (where r_s is the scale radius and r_500c is the radius containing M_500c), using the relation between concentration and mass from Diemer & Joyce (2019).⁴¹ For a lens of mass m and redshift z, the convergence is related to the mass surface density, Σ(m, z), through

$$\begin{eqnarray}&&\kappa (m,z)=\displaystyle \frac{{\rm{\Sigma }}(m,z)}{{{\rm{\Sigma }}}_{\mathrm{cr}}(z)}\equiv {\left(\displaystyle \frac{{c}^{2}}{4\pi G}\displaystyle \frac{{D}_{{\rm{s}}}}{{D}_{{\rm{l}}}(z){D}_{\mathrm{ls}}(z)}\right)}^{-1}{\rm{\Sigma }}(m,z),\end{eqnarray} \tag{ B2 }$$

where D_s, D_l, and D_ls are the angular diameter distances to the last scattering surface, to the cluster, and between the cluster and the last scattering surface, respectively. The expression for the NFW surface density profile is given in Wright & Brainerd (2000).

The two-halo term arises due to the large-scale galaxy-matter power spectrum,

$$\begin{eqnarray}&&{P}_{\mathrm{gm}}^{2{\rm{h}}}(k| m,z)=b(m,z){P}_{{\rm{m}}}(k| z),\end{eqnarray} \tag{ B3 }$$

where ${P}_{{\rm{m}}}(k| z)$ is the linear matter power spectrum and b(m, z) is the halo bias as calculated by Tinker et al. (2010). Note that all calculations of the two-halo term are performed in comoving coordinates. The power spectrum is the Fourier transform of the correlation function,

$$\begin{eqnarray}&&{\xi }_{\mathrm{gm}}(r| m,z)=\displaystyle \frac{1}{2{\pi }^{2}}{\int }_{0}^{\infty }\displaystyle \frac{\sin {kr}}{{kr}}{k}^{2}\,{dk}\,{P}_{\mathrm{gm}}^{2{\rm{h}}}\end{eqnarray} \tag{ B4 }$$

from which the surface density can be calculated as

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{2{\rm{h}}}(R| m,z)=2{\bar{\rho }}_{{\rm{m}}}R{\int }_{0}^{1}\displaystyle \frac{{dx}}{{x}^{2}\sqrt{1-{x}^{2}}}\,{\xi }_{\mathrm{gm}}(R/x| m,z),\end{eqnarray} \tag{ B5 }$$

where ${\bar{\rho }}_{{\rm{m}}}$ is the (time-independent) comoving mean matter density. The two-halo term in Equation (B1) is then

$$\begin{eqnarray}&&{\kappa }_{2{\rm{h}}}(\theta | m,z)={\left(1+z\right)}^{2}\displaystyle \frac{{{\rm{\Sigma }}}_{2{\rm{h}}}(\theta | m,z)}{{{\rm{\Sigma }}}_{{\rm{c}}}},\end{eqnarray} \tag{ B6 }$$

where the additional ${(1+z)}^{2}$ compared to Equation (B2) is due to our use of comoving coordinates in the two-halo term calculation (see, e.g., Dvornik et al. 2018). We then filter both the 1h and 2h components with the Fourier-space filter described in Appendix A to produce a prediction for $\hat{\kappa }(\theta | m,z)$.

In reality, our measurement is the mean over clusters covering a range in masses and redshifts. The halo mass function provides a prediction for the distribution of clusters in mass and redshift, which can be linked to the observable used to construct the sample (namely, richness, $\lambda$) through a mass-observable relation. However, like for any other cluster sample, due to the selection algorithm (in this case, red-sequence overdensities), MaDCoWS is a biased subsample of the underlying cluster population. We therefore follow previous galaxy–galaxy- and cluster-lensing studies (e.g., van den Bosch et al. 2013; Viola et al. 2015; Dvornik et al. 2018; Miyatake et al. 2019) and model the measured lensing signal as

$$\begin{eqnarray}&&\begin{array}{l}\hat{\kappa }(\theta )=\displaystyle \frac{1}{\bar{N}}\displaystyle \int dz\,\displaystyle \frac{d{V}_{{\rm{e}}{\rm{f}}{\rm{f}}}}{dz}\displaystyle \int dm\,\displaystyle \frac{d{n}_{{\rm{h}}}}{dm}\\ \,\times \,\displaystyle \int d\lambda { \mathcal P }({\rm{l}}{\rm{n}}\lambda | \,{\rm{l}}{\rm{n}}m,z){ \mathcal S }(\lambda ,z)\hat{\kappa }(\theta | \,m,z).\end{array}\end{eqnarray} \tag{ B7 }$$

Here $\frac{d{V}_{{\rm{eff}}}}{{dz}}=A\frac{d{N}_{{\rm{cl}}}}{{dz}}$ is the effective volume probed per unit redshift, which accounts for the unknown redshift component in the selection function (and A is a constant that absorbs the unit conversion, but cancels between the numerator and denominator); $\frac{d{n}_{{\rm{h}}}}{{dm}}$ is the halo mass function from Tinker et al. (2010);⁴² ${ \mathcal P }(\mathrm{ln}\lambda | \mathrm{ln}m,z)$ is the conditional probability of richness given cluster mass (defined below), and $\bar{N}$ is the expected number of clusters, which differs from the halo mass function due to the mass-observable relation (see below) and selection function, ${ \mathcal S }(\lambda ,z)$, which is given by

$$\begin{eqnarray}&&\bar{N}=\int {dz}\,\displaystyle \frac{d{V}_{{\rm{eff}}}}{{dz}}\int {dm}\,\displaystyle \frac{d{n}_{{\rm{h}}}}{{dm}}\times \int d\lambda { \mathcal P }({\rm{ln}}\lambda | {\rm{ln}}m,z){ \mathcal S }(\lambda ,z).\end{eqnarray} \tag{ B8 }$$

Consequently, in this framework mean values refer to the mean of the underlying cluster population and are calculated in analogy to Equation (B7) as

$$\begin{eqnarray}&&\begin{array}{l} \langle X \rangle =\displaystyle \frac{1}{\bar{N}}\int dz\,\displaystyle \frac{d{V}_{{\rm{e}}{\rm{f}}{\rm{f}}}}{dz}\int dm\,\displaystyle \frac{d{n}_{{\rm{h}}}}{dm}\\ \,\times \,\int d\lambda { \mathcal P }({\rm{l}}{\rm{n}}\lambda | {\rm{l}}{\rm{n}}m,z){ \mathcal S }(\lambda ,z)\,X.\end{array}\end{eqnarray} \tag{ B9 }$$

For simplicity, we assume that the MaDCoWS sample includes all clusters with richness larger than 20 (and we have discarded all clusters with lower richness), such that

$$\begin{eqnarray}&&{ \mathcal S }(\lambda ,z)={ \mathcal S }(\lambda )={\rm{\Theta }}(\lambda -20),\end{eqnarray} \tag{ B10 }$$

where Θ is the Heaviside step function. In addition, we limit the redshift integral to the range $z\in [0.7,1.8]$. (Note that by using the observed redshift distribution, $\tfrac{{{dN}}_{\mathrm{cl}}}{{dz}}$, instead of the comoving volume element dV/dz in Equation (B7), we are effectively introducing an additional redshift component to the selection function; this choice has no impact on our results.)

Additionally, we assume a linear relation between cluster richness, $\lambda$, and cluster mass, ${M}_{500c}\equiv m$:

$$\begin{eqnarray}&&\mathrm{ln}\lambda =\beta +\mathrm{ln}(m/{10}^{14}\,{M}_{\odot })\equiv \beta +\mathrm{ln}\mu \end{eqnarray} \tag{ B11 }$$

and assume a log-normal conditional probability of $\lambda$ given cluster mass,

$$\begin{eqnarray}&&{ \mathcal P }(\mathrm{ln}\lambda | \mathrm{ln}m,z)=\displaystyle \frac{1}{\sqrt{2\pi }\sigma }\exp \left[-\displaystyle \frac{{\left(\mathrm{ln}\lambda -\mathrm{ln}\mu -\beta \right)}^{2}}{2{\sigma }^{2}}\right].\end{eqnarray} \tag{ B12 }$$

We assume a scatter on the mass-observable relation $\sigma \,\equiv {\sigma }_{\mathrm{ln}\lambda | \mathrm{ln}m}=0.5$, typical of the scatter of different richness definitions (Rozo et al. 2014; Murata et al. 2019). The normalization β is the only free parameter in this model. While its posterior value depends on our assumptions about the selection function and the scatter above, the mean mass of the sample (calculated following Equation (B9)) is robust to these changes as it directly quantifies the amplitude of the lensing signal.

We estimate the best-fit parameters by varying the normalization β and minimizing

$$\begin{eqnarray}&&{\chi }^{2}={\left({\boldsymbol{d}}-{\boldsymbol{m}}\right)}^{T}\cdot {{\boldsymbol{C}}}^{-1}\cdot \left({\boldsymbol{d}}-{\boldsymbol{m}}\right),\end{eqnarray} \tag{ B13 }$$

where ${\boldsymbol{d}}$ and m are the vectors of measured radial profile bins and the model expectations for them, respectively, and C is the covariance matrix of the measurement (see Equation (1)). We use five bins in our fit that encompass a region within 8' of the center of the stamp.

In order to demonstrate the contribution from each term in the model described above, we first fit our measurements with a single NFW profile at z = 1.1, including only the one-halo contribution, with the mass as the sole free parameter. This results in a best-fit mass $\langle {M}_{500c}\rangle =\left(2.1\pm 0.5\right)\times {10}^{14}\,{M}_{\odot }$. We then add the two-halo term, assuming all clusters have the same mass and are located at the same mean redshift, and find $\langle {M}_{500c}\rangle =\left(1.7\pm 0.4\right)\times {10}^{14}\,{M}_{\odot }$. Finally, we implement the full model as described above. We find β = 2.4 ± 0.3, which translates to a population-weighted mean cluster mass $\langle {M}_{500c}\rangle =\left(1.7\pm 0.4\right)\times {10}^{14}\,{M}_{\odot }$, where uncertainties correspond to Δχ² = 1 (i.e., 68.3% credible range). This is the main result of this paper.

For reference, we calculate the mass M_200c within an overdensity of 200 times the critical density for each point in the $({M}_{500c},z)$ grid and perform the integral given by Equation (B9) over the appropriate halo mass function, and find $\langle {M}_{200c}\rangle =(2.5\pm 0.6)\times {10}^{14}\,{M}_{\odot }$.

As discussed above, the best-fit β depends strongly on our modeling assumptions and should not be over-interpreted, but the mean mass changes at most by 5% regardless of our assumptions about the selection function and scaling relation. Specifically, we have varied the scaling relation normalization, β, and intrinsic scatter, σ, by factors of two from the adopted values, as well as selection functions given both by step functions and error functions with midpoints in the range λ = [10, 30], and in all cases find mean masses within 5% of the reported value. The largest shift in mass is introduced by modifying the mass-concentration relation (for reference, the population-weighted mean concentration is ${c}_{500c}=2.5$): increasing (decreasing) the concentration by 20%—which brackets most mass-concentration relations in the literature—increases (decreases) the mean mass by 10%, well within the uncertainties. Redshift uncertainties also do not change our results significantly: applying a systematic shift Δz = ± 0.1 (slightly larger than the 1σ scatter of photometric redshifts compared to spectroscopic redshifts, Gonzalez et al. 2019) only changes the mean mass by less than 3%. Assuming that 5% of the clusters have photometric redshifts that are wrong by Δz = +0.2 (3σ outliers, all in the same direction) has a similar impact.

Our final assessment of modeling systematics regards miscentering. The MaDCoWS algorithm defines the cluster center as the location of the peak amplitude in a smoothed galaxy density map, which can potentially be significantly offset from the center of mass (George et al. 2012; Viola et al. 2015) and results in a reduced central lensing amplitude. However, the physical resolution of our measurements is roughly 1 Mpc, which greatly reduces the impact from miscentering, and we choose not to account for it in our analysis. To test the impact of miscentering, we implement the distribution found by Johnston et al. (2007) for optically selected clusters, exploring typical offsets between the assumed and true centers of up to 0.4 Mpc (roughly $2{r}_{{\rm{s}}}$); the change in mass is at most 5%.