# Masses of Sunyaev-Zel’dovich Galaxy Clusters Detected by The Atacama
Cosmology Telescope: Stacked Lensing Measurements with Subaru HSC Year 3 data

Masato Shirasaki  0000-0002-1706-5797<EMAIL_ADDRESS>National
Astronomical Observatory of Japan (NAOJ), National Institutes of Natural
Science, Mitaka, Tokyo 181-8588, Japan  The Institute of Statistical
Mathematics, Tachikawa, Tokyo 190-8562, Japan  Cristóbal Sifón
0000-0002-8149-1352 Instituto de Física, Pontificia Universidad Católica de
Valparaíso, Casilla 4059, Valparaíso, Chile  Hironao Miyatake
0000-0001-7964-9766 Kobayashi-Maskawa Institute for the Origin of Particles
and the Universe (KMI), Nagoya University, Nagoya, 464-8602, Japan  Institute
for Advanced Research, Nagoya University, Nagoya 464-8601, Japan  Kavli
Institute for the Physics and Mathematics of the Universe (WPI), The
University of Tokyo Institutes for Advanced Study (UTIAS), The University of
Tokyo, Chiba 277-8583, Japan  Erwin Lau  0000-0001-8914-8885 Center for
Astrophysics — Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA
Zhuowen Zhang  Department of Astronomy and Astrophysics, University of
Chicago, 5640 S. Ellis Ave, Chicago, IL, USA, 60637  Fermilab, Kirk & Pine
Road, Batavia, IL, USA, 60510  Neta Bahcall  Department of Astrophysical
Sciences, Peyton Hall, Princeton University, Princeton, NJ USA 08544  Nicholas
Battaglia  0000-0001-5846-0411 Department of Astronomy, Cornell University,
Ithaca, NY 14853, USA  Mark Devlin  0000-0002-3169-9761 Department of Physics
and Astronomy, University of Pennsylvania, 209 South 33rd Street,
Philadelphia, PA, USA  Fermilab, Kirk & Pine Road, Batavia, IL, USA, 60510  Jo
Dunkley  0000-0002-7450-2586 Joseph Henry Laboratories of Physics, Jadwin
Hall, Princeton University, Princeton, NJ, USA 08544  Department of
Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ USA
08544  Arya Farahi  0000-0003-0777-4618 Department of Statistics and Data
Science, The University of Texas at Austin, TX 78712, USA  Matt Hilton
0000-0002-8490-8117 Wits Centre for Astrophysics, School of Physics,
University of the Witwatersrand, Private Bag 3, 2050, Johannesburg, South
Africa  Astrophysics Research Centre, School of Mathematics, Statistics, and
Computer Science, University of KwaZulu-Natal, Westville Campus, Durban 4041,
South Africa  Yen-Ting Lin  Institute of Astronomy and Astrophysics, Academia
Sinica, Taipei 10617, Taiwan  Daisuke Nagai  0000-0002-6766-5942 Department of
Physics, Yale University, New Haven, CT 06520, USA  Suzanne T. Staggs
0000-0002-7020-7301 Joseph Henry Laboratories of Physics, Jadwin Hall,
Princeton University, Princeton, NJ, USA 08544  Tomomi Sunayama
0009-0004-6387-5784 Steward Observatory, University of Arizona, Tucson, AZ
85704, USA  Kobayashi-Maskawa Institute for the Origin of Particles and the
Universe (KMI), Nagoya University, Nagoya, 464-8602, Japan  David Spergel
Simons Foundation, 160 Fifth Avenue NY, NY 10010, USA  Edward J. Wollack
0000-0002-7567-4451 NASA/Goddard Space Flight Center, Greenbelt, MD 20771, USA

###### Abstract

We present a stacked lensing analysis of 96 galaxy clusters selected by the
thermal Sunyaev-Zel’dovich (SZ) effect in maps of the cosmic microwave
background (CMB). We select foreground galaxy clusters with a $5\sigma$-level
SZ threshold in CMB observations from the Atacama Cosmology Telescope, while
we define background source galaxies for the lensing analysis with secure
photometric redshift cuts in Year 3 data of the Subaru Hyper Suprime Cam
survey. We detect the stacked lensing signal in the range of
$0.1<R\,[h^{-1}\mathrm{Mpc}]<100$ in each of three cluster redshift bins,
$0.092<z\leq 0.445$, $0.445<z\leq 0.695$, and $0.695<z\leq 1.180$, with 32
galaxy clusters in each bin. The cumulative signal-to-noise ratios of the
lensing signal are $14.6$, $12.0$, and $6.6$, respectively. Using a halo-based
forward model, we then constrain statistical relationships between the mass
inferred from the SZ observation (i.e. SZ mass) and the total mass derived
from our stacked lensing measurements. At the average SZ mass in the cluster
sample ($2.1-2.4\times 10^{14}\,h^{-1}M_{\odot}$), our likelihood analysis
shows that the average total mass differs from the SZ counterpart by a factor
of $1.3\pm 0.2$, $1.6\pm 0.2$, and $1.6\pm 0.3$ (68%) in the aforementioned
redshift ranges, respectively. Our limits are consistent with previous lensing
measurements, and we find that the cluster modeling choices can introduce a
$1\sigma$-level difference in our parameter inferences.

##  I Introduction

Galaxy clusters are the largest gravitationally bound objects in the universe.
Analytic approaches (e.g. [1, 2]) and numerical simulations (e.g. [3, 4, 5, 6,
7]) predict that the number density of galaxy clusters as a function of their
masses is primarily determined by the variance in linear fluctuations in
cosmic mass density within a Lagrangian comoving scale of
$R=[3M/(4\pi\bar{\rho}_{\mathrm{m0}})]^{1/3}$, where $M$ is the cluster mass
and $\bar{\rho}_{\mathrm{m0}}$ is the average cosmic mass density in comoving
coordinates. Hence, number counts of galaxy clusters offer a unique means of
measuring the formation rate of cosmological structures at their most massive
end, $\sim 10^{14}\,h^{-1}M_{\odot}$, giving physical insight into the mass
variance at the scale of $R=8\,h^{-1}\mathrm{Mpc}$, the quantity referred to
as $\sigma_{8}$ in the literature. This enables one to constrain physics
behind fundamental components in our universe such as dark energy by counting
the number of galaxy clusters (see Ref. [8] for a review).

Cosmological inference with galaxy clusters requires careful assessment of
both selection effects and biases in mass estimates. Among several cluster
finders, observations of thermal Sunyaev Zel’dovich (SZ) effects on the cosmic
microwave background (CMB) [9] have received great attention because the
thermal SZ effect in a single cluster does not depend on the cluster redshift.
Recent observations of the CMB with multiple frequency channels allow
searching for galaxy clusters based on their thermal SZ effects alone. Several
large catalogs of SZ-selected clusters have been constructed [10, 11, 12],
some containing more than $1000$ clusters. In this work, we employ a cluster
catalog from the Atacama Cosmology Telescope (ACT) [12]. Future CMB
observations can further increase the number of the SZ-selected clusters,
enhancing the constraining power of galaxy clusters in fundamental physics
[13]. Several cosmological analyses with the SZ-selected clusters have been
performed [14, 15, 11, 16]; however, their cosmological constraints hinge on
systematic uncertainties in the estimate of cluster masses.

Masses of individual clusters are commonly estimated from X-ray observations
under the assumption of hydrostatic equilibrium (see Ref. [17] for a review).
For several reasons, the assumption behind the X-ray mass estimates can be
invalid in practice. Those include asphericity of the galaxy clusters,
clumpiness in cluster gas density, gas accelerations at cluster outskirts, and
non-thermal pressure components induced by turbulent gas motions (see Ref [18]
and references therein for more details). Recent progress in galaxy imaging
surveys has opened an additional means of estimating cluster masses –
measuring the gravitational lensing effect on background galaxies. The mass
distributed around individual clusters can act as a lens for background
galaxies, introducing a coherent pattern of distortion in the galaxy shapes.
Tangential distortions in the shapes of background galaxies with respect to
the foreground cluster position provide a direct measure of the cluster mass
including invisible dark matter [19].

In this paper, we measure a gravitational-lensing signature induced by the SZ-
selected clusters obtained from the Atacama Cosmology Telescope [12]. We
employ a stacked lensing technique, which is suitable for studying the average
mass density distribution around the clusters of interest [20]. The stacked
lensing measurement allows us to investigate mass density profiles around the
clusters because it combines the lensing effects of individual clusters into a
single profile and thus increases the detection significance by reducing shot
noise in the galaxy-shape measurements. Moreover, spherically-symmetric
modeling of the mass density profile provides a good approximation to the
stacked lensing measurements, because combining several lensing measurements
suppresses any preferences in the orientation of real galaxy clusters given
the statistical isotropy of the universe. For the stacked lensing
measurements, we use the galaxy shape catalog constructed from Year 3 data of
the Subaru Hyper Suprime Cam (HSC) Survey [21]. The HSC has the largest source
number density among its competitors in stage-III lensing surveys, allowing us
to study the stacked lensing measurements of the SZ-selected clusters over a
wide range of cluster redshifts.

We improve upon our previous analysis in Ref. [22] in several ways. The main
improvements on the previous analysis are summarized below.

  1. (i)

We increase the number of the SZ-selected clusters by a factor of 12. Our
cluster sample now consists of 96 objects with a redshift range $0.096<z\leq
1.180$.

  2. (ii)

We estimate statistical errors in our stacked lensing measurements with a set
of realistic mock catalogs of HSC source galaxies and SZ-selected clusters
from ACT observation. Our mock catalogs are constructed from full-sky lensing
simulations and their associated halo catalogs in Ref. [23], including all
sample-variance effects caused by non-linear structure formation.

  3. (iii)

We adopt halo-based forward modeling for our stacked lensing measurements. As
an input, the model requires an estimate of scatter111We mean the scatter is
between clusters in the sample. in the scaling relation between an SZ
observable and the cluster mass; this is the dominant uncertainty in our
modeling. We use the Illustris-TNG simulations [24] to set an informative
prior on the scatter caused by the variation of electron thermal pressure
inside clusters. Also, we include the information on cluster number counts
when inferring our model parameters with the stacked lensing measurements,
helping to break inherent parameter degeneracy in our model.

  4. (iv)

We incorporate possible baryonic effects on cluster mass density profiles with
an efficient emulator of dark-halo statistics developed in Ref. [25]. Hence,
our model includes fully non-linear gravitational effects in the mass
distribution around cluster-sized halos as well as back-reaction effects on
the cluster mass by the presence of stellar and gas components. We marginalize
the modeling uncertainty due to baryonic effects to place constraints on the
cluster masses.

The goal of this paper is to constrain statistical relationships between an SZ
observable and the total cluster mass with our stacked lensing measurements.
To be specific, we focus attention on the mass inferred from the SZ
observation, referred to as SZ mass, and on a possible difference between the
SZ and total cluster mass. Given the completeness of our SZ-selected clusters
as a function of SZ masses, we can infer the total cluster mass from our
stacked lensing measurements by accounting for selection effects in the ACT
observations. Because the estimate of SZ masses relies on a model of cluster
electron pressure calibrated with X-ray observations under the assumption of
hydrostatic equilibrium, our inference of the cluster mass can validate the
conventional hydrostatic-mass assumption in a direct manner.

Figure 1:  ACT SZ-selected clusters in the HSC S19A fields. The colored points
show 96 ACT clusters with the selection threshold of ${\tt
fixed\verb||_{|}SNR}>5$. The blue circles represent the clusters in the
redshift range $0.092<z\leq 0.445$, while orange triangles and the green
squares show clusters at $0.445<z\leq 0.695$ and $0.695<z\leq 1.180$,
respectively. The white regions in each panel highlight the full-depth, full-
color (FDFC) regions in the HSC S19A observations.

Recently, Ref. [26] presented a stacked lensing analysis of the ACT-selected
galaxy clusters using the KiDS-1000 weak lensing catalog. Their analysis
mainly focuses on the calibration of cluster masses at the average cluster
redshift of $z\sim 0.3-0.4$, whereas we aim at studying a possible redshift
evolution of stacked lensing effects around the clusters employing a secure
selection of background galaxies based on HSC photometric redshift
measurements. Note that the authors in Ref. [26] adopted a simple halo model
based on the universal dark-matter density profile [27] to infer the cluster
masses. In this paper, we apply a more flexible model, accounting for baryonic
effects and possible deviations from the universal density profile around
boundaries of galaxy clusters.

The rest of this paper is organized as follows. We describe our data in
Section. II. Next, we summarize details of our stacked lensing measurements in
Section. III. We explain our halo-based model of the stacked lensing
measurements and the setup of our parameter inference in Section IV. Section V
presents the key results, whereas we discuss the limitations of our analysis
in Section VI. Finally, concluding remarks are provided in Section VII. In the
following, $\ln$ represents the natural logarithm and $\log$ is the logarithm
with base 10. Throughout this paper, we adopt $\Lambda$CDM cosmological
parameters below; the average cosmic mass density $\Omega_{\mathrm{m}}=0.279$,
the cosmological constant $\Omega_{\Lambda}=1-\Omega_{\mathrm{m}}=0.721$, the
average baryon density $\Omega_{\mathrm{b}}=0.046$, the present-day Hubble
parameter $H_{0}=100h=70.0\,\mathrm{km}/\mathrm{s}/\mathrm{Mpc}$, the spectral
index of the power spectrum of primordial curvature perturbations
$n_{s}=0.97$, and the square root of the linear mass variance within
$8\,h^{-1}\mathrm{Mpc}$ being $\sigma_{8}=0.82$. Those parameters are
consistent with statistical analyses of the CMB anisotropies in Ref. [28]. Our
analysis requires mock catalogs of the HSC lensing data and the ACT cluster
samples to precisely estimate statistical errors in our lensing measurements.
The mock catalogs in this paper have been constructed from N-body simulations
with the above cosmological model, which is slightly different from the latest
version inferred by the Planck satellite [29].

##  II Data

###  II.1 ACT DR5 clusters

In this paper, we use the cluster sample constructed from ACT observations
from the 2008-2018 observing seasons [12]. The sample was extracted from 98
and 150 GHz observations of a 13,211 $\mathrm{deg}^{2}$ sky, and form part of
the ACT fifth data release (DR5; [30]). A multi-frequency matched filter (e.g.
[31, 32]) was applied to the DR5 maps for cluster searches. When searching for
cluster candidates, the authors in Ref. [12] constructed a signal template
from the Universal Pressure Profile (UPP; [33]) and its associated SZ signal-
to-mass scaling relation. We provide the explicit functional form of the UPP
model in Appendix A. Candidates were then confirmed as galaxy clusters by
cross-matching with external optical and infrared (IR) data sets, including
several spectroscopic surveys [34, 35, 36, 37], as well as optically-selected
cluster samples in Sloan Digital Sky Survey (SDSS), Dark Energy Survey (DES),
and HSC [38, 39, 40]. In addition, some cluster redshifts were estimated by
the zCluster algorithm [41]. The detection significance of each cluster
candidate is characterized by the signal-to-noise ratio at the reference
2.4-arcmin filter scale, referred to as ${\tt fixed\verb|_|SNR}$. In this
paper, we impose a threshold ${\tt fixed\verb|_|SNR}>5$ to have a secure
cluster sample, as explained in Ref. [12]. Further details of the cluster
search and optical/IR follow-up processes are summarized in Ref. [12]. In the
following, we refer to our cluster sample as ACT clusters.

###  II.2 HSC-ACT survey overlap

In the HSC observations acquired from March 2014 to April 2019 (referred to as
S19A), the full-depth full-color (FDFC) footprint has
$474.43\,\mathrm{deg}^{2}$ of overlap with the ACT DR5 cluster search area.
Fig. 1 represents the angular distribution of our cluster sample in five HSC
S19A regions of interest, named XMM, VVDS, GAMA09H, GAMA15H, and WIDE12H. We
then estimate the completeness of the ACT cluster search in terms of the mass
parameter in the UPP model $M_{\mathrm{UPP}}$ (denoted as the SZ mass
hereafter) using Monte Carlo simulations as described in Ref. [12]. It is
beneficial to note that $M_{\mathrm{UPP}}$ is a fiducial mass estimate and we
use it for ease of comparison with previous works. Fig. 2 shows the
completeness as a function of redshift and $M_{\mathrm{UPP}}$ across the HSC
S19A regions222The completeness in this paper is computed with the Monte Carlo
code as of August 2021. Since then, the code has been updated so that it
improves beam convolutions and in-painting objects, introducing small
differences in the completeness. We have confirmed that the latest
completeness changes our model predictions of lensing signals at a $<2\%$
level, but this does not affect our final results significantly.. In the
figure, the dotted line highlights the 90% mass completeness limit, which is
similar to the equivalent limit when averaging over the whole 13,211
$\mathrm{deg}^{2}$ ACT cluster search area. Note that we confirmed that every
cluster in our sample in the HSC S19A regions has complete redshift follow-up.
For our stacked lensing analysis, we divide our cluster sample into three
subsamples by the cluster redshifts so that the number of objects in each
subsample is equal (32). The gray histogram in Fig. 3 shows the redshift
distribution of ACT clusters in the HSC S19A regions with ${\tt
fixed\verb|_|SNR}>5$. We then construct three equal subsamples by imposing
limits $0.092<z\leq 0.445$, $0.445<z\leq 0.695$ and $0.695<z\leq 1.180$, as
highlighted in colored shaded regions in Fig. 3. The average redshifts of the
three samples are given by 0.324, 0.567, and 0.878, respectively. For
reference, the three samples have the average SZ masses in units of
$10^{14}\,h^{-1}M_{\odot}$ of 2.41, 2.39, and 2.09 from lower to higher
redshift bins.

Figure 2:  The completeness for ${\tt fixed\verb||_{|}SNR}>5$ as a function of
redshift, in terms of the mass parameter in the UPP model $M_{\mathrm{UPP}}$.
This is computed for sources in the HSC S19A footprint overlapping the ACT DR5
map. Note that the UPP model in Ref. [33] is adopted to estimate the
completeness. The dotted line shows the 90% completeness limit.

Figure 3:  The redshift distribution of ACT clusters in the HSC S19A fields.
The gray histogram shows the distribution for objects selected by ${\tt
fixed\verb||_{|}SNR}>5$. The colored shaded regions mark the boundaries in
cluster redshifts used to define the three subsamples for our stacked lensing
analysis. We set the boundary so that the number of each subsample is equal.

###  II.3 HSC Y3 source galaxies

HSC is a wide-field optical imager installed at the Subaru Telescope [42, 43].
Under the Subaru Strategic Survey Program (SSP; [44]), HSC conducted a deep
multi-band imaging survey with a target sky coverage of 1,400
$\mathrm{deg}^{2}$, starting in 2014. The 330-night observations in the HSC
SSP program have been completed. In this paper, we use the HSC third-year
galaxy shape catalog (HSC Y3; [21]) based on the HSC-SSP imaging data taken
from March 2014 through April 2019. The HSC Y3 data consist of six distinct
fields (HECTOMAP, GAMA09H, WIDE12H, GAMA15H, XMM, VVDS) and cover
$433.48\,\mathrm{deg}^{2}$ in total. We ignore the HECTOMAP field since it
does not overlap with the ACT DR5 map. The effective number density of source
galaxies marginally depends on the full width at half maximum of $i$-band
Point Spread Function (PSF) reaching $19.9\,\mathrm{arcmin}^{-2}$ on average.
More details of the Y3 shape catalog are found in Ref. [21], but we summarize
the basic information on the HSC Y3 source galaxies below.

The galaxy shapes were estimated on coadded $i$-band images by the re-
Gaussianization PSF correction method [45]. This moment-based method has been
applied to SDSS imaging data and systematic errors have been extensively
characterized [46, 47]. In the HSC Y3 shape catalog, the shapes (ellipticities
to be more precise) of individual galaxies are defined by

$\displaystyle(e_{1},e_{2})=\frac{1-(b/a)^{2}}{1+(b/a)^{2}}(\cos 2\phi,\sin
2\phi),$ (1)

where $\phi$ is position angle and $b/a$ is the minor-to-major axes ratio. Eq.
(1) provides an estimator of lensing distortion [48] and one needs to
introduce a multiplicative correction of the shape to estimate lensing shear
$\gamma$ as [49]

$\displaystyle\hat{\gamma}_{\alpha}=\frac{1}{2{\cal R}}\langle
e_{\alpha}\rangle,$ (2)

where $\alpha=1,2$ are the indices for the galaxy shapes and lensing shear,
$\langle\cdots\rangle$ indicates a weighted average and $\cal R$ is called the
shear responsivity. An inverse variance weight is commonly used for the
average in Eq. (2) and it is given by

$\displaystyle
w_{i}=\frac{1}{e^{2}_{\mathrm{RMS},i}+\sigma^{2}_{\mathrm{e},i}},$ (3)

where $e_{\mathrm{RMS},i}$ and $\sigma_{\mathrm{e},i}$ are the intrinsic shape
noise and the error in the shape measurement for the $i$-th galaxy,
respectively. The measurement error $\sigma_{\mathrm{e}}$ is statistically
estimated from the shape measurements performed on simulated images and the
intrinsic shapes are derived by subtracting the measurement error from the
ellipticity dispersion measured from the real data. The shear responsivity for
a given source population is then computed as

$\displaystyle{\cal
R}=1-\frac{\sum_{i}w_{i}e^{2}_{\mathrm{RMS},i}}{\sum_{i}w_{i}}.$ (4)

To eliminate possible biases in the shear estimation in the re-Gaussianization
method, the galaxy shapes were further calibrated against image simulations
generated with the open source software GalSim [50]. We refer the reader to
Ref. [51] for details of the production of image simulations. The calibration
assumes that Eq. (2) can provide a biased estimator of the true shear
expressed as

$\displaystyle\hat{\gamma}_{\alpha,\mathrm{obs}}=(1+m)\gamma_{\alpha,\mathrm{true}}+c_{\alpha},$
(5)

where $\gamma_{\alpha,\mathrm{true}}$ represents the true shear, $m$ and
$c_{\alpha}$ are correction factors, referred to as the multiplicative and
additive biases, respectively. Extensive use of HSC-like image simulations
allows us to calibrate $m$ and $c_{\alpha}$ as a function of galaxy
properties. For our stacked lensing analysis, we properly include the
correction for both multiplicative and additive biases to produce an unbiased
estimate of lensing shear as described in Section 3.1.2 and 3.1.3 in Ref.
[21].

Five-band photometry in the HSC SSP allows us to estimate redshifts of the HSC
Y3 source galaxies with photometric redshift (photo-z) measurements (see Ref.
[52] for the photo-z statistics in the first-year HSC data). Among several
photo-z methods, we use a neural network-based photo-z conditional density
estimation code ${\tt dnnz}$, as our fiducial choice. The ${\tt dnnz}$
architecture consists of multilayer perceptrons with 5 hidden layers. The
training uses ${\tt cmodel}$ fluxes, unblended convolved fluxes, PSF fluxes,
and galaxy shape information. The construction of the conditional density uses
100 nodes in the output layer, each representing a redshift histogram bin
spanning the range $z=0$ to $7$ (see Nishizawa et al. in prep. for details of
${\tt dnnz}$). To define background source galaxies behind a given ACT
cluster, we apply the P-cut method [53] by using the ${\tt dnnz}$ estimates of
the posterior redshift distribution of each HSC Y3 galaxy. More detail about
the P-cut method is provided in Section III.2.

##  III Stacked lensing measurements

###  III.1 Estimator

Gravitational lensing effects by foreground clusters introduce a coherent
distortion of observed shapes of background source galaxies. In the weak
lensing regime, lensing effects by the cluster’s gravitational potential leave
a tangential distortion in galaxy shapes after one takes the azimuthal average
of the source shapes around the cluster center. For a source galaxy at a
comoving transverse separation $R$ from the lens center, the azimuthal
averaged profile of tangential distortion $\gamma_{t}$ can be expressed as
(see Ref. [19] for a review)

$\displaystyle\gamma_{t}(R)=\frac{\bar{\Sigma}(<R)-\Sigma(R)}{\Sigma_{\mathrm{cr}}(z_{l},z_{s})}\equiv\frac{\Delta\Sigma(R)}{\Sigma_{\mathrm{cr}}(z_{l},z_{s})},$
(6)

where $\Sigma(R)$ is the surface mass density of the lens cluster,
$\bar{\Sigma}(<R)$ is the average projected mass density inside $R$, and
$\Sigma_{\mathrm{cr}}(z_{l},z_{s})$ represents the critical surface mass
density which depends on the redshifts of the source galaxy, $z_{s}$, and the
lens, $z_{l}$. For a single lens-source pair,
$\Sigma_{\mathrm{cr}}(z_{l},z_{s})$ is defined as

$\displaystyle\Sigma_{\mathrm{cr}}(z_{l},z_{s})=\frac{c^{2}}{4\pi
G}\frac{D_{A}(z_{s})}{(1+z_{l})^{2}D_{A}(z_{l})D_{A}(z_{l},z_{s})},$ (7)

where $D_{A}(z_{s})$, $D_{A}(z_{l})$, and $D_{A}(z_{l},z_{s})$ represent the
angular diameter distances of the source and lens from us, and between the
two, respectively. Note that the factor $(1+z_{l})^{-2}$ in Eq. (7) arises
from our use of comoving coordinates [54].

Eq. (6) motivates us to measure the excess surface mass density $\Delta\Sigma$
from real data on the background galaxies’ shapes. In this work, we aim at
measuring the average of $\Delta\Sigma$ over a set of foreground clusters by
stacking tangential components of the shapes of background galaxies, $e_{t}$.
In practice, the estimator of $\Delta\Sigma$ for our cluster sample is given
by (e.g. see Ref. [55] for a similar analysis)

$\displaystyle\widehat{\Delta\Sigma}(R_{i})$ $\displaystyle=$
$\displaystyle\frac{\sum_{\mathrm{ls}\in
R_{i}}w_{\mathrm{ls}}\langle\Sigma^{-1}_{\mathrm{cr}}\rangle_{\mathrm{ls}}\left[e_{t,\mathrm{s}}/\left(2{\cal
R}(R_{i})\right)-c_{t,\mathrm{s}}\right]}{\left[1+K_{\mathrm{sel}}(R_{i})\right]\left[1+K(R_{i})\right]\sum_{\mathrm{ls}\in
R_{i}}w_{\mathrm{ls}}}\Bigg{|}_{R_{i}}$ (8)
$\displaystyle\qquad-(\mathrm{signal\,around\,random\,points})|_{R_{i}},$

where the summation “ls” run over all lens-source pairs that lie in the i-th
radial bin $R_{i}\equiv D_{A}(z_{l})\Delta\theta_{\mathrm{ls}}$;
$D_{A}(z_{l})$ is the comoving angular diameter distance to the l-th ACT lens
cluster at the redshift $z_{l}$, and $\Delta\theta_{\mathrm{ls}}$ is the
angular separation between lens and source in each pair; $e_{t,\mathrm{s}}$ is
the tangential component of ellipticity of the s-th HSC source galaxy;
$c_{t,\mathrm{s}}$ is the additive shear calibration factor as given in the
HSC shape catalog; $\langle\Sigma^{-1}_{\mathrm{cr}}\rangle_{\mathrm{ls}}$ is
the average of the inverse critical surface mass density-weighted with the
photo-z posterior distribution function of each source galaxy, $P_{s}(z_{s})$,
behind the l-th lens cluster:

$\displaystyle\langle\Sigma^{-1}_{\mathrm{cr}}\rangle_{\mathrm{ls}}=\frac{\int_{0}^{z_{\mathrm{max}}}\mathrm{d}z_{s}\,P_{s}(z_{s})\,\Sigma^{-1}_{\mathrm{cr}}(z_{l},z_{s})}{\int_{0}^{z_{\mathrm{max}}}\mathrm{d}z_{s}\,P_{s}(z_{s})},$
(9)

where $z_{\mathrm{max}}=7$ as given by ${\tt dnnz}$. Note that we set
$\Sigma^{-1}_{\mathrm{cr}}=0$ in Eq. (9) when $z_{s}<z_{l}$. The factor
$w_{\mathrm{ls}}$ in Eq. (8) represents the weight, for which we employ an
inverse-variance weighting that is nearly optimal in the shape-noise-dominated
regime [56, 46, 47, 57]:

$\displaystyle
w_{\mathrm{ls}}=w_{\mathrm{l}}w_{\mathrm{s}}\langle\Sigma^{-1}_{\mathrm{cr}}\rangle_{\mathrm{ls}}^{2}$
(10)

where $w_{\mathrm{s}}$ is the weight in the HSC shape catalog as in Eq. (3)
and we adopt $w_{\mathrm{l}}=1$ for the weight of the ACT clusters.

In Eq. (8), we introduce the shear responsivity for the stacked lensing signal
at a given radial bin as

$\displaystyle{\cal R}(R_{i})=1-\frac{\sum_{\mathrm{ls}\in
R_{i}}w_{\mathrm{ls}}e^{2}_{\mathrm{RMS,s}}}{\sum_{\mathrm{ls}\in
R_{i}}w_{\mathrm{ls}}}.$ (11)

Also, we adopt the correction of the multiplicative bias in the shape
measurements as

$\displaystyle 1+K(R_{i})=1+\frac{\sum_{\mathrm{ls}\in
R_{i}}w_{\mathrm{ls}}\,m_{\mathrm{s}}}{\sum_{\mathrm{ls}\in
R_{i}}w_{\mathrm{ls}}},$ (12)

where $m_{\mathrm{s}}$ is the multiplicative bias for the s-th source galaxy.
The factor of $1+K_{\mathrm{sel}}(R_{i})$ accounts for the effect of selection
bias. We set $1+K_{\mathrm{sel}}(R_{i})$ by following the procedure as in
Section 3.1.4 in Ref. [21]. Note that Eq. (12) corrects a 1%-level bias in our
measurements.

Finally, the second term on the right-hand side of Eq. (8) denotes the signal
around a set of random positions, which is measured by replacing lens clusters
with random points. We need to subtract this random signal to correct for
observational systematic effects such as residual systematics in shape
measurements due to an imperfect correction of optical distortions across the
field of view. The number of random points is 100 times as large as that of
lens clusters, where the random catalogs are generated assuming the
hydrostatic equilibrium and analytic halo mass function in Ref. [3]. We should
point out here that this is done with the purpose of reproducing the redshift
distribution because it impacts source selection, but the mass distribution is
irrelevant. In the end, the redshift distribution of the random points is
slightly different from the counterpart of the actual ACT clusters. To correct
the difference in the redshift distribution in our lensing analysis, we set
the weight for the random points (denoted as $w_{\mathrm{r}}$) as

$\displaystyle w_{r}(z)=\frac{{\cal
H}_{\mathrm{cl}}(z;z_{\mathrm{b}},\Delta)}{{\cal
H}_{\mathrm{r}}(z;z_{\mathrm{b}},\Delta)},$ (13)

where ${\cal H}_{\mathrm{cl}}(z;z_{\mathrm{b}},\Delta)$ is the histogram for
the number of the clusters at their redshift range of
$z_{\mathrm{b}}-\Delta/2\leq z<z_{\mathrm{b}}+\Delta/2$, ${\cal
H}_{\mathrm{r}}(z;z_{\mathrm{b}},\Delta)$ is the histogram for the random
points, and each histogram is normalized to unity. To compute the histogram,
we employ 200 linear-spaced bins in the range of $0\leq z<2$. We found that
the signal around random points starts to deviate from zero at $R\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}20\,h^{-1}\mathrm{Mpc}$, which is
caused by a finite area coverage in our lensing analysis.

Figure 4:  Stacked lensing signal, $\Delta\Sigma$ (normalized by the signal
with $p_{\mathrm{cut}}=0.99$ and $\Delta z=0.3$), as a function of lens
redshift threshold $\Delta z$ in the P-cut selection (see Eq. [14] for the
definition). The three panels summarize the results at different cluster
redshift bins of $0.092<z\leq 0.445$, $0.445<z\leq 0.695$ and $0.695<z\leq
1.180$ from left to right. In each panel, different symbols represent the
measurements at different radial bins centered at
$R=0.573\,h^{-1}\mathrm{Mpc}$, $1.234\,h^{-1}\mathrm{Mpc}$ and
$2.658\,h^{-1}\mathrm{Mpc}$. This range of $R$ approximates the virial radii
of a typical cluster-sized halo with $M\sim 10^{14-15}\,h^{-1}M_{\odot}$. The
horizontal dashed line indicates the most stringent selection criteria
($p_{\mathrm{cut}}=0.99$ and $\Delta z=0.3$). The error bars in each panel are
set by the shot noise due to intrinsic scatters in galaxy shapes.

Throughout this paper, we employ a logarithmic binning with 18 radial bins
ranging from $R=0.1\,h^{-1}\mathrm{Mpc}$ to $R=100\,h^{-1}\mathrm{Mpc}$ for
our stacked lensing analysis. For our fiducial choice, we define the center of
each cluster as the position of the brightest cluster galaxy (BCG). In
Appendix B, we perform several systematic tests on our stacked lensing
analyses, including the divergence-free B-mode signal around the lens
clusters, the B-mode signal around random points, and the measurement of the
boost factor that arises from contamination of the source galaxy sample by
galaxies physically associated with lens clusters themselves, since the source
photo-z measurements may be imperfect. We did not find any evidence for such
residual systematic effects in our stacked lensing measurement.

Apart from the null tests, we also examine how sensitive our stacked lensing
analysis is to (i) the definition of cluster centers and (ii) the choice of
photo-z methods (see Appendix B). We found that the magnitude of stacked
lensing signal with the center inferred from the thermal SZ analysis (the SZ
center) tends to become smaller than that using the BCG positions at $R\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}0.3\,h^{-1}\mathrm{Mpc}$ when
stacking clusters in the lowest redshift bin (the difference is at a $\sim
1\sigma$ level). In contrast, the choice of the cluster centers does not
introduce a significant effect in our stacked lensing analysis for the other
two cluster redshift bins. To investigate the photo-z impact, we measure
$\Delta\Sigma$ by using photo-z information estimated with a different photo-z
code, ${\tt mizuki}$, which uses template fitting and physical priors about
galaxy properties [58]. Note that our default photo-z code ${\tt dnnz}$ takes
a data-driven approach using a neural network, while ${\tt mizuki}$ falls into
the category of intuitive photo-z methods with informative priors about galaxy
physics. In this sense, the comparison of $\Delta\Sigma$ with the ${\tt dnnz}$
photo-z and the ${\tt mizuki}$ counterpart enables us to bracket a typical
uncertainty due to the choice of photo-z methods. Imposing the background
selection as described in Section III.2, we did not find any significant
differences in our stacked lensing signals when using the photo-z from ${\tt
dnnz}$ or ${\tt mizuki}$.

###  III.2 Selecting background sources

If member galaxies in a cluster of interest are misidentified as background
galaxies, they will introduce a systematic dilution in the $\Delta\Sigma$
measurement for that cluster. To make a secure selection of background
sources, we adopt the P-cut method as introduced in Ref. [40]. The method
utilizes the photo-z posterior distributions of individual source galaxies and
selects the source galaxy when it satisfies the condition below:

$\displaystyle p_{\mathrm{cut}}<\int_{z_{l}+\Delta
z}^{7.0}\,\mathrm{d}z_{s}\,P_{s}(z_{s}),$ (14)

where $z_{l}$ is the redshift of the l-th cluster, $P_{s}(z_{s})$ is the
photo-z posterior distribution of the s-th source galaxy, $p_{\mathrm{cut}}$
and $\Delta z$ are parameters in the P-cut selection. Ref. [59] have
extensively examined the P-cut selection by varying the parameters for
optically-selected clusters and source galaxies in previous HSC data, finding
that $p_{\mathrm{max}}=0.98$ and $\Delta z=0.2$ are a reasonable choice to
mitigate the dilution effect.

We here investigate the choice of $p_{\mathrm{cut}}=0.98$ and $\Delta z=0.2$
with the HSC Y3 source galaxies and the ACT clusters. We perform the stacked
lensing measurements by using the ${\tt dnnz}$ photo-z information, while
varying the parameters $p_{\mathrm{cut}}$ and $\Delta z$ used in the
background selection. We compare our fiducial $\Delta\Sigma$ signals with
those obtained by assuming $p_{\mathrm{cut}}=0.99$ and $\Delta z=0.3$ to
examine possible systematic effects in the P-cut background selection. Fig. 4
shows $\Delta\Sigma$ in three different radial bins, varying the parameter
$\Delta z$ but fixing $p_{\mathrm{cut}}=0.98$. These are compared to
$\Delta\Sigma$ calculated with $p_{\mathrm{cut}}=0.99$ and $\Delta z=0.3$. As
seen in the figure, we find that the selection with $p_{\mathrm{cut}}=0.98$
and $\Delta z=0.2$ does not introduce any significant dilution effects
compared to the more stringent selection criteria ($p_{\mathrm{cut}}=0.99$ and
$\Delta z=0.3$) at the cluster redshifts of $0.092<z\leq 0.445$ and
$0.445<z\leq 0.695$. In the highest redshift bin of $0.695<z\leq 1.180$, we
see a trend showing possible dilution in our stacked lensing measurement when
adopting $p_{\mathrm{cut}}=0.98$ and $\Delta z=0.2$, while the choice of
$p_{\mathrm{cut}}=0.98$ and $\Delta z=0.3$ provides more converged results.
Given the results of this experiment, we decide to use $\Delta z=0.2$ for the
first and second cluster redshift bins, but $\Delta z=0.3$ is applied to the
$\Delta\Sigma$ measurement at the highest redshift bin.

Because the P-cut method is applied to every single lens-source pair, it is
not suitable to define a single population of background galaxies.
Nevertheless, it is worth computing a representative number density of
background source galaxies by using a single value of the lens redshift.
Following the definition of the effective source number density
$\bar{n}_{\mathrm{gal}}$ in Ref. [60], we find $\bar{n}_{\mathrm{gal}}=10.8$,
$6.48$ and $0.40\,\mathrm{arcmin}^{-2}$ when applying the P-cut selection for
the lens redshifts of $0.324$, $0.567$, and $0.878$, respectively. Note that
the small source number density at the highest cluster redshift bin prevents
us from measuring $\Delta\Sigma$ at $R\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}0.3\,h^{-1}\mathrm{Mpc}$ (there are
no pairs of random points and source galaxies inside this radius after we
impose the P-cut selection).

###  III.3 Statistical uncertainty

Figure 5:  The diagonal components of the covariance matrix for measurements
of the stacked lensing profiles for the ACT clusters. We used 1404 mock
catalogs of the ACT clusters and HSC-Y3 source galaxies to estimate the sample
variance (see text for details). Three panels summarise the results at
different cluster redshift bins ($0.092<z\leq 0.445$, $0.445<z\leq 0.695$, and
$0.695<z\leq 1.180$ from right to left). The colored symbols in each panel
show the full covariance for the lensing profile, while the dashed lines
represent the shape noise contribution alone for comparison.

Precise estimations of the statistical error in $\Delta\Sigma$ measurements
are essential for a robust characterization of ACT clusters. In this paper, we
follow the method given in Refs. [61, 62] to determine the covariance matrix
of statistical errors for our stacked lensing analysis. The covariance matrix
in $\Delta\Sigma$ measurements can be decomposed into two parts:

$\displaystyle\mathrm{Cov}_{ij}=\mathrm{Cov}^{\mathrm{int}}_{ij}+\mathrm{Cov}^{\mathrm{sample}}_{ij},$
(15)

where the indices of $i$ and $j$ denote the i- and j-th radius bins,
$\mathrm{Cov}^{\mathrm{int}}_{ij}$ represents the statistical uncertainty due
to intrinsic scatter in galaxy shapes (so-called shape noise), and
$\mathrm{Cov}^{\mathrm{sample}}_{ij}$ stands for the sample variance arising
from uncorrelated large-scale structures along a line of sight, inherent
variations in mass density distributions inside clusters, and modulation
effects in underlying density fields inside our survey window caused by
density fluctuations at scales longer than the survey window [63, 64].

To estimate the first term of the right-hand side in Eq. (15), we use a random
rotation of the HSC Y3 source galaxy shapes and real ACT clusters. It is worth
noting that we can remove the real lensing signal by rotating the ellipticity
orientation of each source galaxy. We perform the stacked lensing measurements
with randomly-rotated shapes of HSC Y3 galaxies and the ACT clusters in the
same manner as in the actual measurement, including the corrections of
multiplicative bias discussed in Section III.1. We repeat the process of
random rotation and stacked lensing analysis to produce 1000 realizations of
$\Delta\Sigma$, allowing us to estimate the shape noise contribution to the
covariance matrix.

To find the sample covariance of $\mathrm{Cov}^{\mathrm{sample}}_{ij}$ in Eq.
(15), we use a set of mock catalogs of the HSC Y3 source galaxy shapes and ACT
clusters. To construct the mock catalogs in a light-cone volume, we use 108
full-sky, light-cone simulations generated in Ref. [23]. As described in
Appendix C in detail, we populate HSC-like galaxies and ACT-like clusters in
the full-sky, light-cone simulation. Since the HSC-Y3 data still has a small
area coverage (about 433 square degrees), we identify 13 footprints of the
HSC-Y3 data in each of the all-sky, light-cone simulation realizations. We
thus generate $108\times 13=1404$ mock catalogs of the HSC data and the ACT
clusters in total.

We then perform measurements of $\Delta\Sigma$ from each realization using the
same analysis pipelines that are used in the actual measurements of the real
HSC and ACT data. Because the intrinsic scatter and measurement error in
galaxy shapes does not affect the sample covariance of
$\mathrm{Cov}^{\mathrm{sample}}_{ij}$, we use a modified version of
$\Delta\Sigma$ estimator instead of Eq. (8):

$\displaystyle\widehat{\Delta\Sigma}^{\prime}(R_{i})$ $\displaystyle=$
$\displaystyle\frac{\sum_{\mathrm{ls}\in
R_{i}}w_{\mathrm{ls}}\langle\Sigma^{-1}_{\mathrm{cr}}\rangle_{\mathrm{ls}}\gamma^{\mathrm{mock}}_{t,\mathrm{s}}}{\sum_{\mathrm{ls}\in
R_{i}}w_{\mathrm{ls}}}\Bigg{|}_{R_{i}}$ (16)
$\displaystyle\quad-(\mathrm{signal\,around\,random\,points})|_{R_{i}},$

where $\gamma^{\mathrm{mock}}_{t,\mathrm{s}}$ represents the tangential
component of lensing shear for the s-th HSC mock galaxies. Finally, we
estimate the sample covariance for the stacked lensing signal $\Delta\Sigma$
at three cluster redshift bins, from the scatters among the 1404 measurements
of $\widehat{\Delta\Sigma}^{\prime}$.

Fig. 5 shows the diagonal components of the covariance matrix. In each panel,
dashed lines highlight the shape-noise covariance decreasing with radius as
$1/R$. The colored symbols in each panel show the full variance including the
sample variance term, making the variance larger at $R\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}10\,h^{-1}\mathrm{Mpc}$ for the
first and second cluster redshift bins. Note that the shape-noise variance is
dominant in the statistical error of our $\Delta\Sigma$ measurement at the
highest cluster redshift bin because of a secure background selection of the
HSC source galaxies. We also find that the off-diagonal components in the
$\Delta\Sigma$ covariance become prominent at $R\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}3\,h^{-1}\mathrm{Mpc}$ for the
first and second cluster redshift bins, while they are only important at
$R\lower 2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}10\,h^{-1}\mathrm{Mpc}$
for the highest cluster redshift bin.

##  IV model and inference

In this section, we describe a theoretical model to predict our lensing
measurements. Subsection IV.1 provides key ingredients in our halo-based
model, including the definition of halo masses, the scaling relation between
the SZ and halo masses, the halo-model expression of our observables, and the
choice of model parameters. Compared with our previous work [22], we adopt an
informative prior of the scatter in the scaling relation extracted from a
hydrodynamical simulation suite of Illustris TNG [24]. In Subsection IV.2, we
explain how to incoorporate possible baryonic effects with our model by using
a semi-analytic method [65, 66, 67]. Subsection IV.3 gives some physical
intuition about the model-parameter dependence of stacked lensing profiles
$\Delta\Sigma$ based on our model. We then summarize the setup of our
likelihood analysis in Subsection IV.4.

###  IV.1 Halo-based forward model

In this paper, we adopt a forward modeling approach to predict cluster
observables. Our modeling pipeline relies heavily on the Dark Emulator, an
efficient and accurate emulation tool for statistics of dark matter halos
[25]. The Dark Emulator enables us to predict the halo mass function, halo
autocorrelations, and halo-mass cross-correlations as functions of redshift,
halo masses, and cosmological parameters in modest CPU time. As in the Dark
Emulator, we adopt $M\equiv
M_{\mathrm{200b}}=4\pi(R_{\mathrm{200b}})^{3}\bar{\rho}_{\mathrm{m0}}\times
200/3$ for our fiducial halo mass definition, where $R_{\mathrm{200b}}$ is the
spherical halo boundary radius within which the mean mass density is 200 times
the present-day mean mass density. Note that the mass $M$ is not necessarily
equal to the spherical overdensity mass, $M_{\mathrm{200b}}$, when accounting
for baryonic effects as described in Section IV.2. Another spherical
overdensity mass of $M_{\mathrm{500c}}$ is widely adopted for
characterizations of cluster masses in the literature. Within a halo-model
framework provided by the Dark Emulator, we define $M_{\mathrm{500c}}$ as

$\displaystyle M_{\mathrm{500c}}$ $\displaystyle\equiv$
$\displaystyle\frac{4\pi}{3}\,500\,\rho_{\mathrm{crit}}(z)\,R^{3}_{\mathrm{500c}}\,(1+z)^{-3}$
(17) $\displaystyle=$
$\displaystyle\int_{0}^{R_{\mathrm{500c}}}\mathrm{d}r\,4\pi
r^{2}\,\rho_{\mathrm{hm}}(r,M,z,\mathbf{c}),$ (18)

where $\rho_{\mathrm{crit}}(z)$ is the critical density in the universe at a
redshift $z$, $R_{\mathrm{500c}}$ is the halo boundary radius in comoving
coordinates and $\rho_{\mathrm{hm}}(r,M,z,\mathbf{c})$ represents the average
halo mass density profile at the halo-centric radius $r$, as a function of the
redshift $z$, halo mass $M$ and a set of cosmological parameters $\mathbf{c}$,
given by the Dark Emulator. In the following, we omit the cosmological
parameter dependence of halo statistics for simplicity.

Given inputs from the Dark Emulator, a model prediction of our stacked lensing
measurements contains the integral of

$\displaystyle\int\mathrm{d}\log M_{\mathrm{UPP}}\int\mathrm{d}M\,S(\log
M_{\mathrm{UPP}},z)\,{\cal P}(\log M_{\mathrm{UPP}}|M,z)$
$\displaystyle\qquad\qquad\qquad\qquad\times
n_{h}(M,z)\,\Delta\Sigma_{h}(R,M,z),$ (19)

where $S(\log M_{\mathrm{UPP}},z)$ represents the selection function of our
cluster sample, ${\cal P}(\log M_{\mathrm{UPP}}|M,z)$ sets the probability
distribution of the SZ mass, $M_{\mathrm{UPP}}$, for halos with a given mass
$M$ at redshift $z$, $n_{h}(M,z)$ and $\Delta\Sigma_{h}(R,M,z)$ are the halo
mass function and the excess surface mass density profile for a single halo,
respectively. The Dark Emulator can predict two of $n_{h}(M,z)$ and
$\Delta\Sigma_{h}(R,M,z)$ in an efficient way, while the functional form of
${\cal P}(\log M_{\mathrm{UPP}}|M,z)$ must be given in advance. A detailed
description of two key ingredients ${\cal P}(\log M_{\mathrm{UPP}}|M,z)$ and
$\Delta\Sigma_{h}(R,M,z)$ is found below.

#### Scaling relation between the SZ and halo masses

Following our previous work [22], we assume a log-normal distribution for
${\cal P}(\log M_{\mathrm{UPP}}|M)$;

$\displaystyle{\cal P}(\log M_{\mathrm{UPP}}|M)\,\mathrm{d}\log
M_{\mathrm{UPP}}=\frac{\mathrm{d}\log
M_{\mathrm{UPP}}}{\sqrt{2\pi}\sigma_{\log M_{\mathrm{UPP}}}}$
$\displaystyle\qquad\qquad\times\exp\left[-\frac{[\log M_{\mathrm{UPP}}-\log
M^{\prime}_{\mathrm{UPP}}(M)]^{2}}{2\sigma^{2}_{\log
M_{\mathrm{UPP}}}}\right],$ (20)

where $\log M^{\prime}_{\mathrm{UPP}}(M)$ sets the mean relation of $\log
M_{\mathrm{UPP}}$ given halo masses and $\sigma_{\log M_{\mathrm{UPP}}}$
represents the rms in the log-normal distribution. Because the SZ mass
$M_{\mathrm{UPP}}$ is introduced in terms of the spherical overdensity mass
$M_{\mathrm{500c}}$, we parameterize the mean relation of $\log
M^{\prime}_{\mathrm{UPP}}(M)$ as guided by the $M_{\mathrm{500c}}-M$ relation
in Eq. (18). To be specific, we adopt the functional form below;

$\displaystyle\log M^{\prime}_{\mathrm{UPP}}(M)$ $\displaystyle=$
$\displaystyle\log a_{0}+\log M$ (21)
$\displaystyle\qquad\qquad+a_{1}\left(\log M-14.5\right),$

where $M_{\mathrm{UPP}}$ and $M$ are set in units of $M_{\odot}$ in Eq. (21).
For the standard $\Lambda$CDM cosmology, we find that the right-hand side of
Eq. (21) with $0.45<a_{0}<0.6$ and $a_{1}\sim-0.01$ provides a reasonable fit
to the Dark-Emulator prediction of $M_{\mathrm{500c}}-M$ relations at $0<z<1$.

Apart from the mean relation, the scatter of $\sigma_{\log M_{\mathrm{UPP}}}$
has a critical role in comparisons with the model and stacked lensing
measurements. Because of a strong degeneracy between the mean and scatter, it
is difficult to place a meaningful constraint on the scatter in cluster
observable-mass relations with stacked lensing measurements alone [61, 68].
Hence, we introduce an informative prior in the scatter in our model by using
the Illustris TNG simulations [24]. We estimate the covariance matrix of gas
pressure profiles in cluster-sized halos in the simulations by using the
Kernel Localized Linear Regression model [69, 70], which allows us to infer
mass-conditioned estimates of the gas pressure profile. Given the covariance
matrix, we can generate random realizations of the gas pressure profile at a
fixed halo mass and redshift assuming the log-normality of the pressure
profile. More details about our calibration of the covariance model are given
in Appendix D.

To estimate the scatter, $\sigma_{\log M_{\mathrm{UPP}}}$, we take a five-step
procedure as follows;

  * (i)

set a cluster redshift $z$,

  * (ii)

set a cluster mass $\Upsilon=\log(M_{\mathrm{500c}}[h^{-1}M_{\odot}])$ by
drawing a uniform random number between $\Upsilon=13.5$ and 16.0,

  * (iii)

produce a random realization of the gas pressure profile $P_{\mathrm{gas}}(r)$
with the mean profile being the UPP at the redshift $z$ and logarithmic mass
$\Upsilon$ (where $r$ stands for the halo-centric radius and we here assume
$M_{\mathrm{500c}}=M_{\mathrm{UPP}}$ for simplicity),

  * (iv)

perform a least chi-square fitting of $P_{\mathrm{gas}}(r)$ with the UPP
profile varying a single mass parameter $M_{\mathrm{UPP}}$ and,

  * (v)

repeat the steps (ii)-(iv) 100,000 times and then estimate the standard
deviation of $\log M_{\mathrm{UPP,bestfit}}-\log M_{\mathrm{UPP,input}}$ where
$M_{\mathrm{UPP,input}}$ and $M_{\mathrm{UPP,bestfit}}$ are the input SZ mass
parameter given in step (ii) and the best-fit estimate inferred in step (iv),
respectively. We then identify the standard deviation with the scatter,
$\sigma_{\log M_{\mathrm{UPP}}}$.

We also examine a possible redshift dependence of $\sigma_{\log
M_{\mathrm{UPP}}}$ with the above procedure by varying the cluster redshift
between $z=0.05$ and $z=1.50$. In the end, we find that the estimated scatter
from our random realizations of the gas pressure profiles can be described as

$\displaystyle\sigma^{2}_{\mathrm{ref}}(\log
M_{\mathrm{UPP}},z)=\alpha(z)\left\\{\log(M_{\mathrm{UPP}}\,[M_{\odot}])-15\right\\}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\beta(z),$
(22) $\displaystyle\alpha(z)=-1.13\times 10^{-3}\,z+5.73\times 10^{-3},$ (23)
$\displaystyle\beta(z)=-8.79\times 10^{-5}\,z+1.54\times 10^{-2},$ (24)

where Eq. (22) is valid for $M_{\mathrm{UPP}}\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}10^{14}\,M_{\odot}$ and $0.05\leq
z\leq 1.50$. Although we expect that Eq. (22) will provide a reasonable
estimate of $\sigma_{\log M_{\mathrm{UPP}}}$, we introduce a single nuisance
parameter, $C_{0}$, to account for any differences with the simulation-based
model and real observables;

$\displaystyle\sigma_{\log M_{\mathrm{UPP}}}(\log
M_{\mathrm{UPP}},z)=C_{0}\sigma_{\mathrm{ref}}(\log M_{\mathrm{UPP}},z),$ (25)

where $C_{0}$ is our model parameter. In practice, we adopt
$M_{\mathrm{UPP}}=0.5M$ in Eq. (25) to set $\sigma_{\log M_{\mathrm{UPP}}}$,
allowing us to maintain the log-normal distribution of Eq. (20).

Hence we model the mass scaling relation in Eq. (20) by three model parameters
of $a_{0},a_{1}$ and $C_{0}$. We will explore the best-fitting model
parameters that can reproduce both the abundance and lensing profile
measurements simultaneously. Information on the cluster abundance is required
to break a degeneracy between $a_{0}$ and $C_{0}$. Because we perform a
likelihood analysis in each of the three cluster redshift bins independently,
we ignore possible redshift dependence of the mass scaling relation for
simplicity. It is worth noting that our lensing measurement is not designed to
extract the mass dependence of the excess surface mass density $\Delta\Sigma$.
Therefore, the posterior distribution of $a_{1}$ is mostly determined by the
prior information. Nevertheless, there are no strong motivations to fix
$a_{1}$ from an observational point of view at present. To make our analysis
conservative, we decide to vary $a_{1}$ over a wide range.

#### Cluster number count

To break a parameter degeneracy between $a_{0}$ and $C_{0}$, we include
information on the cluster abundance for our parameter inference. Given the
completeness of our cluster search, the cluster abundance in the redshift
range of $z_{\mathrm{min}}\leq z\leq z_{\mathrm{max}}$ can be written as

$\displaystyle N_{\mathrm{cl}}(z_{\mathrm{min}},z_{\mathrm{max}})=4\pi
f_{\mathrm{sky}}\,\int_{z_{\mathrm{min}}}^{z_{\mathrm{max}}}\mathrm{d}z\,\frac{c\,D^{2}_{A}(z)}{H(z)}$
$\displaystyle\times\int\mathrm{d}\log M_{\mathrm{UPP}}\int\mathrm{d}M\,S(\log
M_{\mathrm{UPP}},z)\,{\cal P}(\log M_{\mathrm{UPP}}|M,z)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\times
n_{h}(M,z),$ (26)

where $f_{\mathrm{sky}}=0.0115$ is the area fraction of our cluster search
region, $H(z)$ is the Hubble parameter at a redshift $z$, $D_{A}(z)$ is the
comoving angular diameter distance to the redshift $z$, $S(\log
M_{\mathrm{UPP}},z)$ is the completeness of our ACT cluster catalog as shown
in Fig. 2, ${\cal P}(\log M_{\mathrm{UPP}}|M,z)$ represents the probability
distribution given in Eq. (20), and $n_{h}(M,z)$ is the halo mass function at
the redshift $z$ and halo mass $M$ computed by the Dark Emulator. In Eq. (26),
we set lower and upper limits in the double mass integrals to be
$10^{13.5}\,[h^{-1}M_{\odot}]$ and $10^{16}\,[h^{-1}M_{\odot}]$, respectively.

#### Stacked lensing signal

The stacked lensing profile for halos with mass $M$ and at redshift $z$ probes
the average radial profile of matter distribution around the halos,
$\rho_{\mathrm{hm}}(r,M,z)$. Because a direct output from the Dark Emulator is
the cross-correlation function between the halo distribution and the matter
density fluctuation field $\xi_{\mathrm{hm}}(r,M,z)$, it is useful to express
the average mass density profile in terms of $\xi_{\mathrm{hm}}(r,M,z)$:

$\displaystyle\rho_{\mathrm{hm}}(r,M,z)=\bar{\rho}_{\mathrm{m0}}\left[1+\xi_{\mathrm{hm}}(r,M,z)\right].$
(27)

The cross-correlation is related to the cross-power spectrum
$P_{\mathrm{hm}}(k,M,z)$ via the Fourier transform as

$\displaystyle\xi_{\mathrm{hm}}(r,M,z)=\bar{\rho}_{\mathrm{m0}}\int_{0}^{\infty}\frac{k\mathrm{d}k}{2\pi}P_{\mathrm{hm}}(k,M,z)j_{0}(kr),$
(28)

where $j_{0}(x)$ is the zeroth-order spherical Bessel function. For a given
$P_{\mathrm{hm}}(k)$, the excess surface mass density profile for a single
halo is computed as

$\displaystyle\Delta\Sigma_{h}(R,M,z)=\bar{\rho}_{\mathrm{m0}}\int_{0}^{\infty}\frac{k\mathrm{d}k}{2\pi}P_{\mathrm{hm}}(k,M,z)J_{2}(kR),$
(29)

where $J_{2}(x)$ is the second-order Bessel function.

Taking into account the distribution of halo masses and redshifts for the
clusters as well as the cluster completeness, we can compute a model
prediction for the stacked lensing profile as

$\displaystyle\Delta\Sigma(R;z_{\mathrm{min}},z_{\mathrm{max}})=\frac{4\pi
f_{\mathrm{sky}}}{N_{\mathrm{cl}}(z_{\mathrm{min}},z_{\mathrm{max}})}\int_{z_{\mathrm{min}}}^{z_{\mathrm{max}}}\mathrm{d}z\,\frac{c\,D^{2}_{A}(z)}{H(z)}$
$\displaystyle\times\int\mathrm{d}\log M_{\mathrm{UPP}}\int\mathrm{d}M\,S(\log
M_{\mathrm{UPP}},z)\,{\cal P}(\log M_{\mathrm{UPP}}|M,z)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\times
n_{h}(M,z)\,\Delta\Sigma_{h}(R,M,z).$ (30)

In practice, our choice of a proxy for the cluster center might not match the
true center. To account for this off-centering effect, we follow the model
proposed in Ref. [71] by modifying the halo-matter cross-power spectrum as

$\displaystyle P_{\mathrm{hm}}(k,M,z)$ $\displaystyle\rightarrow$
$\displaystyle{\cal F}_{\mathrm{off}}(k,M,z)\,P_{\mathrm{hm}}(k,M,z),$ (31)
$\displaystyle{\cal F}_{\mathrm{off}}(k,M,z)$ $\displaystyle=$ $\displaystyle
f_{\mathrm{cen}}(M)+[1-f_{\mathrm{cen}}(M)]$ (32)
$\displaystyle\qquad\qquad\qquad\times\tilde{p}_{\mathrm{off}}(k;R_{\mathrm{off}}(z)),$

when computing Eq. (30) in the redshift range of interest. Here
$f_{\mathrm{cen}}$ is the fraction of clusters with no off-centering effects,
while $1-f_{\mathrm{cen}}$ is the fraction of off-centered clusters. Motivated
by a possible mass dependence in $f_{\mathrm{cen}}$ (e.g. see Ref. [61] for
the richness dependence), we adopt a functional form of
$f_{\mathrm{cen}}(M)=f_{0}+f_{1}\ln(M/M_{\mathrm{pivot}})$ where $f_{0}$ and
$f_{1}$ are free parameters to be varied and $M_{\mathrm{pivot}}$ is set to
$3\times 10^{14}\,h^{-1}M_{\odot}$. The function
$p_{\mathrm{off}}(r;R_{\mathrm{off}})$ is the normalized radial profile of
off-centered “central” galaxies with respect to the true center for which we
assume a Gaussian distribution given by
$p_{\mathrm{off}}(r;R_{\mathrm{off}})=\exp(-r^{2}/2/R^{2}_{\mathrm{off}})/[(2\pi)^{3/2}R^{3}_{\mathrm{off}}]$,
where $R_{\mathrm{off}}$ is a parameter to model the typical off-centering
radius. The Fourier transform of $p_{\mathrm{off}}(r;R_{\mathrm{off}})$ is
then given by
$\tilde{p}_{\mathrm{off}}(k;R_{\mathrm{off}})=\exp(-k^{2}R^{2}_{\mathrm{off}}/2)$.
We parametrize the off-centering radius as
$R_{\mathrm{off}}(z)=\alpha_{\mathrm{off}}D_{A}(z)\sigma_{\mathrm{off}}$,
where $\alpha_{\mathrm{off}}$ is a free parameter and $\sigma_{\mathrm{off}}$
is a reference scale for recovered position offsets obtained from source
insertion simulations as done in Ref. [12]. Based on the best estimate
inferred from the simulations in Ref. [12], we set
$\sigma_{\mathrm{off}}=0.28\,\mathrm{arcmin}$ throughout this paper. In
addition to the parameters in the mass scaling relation, we thus have three
free parameters, $f_{0}$, $f_{1}$ and $\alpha_{\mathrm{off}}$, to be varied in
our analysis.

###  IV.2 Baryonic effects

Baryonic feedback can affect the mass density distribution in and around dark
matter halos. In cluster-sized halos, feedback processes from active galactic
nuclei (AGN) would expel intracluster material from the cluster center,
modifying cluster mass density profiles even at Mpc scales (see Ref. [18] for
a recent review). To account for possible baryonic effects on the cluster mass
profiles and to marginalize these effects in our parameter inference, we adopt
a semi-analytic model developed in Refs. [65, 66, 67]. The model is widely
referred to as the baryonification method in the literature; it is designed to
include baryonic effects in gravity-only N-body simulations by shifting
positions of N-body particles with respect to centers of various dark matter
halos. We apply the baryonification method to the outputs of the Dark
Emulator, increasing flexibility in our model of stacked lensing signals.

The baryonification method aims at the transformation of dark-matter-only
(dmo) profiles into modified, dark-matter-baryon (dmb) profiles that include
the effects of dark matter, gas, and stars;

$\displaystyle\rho_{\mathrm{hm+t}}(r,M)\rightarrow\rho_{\mathrm{hm,dmb}}(r,M)$
$\displaystyle\qquad\quad=\rho_{\mathrm{DM}}(r,M)+\rho_{\mathrm{gas}}(r,M)+\rho_{\mathrm{star}}(r,M),$
(33)

where $\rho_{\mathrm{hm+t}}(r,M)$ represents the average mass density profile
at halo mass $M$ predicted by the Dark Emulator, but applying a truncation at
outer regions, $\rho_{\mathrm{hm,dmb}}(r,M)$ consists of dark matter
($\rho_{\mathrm{DM}}$), gas ($\rho_{\mathrm{gas}}$), and stellar density
profile ($\rho_{\mathrm{star}}$). Some truncation needs to be applied to the
original output of the Dark Emulator $\rho_{\mathrm{hm}}$, because the
baryonification assumes that the density profiles of any component in a halo
must be constrained to have a converged enclosed mass in the limit of an
infinite boundary radius [65]. Note that the Dark Emulator provides the
prediction without the common decomposition into one- and two-halo terms,
making the truncation somewhat arbitrary. In this paper, we perform a profile
fitting to the Dark-Emulator output as in Ref. [72] and then find an outer
radius $r_{\mathrm{out}}$ for which the logarithmic slope is steepest. Using
the estimate of the boundary radius $r_{\mathrm{out}}$, we then correct the
Dark-Emulator prediction with a truncation below [73, 74];

$\displaystyle\rho_{\mathrm{hm+t}}(r)=\left\\{\begin{array}[]{ll}\rho_{\mathrm{hm}}(r)&(r\leq
r_{\mathrm{out}})\\\
\rho_{\mathrm{hm}}(r_{\mathrm{out}})\left(\frac{r}{r_{\mathrm{out}}}\right)^{\gamma}\left(\frac{r^{2}_{\mathrm{out}}+\tau^{2}}{r^{2}+\tau^{2}}\right)^{2}&(r>r_{\mathrm{out}})\end{array}\right.,$
(36)

where $\gamma$ is the density slope of $\rho_{\mathrm{hm}}$ at
$r=r_{\mathrm{out}}$ and we set $\tau=3r_{\mathrm{out}}$. More details about
the truncation and the explicit functional forms of $\rho_{\mathrm{DM}}$,
$\rho_{\mathrm{gas}}$ and $\rho_{\mathrm{star}}$ are found in Appendix E.

Given the corrected density profiles, we can compute the enclosed mass profile
as $M_{\chi}(r)\equiv\int_{0}^{r}4\pi s^{2}\mathrm{d}s\,\rho_{\chi}(s)$
($\chi=\mathrm{hm+t},\mathrm{DM},\mathrm{gas},\mathrm{star}$). The
baryonification method then defines a Lagrangian mapping from the initial
(dmo) to the final (dmb) density profile by shifting the location of a
spherical mass shell as

$\displaystyle\frac{r_{f}}{r_{i}}-1$ $\displaystyle=$ $\displaystyle{\cal
A}\left[\left(\frac{M_{i}}{M_{f}}\right)^{n}-1\right],$ (37) $\displaystyle
M_{i}$ $\displaystyle=$ $\displaystyle M_{\mathrm{hm}}(r_{i})$ (38)
$\displaystyle M_{f}$ $\displaystyle=$ $\displaystyle
M_{\mathrm{DM}}(r_{i})+M_{\mathrm{gas}}(r_{f})+M_{\mathrm{star}}(r_{f}),$ (39)

where $r_{i}$ and $r_{f}$ represent the initial and final locations of the
mass shell from the halo center, respectively. In Eq. (37), we have two
parameters, ${\cal A}$ and $n$, governing the adiabatic contraction and
expansion of the mass shells. Throughout this paper, we set ${\cal A}=0.3$ and
$n=2$ as in Ref. [66] (Note that these values can provide a fit to the
simulation results in Ref. [75]). One can solve Eq. (37) for
$\zeta=r_{f}/r_{i}$ iteratively, allowing to have the displacement length
(i.e. $d[r_{i}]=r_{f}-r_{i}$) at a given $r_{i}$. In the end, the halo-matter
cross-correlation with baryonic effects can be expressed as

$\displaystyle\xi_{\mathrm{hm+b}}(r,M,z,\mathbf{p})$ $\displaystyle=$
$\displaystyle\xi_{\mathrm{1h}}(r,M,z,\mathbf{p})+\xi_{\mathrm{2h}}(r,M,z),$
(40) $\displaystyle\xi_{\mathrm{1h}}(r,M,z,\mathrm{p})$ $\displaystyle=$
$\displaystyle\frac{\rho_{\mathrm{hm,dmb}}(r,M,z,\mathbf{p})}{\bar{\rho}_{\mathrm{m0}}},$
(41) $\displaystyle\xi_{\mathrm{2h}}(r,M,z)$ $\displaystyle=$
$\displaystyle\frac{\rho_{\mathrm{hm}}(r,M,z)-\rho_{\mathrm{hm+t}}(r,M,z)}{\bar{\rho}_{\mathrm{m0}}}$
(42) $\displaystyle\qquad\qquad\qquad\qquad\qquad-1$

where $\mathbf{p}$ represents a set of parameters in the baryonification
method. We here assume that an effective two-halo term of
$\rho_{\mathrm{hm}}-\rho_{\mathrm{hm+t}}$ is not affected by baryonic effects.
When including the baryonic effects in the model of $\Delta\Sigma$, we first
perform the Fourier transform of Eq. (40) to have the corresponding halo-
matter power spectrum $P_{\mathrm{hm+b}}$ and then plug in $P_{\mathrm{hm+b}}$
for $P_{\mathrm{hm}}$ in Eq. (29). It is worth noting that the baryonification
method does not introduce any changes in the halo mass function $n_{h}(M)$ as
long as we distinguish each halo by a single halo mass $M$. After the
baryonification is applied, the halo mass $M$ does not mean the spherical halo
mass of $M_{\mathrm{200b}}$ as designed in the Dark Emulator, but should be
regarded as an auxiliary variable to identify individual halos.

The baryonification method introduces various parameters. Because our primary
interest is to measure the masses of cluster-sized halos, it is reasonable to
vary only a small set of the model parameters which are relevant to baryonic
effects on mass density fields at Mpc scales. As examined in Refs. [76, 67],
the baryon mass fraction inside halos is most relevant to baryonic effects on
Mpc-scale matter clustering (also see Refs. [77, 78] for baryonic effects on
cluster mass estimates). Motivated by these previous findings, we allow three
parameters of gas physics to vary in our analysis. These include (i) the mass
dependence $\mu$ in the inner density slope of $\rho_{\mathrm{gas}}$, (ii) the
pivot mass scale $M_{c}$ at which the inner density slope becomes 3/2, (iii)
and the dimensionless outer scale radius of $\rho_{\mathrm{gas}}$ expressed as
$\theta_{\mathrm{ej}}$. See Appendix E for details.

Figure 6:  Parameter dependence of stacked lensing signals $\Delta\Sigma$ for
the ACT clusters at redshifts of $0.092<z\leq 0.445$. The left, middle, and
right panels show $\Delta\Sigma$ when we vary the parameters related to the
mass-scaling relation, off-centering effects, and baryonic effects,
respectively. The top left panel shows the change in $\Delta\Sigma$ caused by
varying the parameter in the mean scaling relation $a_{0}$ and the scatter
$C_{0}$, while the top middle panel summarizes the change by the off-centering
parameters $f_{0}$ and $\alpha_{\mathrm{off}}$ (see Section IV.1 for the
definition of those parameters.). The top right panel represents
$\Delta\Sigma$ for several baryon models based on the baryonification method
as described in Section IV.2. Note that we plot $R\,\Delta\Sigma$ in the top
panels for visualization purposes. The baseline model is shown as the black
solid line in each panel: it does not include any baryonic effects. See
Section IV.3 for the choice of model parameters in this figure. The bottom
panels show the ratio of $\Delta\Sigma$ to our baseline model.

###  IV.3 Parameter dependence of $\Delta\Sigma$

To provide intuition, the dependence of stacked lensing signals on the model
parameters are briefly summarized here. Fig. 6 shows our model of stacked
lensing signal $\Delta\Sigma$ for the ACT clusters in the redshift interval of
$0.092<z\leq 0.445$ as we vary some model parameters. In this figure, we
consider eight cases as follows:

  * (i)

The fiducial model, assuming a mass-scaling relation with $a_{0}=0.5$ and
$a_{1}=0$; scatter in the scaling relation set at $C_{0}=1$; and the off-
centering effects with $f_{0}=0.75$, $f_{1}=0.25$, and
$\alpha_{\mathrm{off}}=1$. This model does not include any baryonic effects;

  * (sl-i)

the model assuming a mass-scaling relation set to a higher $a_{0}$ by 0.1;

  * (sl-ii)

the model with a higher scatter in the mass scaling relation, changing $C_{0}$
by 1.0;

  * (oc-i)

the model with fewer off-centered clusters, with a higher $f_{0}$ by 0.2

  * (oc-ii)

the model with a larger off-centering length with a higher
$\alpha_{\mathrm{off}}$ by 2.0

  * (b-i)

The fiducial model of the mass-scaling relation and off-centering effects, but
including baryonic effects with $\mu=1.0$,
$M_{c}=10^{13.5}\,[h^{-1}M_{\odot}]$, and $\theta_{\mathrm{ej}}=3$;

  * (b-ii)

The fiducial model of the mass-scaling relation and off-centering effects, and
the baryonic effects with $\mu$ increased by 0.5;

  * (b-iii)

The fiducial model of the mass-scaling relation and off-centering effects, and
baryonic effects with $\theta_{\mathrm{ej}}$ increased by 5;

The left panel in Fig. 6 displays the models labeled above as (sl-i) and (sl-
ii); the middle panel shows models (oc-i) and (oc-ii); and models (b-i),
(b-ii), and (b-iii) are summarized in the right panel. Note that the fiducial
model of (i) is shown in every panel by the black solid line for ease of
comparison.

#### Mass-scaling relation

The left panel in Fig. 6 highlights the dependence of $\Delta\Sigma$ on the
mass-scaling relation. In our model, the scaling relation between SZ and total
masses is set by three parameters: $a_{0}$, $a_{1}$, and $C_{0}$. The stacked
lensing signal is most sensitive to two parameters, $a_{0}$ and $C_{0}$.
Because the completeness of our ACT clusters is fixed as shown in Fig. 2, the
average SZ mass in the cluster sample in a given redshift range is not
affected by the choice of our model parameters. Hence, increasing $a_{0}$
indicates that smaller halo masses become more relevant to the average SZ
mass, making the amplitude of $\Delta\Sigma$ smaller. Also, increasing the
scatter in the mass-scaling relation (i.e. a higher $C_{0}$ in our model) can
lead to a smaller $\Delta\Sigma$, because the larger scatter can allow less
massive clusters to contribute to the stacked lensing signal. As shown in the
left panel, there exists a strong degeneracy between the two parameters
$a_{0}$ and $C_{0}$ in our model, which has been mentioned in the literature
[61, 68]. It is worth noting that both $a_{0}$ and $C_{0}$ can affect the
signal not only at $R\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}10\,h^{-1}\mathrm{Mpc}$ but also at
a large separation of $R\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}10\,h^{-1}\mathrm{Mpc}$. This
effect can be explained by the change in the linear halo bias as one varies
the model parameter [68].

#### Off-centering effects

The middle panel in Fig. 6 summarizes off-centering effects in our model of
stacked lensing measurements. Our model characterizes the off-centering
effects by three parameters $f_{0}$, $f_{1}$ and $\alpha_{\mathrm{off}}$; of
these, $f_{0}$ and $\alpha_{\mathrm{off}}$ are most relevant in our case. The
parameter $f_{0}$ is the fraction of clusters with no off-centering at halo
mass of $3\times 10^{14}\,h^{-1}M_{\odot}$. Hence, increasing $f_{0}$ can make
the prediction of $\Delta\Sigma$ at $R\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}1\,h^{-1}\mathrm{Mpc}$ higher. At a
fixed $f_{0}$, we can see how a typical off-centering length affects the
stacked lensing signal by varying the parameter $\alpha_{\mathrm{off}}$.
Increasing $\alpha_{\mathrm{off}}$ means the stacked lensing signal can be
affected by off-centering effects to larger radii.

#### Baryonic effects

The right panel in Fig. 6 represents baryonic effects in the stacked lensing
signal $\Delta\Sigma$ in our model. By comparing the solid line and the blue
dashed line in the figure, we can see the impact of baryonic effects;
including baryons can lead to a more concentrated halo mass profile (see also
Ref. [79]), but it does not affect the signal at large radii which are
governed by gravity alone. We have three parameters in the model of baryonic
effects, enabling the prediction to be highly flexible. For instance,
increasing the parameter $\mu$ at a given $M_{c}$ steepens the inner slope of
the gas density profile, leading to an increase in the mass density
concentration. Also, increasing the parameter $\theta_{\mathrm{ej}}$ indicates
that the gas in a cluster can be distributed more in outer regions beyond the
cluster’s virial radius, allowing the mass density profile to experience some
adiabatic expansion. As shown in the figure, baryonic effects may affect the
mass density profile at Mpc scales, which can not be ignored when infering the
average halo mass from the stacked lensing measurements. Nevertheless, the
information on $\Delta\Sigma$ at large radii can be helpful to constrain the
halo mass even if our model is subject to uncertainty in baryonic physics.

Figure 7:  The observed stacked lensing signals around ACT SZ-selected
clusters and the best-fit model derived from the MCMC analysis. We work with
three different cluster redshift bins (increasing from left to right). In the
upper panels, the colored points represent the observed signal of
$\Delta\Sigma$, while the line corresponds to the best-fit model. The error
bars in each upper panel are estimated with a combination of two mock lensing
analyses (see Section III.3 for details). In the lower panels, we plot the
difference between the observed $\Delta\Sigma$ signal and the model
counterpart. There, we normalize the difference by the diagonal component of
the covariance matrix. The gray region in the lower panels highlights a $\pm
1\sigma$-level difference. For reference, we note the chi-square value
$\chi^{2}$ as defined in Eq. (45), the number of the data elements
$n_{\mathrm{bin}}$, the number of free parameters in the model
$n_{\mathrm{p}}$, and the $p$ value of the $\chi^{2}$ value with respect to
the number of degrees of freedom (i.e. $n_{\mathrm{bin}}-n_{\mathrm{p}}$) in
each upper panel. Our model adopts the Dark Emulator [25] for a baseline
prediction of statistics of cluster-sized halos, and it also takes into
account possible baryonic effects with the prescription of baryonification
[65, 66, 67]. See Section IV.1 and IV.2 for further details on our halo-based
model.

###  IV.4 Likelihood

To constrain the model parameters with the stacked lensing measurements, we
define a log-likelihood function as

$\displaystyle-2\ln{\cal L}(\mbox{\boldmath$\theta$})$ $\displaystyle=$
$\displaystyle\sum_{i,j}d_{i}\mathrm{Cov}^{-1}_{ij}d_{j}+\left(\frac{N_{\mathrm{cl,obs}}-N_{\mathrm{cl}}(\mbox{\boldmath$\theta$})}{2N^{1/2}_{\mathrm{cl,obs}}}\right)^{2},$
(43) $\displaystyle d_{i}$ $\displaystyle\equiv$
$\displaystyle\Delta\Sigma_{\mathrm{obs}}(R_{i})-\Delta\Sigma(R_{i};\mbox{\boldmath$\theta$})$
(44)

where $\Delta\Sigma_{\mathrm{obs}}(R_{i})$ represents the stacked lensing
signal at the i-th radial bin from the cluster center,
$N_{\mathrm{cl,obs}}=32$ is the observed cluster number count,
$\Delta\Sigma(R_{i})$ and $N_{\mathrm{cl}}$ are the cluster number counts
predicted by our halo-based model (see Eq. [26]), and $\theta$ describes a set
of our model parameters. In Eq. (43), we assume a double Poisson error for the
cluster number count to account for possible uncertainties in our model (e.g.
our model ignores a variation of the survey depth in the ACT observation
[12]). Note that the statistical errors on the cluster abundance estimated by
our mock catalogs are consistent with the Poisson expectation with a level of
10%, while the cross covariance between $\Delta\Sigma_{\mathrm{obs}}$ and
$N_{\mathrm{cl}}$ in the mock catalogs is found to be less than 0.1 after we
normalize the cross covariance with the variances of
$\Delta\Sigma_{\mathrm{obs}}$ and $N_{\mathrm{cl}}$. Hence, our likelihood
function primarily uses the information of lensing signals, but it has a
useful contribution from the cluster number counts to break possible
degeneracy among our model parameters.

In the log-likelihood function, we have nine free parameters apart from
cosmological parameters. Six of these parameters set the mass-scaling relation
($a_{0}$, $a_{1}$, and $C_{0}$) and possible off-centering effects in the ACT
cluster sample ($f_{0}$, $f_{1}$ and $\alpha_{\mathrm{off}}$). The other three
parameters control baryonic effects on the stacked lensing signal ($\mu$,
$M_{c}$, and $\theta_{\mathrm{ej}}$). We set flat priors in the ranges of
$0.01\leq a_{0}\leq 0.75$, $-0.2\leq a_{1}\leq 0.2$, $0.2\leq C_{0}\leq 5.0$,
$0\leq f_{0}\leq 1$, $0\leq f_{1}\leq 1$, $0\leq\alpha_{\mathrm{off}}\leq 3$,
$0\leq\mu\leq 2$, $12.0\leq\log(M_{c}\,[h^{-1}M_{\odot}])\leq 15.0$ and
$1\leq\theta_{\mathrm{ej}}\leq 10$ in the likelihood analysis. To compute the
log-likelihood function in our nine-parameter space, we use the Markov Chain
Monte Carlo (MCMC) sampler ${\tt emcee}$ [80]. We set 45 walkers, ran 32000
steps of burn-in, and then 32000 steps to sample the likelihood function. We
confirmed that sampling within and among chains has converged to a level of
$1.5\%$. In this paper, we perform an MCMC analysis for the cluster subsample
at each of the three redshift bins separately.

##  V Results

###  V.1 Parameter inference

Figure 8:  Impact of the cluster abundance on parameter inference. In each
panel, we show the 68% and 95% credible intervals in the two-dimensional
parameter space ($a_{0}$ and $C_{0}$) for three different cluster redshift
bins. The solid colored lines show our fiducial MCMC results including cluster
abundance information, while the dashed ones are the counterparts without
information on the cluster number count. The colored symbol in each panel
highlights the best-fit parameters in our fiducial analysis, and the cross
symbol corresponds to the best-fit values when we do not use the abundance
information in the parameter inference. The strong degeneracy between $a_{0}$
and $C_{0}$ can make the best-fit values off-center in the contours after
marginalization in our high-dimensional parameter space when we work on the
MCMC analysis without the abundance prior.

We here summarize our parameter inference on the mass density profile around
the ACT clusters. Fig. 7 shows the comparison of stacked lensing signals
$\Delta\Sigma$ with our best-fit models. To assess the goodness of fit, we
define a chi-square quantity as

$\displaystyle\chi^{2}$ $\displaystyle=$
$\displaystyle\sum_{i,j}d_{i,\mathrm{bf}}\mathrm{Cov}_{ij}^{-1}d_{j,\mathrm{bf}},$
(45) $\displaystyle d_{i,\mathrm{bf}}$ $\displaystyle=$
$\displaystyle\Delta\Sigma_{\mathrm{obs}}(R_{i})-\Delta\Sigma(R_{i};\mbox{\boldmath$\theta$}_{\mathrm{bf}}),$
(46)

where $\mbox{\boldmath$\theta$}_{\mathrm{bf}}$ represents the best-fit
parameters obtained by our MCMC analysis. We find the chi-square values of
$\chi^{2}=11.62,15.68$ and $13.24$ at the cluster redshifts of $0.092<z\leq
0.445$, $0.445<z\leq 0.695$, and $0.685<z\leq 1.180$, respectively. Our model
consists of nine free parameters and the number of data elements is 18 for the
first two redshift bins and 15 for the highest redshift bin. Hence, our best-
fit model provides a reasonable chi-square value with respect to the number of
degrees of freedom. We expect that an effective number of degrees of freedom
would be larger than the naive expectation of
$n_{\mathrm{bin}}-n_{\mathrm{p}}$ where $n_{\mathrm{bin}}$ and
$n_{\mathrm{p}}$ are the number of data elements and free parameters, because
our measurements are not informative enough to constrain all the free
parameters. For instance, the parameters of the baryonic effects ($\mu$, $\log
M_{c}$, and $\theta_{\mathrm{ej}}$) are not well constrained with our analysis
alone. We provide credible intervals for the nine free parameters in Appendix
F, highlighting that our analysis can put meaningful constraints on only two
parameters, $a_{0}$ and $C_{0}$. Other parameters can be regarded as basically
nuisance parameters in our analysis, accounting for systematic uncertainties
in our modeling of the stacked lensing measurements. As a cross check, we
compute a probability-to-exceed (PTE) for our minimum $\chi^{2}$ by using
$10^{6}$ random data vectors with mean zero and the data covariance. We find
that the PTE values are 0.87, 0.62, and 0.58 in order of increasing redshifts,
showing our best-fit model provides a good fit to the data.

Figure 9:  Similar to Fig. 8. We compare credible intervals when accounting
for the baryonic effects in the model or not. The dashed lines in each panel
represent the 68% and 95% credible intervals when our model does not include
baryonic effects, and the cross symbol shows the best-fit parameters in this
case.

For reference purposes, we compute the cumulative signal-to-noise ratio of our
stacked lensing measurements as

$\displaystyle\left(\frac{\mathrm{S}}{\mathrm{N}}\right)^{2}=\sum_{i,j}\Delta\Sigma_{\mathrm{obs}}(R_{i})\mathrm{Cov}_{ij}^{-1}\Delta\Sigma_{\mathrm{obs}}(R_{j}),$
(47)

where the indices $i,j$ run over the range of
$0.1<R\,[h^{-1}\mathrm{Mpc}]<100$. We find $\mathrm{S/N}=14.6,12.0$ and $6.6$
at the cluster redshifts of $0.092<z\leq 0.445$, $0.445<z\leq 0.695$, and
$0.685<z\leq 1.180$, respectively. The cumulative signal-to-noise ratios
indicate that our measurements can constrain the amplitude of $\Delta\Sigma$
with a level of $6.8\%$, $8.3\%$, and $15\%$ from lower to higher redshift
bins. Our MCMC analysis provides a 68% credible interval of the primary
parameters ($a_{0}$ and $C_{0}$) as

$\displaystyle
a_{0}=\left\\{\begin{array}[]{ll}0.27^{+0.05}_{-0.07}&(0.092<z\leq 0.445)\\\
0.32\pm{0.05}&(0.445<z\leq 0.695)\\\ 0.35^{+0.05}_{-0.06}&(0.695<z\leq
1.180)\end{array}\right.,$ (51)

and

$\displaystyle C_{0}=\left\\{\begin{array}[]{ll}1.3^{+0.9}_{-0.7}&(0.092<z\leq
0.445)\\\ 0.7^{+0.6}_{-0.4}&(0.445<z\leq 0.695)\\\
0.7^{+0.7}_{-0.4}&(0.695<z\leq 1.180)\end{array}\right.,$ (55)

where we have a $15\%$-level constraints on $a_{0}$ with our lensing analysis
in each redshift bin. Note the peak in the posterior distribution of $C_{0}$
is found to be close to the lowest level in the prior ($C_{0}=0.2$) except for
the lowest redshift bin. Hence, our estimate of the 68% credible interval for
$C_{0}$ is meaningful in the lowest redshift bin alone. We also examined a
wide prior of $0.01\leq C_{0}\leq 5.0$ in our likelihood analysis, confirming
that our constraints of lensing masses are not affected by the lower bound of
the $C_{0}$ prior.

Figure 10:  Constraints of the mass bias
$B=M_{\mathrm{500c}}/M_{\mathrm{UPP}}$ from our stacked lensing measurements.
In three left panels, the open circle with error bars represents the 68%
credible interval of $B$ in the different cluster redshift bins (lower to
higher redshift from left to right). Note that we define the mass bias in
logarithmic space as $\log B\equiv\langle\log M_{\mathrm{500c}}|\log
M_{\mathrm{UPP}}\rangle-\log M_{\mathrm{UPP}}$ where $\langle\log
M_{\mathrm{500c}}|\log M_{\mathrm{UPP}}\rangle$ is the average of $\log
M_{\mathrm{500c}}$ at the average $M_{\mathrm{UPP}}$ (see Section V.2 for
details). For comparison, we also show the constraints on the mass bias by
previous lensing analyses in the literature; ACT $\times$ CS82 [81], Planck
$\times$ CLASH [82], Planck $\times$ HSCY1 [83], PSZ2LenS [84], ACT $\times$
HSCY1 [22], Planck $\times$ CCCP [85], SPT $\times$ HST [86], and ACT $\times$
KiDS-1000 [26]. Caution that the limits of Planck $\times$ CLASH, Planck
$\times$ HSCY1, PSZ2LenS, ACT $\times$ HSCY1 and SPT $\times$ HST are set by
the average masses of individual lensing measurements. The limits of ACT
$\times$ CS82, Planck $\times$ CCCP and ACT $\times$ KiDS-1000 should be used
for a fair comparison. For reference purposes, the right panel displays the
constraint on $B$ derived from the auto power spectrum of the Planck thermal
SZ map [87]. Note that some previous studies have used the inverse convention
of quoting $B^{-1}\equiv 1-b$.

Our baseline analysis adopts a loose prior in the cluster abundance as seen in
the second term on the right-hand side of Eq. (43). It is beneficial to see
how important the abundance prior is in our parameter inference. Fig. 8
summarizes the 68% and 95% credible intervals in the two-dimensional plane of
$a_{0}$ versus $C_{0}$ to highlight the impact of the abundance prior. In the
figure, the solid colored lines represent the constraints when we include the
abundance prior, while the dashed lines stand for the corresponding
constraints without the abundance information. For the MCMC analysis made
without the abundance prior, we double the number of MCMC steps to have a
stable parameter inference. We find that the abundance prior is efficient in
rejecting the parameter space at lower $a_{0}$ values in the MCMC analysis.
Note that a lower $a_{0}$ value in our model can predict a higher mean halo
mass at a given SZ mass, providing smaller cluster number counts in each
redshift bin. Hence, information on the cluster abundance can be effective to
break the inherent parameter degeneracy between $a_{0}$ and $C_{0}$ in the
$\Delta\Sigma$ model. We also observe that some parameter regions with large
$a_{0}$ and small $C_{0}$ values are incompatible with the cluster abundance
in the lowest cluster redshift bin of $0.092<z\leq 0.445$, because those
parameters provide an overestimate of the cluster abundance.

We also study the importance of including baryonic effects in the
$\Delta\Sigma$ model in our analysis. Fig. 9 shows the 68% and 95% credible
intervals in the $a_{0}-C_{0}$ plane when we include or exclude baryonic
effects in our modeling. We find that our measurements are not very sensitive
to details of the baryonic effects at present. Nevertheless, we observe a
trend: including baryonic effects increase the credible intervals in the two
lowest cluster redshift bins, implying that it will be more important to
account for baryons in parameter inference as the cumulative signal-to-noise
ratio increases.

###  V.2 Total mass at a given SZ mass

The SZ mass $M_{\mathrm{UPP}}$ is not necessarily equal to the underlying
cluster mass $M_{\mathrm{500c}}$, because the scaling relation between the
electron pressure and the SZ mass is derived under the assumption of
hydrostatic equilibrium. The UPP scaling relation was also derived from X-ray
observations of $z=0.1$ clusters [33]. Also, any unknown systematic errors in
X-ray analyses induce deviations of $M_{\mathrm{UPP}}$ from the true mass of
$M_{\mathrm{500c}}$ even if hydrostatic equilibrium is a good approximation on
cluster scales. We here study the mass bias of $B\equiv
M_{\mathrm{500c}}/M_{\mathrm{UPP}}$ at a fixed SZ mass with our MCMC parameter
inference based on the stacked lensing measurements.

In our halo-based forward model, we can estimate the probability distribution
of $M_{\mathrm{500c}}$ at a given $\log M_{\mathrm{UPP}}$ and redshift $z$ as

$\displaystyle\mathrm{Prob}(M_{\mathrm{500c}}=m|\log M_{\mathrm{UPP}},z)={\cal
N}^{-1}\int\mathrm{d}M\,n_{h}(M,z)$ $\displaystyle\qquad\qquad\times{\cal
P}(\log
M_{\mathrm{UPP}}|M)\,\delta_{\mathrm{D}}\left(f_{\mathrm{500c}}(M,z)-m\right),$
(56) $\displaystyle{\cal N}=\int\mathrm{d}M\,n_{h}(M,z)\,{\cal P}(\log
M_{\mathrm{UPP}}|M),$ (57)

where $\delta_{\mathrm{D}}(x)$ is the one-dimensional Dirac function and we
find the function of $f_{\mathrm{500c}}(M,z)$ as a solution of Eqs. (17) and
(18). Note that the probability density of Eq. (56) can be derived for a given
set of our model parameters. The parameters $a_{0}$, $a_{1}$ and $C_{0}$ set
the function ${\cal P}(\log M_{\mathrm{UPP}}|M)$ and the parameters involved
in baryonic effects ($\mu$, $\log M_{c}$ and $\theta_{\mathrm{ej}}$) affect
the function $f_{\mathrm{500c}}$ by changing the relation between the
spherical overdensity mass $M_{\mathrm{500c}}$ and the baseline halo mass $M$
as given by the Dark Emulator. Using Eq. (56), we then derive the average
logarithmic mass at a given SZ mass as

$\displaystyle\langle\log M_{\mathrm{500c}}|\log
M_{\mathrm{UPP}}\rangle\equiv\int\mathrm{d}M_{\mathrm{500c}}\,\log
M_{\mathrm{500c}}$
$\displaystyle\qquad\qquad\qquad\qquad\times\mathrm{Prob}(M_{\mathrm{500c}}|\log
M_{\mathrm{UPP}},z),$ (58)

where we use the average redshift over the 32 clusters in each subsample to
compute Eq. (58). To be specific, we set $z=0.324,0.567$ and $0.878$ for the
cluster sample in the redshift ranges $0.092<z\leq 0.445$, $0.445<z\leq 0.695$
and $0.695<z\leq 1.180$, respectively.

Fig. 10 summarizes the 68% credible interval of the mass bias333We compute the
ratio of $M_{\mathrm{500c}}/M_{\mathrm{UPP}}$ to be $10^{\langle\log
M_{\mathrm{500c}}|\log M_{\mathrm{UPP}}\rangle-\log M_{\mathrm{UPP}}}$ for our
constraints. inferred from the posterior distribution of our model parameters.
Our stacked lensing measurements are able to constrain the mass bias at the
average SZ mass as

$\displaystyle B=\left\\{\begin{array}[]{ll}1.3\pm 0.2&(0.092<z\leq 0.445)\\\
1.6\pm 0.2&(0.445<z\leq 0.695)\\\ 1.6\pm 0.3&(0.695<z\leq
1.180)\end{array}\right.,$ (62)

so we find marginal evidence of the mass bias ($B\neq 1$) at every redshift
bin in our analysis.

It is worth noting that our limits of Eq. (62) depend on the choice of prior
in our model parameters. As shown in Appendix F, we observe that the mass bias
parameter can correlate with not only $a_{0}$ (the baseline mean scaling
relation) but also $C_{0}$ (the scatter in the scaling relation) in our
posterior distribution. This indicates that any additional information of the
scatter can improve the limits of $B$. Apart from $C_{0}$, $a_{1}$ (describing
the mass dependence in the scaling relation) has a weak correlation with the
derived parameter $B$ in our inference. Future precise stacked-lensing
measurements would require to break the degeneracy between $B$ and $a_{1}$ by
further dividing cluster samples with $M_{\mathrm{UPP}}$ or another mass
proxy.

For comparison, we overlay cluster mass estimates of $M_{\mathrm{500c}}$ from
other weak lensing analyses in Fig. 10. Those include the stacked lensing
analysis of 9 ACT-selected clusters [81], the lensing analysis from the
Cluster Lensing And Supernova survey with Hubble (CLASH) for 21 galaxy
clusters selected by their Planck thermal SZ signal [82], the lensing mass
measurements of 5 Planck-selected clusters using HSC year-1 shear data [83],
the lensing analysis of 35 Planck-selected clusters employing galaxy shapes
from CFHTLenS and RCSLenS [84], the lensing mass measurements of 8 ACT-
selected clusters employing HSC year-1 shear data [22], the scaling relation
between the lensing and the Planck SZ masses derived by the Canadian Cluster
Comparison Project with the help of the Multi Epoch Nearby Cluster Survey
[85], the lensing mass measurements of 30 galaxy clusters from the South Pole
Telescope Sunyaev-Zel’dovich (SPT-SZ) Survey with galaxy shape measurements by
Hubble Space Telescope [86], and the stacked lensing analysis of 157 ACT-
selected clusters using KiDS-1000 observations [26]. In the figure, we compute
the average masses from individual lensing mass measurements in Refs. [82, 83,
84, 22, 86] for visualization purposes. We find that our inference of the mass
bias is consistent with the previous constraints at a $\sim 1\sigma$ level.

Another important constraint on the mass bias has been derived from the power
spectrum analysis of the Planck thermal SZ map [87]. The right panel in Fig.
10 highlights the constraint on the mass bias from this analysis, assuming our
fiducial $\Lambda$CDM cosmology. We caution that Ref. [87] assumes a
deterministic scaling relation between the SZ mass $M_{\mathrm{UPP}}$ and the
spherical overdensity mass $M_{\mathrm{500c}}$, suggesting that a direct
comparison with our results may be inappropriate.

###  V.3 Different choices of cluster forward modeling

We evaluate the impact of our modeling choices by implementing a simpler
forward model with several different assumptions. Here we assume the density
profile to be an NFW profile [27],

$\rho_{\mathrm{1h}}(r)=\frac{\delta_{c}\rho_{c}}{(r/r_{\mathrm{s}})\left(1+r/r_{\mathrm{s}}\right)^{2}}$
(63)

instead of Eq. (IV.2), using the expression for the NFW excess surface density
given by [88]. Here, $\rho_{c}$ is the normalization and $r_{\mathrm{s}}$ is
the scale radius, and

$\delta_{c}=\frac{500}{3}\frac{c_{500}^{3}}{\left[\ln(1+c_{500})-c_{500}/(1+c_{500})\right]}.$
(64)

with $c_{500}\equiv r_{500}/r_{s}$. We assume the mass-concentration relation
of Ref. [89] as implemented in colossus [90]. We rewrite Eq. (21) as

$\log
M^{\prime}_{\mathrm{UPP}}=a_{0}^{\prime}+a_{1}^{\prime}\log\left(\frac{M_{\mathrm{500c}}}{5\times
10^{14}\,M_{\odot}}\right),$ (65)

where for simplicity $a_{1}^{\prime}=1$ while $a_{0}^{\prime}$ is fit to each
of the redshift bins with a uniform prior,
$a_{0}^{\prime}\sim\mathcal{U}(13,16)$. We assume a log-normal scatter (Eq.
IV.1) independent of mass and redshift, $\sigma_{\log
M_{\mathrm{UPP}}}\sim\mathcal{U}(0.01,1)$, in each bin.

The 2-halo term is calculated explicitly as in Ref. [26], i.e., assuming the
cluster power spectrum is biased with respect to the dark matter power
spectrum,

$P_{\mathrm{2h}}(k|M,z)=\mathcal{A}_{\mathrm{2h}}b_{\mathrm{h}}(M)P_{\mathrm{mm}}(k|M,z)$
(66)

with $b_{\mathrm{h}}(M)$ the halo bias prescription of Ref. [91], and using
Eqs. (28)–(29) to calculate the lensing 2-halo term from the power spectrum.
The predicted signal for each redshift bin is calculated using Eq. (IV.1),
adopting the halo mass function of Ref. [3]. In this implementation we
calculate $b_{\mathrm{h}}$, $P_{\mathrm{mm}}$, and $n_{h}$ using the Core
Cosmology Library [92]. We adopt a flat prior
$0.5\leq\mathcal{A}_{\mathrm{2h}}\leq 2.0$ to give the model some freedom over
these choices.

Finally, we assume that a fraction $f_{\mathrm{off}}$ of clusters are off-
centered following a Rayleigh distribution [93] with a characteristic width
$\sigma_{\mathrm{mis}}$, independent of mass and redshift. Consistent with the
fiducial model, we set wide priors for these parameters, making them effective
parameters capturing any small-scale deviations from our model:
$f_{\mathrm{off}}\sim\mathcal{U}(0,1)$ and
$\sigma_{\mathrm{off}}\sim\mathcal{U}(0,1)\,h^{-1}$Mpc.

As in the fiducial models, we fit the model above independently to each of the
redshift bins.

Redshift bin | $a_{0}^{\prime}$ |  $\sigma_{\log M_{\mathrm{UPP}}}$† | $\mathcal{A}_{\mathrm{2h}}$
---|---|---|---
$0.092<z\leq 0.445$ | $13.82_{-0.57}^{+0.56}$ | $\leq 0.48$ | $1.0_{-0.3}^{+0.4}$
$0.445<z\leq 0.695$ | $13.86_{-0.59}^{+0.53}$ | $\leq 0.44$ | $1.1_{-0.4}^{+0.5}$
$0.695<z\leq 1.180$ | $13.86_{-0.60}^{+0.44}$ | $\leq 0.40$ | $1.0_{-0.3}^{+0.5}$
Table 1: Marginalized posterior constraints on the alternative model of Sec.
V.3. †95% upper limits.

We list the best-fit parameters in Table 1. Miscentering parameters are
unconstrained in all cases. The resulting mean mass bias at the average
$M_{\mathrm{UPP}}$ at each redshift bin are

$B^{\mathrm{alt}}=\begin{cases}1.3\pm 0.2&(0.092<z\leq 0.445)\\\ 1.3\pm
0.2&(0.445<z\leq 0.695)\\\ 1.6_{-0.4}^{+0.5}&(0.695<z\leq 1.180),\end{cases}$
(67)

where the minimum $\chi^{2}$ of each redshift bin are 14.2, 20.0, 12.9,
resulting in PTEs of 0.72, 0.33, 0.61. Eq. 67 should be directly compared with
our fiducial limits of Eq. 62. Hence, we find that the mass-bias limits by our
stacked lensing measurements can be affected by choices of cluster modeling
with a $1\sigma$ level at largest. A similar systematic effect has been
reported by our previous analysis [22].

##  VI Limitation

Before concluding, we here summarize the major limitations of our analysis.
The results in Section V depend strongly on our forward model of galaxy
clusters selected by the thermal SZ effect in CMB observations. We note that
our model relies on three assumptions; (i) that the completeness of the SZ-
selected clusters is determined by SZ masses and redshifts alone, (ii) that
the stacked lensing signals around galaxy clusters beyond their virial radii
are determined by total cluster masses, redshifts, and cosmology, and (iii)
that the lensing effect induced by a single cluster is independent of the SZ
effect of that cluster. The validity of these assumptions will be addressed in
future studies. In the following, we briefly discuss what can break the
assumptions and how to improve our halo-based forward modeling of the SZ-
selected clusters.

###  VI.1 Completeness of SZ-selected clusters

The completeness of the list of galaxy clusters in a given SZ observation
plays a critical role in developing a precise prediction of stacked lensing
measurements around the selected clusters. In this paper, we adopt the method
developed in Ref. [12] to estimate the completeness. This is a simulation-
based method assuming the thermal SZ scaling relation in Ref. [33] with a log-
normal scatter and random Gaussian noise in accord with the ACT observation.
Note that the electron-pressure model in Ref. [33] has been constructed with a
cluster sample at $z<0.3$; recent X-ray and thermal SZ observations, however,
have confirmed that it can provide a reasonable fit to the average pressure
profiles in cluster outskirts even for higher redshifts up to $z\sim 1$ [94,
95, 96, 97, 98, 99, 100, 101, 102, 103, 104].

In addition to its average, the scatter in the thermal SZ effect of individual
clusters can induce possible systematic effects in estimating the
completeness. For instance, Ref. [105] developed a semi-analytic model of
spherical-symmetric pressure profiles of intracluster media in the presence of
mass accretion, finding that the scatter in the thermal SZ-mass scaling
relation on cluster scales can be closely related to the amount of assembly-
driven, non-thermal pressure. Because the non-thermal pressure in cluster
outskirts becomes more significant as the clusters experience mass assembly at
a faster accretion rate [106], the thermal SZ observation can preferentially
select galaxy clusters with slower accretion rates (i.e. less non-thermal
pressure support) at a given halo mass. In the end, stacked lensing signals
around the SZ-selected clusters can thus depend on the mass accretion history
of the clusters [72, 107], which our halo-based model does not take into
account. Also, compact radio sources can cause addtional scatter in the
cluster scaling relations. Deep radio-wavelength observations would help to
quantify contaminations from the point sources [108, 109].

To examine the completeness of SZ-selected clusters in detail, it would be
useful to have samples of galaxy clusters selected by weak-lensing
observations of large-scale structures [110, 111, 112]. Lensing-selected
clusters are predominantly sensitive to the underlying mass distributions
around the clusters [113, 114], providing comprehensive information about the
selection function of galaxy clusters. As a matter of fact, the HSC survey is
aimed at constructing a large sample of lensing-selected clusters [115, 116].
Measurements of thermodynamic properties of lensing-selected clusters can give
insight into cluster astrophysics [117].

###  VI.2 Secondary bias effects on galaxy-cluster scales

Our inference of a representative halo mass of the SZ-selected clusters by the
stacked lensing measurements partly relies on a theoretical model of the mass
distribution around the clusters beyond their virial radii. In a halo-model
language (see Ref. [118] for a review of the halo model), the stacked lensing
signal at large radial distances can be explained by the two-halo term,
describing clustering effects of mass density between neighboring clusters
[119]. At the lowest order, the two-halo term is proportional to the linear
biasing factor and the correlation function of mass density fields, which are
merely a function of the cluster masses, redshifts, and cosmological models.

However, numerical simulations have shown that the linear biases of dark
matter halos can depend not only on their masses but also on other properties
defined in their inner density profiles at a given redshift [120, 121]. These
are commonly referred to as secondary bias effects, and can potentially lead
to a biased estimate of the cluster mass from stacked lensing analyses. Given
the correlation between the thermal SZ effect and the mass accretion history,
secondary bias on cluster scales can be important for future lensing analyses
with high precision. Note that a detailed characterization of the secondary
bias of cluster-sized halos requires a careful assessment of the mass assembly
history [122, 123, 107]. A large-volume numerical simulation with a fine-time
sampling is urgently needed to improve our understanding of the secondary bias
on cluster scales.

###  VI.3 Correlated scatter between stacked lensing signals and cluster
observables

Our halo-based model of the stacked lensing signals implicitly assumes that
the SZ mass of a single cluster does not correlate with the inner mass density
profile of the cluster. This is just an assumption and further investigations
into possible correlations between the thermal SZ and lensing signals may be
needed. In fact, numerical simulations have shown that the volume-integrated
Compton $Y$ parameter of a galaxy cluster can correlate with their cluster
lensing mass when one infers those from projected pressure and mass density
fields [124]. The correlation can arise from coherent structures of gas
pressure and mass density along a line of sight. The importance of such
possible correlations in stacked lensing analyses around SZ-selected clusters
is still uncertain, whereas recent studies of optically-selected galaxy
clusters have pointed out that the correlated scatter between the lensing
signal and the cluster richness can provide a biased estimate of the linear
halo bias [125] and the halo mass [126]. To examine the correlated scatter for
SZ-selected clusters, we require a careful analysis based on realistic mock
observations of CMB anisotropies as well as weak lensing shear.

##  VII Conclusion

In this paper, we performed a stacked lensing analysis of galaxy clusters
selected through thermal Sunyaev-Zel’dovich (SZ) effects in observations of
the cosmic microwave background by the Atacama Cosmology Telescope. For the
analysis, we used the galaxy shape catalog constructed from the Subaru Hyper-
Suprime Cam (HSC) Year 3 data. The large number density in the HSC galaxy
catalog allows us to study the stacked lensing signal of galaxy clusters as a
function of cluster redshift. We improved upon our previous analysis [22] by
increasing the number of galaxy clusters, using a precise estimation of
statistical errors in our measurements based on mock galaxy and cluster
catalogs, adopting a simulation-based prior on the scatter in an observable-
mass scaling relation, and accounting for possible baryonic effects on the
stacked lensing signals with a semi-analytic model as motivated by Refs. [65,
66, 67]. We worked with a halo-based forward model to infer statistical
relationships between cluster masses estimated by the thermal SZ analysis
(i.e. the SZ masses) and weak-lensing counterparts. Our model relies on an
accurate emulation tool, the Dark Emulator [25], for a baseline prediction of
statistics of cluster-sized dark matter halos, while modifying the emulator’s
output so that the model includes possible changes in the mass density
distribution due to the existence of stellar and gas components in the
clusters. Our findings in this study are summarized below.

  1. (i)

Imposing a selection of foreground clusters with $>5\sigma$ significance in
the thermal SZ analysis and background source galaxies with secure redshift
cuts, we found that the cumulative signal-to-noise ratio of our stacked
lensing measurements in the range of $0.1<R\,[h^{-1}\mathrm{Mpc}]<100$ reaches
14.6, 12.0 and 6.6 at the cluster redshifts of $0.092<z\leq 0.445$,
$0.445<z\leq 0.695$, and $0.695<z\leq 1.180$, respectively. Our lensing
analysis passed several null tests about B-mode signals and boost factors. We
also confirmed that the measurements are only minimally affected by the choice
of cluster centers and details in photometric redshift estimations.

  2. (ii)

We inferred the statistical relationship between the SZ and lensing masses
with the MCMC analysis including a loose prior about the cluster abundance.
Our best-fit model provided a reasonable chi-squared value with a p-value
larger than 0.01. We found that the abundance prior helps to break the
inherent parameter degeneracy between the mean and the scatter in the mass
scaling relation. Our stacked lensing measurements are able to constrain a
parameter describing the mean scaling relation at a level of 15% at different
cluster redshifts, while the scatter is not well constrained and is mostly
determined by the prior. We also observed that baryonic effects are not
important in our parameter inference at this stage.

  3. (iii)

Our analysis put a constraint of the mass bias parameter
$B=M_{\mathrm{500c}}/M_{\mathrm{UPP}}$, where $M_{\mathrm{500c}}$ is the
spherical overdensity mass with respect to 500 times critical density and
$M_{\mathrm{UPP}}$ is the counterpart estimated by the thermal SZ analysis.
For the average $M_{\mathrm{UPP}}$ of $2.1-2.4\times 10^{14}\,h^{-1}M_{\odot}$
in our cluster selection, we found that $B=1.3\pm 0.2,1.6\pm 0.2$ and $1.6\pm
0.3$ at the redshifts of $0.092<z\leq 0.445$, $0.445<z\leq 0.695$, and
$0.695<z\leq 1.180$, respectively (the errors are set at 68% credible
intervals). The constraints were found to be consistent with ones by previous
lensing analyses at a $1\sigma$ level, whereas the choice of the cluster
forward modeling can shift the mass bias by $\sim 1\sigma$ at largest.

Despite this careful analysis, we realize that the results in this paper
depend on details of our forward modeling approach. We discussed some
limitations of our model in Section VI, prompting further investigation of
possible biases in the thermal SZ-based selection of galaxy clusters and of
secondary bias effects of cluster-sized dark matter halos, as well as detailed
characterization of stacked lensing signals around the thermal SZ-selected
clusters with realistic mock catalogs. The HSC survey will provide their final
release of weak lensing shear data, increasing the survey area by a factor of
$\sim 3$ in a timely manner. A detailed mass calibration of the thermal SZ-
selected clusters by a stacked lensing analysis with the HSC final data can
deepen our understanding of the deficit in cluster masses inferred from
visible and dark elements, unlocking clues about structure formation of the
most massive objects in the universe.

###### Acknowledgements.

Support for ACT was through the U.S. National Science Foundation through
awards AST-0408698, AST-0965625, and AST-1440226 for the ACT project, as well
as awards PHY-0355328, PHY-0855887 and PHY-1214379. Funding was also provided
by Princeton University, the University of Pennsylvania, and a Canada
Foundation for Innovation (CFI) award to UBC. ACT operated in the Parque
Astronómico Atacama in northern Chile under the auspices of the Agencia
Nacional de Investigación y Desarrollo (ANID). The development of multichroic
detectors and lenses was supported by NASA grants NNX13AE56G and NNX14AB58G.
Detector research at NIST was supported by the NIST Innovations in Measurement
Science program. Computing for ACT was performed using the Princeton Research
Computing resources at Princeton University, the National Energy Research
Scientific Computing Center (NERSC), and the Niagara supercomputer 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. We thank the Republic of Chile for
hosting ACT in the northern Atacama, and the local indigenous Licanantay
communities whom we follow in observing and learning from the night sky. The
Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of
Japan and Taiwan, and Princeton University. The HSC instrumentation and
software were developed by the National Astronomical Observatory of Japan
(NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe
(Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research
Organization (KEK), the Academia Sinica Institute for Astronomy and
Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was
contributed by the FIRST program from Japanese Cabinet Office, the Ministry of
Education, Culture, Sports, Science and Technology (MEXT), the Japan Society
for the Promotion of Science (JSPS), Japan Science and Technology Agency
(JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and
Princeton University. This paper makes use of software developed for the Large
Synoptic Survey Telescope. We thank the LSST Project for making their code
available as free software at http://dm.lsst.org. The Pan-STARRS1 Surveys
(PS1) have been made possible through contributions of the Institute for
Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-
Planck Society, and its participating institutes, the Max Planck Institute for
Astronomy, Heidelberg, and the Max Planck Institute for Extraterrestrial
Physics, Garching, The Johns Hopkins University, Durham University, the
University of Edinburgh, Queen’s University Belfast, the Harvard-Smithsonian
Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network
Incorporated, the National Central University of Taiwan, the Space Telescope
Science Institute, the National Aeronautics and Space Administration under
Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA
Science Mission Directorate, the National Science Foundation under Grant No.
AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE) and
the Los Alamos National Laboratory. Based on data collected at the Subaru
Telescope and retrieved from the HSC data archive system, which is operated by
Subaru Telescope and Astronomy Data Center, National Astronomical Observatory
of Japan. This work is supported by MEXT KAKENHI Grant Numbers, 19K14767,
20H05861, 20H01932, 20H05855, 21H05456, 22K21349, 23H00108 and 23K19070. This
work is also supported by Tokai Pathways to Global Excellence (T-GEx), part of
MEXT Strategic Professional Development Program for Young Researchers. CS
acknowledges support from the Agencia Nacional de Investigación y Desarrollo
(ANID) through Basal project FB210003. DN is supported by NSF (AST-2206055 &
2307280) and NASA (80NSSC22K0821 & TM3-24007X) grants. MH acknowledges support
from the National Research Foundation of South Africa (grant no. 137975). A
part of numerical computations was carried out on Cray XC50 at the Center for
Computational Astrophysics in NAOJ.

## References

  * Press and Schechter [1974] W. H. Press and P. Schechter, Formation of galaxies and clusters of galaxies by selfsimilar gravitational condensation, Astrophys. J. 187, 425 (1974).
  * Sheth [1998] R. K. Sheth, An Excursion set model for the distribution of dark matter and dark matter halos, Mon. Not. Roy. Astron. Soc. 300, 1057 (1998), arXiv:astro-ph/9805319 .
  * Tinker  _et al._ [2008] J. L. Tinker, A. V. Kravtsov, A. Klypin, K. Abazajian, M. S. Warren, G. Yepes, S. Gottlober, and D. E. Holz, Toward a halo mass function for precision cosmology: The Limits of universality, Astrophys. J. 688, 709 (2008), arXiv:0803.2706 [astro-ph] .
  * Bhattacharya  _et al._ [2011] S. Bhattacharya, K. Heitmann, M. White, Z. Lukić, C. Wagner, and S. Habib, Mass Function Predictions Beyond $\Lambda$CDM, Astrophys. J. 732, 122 (2011), arXiv:1005.2239 [astro-ph.CO] .
  * Watson  _et al._ [2013] W. A. Watson, I. T. Iliev, A. D’Aloisio, A. Knebe, P. R. Shapiro, and G. Yepes, The halo mass function through the cosmic ages, Mon. Not. Roy. Astron. Soc. 433, 1230 (2013), arXiv:1212.0095 [astro-ph.CO] .
  * Despali  _et al._ [2016] G. Despali, C. Giocoli, R. E. Angulo, G. Tormen, R. K. Sheth, G. Baso, and L. Moscardini, The universality of the virial halo mass function and models for non-universality of other halo definitions, Mon. Not. Roy. Astron. Soc. 456, 2486 (2016), arXiv:1507.05627 [astro-ph.CO] .
  * Shirasaki  _et al._ [2021] M. Shirasaki, T. Ishiyama, and S. Ando, Virial Halo Mass Function in the Planck Cosmology, Astrophys. J. 922, 89 (2021), arXiv:2108.11038 [astro-ph.CO] .
  * Allen  _et al._ [2011] S. W. Allen, A. E. Evrard, and A. B. Mantz, Cosmological Parameters from Observations of Galaxy Clusters, Ann. Rev. Astron. Astrophys. 49, 409 (2011), arXiv:1103.4829 [astro-ph.CO] .
  * Sunyaev and Zeldovich [1972] R. A. Sunyaev and Y. B. Zeldovich, The Observations of relic radiation as a test of the nature of X-Ray radiation from the clusters of galaxies, Comments Astrophys. Space Phys. 4, 173 (1972).
  * Ade  _et al._ [2016a] P. A. R. Ade _et al._ (Planck), Planck 2015 results. XXVII. The Second Planck Catalogue of Sunyaev-Zeldovich Sources, Astron. Astrophys. 594, A27 (2016a), arXiv:1502.01598 [astro-ph.CO] .
  * Bocquet  _et al._ [2019] S. Bocquet _et al._ (SPT), Cluster Cosmology Constraints from the 2500 deg2 SPT-SZ Survey: Inclusion of Weak Gravitational Lensing Data from Magellan and the Hubble Space Telescope, Astrophys. J. 878, 55 (2019), arXiv:1812.01679 [astro-ph.CO] .
  * Hilton  _et al._ [2021] M. Hilton _et al._ (ACT, DES), The Atacama Cosmology Telescope: A Catalog of $>$4000 Sunyaev–Zel’dovich Galaxy Clusters, Astrophys. J. Suppl. 253, 3 (2021), arXiv:2009.11043 [astro-ph.CO] .
  * Madhavacheril  _et al._ [2017] M. S. Madhavacheril, N. Battaglia, and H. Miyatake, Fundamental physics from future weak-lensing calibrated Sunyaev-Zel’dovich galaxy cluster counts, Phys. Rev. D 96, 103525 (2017), arXiv:1708.07502 [astro-ph.CO] .
  * Hasselfield  _et al._ [2013] M. Hasselfield _et al._ , The Atacama Cosmology Telescope: Sunyaev-Zel’dovich selected galaxyclusters at 148 GHz from three seasons of data, JCAP 07, 008, arXiv:1301.0816 [astro-ph.CO] .
  * Ade  _et al._ [2016b] P. A. R. Ade _et al._ (Planck), Planck 2015 results. XXIV. Cosmology from Sunyaev-Zeldovich cluster counts, Astron. Astrophys. 594, A24 (2016b), arXiv:1502.01597 [astro-ph.CO] .
  * Salvati  _et al._ [2022] L. Salvati _et al._ , Combining Planck and SPT Cluster Catalogs: Cosmological Analysis and Impact on the Planck Scaling Relation Calibration, Astrophys. J. 934, 129 (2022), arXiv:2112.03606 [astro-ph.CO] .
  * Rosati  _et al._ [2002] P. Rosati, S. Borgani, and C. Norman, The evolution of x-ray clusters of galaxies, Ann. Rev. Astron. Astrophys. 40, 539 (2002), arXiv:astro-ph/0209035 .
  * Pratt  _et al._ [2019] G. W. Pratt, M. Arnaud, A. Biviano, D. Eckert, S. Ettori, D. Nagai, N. Okabe, and T. H. Reiprich, The galaxy cluster mass scale and its impact on cosmological constraints from the cluster population, Space Sci. Rev. 215, 25 (2019), arXiv:1902.10837 [astro-ph.CO] .
  * Umetsu [2020] K. Umetsu, Cluster–galaxy weak lensing, Astron. Astrophys. Rev. 28, 7 (2020), arXiv:2007.00506 [astro-ph.CO] .
  * Okabe  _et al._ [2010] N. Okabe, M. Takada, K. Umetsu, T. Futamase, and G. P. Smith, LoCuSS: Subaru Weak Lensing Study of 30 Galaxy Clusters, Publ. Astron. Soc. Jap. 62, 811 (2010), arXiv:0903.1103 [astro-ph.CO] .
  * Li  _et al._ [2022] X. Li _et al._ , The three-year shear catalog of the Subaru Hyper Suprime-Cam SSP Survey, Publ. Astron. Soc. Jap. 74, 421 (2022), arXiv:2107.00136 [astro-ph.CO] .
  * Miyatake  _et al._ [2019] H. Miyatake _et al._ , Weak-lensing Mass Calibration of ACTPol Sunyaev–Zel’dovich Clusters with the Hyper Suprime-Cam Survey, Astrophys. J. 875, 63 (2019), arXiv:1804.05873 [astro-ph.CO] .
  * Takahashi  _et al._ [2017] R. Takahashi, T. Hamana, M. Shirasaki, T. Namikawa, T. Nishimichi, K. Osato, and K. Shiroyama, Full-sky Gravitational Lensing Simulation for Large-area Galaxy Surveys and Cosmic Microwave Background Experiments, Astrophys. J. 850, 24 (2017), arXiv:1706.01472 [astro-ph.CO] .
  * Springel  _et al._ [2018] V. Springel _et al._ , First results from the IllustrisTNG simulations: matter and galaxy clustering, Mon. Not. Roy. Astron. Soc. 475, 676 (2018), arXiv:1707.03397 [astro-ph.GA] .
  * Nishimichi  _et al._ [2019] T. Nishimichi _et al._ , Dark Quest. I. Fast and Accurate Emulation of Halo Clustering Statistics and Its Application to Galaxy Clustering, Astrophys. J. 884, 29 (2019), arXiv:1811.09504 [astro-ph.CO] .
  * Robertson  _et al._ [2024] N. C. Robertson _et al._ , ACT-DR5 Sunyaev-Zel’dovich clusters: Weak lensing mass calibration with KiDS, Astron. Astrophys. 681, A87 (2024), arXiv:2304.10219 [astro-ph.CO] .
  * Navarro  _et al._ [1997] J. F. Navarro, C. S. Frenk, and S. D. M. White, A Universal density profile from hierarchical clustering, Astrophys. J. 490, 493 (1997), arXiv:astro-ph/9611107 .
  * Hinshaw  _et al._ [2013] G. Hinshaw _et al._ (WMAP), Nine-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Parameter Results, Astrophys. J. Suppl. 208, 19 (2013), arXiv:1212.5226 [astro-ph.CO] .
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Naess  _et al._ [2020] S. Naess _et al._ , The Atacama Cosmology Telescope: arcminute-resolution maps of 18 000 square degrees of the microwave sky from ACT 2008–2018 data combined with Planck, JCAP 12, 046, arXiv:2007.07290 [astro-ph.IM] .
  * Melin  _et al._ [2006] J.-B. Melin, J. G. Bartlett, and J. Delabrouille, Catalog extraction in sz cluster surveys: a matched filter approach, Astron. Astrophys. 459, 341 (2006), arXiv:astro-ph/0602424 .
  * Williamson  _et al._ [2011] R. Williamson _et al._ , An SZ-selected sample of the most massive galaxy clusters in the 2500-square-degree South Pole Telescope survey, Astrophys. J. 738, 139 (2011), arXiv:1101.1290 [astro-ph.CO] .
  * Arnaud  _et al._ [2010] M. Arnaud, G. W. Pratt, R. Piffaretti, H. Boehringer, J. H. Croston, and E. Pointecouteau, The universal galaxy cluster pressure profile from a representative sample of nearby systems (REXCESS) and the Y_SZ-M_500 relation, Astron. Astrophys. 517, A92 (2010), arXiv:0910.1234 [astro-ph.CO] .
  * Blake  _et al._ [2016] C. Blake _et al._ , The 2-degree Field Lensing Survey: design and clustering measurements, Mon. Not. Roy. Astron. Soc. 462, 4240 (2016), arXiv:1608.02668 [astro-ph.CO] .
  * Childress  _et al._ [2017] M. J. Childress _et al._ (DES), OzDES multifibre spectroscopy for the Dark Energy Survey: 3-yr results and first data release, Mon. Not. Roy. Astron. Soc. 472, 273 (2017), arXiv:1708.04526 [astro-ph.CO] .
  * Ahumada  _et al._ [2020] R. Ahumada _et al._ (SDSS-IV), The 16th Data Release of the Sloan Digital Sky Surveys: First Release from the APOGEE-2 Southern Survey and Full Release of eBOSS Spectra, Astrophys. J. Suppl. 249, 3 (2020), arXiv:1912.02905 [astro-ph.GA] .
  * Scodeggio  _et al._ [2018] M. Scodeggio _et al._ , The VIMOS Public Extragalactic Redshift Survey (VIPERS). Full spectroscopic data and auxiliary information release (PDR-2), Astron. Astrophys. 609, A84 (2018), arXiv:1611.07048 [astro-ph.GA] .
  * Rykoff  _et al._ [2014] E. S. Rykoff _et al._ (SDSS), redMaPPer I: Algorithm and SDSS DR8 Catalog, Astrophys. J. 785, 104 (2014), arXiv:1303.3562 [astro-ph.CO] .
  * Rykoff  _et al._ [2016] E. S. Rykoff _et al._ (DES), The redMaPPer Galaxy Cluster Catalog From DES Science Verification Data, Astrophys. J. Suppl. 224, 1 (2016), arXiv:1601.00621 [astro-ph.CO] .
  * Oguri  _et al._ [2018] M. Oguri _et al._ , An optically-selected cluster catalog at redshift 0.1$<$ z$<$ 1.1 from the Hyper Suprime-Cam Subaru Strategic Program S16A data, Publ. Astron. Soc. Jap. 70, S20 (2018), arXiv:1701.00818 [astro-ph.CO] .
  * Hilton  _et al._ [2018] M. Hilton _et al._ (ACT), The Atacama Cosmology Telescope: The Two-Season ACTPol Sunyaev-Zel’dovich Effect Selected Cluster Catalog, Astrophys. J. Suppl. 235, 20 (2018), arXiv:1709.05600 [astro-ph.CO] .
  * Furusawa  _et al._ [2018] H. Furusawa, M. Koike, T. Takata, Y. Okura, H. Miyatake, R. H. Lupton, S. Bickerton, P. A. Price, J. Bosch, N. Yasuda, S. Mineo, Y. Yamada, S. Miyazaki, F. Nakata, S. Koshida, Y. Komiyama, Y. Utsumi, S. Kawanomoto, E. Jeschke, J. Noumaru, K. Schubert, I. Iwata, F. Finet, T. Fujiyoshi, A. Tajitsu, T. Terai, and C.-H. Lee, The on-site quality-assurance system for Hyper Suprime-Cam: OSQAH, Publ. Astron. Soc. Jap. 70, S3 (2018).
  * Miyazaki  _et al._ [2018] S. Miyazaki, Y. Komiyama, S. Kawanomoto, Y. Doi, H. Furusawa, T. Hamana, Y. Hayashi, H. Ikeda, Y. Kamata, H. Karoji, M. Koike, T. Kurakami, S. Miyama, T. Morokuma, F. Nakata, K. Namikawa, H. Nakaya, K. Nariai, Y. Obuchi, Y. Oishi, N. Okada, Y. Okura, P. Tait, T. Takata, Y. Tanaka, M. Tanaka, T. Terai, D. Tomono, F. Uraguchi, T. Usuda, Y. Utsumi, Y. Yamada, H. Yamanoi, H. Aihara, H. Fujimori, S. Mineo, H. Miyatake, M. Oguri, T. Uchida, M. M. Tanaka, N. Yasuda, M. Takada, H. Murayama, A. J. Nishizawa, N. Sugiyama, M. Chiba, T. Futamase, S.-Y. Wang, H.-Y. Chen, P. T. P. Ho, E. J. Y. Liaw, C.-F. Chiu, C.-L. Ho, T.-C. Lai, Y.-C. Lee, D.-Z. Jeng, S. Iwamura, R. Armstrong, S. Bickerton, J. Bosch, J. E. Gunn, R. H. Lupton, C. Loomis, P. Price, S. Smith, M. A. Strauss, E. L. Turner, H. Suzuki, Y. Miyazaki, M. Muramatsu, K. Yamamoto, M. Endo, Y. Ezaki, N. Ito, N. Kawaguchi, S. Sofuku, T. Taniike, K. Akutsu, N. Dojo, K. Kasumi, T. Matsuda, K. Imoto, Y. Miwa, M. Suzuki, K. Takeshi, and H. Yokota, Hyper Suprime-Cam: System design and verification of image quality, Publ. Astron. Soc. Jap. 70, S1 (2018).
  * Aihara  _et al._ [2018] H. Aihara _et al._ , The Hyper Suprime-Cam SSP Survey: Overview and Survey Design, Publ. Astron. Soc. Jap. 70, S4 (2018), arXiv:1704.05858 [astro-ph.IM] .
  * Hirata and Seljak [2003] C. M. Hirata and U. Seljak, Shear calibration biases in weak lensing surveys, Mon. Not. Roy. Astron. Soc. 343, 459 (2003), arXiv:astro-ph/0301054 .
  * Mandelbaum  _et al._ [2005] R. Mandelbaum, C. M. Hirata, U. Seljak, J. Guzik, N. Padmanabhan, C. Blake, M. R. Blanton, R. Lupton, and J. Brinkmann, Systematic errors in weak lensing: Application to SDSS galaxy-galaxy weak lensing, Mon. Not. Roy. Astron. Soc. 361, 1287 (2005), arXiv:astro-ph/0501201 .
  * Mandelbaum  _et al._ [2013] R. Mandelbaum, A. Slosar, T. Baldauf, U. Seljak, C. M. Hirata, R. Nakajima, R. Reyes, and R. E. Smith, Cosmological parameter constraints from galaxy-galaxy lensing and galaxy clustering with the SDSS DR7, Mon. Not. Roy. Astron. Soc. 432, 1544 (2013), arXiv:1207.1120 [astro-ph.CO] .
  * Miralda-Escude [1991] J. Miralda-Escude, The Correlation Function of Galaxy Ellipticities Produced by Gravitational Lensing, Astrophys. J. 380, 1 (1991).
  * Bernstein and Jarvis [2002] G. M. Bernstein and M. Jarvis, Shapes and shears, stars and smears: optimal measurements for weak lensing, Astron. J. 123, 583 (2002), arXiv:astro-ph/0107431 .
  * Rowe  _et al._ [2015] B. T. P. Rowe, M. Jarvis, R. Mandelbaum, G. M. Bernstein, J. Bosch, M. Simet, J. E. Meyers, T. Kacprzak, R. Nakajima, J. Zuntz, H. Miyatake, J. P. Dietrich, R. Armstrong, P. Melchior, and M. S. S. Gill, GALSIM: The modular galaxy image simulation toolkit, Astronomy and Computing 10, 121 (2015), arXiv:1407.7676 [astro-ph.IM] .
  * Mandelbaum  _et al._ [2018] R. Mandelbaum _et al._ , Weak lensing shear calibration with simulations of the HSC survey, Mon. Not. Roy. Astron. Soc. 481, 3170 (2018), arXiv:1710.00885 [astro-ph.CO] .
  * Tanaka  _et al._ [2018] M. Tanaka, J. Coupon, B.-C. Hsieh, S. Mineo, A. J. Nishizawa, J. Speagle, H. Furusawa, S. Miyazaki, and H. Murayama, Photometric Redshifts for Hyper Suprime-Cam Subaru Strategic Program Data Release 1, Publ. Astron. Soc. Jap. 70, S9 (2018), arXiv:1704.05988 [astro-ph.GA] .
  * Oguri [2014] M. Oguri, A cluster finding algorithm based on the multiband identification of red sequence galaxies, Mon. Not. Roy. Astron. Soc. 444, 147 (2014), arXiv:1407.4693 [astro-ph.CO] .
  * Mandelbaum  _et al._ [2006] R. Mandelbaum, U. Seljak, R. J. Cool, M. Blanton, C. M. Hirata, and J. Brinkmann, Density profiles of galaxy groups and clusters from SDSS galaxy-galaxy weak lensing, Mon. Not. Roy. Astron. Soc. 372, 758 (2006), arXiv:astro-ph/0605476 .
  * Miyatake  _et al._ [2022a] H. Miyatake _et al._ , Cosmological inference from an emulator based halo model. II. Joint analysis of galaxy-galaxy weak lensing and galaxy clustering from HSC-Y1 and SDSS, Phys. Rev. D 106, 083520 (2022a), arXiv:2111.02419 [astro-ph.CO] .
  * Sheldon  _et al._ [2004] E. S. Sheldon _et al._ (SDSS), The Galaxy - mass correlation function measured from weak lensing in the SDSS, Astron. J. 127, 2544 (2004), arXiv:astro-ph/0312036 .
  * Shirasaki and Takada [2018] M. Shirasaki and M. Takada, Stacked lensing estimators and their covariance matrices: excess surface mass density versus lensing shear, Mon. Not. Roy. Astron. Soc. 478, 4277 (2018), arXiv:1802.09696 [astro-ph.CO] .
  * Tanaka [2015] M. Tanaka, Photometric Redshift with Bayesian Priors on Physical Properties of Galaxies, Astrophys. J. 801, 20 (2015), arXiv:1501.02047 [astro-ph.GA] .
  * Medezinski  _et al._ [2018] E. Medezinski _et al._ , Source Selection for Cluster Weak Lensing Measurements in the Hyper Suprime-Cam Survey, Publ. Astron. Soc. Jap. 70, 30 (2018), arXiv:1706.00427 [astro-ph.CO] .
  * Chang  _et al._ [2013] C. Chang, M. Jarvis, B. Jain, S. M. Kahn, D. Kirkby, A. Connolly, S. Krughoff, E. Peng, and J. R. Peterson, The Effective Number Density of Galaxies for Weak Lensing Measurements in the LSST Project, Mon. Not. Roy. Astron. Soc. 434, 2121 (2013), arXiv:1305.0793 [astro-ph.CO] .
  * Murata  _et al._ [2018] R. Murata, T. Nishimichi, M. Takada, H. Miyatake, M. Shirasaki, S. More, R. Takahashi, and K. Osato, Constraints on the mass-richness relation from the abundance and weak lensing of SDSS clusters, Astrophys. J. 854, 120 (2018), arXiv:1707.01907 [astro-ph.CO] .
  * Miyatake  _et al._ [2022b] H. Miyatake, Y. Kobayashi, M. Takada, T. Nishimichi, M. Shirasaki, S. Sugiyama, R. Takahashi, K. Osato, S. More, and Y. Park, Cosmological inference from an emulator based halo model. I. Validation tests with HSC and SDSS mock catalogs, Phys. Rev. D 106, 083519 (2022b), arXiv:2101.00113 [astro-ph.CO] .
  * Takada and Hu [2013] M. Takada and W. Hu, Power Spectrum Super-Sample Covariance, Phys. Rev. D 87, 123504 (2013), arXiv:1302.6994 [astro-ph.CO] .
  * Takahashi  _et al._ [2019] R. Takahashi, T. Nishimichi, M. Takada, M. Shirasaki, and K. Shiroyama, Covariances for cosmic shear and galaxy–galaxy lensing in the response approach, Mon. Not. Roy. Astron. Soc. 482, 4253 (2019), arXiv:1805.11629 [astro-ph.CO] .
  * Schneider and Teyssier [2015] A. Schneider and R. Teyssier, A new method to quantify the effects of baryons on the matter power spectrum, JCAP 12, 049, arXiv:1510.06034 [astro-ph.CO] .
  * Schneider  _et al._ [2019] A. Schneider, R. Teyssier, J. Stadel, N. E. Chisari, A. M. C. Le Brun, A. Amara, and A. Refregier, Quantifying baryon effects on the matter power spectrum and the weak lensing shear correlation, JCAP 03, 020, arXiv:1810.08629 [astro-ph.CO] .
  * Giri and Schneider [2021] S. K. Giri and A. Schneider, Emulation of baryonic effects on the matter power spectrum and constraints from galaxy cluster data, JCAP 12 (12), 046, arXiv:2108.08863 [astro-ph.CO] .
  * Wu  _et al._ [2021] H.-Y. Wu, D. H. Weinberg, A. N. Salcedo, and B. D. Wibking, Cosmology with galaxy cluster weak lensing: statistical limits and experimental design, Astrophys. J. 910, 28 (2021), arXiv:2012.01956 [astro-ph.CO] .
  * Farahi  _et al._ [2018] A. Farahi, A. E. Evrard, I. McCarthy, D. J. Barnes, and S. T. Kay, Localized massive halo properties in BAHAMAS and MACSIS simulations: scalings, lognormality, and covariance, Mon. Not. Roy. Astron. Soc. 478, 2618 (2018), arXiv:1711.04922 [astro-ph.CO] .
  * Anbajagane  _et al._ [2020] D. Anbajagane, A. E. Evrard, A. Farahi, D. J. Barnes, K. Dolag, I. G. McCarthy, D. Nelson, and A. Pillepich, Stellar property statistics of massive haloes from cosmological hydrodynamics simulations: common kernel shapes, Mon. Not. Roy. Astron. Soc. 495, 686 (2020), arXiv:2001.02283 [astro-ph.GA] .
  * Oguri and Takada [2011] M. Oguri and M. Takada, Combining cluster observables and stacked weak lensing to probe dark energy: Self-calibration of systematic uncertainties, Phys. Rev. D 83, 023008 (2011), arXiv:1010.0744 [astro-ph.CO] .
  * Diemer and Kravtsov [2014] B. Diemer and A. V. Kravtsov, Dependence of the outer density profiles of halos on their mass accretion rate, Astrophys. J. 789, 1 (2014), arXiv:1401.1216 [astro-ph.CO] .
  * Baltz  _et al._ [2009] E. A. Baltz, P. Marshall, and M. Oguri, Analytic models of plausible gravitational lens potentials, JCAP 01, 015, arXiv:0705.0682 [astro-ph] .