# Weak-Lensing Shear-Selected Galaxy Clusters from the Hyper Suprime-Cam
Subaru Strategic Program:
I. Cluster Catalog, Selection Function and Mass–Observable Relation

Kai-Feng Chen1,2,⋆ 0000-0002-3839-0230 I-Non Chiu3 0000-0002-5819-6566
Masamune Oguri4,5 0000-0003-3484-399X Yen-Ting Lin6 0000-0001-7146-4687
Hironao Miyatake7,8,9 0000-0001-7964-9766 Satoshi Miyazaki10
0000-0002-1962-904X Surhud More11,9 0000-0002-2986-2371 Takashi Hamana12
Markus M. Rau13,14 0000-0003-3709-1324 Tomomi Sunayama15,7 0009-0004-6387-5784
Sunao Sugiyama16,9 0000-0003-1153-6735 Masahiro Takada9 0000-0002-5578-6472
1MIT Kavli Institute, Massachusetts Institute of Technology, Cambridge, MA
02139, USA  2Department of Physics, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA  3Department of Physics, National Cheng Kung
University, 70101 Tainan, Taiwan  4Center for Frontier Science, Chiba
University, 1-33 Yayoi-cho, Inage-ku, Chiba 263-8522, Japan  5Department of
Physics, Graduate School of Science, Chiba University, 1-33 Yayoi-Cho, Inage-
Ku, Chiba 263-8522, Japan  6Institute of Astronomy and Astrophysics Academia
Sinica, Taipei 106216, Taiwan  7Kobayashi-Maskawa Institute for the Origin of
Particles and the Universe (KMI), Nagoya University, Nagoya, 464-8602, Japan
8Institute for Advanced Research, Nagoya University, Nagoya 464-8601, Japan
9Kavli 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  10Subaru Telescope, National Astronomical
Observatory of Japan, 650 North Aohoku Place Hilo, HI 96720, USA  11The Inter-
University Centre for Astronomy and Astrophysics, Post bag 4, Ganeshkhind,
Pune 411007, India  12National Astronomical Observatory of Japan, National
Institutes of Natural Sciences, Mitaka, Tokyo 181-8588, Japan  13High Energy
Physics Division, Argonne National Laboratory, Lemont, IL 60439, USA
14McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon
University, 5000 Forbes Avenue, Pittsburgh, PA 15213, USA  15Department of
Astronomy and Steward Observatory, University of Arizona, 933 North Cherry
Avenue, Tucson, AZ 85719, USA  16Center for Particle Cosmology, Department of
Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA

###### Abstract

We present the first step towards deriving cosmological constraints through
the abundances of galaxy clusters selected in a $510\,\mathrm{deg}^{2}$ weak-
lensing aperture mass map, constructed with the Year-Three shear catalog from
the Hyper Suprime-Cam Subaru Strategic Program. We adopt a conservative source
galaxy selection to construct a sample of $129$ weak-lensing peaks with a
signal-to-noise ratio above $4.7$. We use semi-analytical injection
simulations to derive the selection function and the mass–observable relation
of our sample. These results take into account complicated uncertainties
associated with weak-lensing measurements, such as the non-uniform survey
depth and the complex survey geometry, projection effects from uncorrelated
large-scale structures, and the intrinsic alignment of source galaxies. We
also propose a novel modeling framework to make parts of the mass–observable
relation insensitive to assumed cosmological parameters. Such a framework not
only offers a great computational advantage to cosmological studies, but can
also benefit future astrophysical studies using shear-selected clusters. Our
results are an important step towards utilizing these cluster samples that are
constructed nearly independent of any baryonic assumptions in upcoming deep-
and-wide lensing surveys from the Vera Rubin Observatory, Euclid, and the
Nancy Grace Roman Space Telescope.

††thanks<EMAIL_ADDRESS>
##  1 Introduction

Observations of large-scale structures (LSS) probe the cosmic expansion
history and provide stringent constraints on the matter-energy content of the
universe. Over the last few decades, the standard $\Lambda$ Cold Dark Matter
($\Lambda$CDM) paradigm has achieved great phenomenological success in
describing these observations (e.g., Frieman et al., 2008; Weinberg et al.,
2013). However, the increasingly precise cosmological measurements have raised
the question of whether observations from different tracers across a wide
range of cosmic time can be consistently explained by the same set of
parameters in the $\Lambda$CDM model (see Abdalla et al., 2022, for a recent
summary). While these discrepancies among observations can be an exciting
indicator for new physics beyond $\Lambda$CDM, they are equally likely to be
caused by unaccounted-for systematic effects or inadequate uncertainty
modeling. Accurately characterizing observational data with a comprehensive
model of known physical and systematic effects is therefore an important task
for cosmologists in this decade.

An example of these tensions is in the measurements of
$S_{8}\coloneqq\sigma_{8}\sqrt{\Omega_{\mathrm{m}}/0.3}$, the present-day
amplitude of the matter fluctuations $\sigma_{8}$ scaled by the square root of
the total matter density $\Omega_{\mathrm{m}}$. Early-universe observations
from the cosmic microwave background (CMB, Planck Collaboration et al., 2020;
Aiola et al., 2020) suggest a higher $S_{8}$ than those inferred from the late
universe with cosmic shear (e.g. Asgari et al., 2021; Amon et al., 2022; Secco
et al., 2022; Li et al., 2023; Dalal et al., 2023), galaxy clustering and
galaxy-galaxy lensing ([e.g.][]DESY3_redMaGiC2x2:2022; Porredon et al., 2022),
cluster abundance (e.g. Planck Collaboration et al., 2016; Bocquet et al.,
2019; Chiu et al., 2023), and a combination of these probes (e.g. To et al.,
2021; Heymans et al., 2021; Abbott et al., 2022; More et al., 2023; Miyatake
et al., 2023).

In particular, observations utilizing the effect of weak gravitational lensing
(hereinafter weak lensing; see Mandelbaum, 2018, for a comprehensive review)
offer some of the strongest constraints on $S_{8}$ from the late universe.
These observations measure small but coherent shape distortions, often
referred to as shears, of distant galaxies (source) by intervening foreground
structures (lens), giving us a direct probe of matter density fluctuations.
The first cosmological application of weak lensing comes from the detection of
cosmic shear variance (Bacon et al., 2000; Kaiser et al., 2000; Van Waerbeke
et al., 2000; Wittman et al., 2000), and subsequent cosmological constraints
from cosmic shear power spectra and correlation functions have shown the power
of weak lensing with just the Gaussian summary statistics (Kilbinger, 2015).
Meanwhile, it is of increasing interest to extract additional cosmological
information from higher-order statistics in weak-lensing observations.
Examples of these summary statistics include number counts of peaks (Dietrich
& Hartlap, 2010; Hamana et al., 2015; Liu et al., 2015b, a; Kacprzak et al.,
2016; Shan et al., 2018; Martinet et al., 2018; Harnois-Déraps et al., 2021;
Zürcher et al., 2022; Liu et al., 2023; Marques et al., 2023) or higher
moments from weak-lensing mass maps (Van Waerbeke et al., 2013; Petri et al.,
2015; Peel et al., 2018; Chang et al., 2018; Gatti et al., 2020, 2022;
Anbajagane et al., 2023), three-point correlation functions or bispectra
(Takada & Jain, 2003, 2004; Dodelson & Zhang, 2005; Vafaei et al., 2010; Bergé
et al., 2010; Semboloni et al., 2011; Fu et al., 2014), Minkowski functionals
(Kratochvil et al., 2012; Petri et al., 2013; Vicinanza et al., 2019; Parroni
et al., 2020), density split statistics (Gruen et al., 2018; Friedrich et al.,
2018), and direct field-level inference (Jeffrey et al., 2020; Fluri et al.,
2022). Obtaining accurate and competitive cosmological constraints from these
higher-order statistics will serve as an important consistency test among
weak-lensing probes and help break parameter degeneracies.

In a series of two papers (this work and Chiu et al., 2024), we focus on
constraining cosmology with high signal-to-noise ratio peaks detected on the
weak-lensing aperture mass maps (Schneider, 1996; Oguri et al., 2021). These
peaks are associated with massive dark matter halos where clusters of galaxies
reside. The number counts of galaxy clusters, modeled by the halo mass
function, are sensitive to both the geometry and the structure formation
history in our universe (e.g., Allen et al., 2011). Tight cosmological
constraints have been obtained with cluster samples selected with X-ray
signals (Vikhlinin et al., 2009; Mantz et al., 2010, 2014; Chiu et al., 2023;
Ghirardini et al., 2024) from the thermal bremsstrahlung emission, at
millimeter wavelengths (Planck Collaboration et al., 2016; de Haan et al.,
2016; Bocquet et al., 2019) by the thermal Sunyaev–Zel’dovich effect (Sunyaev
& Zeldovich, 1972), and in the optical (Rozo et al., 2010; Abbott et al.,
2020; Costanzi et al., 2021; Lesci et al., 2022a, b; Sunayama et al., 2023)
via overdensities of galaxies. However, all of these samples rely on baryonic
tracers of clusters and require complex astrophysical modeling to account for
the selection effect and the mass–observable relation. Incorrect or
insufficient modeling of these baryonic effects could lead to biases in the
cosmological constraints (Salvati et al., 2020; Grandis et al., 2021). Cluster
samples constructed from weak-lensing maps (Wittman et al., 2001; Schirmer et
al., 2007; Miyazaki et al., 2007; Shan et al., 2012; Miyazaki et al., 2015,
2018; Oguri et al., 2021), on the other hand, allow us to determine the
selection effect based purely on the theory of gravity, offering a direct link
between the weak-lensing observable and the underlying halo mass. Shear-
selected clusters, therefore, provide a powerful cosmology probe that
complements both cosmic shear power spectrum and traditional cluster
cosmology.

Despite the direct link between the weak-lensing observable and the underlying
projected density field, it is still essential to construct a comprehensive
modeling framework that takes into account all physical and systematic effects
associated with weak-lensing observations. In the past, weak-lensing peak
counts were often modeled analytically by using halo models on a Gaussian
random field (Fan et al., 2010; Shan et al., 2018) or semi-analytically by
injecting synthetic halo profiles into N-body simulations (Lin & Kilbinger,
2015). Important systematic effects such as dilution by cluster member
galaxies (Medezinski et al., 2018; Hamana et al., 2020; Oguri et al., 2021),
intrinsic alignment of source galaxies (Kacprzak et al., 2016; Harnois-Déraps
et al., 2021; Zhang et al., 2022), and other baryon physics (Osato et al.,
2015; Weiss et al., 2019; Coulton et al., 2020; Osato et al., 2021; Lee et
al., 2023) are hard to model and can only be accounted for using nuisance
parameters with strong priors. On the other hand, ray-tracing simulations have
also been carried out to study weak-lensing peak counts (Dietrich & Hartlap,
2010; Liu et al., 2015a; Kacprzak et al., 2016). However, these simulations
are generally expensive to compute and can only be carried out on a limited
grid of cosmological parameters.

In this work, we introduce a novel semi-analytical framework to model weak-
lensing peak counts. Instead of painting halos into N-body simulations, we
inject halos into the observed weak-lensing mass maps. This is achieved by
adding the lensing signals from synthetic halos onto the shear of real source
galaxies in the shape catalog. As the shape of these source galaxies already
contains realistic observational noise and real-world weak-lensing systematic
uncertainties, this allows us to accurately characterize the measurement
uncertainties on the weak-lensing observable. Meanwhile, by performing
injection simulations across the survey footprint, we also take into account
non-uniform imaging depth and the complex geometry of the survey due to bright
star masks and artifacts. Moreover, we mitigate the contamination from cluster
member galaxies by employing a stringent source selection (Oguri et al., 2021)
while still maintaining a high source density. This is only possible thanks to
the incredible imaging depth provided by the Subaru Hyper Suprime-Cam (HSC)
survey (Aihara et al., 2018). While these injection simulations are still
expensive to perform, in this work we also introduce a novel parametrization
method and choose halo properties that will make the derived selection
function and scaling relation independent of the underlying cosmology. The
framework and results developed in this paper will be incorporated into Chiu
et al. (2024) to obtain cosmological constraints.

This paper is organized as follows. In Sec. 2, we introduce the weak lensing
data, the construction of aperture mass maps, and the cluster catalog. Our
modeling framework and the details of the injection simulations are discussed
in Sec. 3. Results of the injection simulations, and the derived selection
function and scaling relation are presented in Sec. 4. In particular,
important validation tests of our modeling framework are shown in Sec. 4.3.
Conclusions are given in Sec. 5.

##  2 Data

###  2.1 Weak-lensing shape catalog

The shape catalog from the HSC-SSP S19A internal data release (Li et al.,
2022, hereinafter the HSC-Y3 data) contains roughly 36 million galaxies across
six fields (XMM, GAMA09H, WIDE12H, GAMA15H, VVDS, HECTOMAP)111The complete HSC
survey consists of only three distinct patches. However, as of year three, the
survey footprint is still fragmented. with an averaged
$19.9~{}\mathrm{arcmin}^{-2}$ effective galaxy number density. The galaxy
shapes are measured with the re-Gaussianization method (Hirata & Seljak, 2003)
and the shear estimation bias is calibrated using image simulations
(Mandelbaum et al., 2018a). A series of null tests performed in Li et al.
(2022) has demonstrated that various residual systematics are well under
controlled, making this a catalog suitable for accurate weak-lensing
measurements.

To further improve the accuracy of the weak-lensing observable associated with
galaxy clusters, we employ a conservative selection of source galaxies
depending on their photometric redshift. It is known that cluster member
galaxies dilute weak-lensing signals as their shapes are not distorted by the
cluster itself, thus reducing the average ellipticity around the galaxy
cluster (Medezinski et al., 2018; Hamana et al., 2020; Oguri et al., 2021).
Oguri et al. (2021) showed that a redshift cut on the source sample can
effectively mitigate the dilution effect for clusters selected in the weak-
lensing mass maps. Following their work, we therefore require all source
galaxies in our sample to have

$\int_{0.7}^{\infty}P(z)\,\mathrm{d}z\geq 0.95.$ (1)

Here, the probability density function $P(z)$ for each source galaxy is
obtained with an empirical fitting method from the Direct Empirical
Photometric code (DEmP; Hsieh & Yee, 2014; Tanaka et al., 2018; Nishizawa et
al., 2020). With this conservative cut, the total number of source galaxies
reduces to roughly 16.9 million galaxies with an averaged
$9.4\,\mathrm{arcmin}^{-2}$ effective galaxy number density. The combined
probability density function of source galaxies’ redshift before and after
this selection is shown as the dashed gray and solid green curve in the left
panel of Fig. 1.

1.175 fig1.pdf

Figure 1: Left: Combined probability density function of source galaxies’
photometric redshift in the HSC-Y3 shape catalog before (dashed grey) and
after (solid green) the selection criterion given in eq. (1). The redshift
distribution of weak-lensing peaks with signal-to-noise ratio larger than 4.7
is shown in the blue histogram. The redshift of the weak-lensing peaks is
obtained from cross-matching with optical cluster catalogs as described in
Sec. 2.3. Right: The aperture mass filters for convergence (red) and
tangential shear (blue) adopted in this work.

1.175 s19a-mass_map.png

Figure 2: Mass maps derived from the HSC-Y3 data with the truncated isothermal
filter given in Fig. 1. Weak-lensing peaks with a signal-to-noise ratio larger
than 4.7 are circled in green. These peaks are matched to optical clusters in
the following order: CAMIRA cluster with a high $P(N_{\mathrm{mem}},z|\nu)$
(purple square); low-redshift clusters found in the SDSS (brown hexagon); any
other optical clusters with $N_{\mathrm{mem}}\geq 20$ (blue diamond).

###  2.2 Construction of Mass Maps

To identify signals from galaxy clusters in our weak-lensing shape catalog, we
create aperture mass maps $M_{\kappa}(\boldsymbol{\mathbf{\theta}})$
(Schneider, 1996) defined as

$M_{\kappa}(\boldsymbol{\mathbf{\theta}})=\int\kappa(\boldsymbol{\mathbf{\theta^{\prime}}})U(|\boldsymbol{\mathbf{\theta^{\prime}}}-\boldsymbol{\mathbf{\theta}}|)\,\mathrm{d}\boldsymbol{\mathbf{\theta^{\prime}}},$
(2)

where $\kappa(\boldsymbol{\mathbf{\theta}})$ is a map of the weak-lensing
convergence and $U(\theta)$ is a filter that is chosen to maximize the signals
from galaxy clusters. We require the filter to be compensated, i.e.,

$\int U(\theta)\theta\,\mathrm{d}\theta=0,$ (3)

so that the aperture mass maps can be equivalently expressed in terms of the
tangential shear maps (Kaiser & Squires, 1993)

$\displaystyle M_{\kappa}(\boldsymbol{\mathbf{\theta}})$
$\displaystyle=\int\kappa(\boldsymbol{\mathbf{\theta^{\prime}}})U(|\boldsymbol{\mathbf{\theta^{\prime}}}-\boldsymbol{\mathbf{\theta}}|)\,\mathrm{d}\boldsymbol{\mathbf{\theta^{\prime}}}$
(4)
$\displaystyle=\int\gamma_{+}(\boldsymbol{\mathbf{\theta^{\prime}}};\theta)Q(|\boldsymbol{\mathbf{\theta^{\prime}}}-\boldsymbol{\mathbf{\theta}}|)\,\mathrm{d}\boldsymbol{\mathbf{\theta^{\prime}}},$

where
$\gamma_{+}(\boldsymbol{\mathbf{\theta^{\prime}}};\boldsymbol{\mathbf{\theta}})$
is the tangential shear at position $\boldsymbol{\mathbf{\theta^{\prime}}}$
with respect to $\boldsymbol{\mathbf{\theta}}$ and the filter $Q(\theta)$ is
related to $U(\theta)$ as (Kaiser et al., 1994)

$Q(\theta)=\frac{2}{\theta^{2}}\int_{0}^{\theta}U(\theta^{\prime})\theta^{\prime}\,\mathrm{d}\theta^{\prime}-U(\theta).$
(5)

In this work, we adopt a truncated isothermal profile (Schneider, 1996)

$U(\theta)=\begin{cases}1&\left(\theta\leq x_{1}\theta_{R}\right)\\\
\displaystyle\frac{1}{1-c}\left(\frac{x_{1}\theta_{R}}{\sqrt{\left(\theta-
x_{1}\theta_{R}\right)^{2}+\left(x_{1}\theta_{R}\right)^{2}}}-c\right)&\left(x_{1}\theta_{R}\leq\theta\leq
x_{2}\theta_{R}\right)\\\
\displaystyle\frac{b}{\theta_{R}^{3}}\left(\theta_{R}-\theta\right)^{2}\left(\theta-\alpha\theta_{R}\right)&\left(x_{2}\theta_{R}\leq\theta\leq\theta_{R}\right)\\\
0&\left(\theta_{R}\leq\theta\right)\end{cases}$ (6)

where the constants $\alpha,b,c$ are chosen so that the filter $U$ and its
first derivative are both continuous at $\theta=x_{2}\theta_{R}$ and the
normalization condition in Eq. (3) is satisfied. The overall shape of such a
filter is designed to roughly match the surface density distribution of a
galaxy cluster. Meanwhile, a constant filter for the convergence at
$\theta\leq x_{1}\theta_{R}$ makes us insensitive to the potential non-linear
and baryonic systematic uncertainties in the innermost region of a cluster,
and a smoothly truncated filter reduces possible discontinuities in our
aperture mass maps from the complex survey geometry. Following Oguri et al.
(2021), we adopt $x_{1}=0.121$, $x_{2}=0.36$, and
$\theta_{R}=16.6\,\mathrm{arcmin}$, which has been shown to be the most
effective in mitigating the dilution effect from member galaxies and in
suppressing the noise from uncorrelated large-scale structures. The resulting
filter $U(\theta)$ and $Q(\theta)$ is shown in the right panel of Fig. 1.

Meanwhile, we define the noise of the aperture mass map at each location
$\boldsymbol{\mathbf{\theta}}$ to be the standard deviation of
$M_{\kappa}(\boldsymbol{\mathbf{\theta}};\phi_{1},\phi_{2},\ldots,\phi_{n})$
where $\phi_{1},\phi_{2},\ldots,\phi_{n}$ are random variables drawn from a
uniform distribution between $[0,2\pi)$ denoting the position angle of source
galaxies. We thus have

$\displaystyle\sigma^{2}\left(\boldsymbol{\mathbf{\theta}}_{0}\right)=$
$\displaystyle\frac{1}{2\pi}\int_{0}^{2\pi}d\phi_{1}\cdots
d\phi_{n}M_{\kappa}^{2}(\boldsymbol{\mathbf{\theta}};\phi_{1},\ldots,\phi_{n})$
(7) $\displaystyle-\frac{1}{2\pi}\left[\int_{0}^{2\pi}d\phi_{1}\cdots
d\phi_{n}M_{\kappa}(\boldsymbol{\mathbf{\theta}};\phi_{1},\ldots,\phi_{n})\right]^{2}.$

In practice, we draw $500$ realizations of
$M_{\kappa}(\boldsymbol{\mathbf{\theta}};\phi_{1},\allowbreak\phi_{2},\ldots,\phi_{n})$
by randomly rotating galaxies and take the standard deviation of them as the
noise. The signal-to-noise ratio map $\nu(\boldsymbol{\mathbf{\theta}})$ is
then defined as

$\nu(\boldsymbol{\mathbf{\theta}})=\frac{M_{\kappa}(\boldsymbol{\mathbf{\theta}})}{\sigma(\boldsymbol{\mathbf{\theta}})}$
(8)

Our weak-lensing maps are sampled on a two-dimensional rectangular grid in the
$\boldsymbol{\mathbf{\theta}}$-coordinate with a pixel size of $0.25$ arcmin.
Fig. 2 shows six signal-to-noise maps derived from the truncated isothermal
filter. For detailed procedure of calculating these weak-lensing maps, we
refer the reader to Sec. 3.4 of Oguri et al. (2021).

###  2.3 Cluster catalogs

From the weak-lensing maps, we select peaks given a signal-to-noise ratio
threshold $\nu_{\mathrm{min}}$. Peaks are identified as pixels that are higher
in value than all eight pixels around them. We discard peaks if they are
within a $5$ arcmin radius of any other stronger peaks to avoid double
counting a single structure. We also discard peaks if they are situated on the
boundaries of the map. To check that the peaks we detected on the aperture
mass maps are indeed associated with massive clusters, we cross-match the
weak-lensing peaks with various optical cluster catalogs, with a primary focus
on the CAMIRA cluster sample (Oguri et al., 2018) derived from the HSC S21A
photometric data. We also check the cross-match with the redMaPPer (Rykoff et
al., 2014) and the WHL (Wen et al., 2012; Wen & Han, 2015) cluster samples
constructed from the overlapping Sloan Digital Sky Survey Data Release 8
(SDSS, Aihara et al., 2011), mainly to inspect low redshift clusters at
$z<0.1$, a regime where the CAMIRA sample does not cover.

For each weak-lensing peaks, we search for optical clusters within a $10$
arcmin radius and keep clusters that are within $1.5\,h^{-1}\mathrm{Mpc}$
using the redshift of the matched optical clusters. As each weak-lensing peak
could be matched to multiple optical clusters, we develop the following
framework to determine the best-matched optical counterpart. First, we search
for potential counterparts within the CAMIRA catalog using the above-mentioned
criteria based on both angular and physical distances. For each of the cross-
matched CAMIRA candidate of the observed weak-lensing peak with a signal-to-
noise ratio $\nu$, we adopt the richness-to-mass relation from Murata et al.
(2019) to estimate the probability of observing the counterpart with a
richness of $N_{\mathrm{mem}}$ at redshift $z_{\mathrm{cl}}$, denoting as
$P(N_{\mathrm{mem}},z_{\mathrm{cl}}|\nu)$ (See Appendix A for more details).
We define a cross-matched CAMIRA candidate as a high-confidence counterpart if
the richness and redshift of that cluster fall within the $3\sigma$ confidence
interval around the peak of the distribution
$P(N_{\mathrm{mem}},z_{\mathrm{cl}}|\nu)$. The best-matched optical
counterpart is chosen as the most likely among all the high-confidence
counterparts. If a peak does not admit a high-confidence counterpart, we then
search for the nearest low-$z$ counterpart in the redMaPPer and WHL catalogs
that has $z\leq 0.2$ and $N_{\mathrm{mem}}\geq 20$. Lastly, a peak is said to
have any optical counterpart if it is matched to any optical clusters with
$N_{\mathrm{mem}}\geq 15$ in either CAMIRA, redMaPPer or WHL cluster catalogs.

1.175 counterpart.pdf

Figure 3: Number of peaks as a function of the selection threshold
$\nu_{\mathrm{min}}$. The percentage of peaks with (high-confidence) optical
counterparts is shown in the (dashed purple) solid blue curve.

Fig. 2 shows an example of a weak-lensing peak sample constructed with
$\nu_{\mathrm{min}}=4.7$. Under this signal-to-noise ratio threshold, there
are in total $129$ peaks identified in the entire HSC-Y3 footprint. Among
which, $112$ are found to have high-confidence optical counterparts within the
CAMIRA catalog; $2$ peaks do not have high-confidence CAMIRA matches but admit
low-$z$ counterparts in SDSS, and $12$ are matched to optical counterparts
with $N_{\mathrm{mem}}\geq 15$ and are within $1.5~{}h^{-1}\mathrm{Mpc}$ to
the peak. While there are $3$ peaks that do not have any optical counterpart
under the distance and richness criteria given above, two of them (located at
$12\mathrm{h}39\mathrm{m}47.0645\mathrm{s},-00\mathrm{d}33\mathrm{m}52.8076\mathrm{s}$
and
$22\mathrm{h}24\mathrm{m}51.0696\mathrm{s}+03\mathrm{d}52\mathrm{m}47.9063\mathrm{s}$)
are surrounded by multiple lower-richness CAMIRA clusters and the remaining
one
($02\mathrm{h}33\mathrm{m}40.6278\mathrm{s}-03\mathrm{d}56\mathrm{m}00.1856\mathrm{s}$)
is near optical clusters identified by the Canada-France-Hawaii Telescope at
$z\sim 0.5$ (Durret et al., 2011; Sarron et al., 2018). Therefore, these peaks
are likely to be the results of chance projection of multiple halos along the
line-of-sight. Such a scenario is taken into account in our modeling framework
given in Sec. 3. A list of weak-lensing peaks in the HSC-Y3 footprint with
$\nu\geq 4.7$ and their best-mateched optical counterpart is given in Table 3.

The redshift distribution of the best-matched optical counterpart for weak-
lensing peaks with $\nu\geq 4.7$ is given as the blue histogram in the left
panel of Fig. 1. We see that all but one weak-lensing peaks are located at
$z\leq 0.7$. As we only utilize source galaxies with significant probability
of being above redshift $0.7$, combining with the fact that we exclude signal
from the innermost region of the halo, the weak-lensing observable we obtained
is safe from the dilution of cluster member galaxies.

Fig. 3 shows the number of peaks as a function of the selection threshold
$\nu_{\mathrm{min}}$ and the percentage of peaks with (high confidence)
optical counterparts. We observe that $96\%\,(100\%)$ of peaks with $\nu\geq
4.0\,(5.0)$ admit optical counterpart, indicating these weak-lensing peaks to
be an effective way to select galaxy clusters without any baryonic
assumptions. In Part II of our analysis (Chiu et al., 2024), we will utilize
these peak samples constructed with different selection threshold
$\nu_{\mathrm{min}}$ to test the robustness of our cosmology constraints.

We note that while the optical counterparts identified in this section can
provide extra information such as redshift and optical richness for the weak-
lensing peaks, in this work we merely utilize them to better understand our
samples. For our cosmological analysis (Chiu et al., 2024), we will only rely
on the weak-lensing observable and do not employ these optical information to
derive cosmological constraints.

ccccccccccc

List of weak-lensing peaks in the HSC-Y3 footprint with $\nu\geq 4.7$ and
their best-matched optical counterpart

Peak ID & RA Dec Peak $\mathcal{S}/\mathcal{N}$ $\textrm{RA}_{\mathrm{cl}}$
$\textrm{Dec}_{\mathrm{cl}}$ Redshift Richness Separation Optical Catalog
(deg) (deg) (deg) (deg) ($h^{-1}\mathrm{Mpc}$)

0 & 197.865421 -1.320876 10.38 197.886195 -1.332920 0.18 166.43 0.18 CAMIRA
1 143.800130 +0.813507 10.18 143.801208 +0.825611 0.35 68.73 0.15 CAMIRA
2 354.413644 +0.252496 9.41 354.415554 +0.271375 0.25 111.13 0.19 CAMIRA
3 145.118208 +2.452905 8.48 145.102438 +2.477641 0.16 65.08 0.20 CAMIRA
4 180.079850 +3.345811 8.05 180.105635 +3.347094 0.14 114.63 0.16 CAMIRA
5 140.561555 +3.804756 7.61 140.545649 +3.778190 0.25 69.27 0.30 CAMIRA
6 130.592751 +1.634875 7.34 130.591163 +1.640656 0.42 77.22 0.08 CAMIRA
7 37.928678 -4.887543 7.22 37.921574 -4.882613 0.19 115.79 0.07 CAMIRA
8 336.039295 +0.323747 7.05 336.040891 +0.325988 0.18 50.09 0.02 CAMIRA
9 146.207093 +2.726571 7.02 146.196298 +2.770387 0.22 84.54 0.40 CAMIRA
10 188.778711 -0.827654 7.01 188.778610 -0.859100 0.11 28.23 0.16 WHL
11 221.087980 +0.187586 6.99 221.097548 +0.214565 0.26 36.61 0.29 CAMIRA
12 143.632430 -0.379970 6.87 143.630004 -0.385442 0.34 46.96 0.07 CAMIRA
13 151.196876 -0.685035 6.76 151.215018 -0.662458 0.18 42.59 0.22 CAMIRA
14 179.027305 -0.338319 6.70 179.045036 -0.350253 0.25 83.20 0.21 CAMIRA
15 134.467340 +3.207578 6.65 134.474796 +3.176485 0.19 81.39 0.25 CAMIRA
16 179.703694 -1.081777 6.50 179.670749 -1.066721 0.15 34.72 0.24 CAMIRA
17 213.761632 -0.439816 6.49 213.651488 -0.406581 0.15 87.57 0.75 CAMIRA
18 156.740791 +0.080817 6.49 156.794729 +0.006912 0.34 30.61 1.10 redMaPPer
19 235.482589 +42.819316 6.45 235.487751 +42.821906 0.24 48.72 0.04 CAMIRA
20 129.811107 +2.449030 6.27 129.747767 +2.483233 0.35 24.28 0.90 CAMIRA
… … … … … … … … … …

Table 3 is published in its entirety in the machine-readable
format222https://github.com/inonchiu/hsc_shear_selected_clusters. A portion is
shown here for guidance regarding its form and content.

##  3 Modeling Framework and Injection Simulations

One of the biggest challenges in cluster cosmology is to understand and
quantify the mass-observable relation. For weak-lensing shear-selected
clusters, the most direct observable is the signal-to-noise ratio $\nu$
associated with the peak detected on the mass maps. Other observables, such as
cluster redshifts, could be obtained through matching peaks with other cluster
catalogs with available redshift information (Miyazaki et al., 2018; Oguri et
al., 2021) or through a dedicated confirmation tool (e.g., Klein et al., 2018,
2022). In this work, we only consider the cluster number counts as a function
of the signal-to-noise ratio $\nu$.

The number count of weak-lensing peaks as a function of their signal-to-noise
ratio $N(\nu)$ is related to the underlying halo mass function through

$\displaystyle\frac{\mathrm{d}N(\nu|\boldsymbol{\mathbf{p}})}{\mathrm{d}\nu}=\Omega_{\mathrm{survey}}\iint\,\mathrm{d}M\mathrm{d}z$
$\displaystyle\Bigg{[}\frac{\mathrm{d}n(M|z,\boldsymbol{\mathbf{p}})}{\mathrm{d}M}\frac{\mathrm{d}V(z|\boldsymbol{\mathbf{p}})}{\mathrm{d}z}$
(9) $\displaystyle\times P(\nu|M,z,\boldsymbol{\mathbf{p}})\Bigg{]},$

where $\Omega_{\mathrm{survey}}$ is the survey area,
${\mathrm{d}n(M|z,\mathbf{p})}/{\mathrm{d}M}$ is the halo mass function at a
given redshift $z$, $\mathrm{d}V(z|\boldsymbol{\mathbf{p}})/\mathrm{d}z$ is
the differential comoving volume element, and
$P(\nu|M,z,\boldsymbol{\mathbf{p}})$ is the probability of observing $\nu$
given the cluster halo mass $M$ at the redshift $z$. We note that all of these
terms are sensitive to cosmology, which is incorporated in the parameter
vector $\boldsymbol{\mathbf{p}}$. A correct prediction of the peak number
count therefore heavily depends on the probability distribution
$P(\nu|M,z,\boldsymbol{\mathbf{p}})$, widely referred to as the
mass–observable relation. It is precisely the goal of this paper to model the
mass–observable relation for these weak-lensing peaks.

A widely adopted analytical model for weak-lensing peak counts was proposed in
Fan et al. (2010) and has been applied to constrain cosmology in Shan et al.
(2012, 2018) and Liu et al. (2023). However, such a model heavily relies on
the assumption that various sources of scatter on the weak-lensing observable
are Gaussian in nature and cannot probe the full complexity of these
systematic uncertainties. Lin & Kilbinger (2015) proposed a semi-analytical
model by injecting halo into N-body simulation to quantify the full impact of
projection from LSS, intrinsic alignment, and survey masks. In this work, we
adopt a more sophisticated semi-analytical model to quantify the
mass–observable relation. Instead of painting halos into N-body simulations,
we inject the lensing signals of halos into the observed weak-lensing mass
maps to determine $P(\nu|M,z,\boldsymbol{\mathbf{p}})$. This is achieved by
injecting the synthetic shear distortions into the ellipticity of the source
galaxies observed in the HSC-Y3 shape catalog.

By injecting the signals of synthetic clusters into the observed shape
catalogs, the resulting weak-lensing observable naturally accounts for the
observed noise and the systematic uncertainties inherent in the real universe.
This is because the shape distortion of these source galaxies already
encompasses contributions such as shape noise, the projection from
uncorrelated LSS, and the intrinsic alignment of source galaxies. Meanwhile,
by repeatedly injecting halos at various positions in our survey footprint, we
are able to characterize not only the complicated geometry of the survey
footprint, including the masking of bright stars, but also the variation in
imaging depth. In addition, our approach accounts for any numerical effects
associated with mass-map making in the real analysis, such as the flat-sky
projection and pixelization.

However, repeating these injection simulations for all possible cosmological
parameters $\boldsymbol{\mathbf{p}}$ to quantify
$P(\nu|M,z,\boldsymbol{\mathbf{p}})$ is computationally challenging. One
possible approach is to perform simulations on a grid of cosmological models
and leverage an emulator to interpolate among different cosmologies (e.g. Liu
et al., 2015a; Kacprzak et al., 2016; Nishimichi et al., 2019). Such a method
is still computationally challenging to scale up with an increasing amount of
cosmological and nuisance parameters that are necessary for cluster cosmology
analyses.

In this work, we propose a novel approach to dissect the selection function
and the mass–observable relation in a physical way and factor out the
cosmological dependence in the part that requires the computationally
intensive injection simulations. To be more specific, we re-write
$P(\nu|M,z,\boldsymbol{\mathbf{p}})$ into

$\displaystyle\,P(\nu|M,z,\boldsymbol{\mathbf{p}})$ (10) $\displaystyle=$
$\displaystyle\iint\mathrm{d}\hat{M_{\kappa}}\mathrm{d}\theta_{\mathrm{s}}P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})P(\hat{M_{\kappa}},\theta_{\mathrm{s}}|M,z,\boldsymbol{\mathbf{p}})\,.$

Here, $\hat{M_{\kappa}}$ and $\theta_{\mathrm{s}}$ are parameters representing
the analytic aperture mass and the characteristic angular scale of the
cluster, which will be precisely defined later in Sec. 3.3. These parameters
capture the shape of the lensing profile. This way, we divide the
mass–observable relation $P(\nu|M,z,\boldsymbol{\mathbf{p}})$ into two pieces:
(1) $P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$ is the re-parametrized
scaling relation that connects the shape of a lensing profile to the final
observable $\nu$; and (2)
$P(\hat{M_{\kappa}},\theta_{\mathrm{s}}|M,z,\boldsymbol{\mathbf{p}})$, which
converts the cluster’s physical properties $\left(M,z\right)$ to the shape of
the lensing profile parameterized by
$\left(\hat{M_{\kappa}},\theta_{\mathrm{s}}\right)$. The first term,
$P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$, accounts for the measurement
uncertainties including shape noise, cosmic shear, shape measurement bias, the
variation in survey depth, etc. This needs to be determined through
computationally expensive injection simulations. However, given a lensing
profile determined by $\left(\hat{M_{\kappa}},\theta_{\mathrm{s}}\right)$, the
process of injecting a synthetic cluster signal, re-detecting it on the
resulting mass maps, and quantifying the relation
$P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$ is insensitive to cosmological
parameters and can be computed in advance. The cosmological dependency goes
entirely into the conversion
$P(\hat{M_{\kappa}},\theta_{\mathrm{s}}|M,z,\boldsymbol{\mathbf{p}})$. This
term also captures the intrinsic uncertainties, such as the intrinsic scatter
in the weak-lensing observable due to the diversity of halo profiles (e.g.,
halo concentration and triaxiality) and the weak-lensing mass bias arising
from the inaccurate assumption about the halo profiles. The modeling of
$P(\hat{M_{\kappa}},\theta_{\mathrm{s}}|M,z,\boldsymbol{\mathbf{p}})$ will be
examined in-depth in our cosmological analysis (Chiu et al., 2024).

In this paper, we focus on deriving
$P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$. In Sec. 3.1, we discuss the
mock cluster samples used in the injection simulations. Details of the
injection simulations are given in Sec. 3.2. In Sec. 3.3, we introduce the
parameters $\left(\hat{M_{\kappa}},\theta_{\mathrm{s}}\right)$ and provide an
analytical argument on how they determine the lensing profile, thus ensuring
the cosmological independence of
$P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$.

###  3.1 The Mock Cluster Samples

We generate mock cluster catalogs analytically and inject their weak-lensing
signals into the observed shape catalogs to quantify the mass–observable
relation $P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$. With the goal of
sampling the parameter space of
$\left(\hat{M_{\kappa}},\theta_{\mathrm{s}}\right)$ as uniform and as complete
as possible, we adopt the following sampling strategy: We first create a fine
uniform grid on the mass $M_{200\mathrm{c}}$ and redshift $z_{\mathrm{cl}}$
space. Here, $M_{200\mathrm{c}}$ stands for the mass enclosed within the
radius of $r_{200\mathrm{c}}$, where the average density inside is $200$ times
the critical density of the universe. For each point on the
$M_{200\mathrm{c}}$–$z_{\mathrm{cl}}$ plane, we create our cluster sample with
a uniform sampling on the characteristic angular size $\theta_{\mathrm{s}}$ of
the cluster. Here, $\theta_{\mathrm{s}}$ is linked to the halo concentration
parameter in a cosmology-dependent way:

$\theta_{\mathrm{s}}\coloneqq\frac{r_{s}}{D_{A}(z_{\mathrm{cl}})};\quad\textrm{where
}r_{s}\coloneqq\frac{r_{200\mathrm{c}}}{c_{200\mathrm{c}}},$ (11)

and $D_{\mathrm{A}}$ is the angular diameter distance. Note that a fixed set
of $(M_{200\mathrm{c}},z_{\mathrm{cl}},\theta_{\mathrm{s}})$ uniquely
determines the halo concentration $c_{200\mathrm{c}}$ which in turn fixes the
shape of the halo density profile if we assume the spherical Navarro, Frenk, &
White 1997 (hereinafter NFW) model. It is worth mentioning that a uniform
sampling in the space of $M_{200\mathrm{c}},z_{\mathrm{cl}}$, and
$\theta_{\mathrm{s}}$ does not imply a uniform sampling in $\hat{M}_{\kappa}$
and $\theta_{\mathrm{s}}$. However, the sampling range in
$(M_{200c},z_{\mathrm{cl}},\theta_{\mathrm{s}})$ used to generate mock
clusters is large enough to cover the the range of
$\left(\hat{M}_{\kappa},\theta_{\mathrm{s}}\right)$ of interest.

In practice, we sample 2554880 clusters for each connected field in our
survey. These clusters are sampled uniformly in the log-parameter space with
the boundaries of $5\times 10^{12}\,M_{\odot}h^{-1}\leq M_{200\mathrm{c}}\leq
5\times 10^{16}\,M_{\odot}h^{-1}$, $0.02\leq z_{\mathrm{cl}}\leq 2.00$, and
$10^{-2}\,\mathrm{arcmin}\leq\theta_{\mathrm{s}}\leq 3\times
10\,\mathrm{arcmin}$.

1.175 injection_before_after_diff.pdf

Figure 4: Examples of injecting cluster lensing signal into the HSC-Y1 XMM
field shape catalog following the procedure outlined in Sec. 3.2. The left and
the right panel each shows two cases that are accepted or rejected by step (v)
in Sec. 3.2 in order to remove cases where we selected peaks that are not
mainly associated with the injected halos. The dashed circles in these plots
indicate the $5$ arc minute range around
$\boldsymbol{\mathbf{x}}_{\mathbf{diff}}$ above which we reject a detection.

###  3.2 Injection Simulations

To quantify the mass–observable relation unique to each field in our survey,
we inject the weak-lensing signals of mock clusters into the observed shape
catalog. For each of the connect fields in our survey and each of the mock
clusters, we perform an injection in the following steps:

  1. (i)

We randomly assign a position $\boldsymbol{\mathbf{\theta}}$ in the minimum
bounding rectangle of a given field. As a result, mock clusters are injected
not only into the region where we observe source galaxies but also into masked
areas or areas close to the survey boundary. This is because it is possible to
have a cluster where its center is masked, but the outskirts of its lensing
profile could still create an off-centered peak in the mass map due to the
presence of shape noise and cosmic shear.

  2. (ii)

For each randomly assigned sky position $\boldsymbol{\mathbf{\theta}}$, we cut
out a $4\theta_{\mathrm{out}}\times 4\theta_{\mathrm{out}}$ square centered at
$\boldsymbol{\mathbf{\theta}}$. $\theta_{\mathrm{out}}$ denotes the radius of
the filter we choose to create the weak-lensing mass maps. Here, we adopt the
flat-sky approximation and use the tangent plane projection as described in
Oguri et al. (2021). Then, for every source galaxy in the cut-out, we assign
it a redshift value $z_{\mathrm{gal}}$ based on a random draw from its
photometric redshift distribution $P(z_{\mathrm{gal}})$. With this, we can
calculate the distortion due to the injected cluster on every source galaxy
within the patch.

  3. (iii)

To calculate the lensing signal, we assume that the mass density profile of
the mock cluster follows a spherical NFW profile.

$\displaystyle\rho_{\mathrm{NFW}}\left(r\right)$
$\displaystyle=\frac{\rho_{\mathrm{s}}}{\left(r/r_{\mathrm{s}}\right)\left(1+r/r_{\mathrm{s}}\right)^{2}}$
(12) $\displaystyle\rho_{s}$
$\displaystyle\coloneqq\frac{200\rho_{\mathrm{cri}}\left(z_{\mathrm{cl}}\right)}{3}\frac{c_{200\mathrm{c}}^{3}}{\ln(1+c_{200\mathrm{c}})-c_{200\mathrm{c}}/(1+c_{200\mathrm{c}})}\,.$

The characteristic radius $r_{s}$ and the concentration parameter
$c_{200\mathrm{c}}$ is derived from the given set of
$(M_{200\mathrm{c}},z_{\mathrm{cl}},\theta_{\mathrm{s}})$ and cosmology. It is
then straightforward to calculate the shear and the convergence of a source
galaxy at redshift $z_{\mathrm{gal}}$ and the sky position
$\boldsymbol{\mathbf{\theta}}_{\mathrm{gal}}$:

$\displaystyle\gamma_{\mathrm{NFW}}\left(\theta_{\mathrm{gal}}\right)$
$\displaystyle=\frac{\bar{\Sigma}(<\theta_{\mathrm{gal}})-\Sigma(\theta_{\mathrm{gal}})}{\Sigma_{\mathrm{cri}}(z_{\mathrm{cl}},z_{\mathrm{gal}})}$
(13) $\displaystyle\kappa_{\mathrm{NFW}}\left(\theta_{\mathrm{gal}}\right)$
$\displaystyle=\frac{\Sigma(\theta_{\mathrm{gal}})}{\Sigma_{\mathrm{cri}}(z_{\mathrm{cl}},z_{\mathrm{gal}})}\,,$

where the projected surface density $\Sigma(\theta)$ is defined as,

$\Sigma\left(\theta\right)\coloneqq\int\mathrm{d}x_{3}\,\rho_{\mathrm{NFW}}\left(\sqrt{x_{3}^{2}+\Big{(}\theta
D_{\mathrm{A}}\left(z_{\mathrm{cl}}\right)\Big{)}^{2}}\right),$ (14)

and the critical surface density
$\Sigma_{\mathrm{cri}}^{-1}(z_{\mathrm{cl}},z_{\mathrm{gal}})$ is defined as

$\displaystyle\begin{cases}0&\quad\text{if }z_{\mathrm{gal}}\leq
z_{\mathrm{cl}},\\\ &\\\ \displaystyle\frac{4\pi
G}{c^{2}}\frac{D_{A}(z_{\mathrm{cl}})D_{A}(z_{\mathrm{cl}},z_{\mathrm{gal}})}{D_{A}(z_{\mathrm{gal}})}&\quad\text{otherwise}.\end{cases}$
(15)

The shear on the source galaxy due to the injected halo is expressed in terms
of shape measurement $(e_{1},e_{2})$ as

$e^{\mathrm{clu}}_{\mathrm{gal},i}=2\mathcal{R}_{\mathrm{gal}}\left[(1+m_{\mathrm{gal}})\frac{\gamma_{\mathrm{NFW},i}(\theta_{\mathrm{gal}})}{1-\kappa_{\mathrm{NFW}}(\theta_{\mathrm{gal}})}\right],$
(16)

where $\mathcal{R}_{\mathrm{gal}}$ is the per-galaxy shear responsivity and
$m_{\mathrm{gal}}$ is the per-galaxy calibration bias (Mandelbaum et al.,
2018b). The resulting shear is then added into the observed ellipticity of the
source galaxy. In the weak-lensing limit, the ellipticities of the source
galaxies before ($e^{\mathrm{before}}_{\mathrm{gal},i}$) and after
($e^{\mathrm{after}}_{\mathrm{gal},i}$) the injection is related to each other
as

$e^{\mathrm{after}}_{\mathrm{gal},i}=e^{\mathrm{before}}_{\mathrm{gal},i}+e^{\mathrm{clu}}_{\mathrm{gal},i}\,.$
(17)

We note that while we model the halo lensing signal by assuming a spherical
NFW profile here, we will model the deviation from this assumption through a
weak-lensing-mass-to-mass scaling relation in Chiu et al. (2024).

  4. (iv)

The shape catalog with the injected weak-lensing signal of the mock cluster is
then passed into the pipeline described in Sec. 2.2 to generate the signal-to-
noise ratio map. Following the procedure of constructing the observed peak
catalog, we first identify all the positive peaks in the resulting map. We
then discard all the peaks that are within a $5$ arcmin radius of any other
stronger peaks. The remaining peak that is the closest to the injection center
is selected. Therefore, for each of the halo we injected, we obtain two
observable, the peak signal-to-noise ratio $\nu_{\mathrm{inj}}$ and its
location $\boldsymbol{\mathbf{x}}_{\mathrm{inj}}$.

  5. (v)

Lastly, we need to determine whether the detected peak is mainly associated
with the injected halo instead of any existing structures on the weak-lensing
map. To address this, we introduce the difference map defined as

$\nu_{\mathrm{diff}}(\boldsymbol{\mathbf{\theta}})\coloneqq\frac{M^{\mathrm{after}}_{\kappa}(\boldsymbol{\mathbf{\theta}})-M^{\mathrm{before}}_{\kappa}(\boldsymbol{\mathbf{\theta}})}{\sigma(\boldsymbol{\mathbf{\theta}})}.$
(18)

We select the closest peak on the difference map in the same manner as
described in step (iv). The resulting peak on the difference map at the
location $\boldsymbol{\mathbf{x}}_{\mathrm{diff}}$ with the signal-to-noise
ratio $\nu_{\mathrm{diff}}$ serves as a proxy for the expected signal-to-noise
of the injected halo in the absence of cosmic shear and shape noise. With
$\left(\boldsymbol{\mathbf{x}}_{\mathrm{inj}},\nu_{\mathrm{inj}}\right)$ and
$\left(\boldsymbol{\mathbf{x}}_{\mathrm{diff}},\nu_{\mathrm{diff}}\right)$ of
the peaks detected in the post-injection and difference maps respectively, we
consider the injected halo to be detected if

$|\boldsymbol{\mathbf{x}}_{\mathrm{inj}}-\boldsymbol{\mathbf{x}}_{\mathrm{diff}}|\leq
5^{\prime};\quad|\nu_{\mathrm{inj}}-\nu_{\mathrm{diff}}|\leq\Delta_{s}\,,$
(19)

where $\Delta_{s}$ is a nuisance parameter that is not fixed a priori and
needs to be calibrated by the data.

Fig. 4 shows four examples in our injection simulation following the procedure
outlined above. For each example, we show the weak-lensing map before and
after we inject the mock cluster lensing signal as described in step (iii)
together with the difference map introduced in step (v). The right panel
demonstrates two cases that are rejected in step (v) to avoid detecting peaks
that are not mainly associated with the injected halos. The top row shows a
rejected case as it fails to meet the
$|\boldsymbol{\mathbf{x}}_{\mathrm{inj}}-\boldsymbol{\mathbf{x}}_{\mathrm{diff}}|\leq
5^{\prime}$ criterion. The bottom row shows another case that is rejected as
it has $|\nu_{\mathrm{inj}}-\nu_{\mathrm{diff}}|>\Delta_{s}$ if we choose
$\Delta_{s}=3.0$. This does not imply that the detected peak here is
completely unassociated with the injected halo. The same case would be
considered as a detection if we allow for a larger ranging of scattering and
choose $\Delta_{s}=4.0$ for instance. Therefore, the allowed range of
scattering $\Delta_{s}$ cannot be determined through injection simulation
alone and must be calibrated by leveraging external information. In Sec. 4.3,
we discuss how $\Delta_{s}$ affects the number count of weak lensing peaks
predicted by our framework and leverage optical counterparts to inform
$\Delta_{s}$. In our cosmological analysis (Chiu et al., 2024), the scatter
parameter $\Delta_{s}$ will be self-calibrated by the number count $N(\nu)$ of
the shear-selected clusters.

Table 1: Different sets of reference cosmological parameters adopted in this
study.

Vanilla | $h=0.7\,\Omega_{m}=0.3,\,\Omega_{b}=0.05,\,\Omega_{k}=0,$
---|---
| $\sigma_{8}=0.80,\,n_{s}=0.95,\,w_{0}=-1.0$
Exotic | $h=0.7,\,\Omega_{m}=0.2,\,\Omega_{b}=0.05,\,\Omega_{k}=0,$
| $\sigma_{8}=0.65,\,n_{s}=1.90,\,w_{0}=-1.5$
Plancka | $h=0.673,\,\Omega_{m}=0.315,\,\Omega_{b}=0.049,\,\Omega_{k}=0,$
| $\sigma_{8}=0.812,\,n_{s}=0.9649,\,w_{0}=-1.0$

1Planck Collaboration et al. (2020)

We note that in calculating the lensing signal in step (iii), we need to
specify a certain cosmology. In Sec. 3.3, we show that we can parameterize the
mass–observable relation in a way that is independent of the choice
cosmological parameters. To verify this, we perform the injection simulations
outlined here for three different sets of cosmological parameters. These
different sets of parameters are summarized in Table 1 and the result is given
in Sec. 4.3.

We also note that as we adopt a fully blinded framework for our cosmological
analysis by working with three copies of blinded shape catalogs. These
injection simulations are carried out for each of the blinded catalog. We
utilize unblinded shape catalog in a small subfield (the XMM field) from the
previous HSC S16A internal data release (Mandelbaum et al., 2018a, hereinafter
the HSC-Y1 data) to validate our pipeline and our modeling framework. All
figures shown in this section and the next are based on injection simulations
performed on the HSC-Y1 XMM field. For details of the HSC blinding strategy,
we refer the reader to Dalal et al. (2023); Li et al. (2023) and Part II of
this series (Chiu et al., 2024).

To summarize, by repeatedly performing the injection simulations, we can
obtain a mapping between any cluster properties to their observed signal-to-
noise ratio on the weak lensing maps. This allows us to calculate the
selection effect as a function of these cluster properties and also derive the
mass–observable relation for our cluster sample. In the following, we
demonstrate that we can choose two particular parameters to completely
describe the cluster lensing profile and parameterize the mass–observable
relation in a way that is independent of cosmological parameters.

1.175 single_inject.pdf

Figure 5: Examples of injecting a single mock halo at various locations on the
weak lensing mass map derived from the HSC-Y1 data. The singal-to-noise ratio
$\nu_{\mathrm{inj}}$ after each injection is given in the middle panel. Here,
not detected indicates a peak is not found on the map after we inject the
halo, most likely due to masking; and rejected indicates a peak is found but
does not satisfy the requirement given in step (v) in Sec. 3.2 for us to
determine that it is associated with the injected halo. The rightmost two
panels show the distributions of peak signal-to-noise ratio $\nu$ and the
distance between the observed peak and the halo center from 5000 injections of
this same halo. As a comparison, the distribution of peak signal-to-noise
ratio predicted from an uniform shape noise is indicated as the dashed blue
curve.

###  3.3 Re-parametrization of Scaling Relation

The weak-lensing profile of a given cluster involves a calculation that is
sensitive to cosmology. However, once we fix the shape of a lensing profile in
angular space, the conversion to the shear of the source galaxies, the
construction of weak-lensing maps, and the detection of peaks are all
independent of cosmological parameters. Thus, the key to parameterize the
selection function and the mass–observable relation into a cosmology-
independent form is to find parameters that describe the shape of the cluster
lensing profile on the weak-lensing map.

Following Sec.2.2, the lensing signal around a galaxy cluster placed at the
origin of a weak-lensing map is given by

$\displaystyle M_{\kappa}(\boldsymbol{\mathbf{\theta}})$
$\displaystyle=\iint\mathrm{d}\boldsymbol{\mathbf{\theta^{\prime}}}\kappa(\boldsymbol{\mathbf{\theta^{\prime}}})U(|\boldsymbol{\mathbf{\theta^{\prime}}}-\boldsymbol{\mathbf{\theta}}|)$
(20)
$\displaystyle=\iint\mathrm{d}\boldsymbol{\mathbf{\theta^{\prime}}}\frac{\Sigma(\boldsymbol{\mathbf{\theta}}^{\prime})}{\Sigma_{\textrm{cri}}(\boldsymbol{\mathbf{\theta}}^{\prime})}U(|\boldsymbol{\mathbf{\theta^{\prime}}}-\boldsymbol{\mathbf{\theta}}|).$

Although we do not know the true redshift of the source galaxies to determine
$\Sigma_{\textrm{cri}}$ and the true halo density profile to calculate
$\Sigma$, we can estimate $M_{\kappa}(\boldsymbol{\mathbf{\theta}})$ using the
ensemble average of
$\Sigma_{\mathrm{cri}}^{-1}(z_{\mathrm{cl}},z_{\textrm{gal}})$ and adopt the
lensing signal from an NFW-like halo. The lensing signal from an NFW-like halo
on the weak-lensing map is

$M_{\kappa}(\boldsymbol{\mathbf{\theta}})\approx
2\pi\rho_{s}r_{s}\left\langle\Sigma_{\mathrm{c}}^{-1}(z_{\mathrm{cl}})\right\rangle\int_{0}^{\theta_{\mathrm{out}}}U(\theta-\theta^{\prime})f\left(\frac{D_{A}(z_{\mathrm{cl}})\theta^{\prime}}{r_{s}}\right)\theta^{\prime}\mathrm{d}\theta^{\prime}$
(21)

where $f(x)$, given by Wright & Brainerd (2000), is

$f(x)\coloneqq\begin{cases}\displaystyle\frac{2}{x^{2}-1}\left[1-\frac{2}{\sqrt{1-x^{2}}}\operatorname{arctanh}\sqrt{\frac{1-x}{1+x}}\right]\quad
x<1,\\\ \displaystyle\frac{2}{3}\quad x=1,\\\
\displaystyle\frac{2}{x^{2}-1}\left[1-\frac{2}{\sqrt{x^{2}-1}}\operatorname{arctanh}\sqrt{\frac{x-1}{x+1}}\right]\quad
x>1.\end{cases}$ (22)

Because of the shape of our chosen filter, it was shown that the aperture mass
is less sensitive to the change of halo density profiles (Chen et al., 2020).
Together with a realistic estimate of
$\langle\Sigma_{\mathrm{cri}}^{-1}\rangle$ using source galaxies’ photometric
redshift distributions, this gives us a faithful representation of a cluster’s
signature on the weak-lensing maps. If we perform a change of variable to
$x\coloneqq D_{A}\theta^{\prime}/r_{s}$ and use $\theta_{s}=r_{s}/D_{A}$, we
get

$\displaystyle M_{\kappa}(\boldsymbol{\mathbf{\theta}})\approx$ $\displaystyle
2\pi\rho_{s}r^{3}_{s}D_{A}^{-2}\langle\Sigma_{\mathrm{cri}}^{-1}\rangle\int_{0}^{x_{\mathrm{out}}\coloneqq\theta_{\mathrm{out}}/\theta_{\mathrm{s}}}U\left(\theta-x\theta_{\mathrm{s}}\right)f(x)x\mathrm{d}x$
(23) $\displaystyle\eqqcolon$
$\displaystyle\hat{M_{\kappa}}\frac{\displaystyle\int_{0}^{x_{\mathrm{out}}}U\left(\theta-x\theta_{\mathrm{s}}\right)f(x)x\mathrm{d}x}{\displaystyle\int_{0}^{x_{\mathrm{out}}}U\left(x\theta_{\mathrm{s}}\right)f(x)x\mathrm{d}x}.$

Here, we define $\hat{M_{\kappa}}$ to be the peak value of the aperture mass
profile for an NFW-like halo

$\displaystyle\hat{M_{\kappa}}(M,z,\theta_{\mathrm{s}})$
$\displaystyle\coloneqq\left\langle\iint\mathrm{d}\boldsymbol{\mathbf{\theta^{\prime}}}\kappa_{\mathrm{NFW}}(\boldsymbol{\mathbf{\theta^{\prime}}})U(|\boldsymbol{\mathbf{\theta^{\prime}}}|)\right\rangle$
(24)
$\displaystyle=2\pi\rho_{s}r^{3}_{s}D_{A}^{-2}\langle\Sigma_{\mathrm{cri}}^{-1}\rangle\int_{0}^{x_{\mathrm{out}}\coloneqq\theta_{\mathrm{out}}/\theta_{\mathrm{s}}}U\left(x\cdot\theta_{\mathrm{s}}\right)f(x)x\mathrm{d}x,$

in which $\rho_{s},r_{s}$ are obtained from $M,z,\theta_{\mathrm{s}}$ assuming
a NFW-like halo, $D_{A}$ is the angular diameter distance at the halo’s
redshift, and $\left\langle\Sigma_{\mathrm{cri}}^{-1}\right\rangle$ is
obtained by an ensemble average of $\Sigma_{\mathrm{cri}}^{-1}$ using the
source galaxy redshift distribution averaged across the entire field where we
inject the halo into.

We thus see that the peak value of the estimated aperture mass profile
$\hat{M_{\kappa}}$ and the characteristic angular scale $\theta_{\mathrm{s}}$
completely fix the lensing profile for an NFW-like halo. This indicates that
the parameter set $\left(\hat{M_{\kappa}},\theta_{\mathrm{s}}\right)$ is the
sufficient statistic for the distribution $P(\nu)$ and motivates us to
postulate that $P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$ is insensitive to
cosmology. Our claim is validated numerically in Sec. 4.3.

##  4 Results

The injection simulations performed in Sec. 3 allow us to obtain a mapping
between the lensing properties
$\left(\hat{M_{\kappa}},\theta_{\mathrm{s}}\right)$ of a cluster to its
observed signal-to-noise ratio $\nu$. In this section, we show results of
these injections in Sec. 4.1 and the corresponding selection function and
mass–observable relation in Sec. 4.2. In Sec. 4.3, we perform validation tests
to examine the predicted number counts of weak lensing peaks. As we adopt a
fully blinded framework for our cosmological analysis by working with three
copies of blinded HSC-Y3 shape catalogs, the results shown in this section are
based on the data from a small subfield (the XMM field) in the HSC-Y1 shape
catalog.

1.175 selfunc.pdf

Figure 6: Completeness function and scaling relation derived from injecting
clusters into the HSC-Y1 XMM field. From left to right: the completeness
function
$\mathcal{C}_{\nu_{\mathrm{min}}}(\hat{M_{\kappa}},\theta_{\mathrm{s}})$; the
distribution $P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$ at three different
points on the $\hat{M_{\kappa}}$–$\theta_{\mathrm{s}}$ plane (indicated as
crosses in the leftmost panel); and the marginalized $\hat{M_{\kappa}}$–$\nu$
scaling relation, with each dot indicating the result from one injection. The
shaded regions in the rightmost panel represent the 68.45% confidence interval
of the distribution $P(\nu|\hat{M_{\kappa}})$ as a function of
$\hat{M_{\kappa}}$ derived from injection simulations (blue) and from
analytical calculations assuming only the shape noise (green).

###  4.1 Results from Injection Simulations

Fig 5 shows the results of injecting the same cluster into different locations
on the weak-lensing mass map. An example of the weak-lensing mass maps before
and after the injection is presented in the leftmost panel. Meanwhile, the
central panel shows that the same cluster is detected with various signal-to-
noise ratios $\nu$ at different locations, revealing significant scatter due
to the variations in both the signal from large-scale structures and the shape
noise observed across the survey footprint. By performing the injections
repeatedly on different locations on the maps, we are able to fully capture
the measurement uncertainties in $\nu$ directly from the data without any
assumption about the nature of these uncertainties. In some cases, the cluster
is not detected due to the survey masking; while in other cases, a peak is
found but is rejected due to the criterion in eq. 19. This means that either
the peak has large angular offset to the injected center ($>5\arcmin$), or the
signal-to-noise ratio $\nu$ of the peak is significantly scattered beyond the
range we consider ($|\nu_{\mathrm{inj}}-\nu_{\mathrm{diff}}|>\Delta_{s}$).

In the rightmost panels in Fig 5, we show the distribution of the signal-to-
noise ratio (top panel) and the offset between the detected peak and the
injection location (bottom panel). In the top panel, we also show the
distribution of $\nu$ assuming only the presence of an uniform and Gaussian
shape noise (blue dashed line) as a comparison. This is obtained by a normal
distribution centered at $\hat{M_{\kappa}}(M,z,c)/\sigma$ with $\sigma$
evaluated as (van Waerbeke, 2000; Mandelbaum et al., 2018b)

$\sigma^{2}=\frac{\sigma^{2}_{\epsilon}}{2n_{\mathrm{gal}}}\int_{0}^{\infty}2\pi|Q(\theta)|^{2}~{}\theta\mathrm{d}\theta,$
(25)

with $\sigma_{\epsilon}=0.4$ and $n_{\mathrm{gal}}$ estimated from the average
galaxy number density over the entire HSC-Y1 XMM subfield. While the mean of
the two distributions agrees with one another, the signal-to-noise values
predicted by the injection simulations exhibit larger scatter. This suggests
there are non-negligible contributions to the signal-to-noise ratio $\nu$ from
large-scale structures, intrinsic alignments, and the non-uniform survey
depth. In the bottom panel, the miscentering distribution shows that the
angular offset between the detected peaks and the true centers is mostly at a
level of $\lesssim 2\arcmin$. This is expected as the inner radius of the
filter we adopted is roughly $2\arcmin$ as shown in the right panel of Fig. 1.
In addition, we can obtain the probability of detecting this halo in this
subfield at a given signal-to-noise ratio threshold $\nu_{\mathrm{min}}$ by
integrating the distribution $P(\nu)$,

$\mathcal{C}_{\nu_{\mathrm{min}}}\coloneqq\int_{\nu_{\mathrm{min}}}^{\infty}P(\nu)\,\mathrm{d}\nu.$
(26)

For this particular halo, we obtain that it has a $\sim 33\%$ chance to be
detected in a shear-selected sample constructed with $\nu_{\mathrm{min}}=4.7$.
We note that this number depends on the area where we inject halos into. As we
inject into an area that is larger than the actual survey footprint, this
probability represents the detectability of the cluster multiplied by the
ratio of the survey footprint to the area where we inject halos into.

###  4.2 Completeness Function and Scaling Relation

With the results of the injection simulation for each cluster in our mock
catalog, we now obtain the scaling relation
$P(\nu|\hat{M}_{\kappa},\theta_{\mathrm{s}})$ for the weak-lensing observable.
Here, we also define the completeness function for a given selection threshold
$\nu_{\mathrm{min}}$ as

$\mathcal{C}_{\nu_{\mathrm{min}}}(\hat{M_{\kappa}},\theta_{\mathrm{s}})\coloneqq\int_{\nu_{\mathrm{min}}}^{\infty}P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})\mathrm{d}\nu$
(27)

In practice, $P(\nu|\hat{M}_{\kappa},\theta_{\mathrm{s}})$ is evaluated by
binning the mock clusters in the logarithmic space of
$\left(\hat{M_{\kappa}},\theta_{\mathrm{s}}\right)$ and calculate the
probability distribution of the signal-to-noise ratio $\nu$ at each bin.

Fig. 6 shows the selection function and scaling relation derived from
injection simulations performed over the HSC-Y1 XMM field with the selection
$\nu_{\mathrm{min}}=4.7$. The leftmost panel shows the selection function
$\mathcal{C}(\hat{M_{\kappa}},\theta_{\mathrm{s}})$; the middle panel shows
the scaling relation $P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$ at three
different points on the $\hat{M_{\kappa}}$–$\theta_{\mathrm{s}}$ plane
(indicated as crosses in the leftmost panel); and the rightmost panel shows
the scatter of $\nu$ as a function of $\hat{M_{\kappa}}$, with each dot
indicating the result from one injection. In the leftmost panel, we observe
that the completeness function does not reach unity even for the most massive
halos. As discussed before, this is because we inject halos into an area that
is larger than the actual survey footprint. The completeness function
saturates at the effective coverage of the survey footprint with respect to
the injection area.

Similar to Fig. 5, the dashed histograms in the middle panel show the
distributions of $\nu$ at different
$\left(\hat{M_{\kappa}},\theta_{\mathrm{s}}\right)$ assuming only the effect
of an uniform shape noise. We observe that the scatter in $\nu$ from the
injection-based results is larger than those of the shape-noise-only
distributions, regardless of the quantities
$\left(\hat{M_{\kappa}},\theta_{\mathrm{s}}\right)$. This is more clearly seen
in the rightmost panel, where we show the 68.45% confidence interval of the
distribution $P(\nu|\hat{M_{\kappa}})$ as a function of $\hat{M_{\kappa}}$
derived from injection simulations (blue) and from analytical calculations
assuming only the shape noise (green). At the high-$\hat{M_{\kappa}}$ end,
halos that have been partly obscured by the survey boundaries or masks could
still produce significant peaks to be detected, allowing more down-scatter; at
the low-$\hat{M_{\kappa}}$ end, significantly stronger up-scatter is predicted
by the injection simulations due to the presence of large-scale structures.

###  4.3 Validations

In this work, we derive the scaling relation
$P\left(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}}\right)$ directly from data to
account for measurement uncertainties on the weak-lensing observable $\nu$. In
what follows, we perform several validation tests on our modeling framework by
comparing cluster number counts from both data and simulation.

Following eq. (9) and (10), the number counts of shear-selected clusters as a
function of the signal-to-noise ratio $\nu$ are calculated as

$\displaystyle\frac{\mathrm{d}N(\nu|\boldsymbol{\mathbf{p}})}{\mathrm{d}\nu}=\Omega_{\mathrm{survey}}\iint\mathrm{d}M\mathrm{d}z\Bigg{[}\frac{\mathrm{d}n(M|z,\boldsymbol{\mathbf{p}})}{\mathrm{d}M}\frac{\mathrm{d}V(z|\boldsymbol{\mathbf{p}})}{\mathrm{d}z}$
(28)
$\displaystyle\times\iint\mathrm{d}\hat{M_{\kappa}}\mathrm{d}\theta_{\mathrm{s}}P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})P(\hat{M_{\kappa}},\theta_{\mathrm{s}}|M,z,\boldsymbol{\mathbf{p}})\Bigg{]}\,,$

where $\mathrm{d}n/\mathrm{d}M$ is modeled by the halo mass function given in
Bocquet et al. (2016) and $P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$ is
obtained by inject simulations performed over a $46.37\,\deg^{2}$ area that
includes the HSC-Y1 XMM footprint. While the actual survey footprint is
smaller, since $P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$ already contains
information about the effective coverage of the survey footprint in the
injection area, $\Omega_{\mathrm{survey}}$ is chosen accordingly to be
$46.37\,\deg^{2}$. The cosmological parameters $\boldsymbol{\mathbf{p}}$ are
fixed to the values measured by Planck Collaboration et al. (2020), which is
given in Table 1.

The remaining component yet to be specified is
$P(\hat{M_{\kappa}},\theta_{\mathrm{s}}|M,z,\boldsymbol{\mathbf{p}})$, the
cosmology-dependent conversion from halos’ physical properties $(M,z)$ to
their lensing properties $\left(M_{\kappa},\theta_{\mathrm{s}}\right)$. This
term will be comprehensively modeled in our cosmological analysis (Chiu et
al., 2024) to self-consistently take into account the intrinsic uncertainties
associated with the halo. For the purpose of validation, we adopt a simplified
modeling framework and write

$P(\hat{M_{\kappa}},\theta_{\mathrm{s}}|M,z)=\\\
\int\mathrm{d}M_{\mathrm{WL}}\,P(M_{\mathrm{WL}}|M,z)P(\hat{M_{\kappa}},\theta_{\mathrm{s}}|M_{\mathrm{WL}},z),$
(29)

where we introduce the distribution $P(M_{\mathrm{WL}}|M,z)$ to capture the
bias and uncertainties arose from the inaccurate assumption about the halo
density profile due to the presence of, e.g., the halo triaxiality (e.g.,
Becker & Kravtsov, 2011; Hamana et al., 2012). Here, we assume a redshift-
independent model for $P(M_{\mathrm{WL}}|M,z)$ where we consider $\ln
M_{\mathrm{WL}}\sim\mathcal{N}(\ln M,0.18^{2})$ (Chen et al., 2020).
Meanwhile, we have

$P(\hat{M_{\kappa}},\theta_{\mathrm{s}}|M_{\mathrm{WL}},z)=P(\hat{M_{\kappa}}|\theta_{\mathrm{s}},M_{\mathrm{WL}},z)P(\theta_{\mathrm{s}}|M_{\mathrm{WL}},z),$
(30)

where $P(\hat{M_{\kappa}}|\theta_{\mathrm{s}},M_{\mathrm{WL}},z)$ is a Dirac
delta function centered at
$\hat{M_{\kappa}}(M_{\mathrm{WL}},z,\theta_{\mathrm{s}})$ given by Eq. (24),
and $P(\theta_{\mathrm{s}}|M_{\mathrm{WL}},z)$ links the angular size of a
halo to its weak lensing mass $M_{\mathrm{WL}}$at a given redshift, which we
model through the mass–concentration relation. Here, we assume that the
mass–concentration relation follows a log-normal distribution around the mean
value predicted by Diemer & Joyce (2019) with a scatter $\sigma_{\ln c}=0.30$.

1.175 validation.pdf

Figure 7: Validation test on the invariance of the scaling relation
$P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$ derived under different
reference cosmologies. The three black lines (solid, dashed, dash-dotted) show
the number counts predicted for the HSC-Y1 XMM field using
$P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$ derived under three different
sets of cosmological parameters while the dotted blue line is the prediction
considering only an uniform shape noise. These results are compared to the
number counts in a mock cluster samples (orange circles) and the number counts
of real observed peak (green squares) in this field. The mock sample is
constructed through drawing a Poisson realization of the halo mass function
and directly inject them into the HSC-Y1 shape catalog to determine which of
them can be observed (see Sec. 4.3.2). For visualization purpose, the green
data points are shifted to the right by 0.05 to avoid overlap.

####  4.3.1 Cosmological Independence

One of the key results in this work is that the scaling relation
$P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$, which accounts for the
measurement uncertainties, is independent of the assumed cosmology in which it
is derived. Since a rigorous model for the measurement uncertainties requires
computationally intensive injection simulations, this allows us to compute
$P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$ in advance and apply them to the
cosmological analysis in Chiu et al. (2024).

To validate this, we repeat the injection simulations in Sec. 3.2 under three
different sets of cosmological parameters $\boldsymbol{\mathbf{p}}_{i}$ and
use the resulting
$P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}},\boldsymbol{\mathbf{p}}_{i})$ to
predict number counts

$\displaystyle\frac{\mathrm{d}N_{i}(\nu)}{\mathrm{d}\nu}=\Omega_{\mathrm{survey}}\iint\mathrm{d}M\mathrm{d}z\Bigg{[}\frac{\mathrm{d}n(M|z,\boldsymbol{\mathbf{p}})}{\mathrm{d}M}\frac{\mathrm{d}V(z|\boldsymbol{\mathbf{p}})}{\mathrm{d}z}$
(31)
$\displaystyle\times\iint\mathrm{d}\hat{M_{\kappa}}\mathrm{d}\theta_{\mathrm{s}}P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}},\boldsymbol{\mathbf{p}}_{i})P(\hat{M_{\kappa}},\theta_{\mathrm{s}}|M,z,\boldsymbol{\mathbf{p}})\Bigg{]}\,.$

Here, $\boldsymbol{\mathbf{p}}$ is fixed to the Planck cosmology to compute
terms that are known to depend on cosmology, while
$\boldsymbol{\mathbf{p}}_{i}$ represents the three sets of cosmological
parameters given in Table 1 that is used to compute
$P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$. The number counts predicted by
the three different
$P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}},\boldsymbol{\mathbf{p}}_{i})$ are
shown in Fig. 7. As can be seen, the resulting cluster number counts
$\mathrm{d}N_{i}/\mathrm{d}\nu$ are in excellent agreement with one another,
strongly demonstrating that the scaling relation
$P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$ is insensitive to the underlying
cosmology. This offers great computational advantages as it allows us to
compute $P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$ in advance and we will
only need to compute the other cosmology-dependent terms in eq.
(LABEL:eq:number_count_model) in the cosmological analysis.

####  4.3.2 Consistency in the Modeling Framework

To ensure the self-consistency and the correctness of our modeling framework,
we further compare the predicted number counts from eq.
(LABEL:eq:number_count_model) with the number counts obtained from a mock
observation and from the real observed catalog.

First, we test whether our modeling framework can recover the number counts of
a mock sample generated using a set of known cosmological and nuisance
parameters. The mock catalog is created as a Poisson realization of the halo
mass function (Bocquet et al., 2016) up to redshift 2 over an area of
$\Omega_{\mathrm{survey}}=46.37,\deg^{2}$. Each mock cluster is assigned a
weak-lensing mass, $M_{\mathrm{WL}}$, with respect to its true mass following
a log-normal distribution, $\ln M_{\mathrm{WL}}\sim\mathcal{N}(\ln
M,0.18^{2})$. Similarly, the concentration of each mock cluster is randomly
drawn from a log-normal distribution with a mean value predicted by Diemer &
Joyce (2019) and a scatter of $\sigma_{\ln c}=0.30$. Following steps (ii)-(iv)
in Sec.,3.2, these mock clusters are then injected into the XMM subfield of
the HSC-Y1 shape catalog, with the detected clusters forming a mock sample of
shear-selected peaks.

The number counts of the mock sample is shown as orange circles in Fig. 7. The
excenllent agreement between the number counts from a mock observations and
those predicted by eq. (LABEL:eq:number_count_model) indicates that the
modeling framework we develop is self-consistent. Moreover, we compare these
results directly with the number counts of shear-selected clusters observed in
this field (green squares). The good agreement (at a level of $\lesssim
1\sigma$) seen in Fig. 7 suggests that the modeling framework provides a good
description of the data.

1.175 dNdRichness-validation.pdf

Figure 8: Validation test on the shape of our predicted number counts of weak
lensing peaks’ optical counterparts as a function of their optical richness
given the scaling relation $P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$
derived with different allowed range of scattering $\Delta_{s}$ (see step (v)
in Sec. 3.2). These predictions are compared with the number counts of the
closest CAMIRA counterpart for weak lensing peaks found in the HSC-Y1 XMM
field with a signal-to-noise ratio larger than 4.7 (green squares). A
statistically more constraining result using the entire HSC-Y3 field obtained
after unblinding the cosmological analysis in Chiu et al. (2024) is shown in
Fig. 9 in Appendix A.

####  4.3.3 Comparison with Cross-matched Optical Clusters

Lastly, we compare the richness distributions between the observed shear-
selected clusters’ optical counterparts and the theoretical prediction using a
calibrated richness–mass relation. Based on the optical richness
$N_{\mathrm{mem}}$, this comparison provides not only an end-to-end test on
the modeling framework but also a constraint on the scattering parameter
$\Delta_{\mathrm{s}}$ used in the injection simulations (see Sec. 3.2).

Similar to eq. (LABEL:eq:number_count_model), we can write the number counts
of optical counterparts as a function of their richness as

$\frac{\mathrm{d}N\left(N_{\mathrm{mem}}\right)}{\mathrm{d}N_{\mathrm{mem}}}=\iiint\Theta\left(\nu-\nu_{\mathrm{min}}\right)P\left(N_{\mathrm{mem}},\nu|M,z\right)\times\\\
\frac{\mathrm{d}N\left(M,z|\boldsymbol{\mathbf{p}}\right)}{\mathrm{d}M\mathrm{d}z}\mathrm{d}\nu\mathrm{d}M\mathrm{d}z\,,$
(32)

where $\Theta\left(\nu-\nu_{\mathrm{min}}\right)$ is the Heaviside step
function that denotes the selection threshold, and the probability
$P\left(N_{\mathrm{mem}},\nu|M,z\right)$ describes the joint distribution of
$N_{\mathrm{mem}}$ and $\nu$ at a given cluster mass $M$ and redshift $z$. To
first-order approximation, we assume that the optical richness
$N_{\mathrm{mem}}$ to be independent of the weak-lensing signal-to-noise ratio
$\nu$,

$P\left(N_{\mathrm{mem}},\nu|M,z\right)=P\left(N_{\mathrm{mem}}|M,z\right)P\left(\nu|M,z\right)\,.$
(33)

We make use of the richness–mass scaling relation from Murata et al. (2019) to
evaluate $P\left(N_{\mathrm{mem}}|M,z\right)$. Together with the scaling
relation $P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$ derived from the
injection simulations, we can compute the predicted richness distribution of
the shear-selected clusters. When evaluating eq. (32), we note that we
restrict the redshift integral to be $z\geq 0.1$, as we compare the richness
distribution with that of the optical counterparts from CAMIRA which only
covers $z\geq 0.1$.

On the other hand, we cross-match the weak-lensing peaks observed in the
HSC-Y1 XMM field with optical clusters that have $N_{\mathrm{mem}}\geq 15$ in
the CAMIRA catalog. An optical counterpart can be identified within a 5 radius
around all peaks with $\nu\geq 4.7$. This gives the shear-selected clusters an
observed richness distribution that we can compare with. Note that in this
comparison the counterparts are obtained by a simple matching in the sky
coordinate and do not utilize the information of richness and redshift of the
optical clusters as done in Section 2.3 and Appendix A. By doing so, the
optical counterparts are determined without any prior knowledge of their
richness distribution, enabling an independent assessment on the resulting
weak-lensing scaling relation and selection function by comparing the richness
distribution of the counterparts.

Fig. 8 shows the comparison between the observed richness distribution and the
theoretical predictions using the scaling relation
$P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$ derived with different allowed
range of scattering $\Delta_{s}$ (see step (v) in Sec. 3.2). This plot is
generated with the cosmological parameters fixed to those measured by Planck.
The most distinct feature of changing the parameter $\Delta_{s}$ is revealed
by the overall shape of $\mathrm{d}N/\mathrm{d}N_{\mathrm{mem}}$: The smaller
(larger) value of $\Delta_{s}$ results in a decreasing (increasing) number of
clusters at the low-richness end. This is expected, given that the more low-
mass systems would be up-scattered into the shear-selected sample with an
increasing $\Delta_{s}$. We see that the observed richness distribution in the
HSC-Y1 XMM field (green squares) demonstrates agreement with the richness
distributions predicted with $\Delta_{s}\sim 4.5$. A statistically more
constraining result using the entire HSC-Y3 field obtained after unblinding
the cosmological analysis in Chiu et al. (2024) is shown in Fig. 9 in Appendix
A.

We stress that the normalization and shape of
$\mathrm{d}N/\mathrm{d}N_{\mathrm{mem}}$ are also affected by the cosmological
parameters (e.g., $\Omega_{\mathrm{m}}$ and $\sigma_{8}$) and the scatter in
both the richness–mass scaling relation and the weak-lensing mass bias.
Therefore, we must simultaneously constrain the parameter $\Delta_{s}$
together with all the cosmological and nuisance parameters in a self-
calibrating way. Moreover, because the optical counterparts shown in Fig. 8
are determined by only the matching in the sky coordinate, we expect some
falsely matched counterparts at the low-richness end compared to that obtained
from the more advanced matching including both the cluster redshift and
richness (See Appendix A). It is feasible but challenging to include the
richness distribution of the correctly matched counterparts into the data
vector and self-consistently constrain the parameter $\Delta_{s}$. We
therefore do not attempt to determine $\Delta_{s}$ here and leave it to future
work. In Chiu et al. (2024), $\Delta_{s}$ will be constrained together with
other parameters based on the number counts $N\left(\nu\right)$ only.

##  5 Conclusions

We have presented the construction of a shear-selected cluster sample through
identifying peaks on weak-lensing aperture mass maps derived from the HSC-Y3
shear catalog. Thanks to the exquisite depth of the HSC survey, we have
selected source galaxies far in the background (with significant probability
of being above redshift $0.7$) when deriving our weak-lensing maps to reduce
the dilution of lensing signal from cluster member galaxies. We have also
chosen a filter that minimizes contributions from the center of the cluster
and from correlated large-scale structures. These choices have made our
cluster weak-lensing observable less sensitive to various bias associated with
weak-lensing peaks known in the literature. Over $\sim 510\,\mathrm{deg}^{2}$,
we have identified $129$ weak-lensing peaks with signal-to-noise ratio
$\nu\geq 4.7$. We have cross-matched these peaks to optical cluster catalogs
and confirmed that $\geq 99\%$ of peaks are associated with galaxy clusters.

To utilize this cluster sample for cosmological studies, we have performed
semi-analytical simulations to comprehensively model the mass–observable
relation and the selection effect on this cluster sample. We have injected
weak-lensing signal of NFW-like halos into the observed shape catalog to
capture measurement uncertainties from our actual survey. This method allows
us to simulate real-world weak lensing systematic effects such as the line-of-
sight large-scale structure, intrinsic alignment of source galaxies, and the
miss-centering of weak-lensing peaks. It also enables us to accurately
characterize the impact of the complex survey geometry and non-uniform survey
depth. For each of the fields in the HSC survey, we have injected halos across
a wide range of physical properties to derive the mass–observable relation as
a function these cluster properties. We have repeated this process for three
doubly-blinded shape catalogs.

As these semi-analytical simulations are still computationally intensive to
perform, we have developed a novel parametrization scheme to arrange the weak-
lensing scaling relation in terms of cluster lensing properties. If one has
chosen some parameters that can capture the shape of the cluster’s weak-
lensing profile in the angular space, the mapping between weak-lensing profile
to the weak-lensing observable will only be affected by observational
uncertainties that are insensitive to cosmological parameters. This allows us
to compute the weak-lensing scaling relation due to complex observational
uncertainties in advance and greatly enhances our ability to characterize the
shear-selected cluster sample in a rigorous but computationally feasible way.
In this work, we have chosen to parameterize the selection function and the
mass–observable relation with the analytic cluster aperture mass
$\hat{M_{\kappa}}$ and the characteristic angular scale of the cluster
$\theta_{\mathrm{s}}$. An analytical argument have been given in Sec. 3.3 to
illustrate why these two parameters are sufficient to describe the cluster
lensing profile. The cosmology-dependent conversion between cluster physical
properties such as mass and redshift to these lensing properties will be
presented in Chiu et al. (2024).

We have validated the weak-lensing scaling relation derived from our injection
scheme is indeed insensitive to changes of cosmological parameters. We have
derived three sets of scaling relations under three reference cosmologies
(Table 1). The resulting weak-lensing peak number counts derived from these
scaling relations have been shown to agree with each other excellently (Fig.
7). We have also compared these model predictions with the number counts of a
mock sample obtained by directly injecting a realization of halo catalog in
the weak-lensing mass maps. These tests have demonstrated the self-consistency
of our modeling framework. The correctness of our model has been further
verified by comparing our predictions with observed number counts of shear-
selected clusters as a function of both weak-lensing signal-to-noise ratio and
richness of the optical counterparts. These validation tests have been carried
out using a small subfield of the HSC-Y1 shear catalog so that we still
respect the blinded analysis.

Our results enable the cosmological study which will be presented in Chiu et
al. (2024). In Chiu et al. (2024), we will discuss our cosmological pipeline
that models uncertainties arisen from photometric redshift, and the deviation
of our analytical description of the halo profiles from true halo properties
such as halo triaxiality. The cosmological constraints derived from these
shear-selected clusters will shed light on the $S_{8}$ tension by serving as a
consistent check with other cosmological constraints obtained from weak
lensing. Looking forward, the methodology developed in this paper can be
applied to weak-lensing data from Stage-IV surveys (Ivezić et al., 2019;
Laureijs et al., 2011; Spergel et al., 2015) to serve as a strong complement
to cosmic shear as we further explore the non-Gaussian information from these
weak lensing data.

## Acknowledgements

K.-F.C. acknowledges support from the Taiwan Think Global Education Trust
Scholarship and the Taiwan Ministry of Education’s Government Scholarship to
Study Abroad. K.-F.C. thanks the computing resources provided by the Academia
Sinica Institute for Astronomy and Astrophysics in Taiwan which are maintained
by Dr. Bau-Ching Hsieh and others. K.-F.C. thanks the hospitality of Adrian
Liu at McGill University where this work is carried out in part. I-Non Chiu
acknowledges the support from the National Science and Technology Council in
Taiwan (Grant NSTC 111-2112-M-006-037-MY3) and the computing resources
provided by the National Center for High-Performance Computing (NCHC) in
Taiwan. YTL acknowledges supports by the grant NSTC 112-2112-M-001-061. This
work was supported by JSPS KAKENHI Grant Numbers JP20H05856, JP22H01260,
JP22K21349.

The Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) is led by the
astronomical communities of Japan and Taiwan, and Princeton University. The
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.
The survey was made possible by funding contributed by 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.

The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been
made possible through contributions by 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, the 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 grant No.
AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the
Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

This paper is based on data collected at the Subaru Telescope and retrieved
from the HSC data archive system, which is operated by the Subaru Telescope
and Astronomy Data Center (ADC) at NAOJ. Data analysis was in part carried out
with the cooperation of Center for Computational Astrophysics (CfCA), NAOJ. We
are honored and grateful for the opportunity of observing the Universe from
Maunakea, which has the cultural, historical and natural significance in
Hawaii.

This paper makes use of software developed for Vera C. Rubin Observatory. We
thank the Rubin Observatory for making their code available as free software
at http://pipelines.lsst.io/. Numerical calculations are performed through the
Python packages numpy (Oliphant, 2006–) and scipy (Virtanen et al., 2020).
Many of the cosmological and astrophysical calculations in this work rely on
the routines wrapped up in colossus (Diemer, 2018) and astropy (Astropy
Collaboration et al., 2022). Plots are made available thanks to matplotlib
(Hunter, 2007) and seaborn (Waskom, 2021).

## References

  * Abbott et al. (2020) Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2020, Phys. Rev. D, 102, 023509, doi: 10.1103/PhysRevD.102.023509
  * Abbott et al. (2022) —. 2022, Phys. Rev. D, 105, 023520, doi: 10.1103/PhysRevD.105.023520
  * Abdalla et al. (2022) Abdalla, E., Abellán, G. F., Aboubrahim, A., et al. 2022, Journal of High Energy Astrophysics, 34, 49, doi: 10.1016/j.jheap.2022.04.002
  * Aihara et al. (2011) Aihara, H., Allende Prieto, C., An, D., et al. 2011, ApJS, 193, 29, doi: 10.1088/0067-0049/193/2/29
  * Aihara et al. (2018) Aihara, H., Arimoto, N., Armstrong, R., et al. 2018, PASJ, 70, S4, doi: 10.1093/pasj/psx066
  * Aiola et al. (2020) Aiola, S., Calabrese, E., Maurin, L., et al. 2020, J. Cosmology Astropart. Phys, 2020, 047, doi: 10.1088/1475-7516/2020/12/047
  * Allen et al. (2011) Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409, doi: 10.1146/annurev-astro-081710-102514
  * Amon et al. (2022) Amon, A., Gruen, D., Troxel, M. A., et al. 2022, Phys. Rev. D, 105, 023514, doi: 10.1103/PhysRevD.105.023514
  * Anbajagane et al. (2023) Anbajagane, D., Chang, C., Banerjee, A., et al. 2023, arXiv e-prints, arXiv:2308.03863, doi: 10.48550/arXiv.2308.03863
  * Asgari et al. (2021) Asgari, M., Lin, C.-A., Joachimi, B., et al. 2021, A&A, 645, A104, doi: 10.1051/0004-6361/202039070
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Bacon et al. (2000) Bacon, D. J., Refregier, A. R., & Ellis, R. S. 2000, MNRAS, 318, 625, doi: 10.1046/j.1365-8711.2000.03851.x
  * Becker & Kravtsov (2011) Becker, M. R., & Kravtsov, A. V. 2011, ApJ, 740, 25, doi: 10.1088/0004-637X/740/1/25
  * Bergé et al. (2010) Bergé, J., Amara, A., & Réfrégier, A. 2010, ApJ, 712, 992, doi: 10.1088/0004-637X/712/2/992
  * Bocquet et al. (2016) Bocquet, S., Saro, A., Dolag, K., & Mohr, J. J. 2016, MNRAS, 456, 2361, doi: 10.1093/mnras/stv2657
  * Bocquet et al. (2019) Bocquet, S., Dietrich, J. P., Schrabback, T., et al. 2019, ApJ, 878, 55, doi: 10.3847/1538-4357/ab1f10
  * Chang et al. (2018) Chang, C., Pujol, A., Mawdsley, B., et al. 2018, MNRAS, 475, 3165, doi: 10.1093/mnras/stx3363
  * Chen et al. (2020) Chen, K.-F., Oguri, M., Lin, Y.-T., & Miyazaki, S. 2020, ApJ, 891, 139, doi: 10.3847/1538-4357/ab74d3
  * Chiu et al. (2024) Chiu, I. N., Chen, K.-F., Oguri, M., & others. 2024, arXiv e-prints.  https://arxiv.org/abs/xxxx.xxxxxx
  * Chiu et al. (2023) Chiu, I. N., Klein, M., Mohr, J., & Bocquet, S. 2023, MNRAS, 522, 1601, doi: 10.1093/mnras/stad957
  * Costanzi et al. (2021) Costanzi, M., Saro, A., Bocquet, S., et al. 2021, Phys. Rev. D, 103, 043522, doi: 10.1103/PhysRevD.103.043522
  * Coulton et al. (2020) Coulton, W. R., Liu, J., McCarthy, I. G., & Osato, K. 2020, MNRAS, 495, 2531, doi: 10.1093/mnras/staa1098
  * Dalal et al. (2023) Dalal, R., Li, X., Nicola, A., et al. 2023, arXiv e-prints, arXiv:2304.00701, doi: 10.48550/arXiv.2304.00701
  * de Haan et al. (2016) de Haan, T., Benson, B. A., Bleem, L. E., et al. 2016, ApJ, 832, 95, doi: 10.3847/0004-637X/832/1/95
  * Diemer (2018) Diemer, B. 2018, ApJS, 239, 35, doi: 10.3847/1538-4365/aaee8c
  * Diemer & Joyce (2019) Diemer, B., & Joyce, M. 2019, ApJ, 871, 168, doi: 10.3847/1538-4357/aafad6
  * Dietrich & Hartlap (2010) Dietrich, J. P., & Hartlap, J. 2010, MNRAS, 402, 1049, doi: 10.1111/j.1365-2966.2009.15948.x
  * Dodelson & Zhang (2005) Dodelson, S., & Zhang, P. 2005, Phys. Rev. D, 72, 083001, doi: 10.1103/PhysRevD.72.083001
  * Durret et al. (2011) Durret, F., Adami, C., Cappi, A., et al. 2011, A&A, 535, A65, doi: 10.1051/0004-6361/201116985
  * Fan et al. (2010) Fan, Z., Shan, H., & Liu, J. 2010, ApJ, 719, 1408, doi: 10.1088/0004-637X/719/2/1408
  * Fluri et al. (2022) Fluri, J., Kacprzak, T., Lucchi, A., et al. 2022, Phys. Rev. D, 105, 083518, doi: 10.1103/PhysRevD.105.083518
  * Friedrich et al. (2018) Friedrich, O., Gruen, D., DeRose, J., et al. 2018, Phys. Rev. D, 98, 023508, doi: 10.1103/PhysRevD.98.023508
  * Frieman et al. (2008) Frieman, J. A., Turner, M. S., & Huterer, D. 2008, ARA&A, 46, 385, doi: 10.1146/annurev.astro.46.060407.145243
  * Fu et al. (2014) Fu, L., Kilbinger, M., Erben, T., et al. 2014, MNRAS, 441, 2725, doi: 10.1093/mnras/stu754
  * Gatti et al. (2020) Gatti, M., Chang, C., Friedrich, O., et al. 2020, MNRAS, 498, 4060, doi: 10.1093/mnras/staa2680
  * Gatti et al. (2022) Gatti, M., Jain, B., Chang, C., et al. 2022, Phys. Rev. D, 106, 083509, doi: 10.1103/PhysRevD.106.083509
  * Ghirardini et al. (2024) Ghirardini, V., Bulbul, E., Artis, E., et al. 2024, arXiv e-prints, arXiv:2402.08458, doi: 10.48550/arXiv.2402.08458
  * Grandis et al. (2021) Grandis, S., Mohr, J. J., Costanzi, M., et al. 2021, MNRAS, 504, 1253, doi: 10.1093/mnras/stab869
  * Gruen et al. (2018) Gruen, D., Friedrich, O., Krause, E., et al. 2018, Phys. Rev. D, 98, 023507, doi: 10.1103/PhysRevD.98.023507
  * Hamana et al. (2012) Hamana, T., Oguri, M., Shirasaki, M., & Sato, M. 2012, MNRAS, 425, 2287, doi: 10.1111/j.1365-2966.2012.21582.x
  * Hamana et al. (2015) Hamana, T., Sakurai, J., Koike, M., & Miller, L. 2015, PASJ, 67, 34, doi: 10.1093/pasj/psv034
  * Hamana et al. (2020) Hamana, T., Shirasaki, M., & Lin, Y.-T. 2020, PASJ, 72, 78, doi: 10.1093/pasj/psaa068
  * Harnois-Déraps et al. (2021) Harnois-Déraps, J., Martinet, N., Castro, T., et al. 2021, MNRAS, 506, 1623, doi: 10.1093/mnras/stab1623
  * Heymans et al. (2021) Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140, doi: 10.1051/0004-6361/202039063
  * Hirata & Seljak (2003) Hirata, C., & Seljak, U. 2003, MNRAS, 343, 459, doi: 10.1046/j.1365-8711.2003.06683.x
  * Hsieh & Yee (2014) Hsieh, B. C., & Yee, H. K. C. 2014, ApJ, 792, 102, doi: 10.1088/0004-637X/792/2/102
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Jeffrey et al. (2020) Jeffrey, N., Lanusse, F., Lahav, O., & Starck, J.-L. 2020, MNRAS, 492, 5023, doi: 10.1093/mnras/staa127
  * Kacprzak et al. (2016) Kacprzak, T., Kirk, D., Friedrich, O., et al. 2016, MNRAS, 463, 3653, doi: 10.1093/mnras/stw2070
  * Kaiser & Squires (1993) Kaiser, N., & Squires, G. 1993, ApJ, 404, 441, doi: 10.1086/172297
  * Kaiser et al. (1994) Kaiser, N., Squires, G., Fahlman, G., & Woods, D. 1994, in 14th Moriond Astrophysics Meeting: Clusters of Galaxies (Gif-sur-Yvette: Edition Frontieres), 269–282
  * Kaiser et al. (2000) Kaiser, N., Wilson, G., & Luppino, G. A. 2000, arXiv e-prints, astro, doi: 10.48550/arXiv.astro-ph/0003338
  * Kilbinger (2015) Kilbinger, M. 2015, Reports on Progress in Physics, 78, 086901, doi: 10.1088/0034-4885/78/8/086901
  * Klein et al. (2018) Klein, M., Mohr, J. J., Desai, S., et al. 2018, MNRAS, 474, 3324, doi: 10.1093/mnras/stx2929
  * Klein et al. (2022) Klein, M., Oguri, M., Mohr, J. J., et al. 2022, A&A, 661, A4, doi: 10.1051/0004-6361/202141123
  * Kratochvil et al. (2012) Kratochvil, J. M., Lim, E. A., Wang, S., et al. 2012, Phys. Rev. D, 85, 103513, doi: 10.1103/PhysRevD.85.103513
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints, arXiv:1110.3193, doi: 10.48550/arXiv.1110.3193
  * Lee et al. (2023) Lee, M. E., Lu, T., Haiman, Z., Liu, J., & Osato, K. 2023, MNRAS, 519, 573, doi: 10.1093/mnras/stac3592
  * Lesci et al. (2022a) Lesci, G. F., Marulli, F., Moscardini, L., et al. 2022a, A&A, 659, A88, doi: 10.1051/0004-6361/202040194
  * Lesci et al. (2022b) Lesci, G. F., Nanni, L., Marulli, F., et al. 2022b, A&A, 665, A100, doi: 10.1051/0004-6361/202243538
  * Li et al. (2022) Li, X., Miyatake, H., Luo, W., et al. 2022, PASJ, 74, 421, doi: 10.1093/pasj/psac006
  * Li et al. (2023) Li, X., Zhang, T., Sugiyama, S., et al. 2023, arXiv e-prints, arXiv:2304.00702, doi: 10.48550/arXiv.2304.00702
  * Lin & Kilbinger (2015) Lin, C.-A., & Kilbinger, M. 2015, A&A, 576, A24, doi: 10.1051/0004-6361/201425188
  * Liu et al. (2015a) Liu, J., Petri, A., Haiman, Z., et al. 2015a, Phys. Rev. D, 91, 063507, doi: 10.1103/PhysRevD.91.063507
  * Liu et al. (2023) Liu, X., Yuan, S., Pan, C., et al. 2023, MNRAS, 519, 594, doi: 10.1093/mnras/stac2971
  * Liu et al. (2015b) Liu, X., Pan, C., Li, R., et al. 2015b, MNRAS, 450, 2888, doi: 10.1093/mnras/stv784
  * Mandelbaum (2018) Mandelbaum, R. 2018, ARA&A, 56, 393, doi: 10.1146/annurev-astro-081817-051928
  * Mandelbaum et al. (2018a) Mandelbaum, R., Lanusse, F., Leauthaud, A., et al. 2018a, MNRAS, 481, 3170, doi: 10.1093/mnras/sty2420
  * Mandelbaum et al. (2018b) Mandelbaum, R., Miyatake, H., Hamana, T., et al. 2018b, PASJ, 70, S25, doi: 10.1093/pasj/psx130
  * Mantz et al. (2010) Mantz, A., Allen, S. W., Rapetti, D., & Ebeling, H. 2010, MNRAS, 406, 1759, doi: 10.1111/j.1365-2966.2010.16992.x
  * Mantz et al. (2014) Mantz, A. B., Allen, S. W., Morris, R. G., et al. 2014, MNRAS, 440, 2077, doi: 10.1093/mnras/stu368
  * Marques et al. (2023) Marques, G. A., Liu, J., Shirasaki, M., et al. 2023, arXiv e-prints, arXiv:2308.10866, doi: 10.48550/arXiv.2308.10866
  * Martinet et al. (2018) Martinet, N., Schneider, P., Hildebrandt, H., et al. 2018, MNRAS, 474, 712, doi: 10.1093/mnras/stx2793
  * Medezinski et al. (2018) Medezinski, E., Battaglia, N., Umetsu, K., et al. 2018, PASJ, 70, S28, doi: 10.1093/pasj/psx128
  * Miyatake et al. (2023) Miyatake, H., Sugiyama, S., Takada, M., et al. 2023, arXiv e-prints, arXiv:2304.00704, doi: 10.48550/arXiv.2304.00704
  * Miyazaki et al. (2007) Miyazaki, S., Hamana, T., Ellis, R. S., et al. 2007, ApJ, 669, 714, doi: 10.1086/521621
  * Miyazaki et al. (2015) Miyazaki, S., Oguri, M., Hamana, T., et al. 2015, ApJ, 807, 22, doi: 10.1088/0004-637X/807/1/22
  * Miyazaki et al. (2018) —. 2018, PASJ, 70, S27, doi: 10.1093/pasj/psx120
  * More et al. (2023) More, S., Sugiyama, S., Miyatake, H., et al. 2023, arXiv e-prints, arXiv:2304.00703, doi: 10.48550/arXiv.2304.00703
  * Murata et al. (2019) Murata, R., Oguri, M., Nishimichi, T., et al. 2019, PASJ, 71, 107, doi: 10.1093/pasj/psz092
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493, doi: 10.1086/304888
  * Nishimichi et al. (2019) Nishimichi, T., Takada, M., Takahashi, R., et al. 2019, ApJ, 884, 29, doi: 10.3847/1538-4357/ab3719
  * Nishizawa et al. (2020) Nishizawa, A. J., Hsieh, B.-C., Tanaka, M., & Takata, T. 2020, arXiv e-prints, arXiv:2003.01511, doi: 10.48550/arXiv.2003.01511
  * Oguri et al. (2018) Oguri, M., Lin, Y.-T., Lin, S.-C., et al. 2018, PASJ, 70, S20, doi: 10.1093/pasj/psx042
  * Oguri et al. (2021) Oguri, M., Miyazaki, S., Li, X., et al. 2021, PASJ, 73, 817, doi: 10.1093/pasj/psab047
  * Oliphant (2006–) Oliphant, T. 2006–, NumPy: A guide to NumPy, USA: Trelgol Publishing.  http://www.numpy.org/
  * Osato et al. (2021) Osato, K., Liu, J., & Haiman, Z. 2021, MNRAS, 502, 5593, doi: 10.1093/mnras/stab395
  * Osato et al. (2015) Osato, K., Shirasaki, M., & Yoshida, N. 2015, ApJ, 806, 186, doi: 10.1088/0004-637X/806/2/186
  * Pandey et al. (2022) Pandey, S., Krause, E., DeRose, J., et al. 2022, Phys. Rev. D, 106, 043520, doi: 10.1103/PhysRevD.106.043520
  * Parroni et al. (2020) Parroni, C., Cardone, V. F., Maoli, R., & Scaramella, R. 2020, A&A, 633, A71, doi: 10.1051/0004-6361/201935988
  * Peel et al. (2018) Peel, A., Pettorino, V., Giocoli, C., Starck, J.-L., & Baldi, M. 2018, A&A, 619, A38, doi: 10.1051/0004-6361/201833481
  * Petri et al. (2013) Petri, A., Haiman, Z., Hui, L., May, M., & Kratochvil, J. M. 2013, Phys. Rev. D, 88, 123002, doi: 10.1103/PhysRevD.88.123002
  * Petri et al. (2015) Petri, A., Liu, J., Haiman, Z., et al. 2015, Phys. Rev. D, 91, 103511, doi: 10.1103/PhysRevD.91.103511
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A24, doi: 10.1051/0004-6361/201525833
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Porredon et al. (2022) Porredon, A., Crocce, M., Elvin-Poole, J., et al. 2022, Phys. Rev. D, 106, 103530, doi: 10.1103/PhysRevD.106.103530
  * Rozo et al. (2010) Rozo, E., Wechsler, R. H., Rykoff, E. S., et al. 2010, ApJ, 708, 645, doi: 10.1088/0004-637X/708/1/645
  * Rykoff et al. (2014) Rykoff, E. S., Rozo, E., Busha, M. T., et al. 2014, ApJ, 785, 104, doi: 10.1088/0004-637X/785/2/104
  * Salvati et al. (2020) Salvati, L., Douspis, M., & Aghanim, N. 2020, A&A, 643, A20, doi: 10.1051/0004-6361/202038465
  * Sarron et al. (2018) Sarron, F., Martinet, N., Durret, F., & Adami, C. 2018, A&A, 613, A67, doi: 10.1051/0004-6361/201731981
  * Schirmer et al. (2007) Schirmer, M., Erben, T., Hetterscheidt, M., & Schneider, P. 2007, A&A, 462, 875, doi: 10.1051/0004-6361:20065955
  * Schneider (1996) Schneider, P. 1996, MNRAS, 283, 837, doi: 10.1093/mnras/283.3.837
  * Secco et al. (2022) Secco, L. F., Samuroff, S., Krause, E., et al. 2022, Phys. Rev. D, 105, 023515, doi: 10.1103/PhysRevD.105.023515
  * Semboloni et al. (2011) Semboloni, E., Schrabback, T., van Waerbeke, L., et al. 2011, MNRAS, 410, 143, doi: 10.1111/j.1365-2966.2010.17430.x
  * Shan et al. (2012) Shan, H., Kneib, J.-P., Tao, C., et al. 2012, ApJ, 748, 56, doi: 10.1088/0004-637X/748/1/56
  * Shan et al. (2018) Shan, H., Liu, X., Hildebrandt, H., et al. 2018, MNRAS, 474, 1116, doi: 10.1093/mnras/stx2837
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv e-prints, arXiv:1503.03757, doi: 10.48550/arXiv.1503.03757
  * Sunayama et al. (2023) Sunayama, T., Miyatake, H., Sugiyama, S., et al. 2023, arXiv e-prints, arXiv:2309.13025, doi: 10.48550/arXiv.2309.13025
  * Sunyaev & Zeldovich (1972) Sunyaev, R. A., & Zeldovich, Y. B. 1972, Comments on Astrophysics and Space Physics, 4, 173
  * Takada & Jain (2003) Takada, M., & Jain, B. 2003, MNRAS, 344, 857, doi: 10.1046/j.1365-8711.2003.06868.x
  * Takada & Jain (2004) —. 2004, MNRAS, 348, 897, doi: 10.1111/j.1365-2966.2004.07410.x
  * Tanaka et al. (2018) Tanaka, M., Coupon, J., Hsieh, B.-C., et al. 2018, PASJ, 70, S9, doi: 10.1093/pasj/psx077
  * To et al. (2021) To, C., Krause, E., Rozo, E., et al. 2021, Phys. Rev. Lett., 126, 141301, doi: 10.1103/PhysRevLett.126.141301
  * Vafaei et al. (2010) Vafaei, S., Lu, T., van Waerbeke, L., et al. 2010, Astroparticle Physics, 32, 340, doi: 10.1016/j.astropartphys.2009.10.003
  * van Waerbeke (2000) van Waerbeke, L. 2000, MNRAS, 313, 524, doi: 10.1046/j.1365-8711.2000.03259.x
  * Van Waerbeke et al. (2000) Van Waerbeke, L., Mellier, Y., Erben, T., et al. 2000, A&A, 358, 30, doi: 10.48550/arXiv.astro-ph/0002500
  * Van Waerbeke et al. (2013) Van Waerbeke, L., Benjamin, J., Erben, T., et al. 2013, MNRAS, 433, 3373, doi: 10.1093/mnras/stt971
  * Vicinanza et al. (2019) Vicinanza, M., Cardone, V. F., Maoli, R., et al. 2019, Phys. Rev. D, 99, 043534, doi: 10.1103/PhysRevD.99.043534
  * Vikhlinin et al. (2009) Vikhlinin, A., Kravtsov, A. V., Burenin, R. A., et al. 2009, ApJ, 692, 1060, doi: 10.1088/0004-637X/692/2/1060
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, doi: 10.1038/s41592-019-0686-2
  * Waskom (2021) Waskom, M. L. 2021, Journal of Open Source Software, 6, 3021, doi: 10.21105/joss.03021
  * Weinberg et al. (2013) Weinberg, D. H., Mortonson, M. J., Eisenstein, D. J., et al. 2013, Phys. Rep., 530, 87, doi: 10.1016/j.physrep.2013.05.001
  * Weiss et al. (2019) Weiss, A. J., Schneider, A., Sgier, R., et al. 2019, J. Cosmology Astropart. Phys, 2019, 011, doi: 10.1088/1475-7516/2019/10/011
  * Wen & Han (2015) Wen, Z. L., & Han, J. L. 2015, ApJ, 807, 178, doi: 10.1088/0004-637X/807/2/178
  * Wen et al. (2012) Wen, Z. L., Han, J. L., & Liu, F. S. 2012, ApJS, 199, 34, doi: 10.1088/0067-0049/199/2/34
  * Wittman et al. (2001) Wittman, D., Tyson, J. A., Margoniner, V. E., Cohen, J. G., & Dell’Antonio, I. P. 2001, ApJ, 557, L89, doi: 10.1086/323173
  * Wittman et al. (2000) Wittman, D. M., Tyson, J. A., Kirkman, D., Dell’Antonio, I., & Bernstein, G. 2000, Nature, 405, 143, doi: 10.1038/35012001
  * Wright & Brainerd (2000) Wright, C. O., & Brainerd, T. G. 2000, ApJ, 534, 34, doi: 10.1086/308744
  * Zhang et al. (2022) Zhang, T., Liu, X., Wei, C., et al. 2022, ApJ, 940, 96, doi: 10.3847/1538-4357/ac9a4c
  * Zürcher et al. (2022) Zürcher, D., Fluri, J., Sgier, R., et al. 2022, MNRAS, 511, 2075, doi: 10.1093/mnras/stac078

##  Appendix A Confidence level of Cross-matched Optical Counterparts

To determine the best-matched optical counterparts in the CAMIRA catalog, we
develop the following framework to estimate the likelihood
$P(N_{\mathrm{mem}},z_{\mathrm{cl}}|\nu)$ and choose the best-match
counterpart to be the most likely.

$\displaystyle P(N_{\mathrm{mem}},z_{\mathrm{cl}}|\nu)$
$\displaystyle=\frac{1}{P(\nu)}\int
P(N_{\mathrm{mem}},z_{\mathrm{cl}},\nu,M)\,\mathrm{d}M$ (34)
$\displaystyle=\frac{1}{P(\nu)}\int
P(N_{\mathrm{mem}},\nu,|M,z_{\mathrm{cl}})P(M,z_{\mathrm{cl}})\,\mathrm{d}M$

To first-order approximation, we assume that the optical richness
$N_{\mathrm{mem}}$ to be independent of the weak-lensing signal-to-noise ratio
$\nu$ and write

$P\left(N_{\mathrm{mem}},\nu|M,z\right)=P\left(N_{\mathrm{mem}}|M,z\right)P\left(\nu|M,z\right).$
(35)

We make use of the richness–mass scaling relation from Murata et al. (2019) to
evaluate $P\left(N_{\mathrm{mem}}|M,z\right)$ as

$P\left(N_{\mathrm{mem}}|M,z\right)=\frac{1}{N_{\mathrm{mem}}\sqrt{2\pi}\sigma_{\ln
N\mid M,z}}\exp\left(-\frac{x^{2}(N_{\mathrm{mem}},M,z)}{2\sigma_{\ln N\mid
M,z}^{2}}\right),$ (36)

where $x(N,M,z)$ is given as

$x(N,M,z)\coloneqq\ln
N_{\mathrm{mem}}-\left[A+B\ln\left(\frac{M}{M_{\text{pivot
}}}\right)+B_{z}\ln\left(\frac{1+z}{1+z_{\text{pivot
}}}\right)+C_{z}\left[\ln\left(\frac{1+z}{1+z_{\text{pivot
}}}\right)\right]^{2}\right].$ (37)

The redshift-dependent scatter $\sigma_{\ln N\mid M,z}$ is modeled as

$\sigma_{\ln N\mid
M,z}\coloneqq\sigma_{0}+q\ln\left(\frac{M}{M_{\mathrm{pivot}}}\right)+q_{z}\ln\left(\frac{1+z}{1+z_{\mathrm{pivot}}}\right)+p_{z}\left[\ln\left(\frac{1+z}{1+z_{\mathrm{pivot}}}\right)\right]^{2}.$
(38)

Murata et al. (2019) obtained constraints on the parameter
$(A,B,B_{z},C_{z},\sigma_{0},q,q_{z},p_{z})$ by fitting to the HSC-Y1 shear
data assuming a pivot mass $M_{\mathrm{pivot}}=3\times
10^{14}\,h^{-1}M_{\odot}$ and redshift $z_{\mathrm{pivot}}=0.6$. Here, we
adopt the constraint derived under the Planck Collaboration et al. (2020)
cosmology
$(A,B,B_{z},C_{z},\sigma_{0},q,q_{z},p_{z})=(3.36,0.83,-0.20,3.51,0.19,-0.02,0.23,1.26)$.
We note that this constraint is obtained for the optical clusters in the
CAMIRA catalog with $N_{\mathrm{mem}}\geq 15$ and $0.1\leq z_{\mathrm{cl}}\leq
1.0$. Therefore, we only cross-match to CAMIRA clusters within these ranges.
For the other terms in eq. (34), $P\left(\nu|M,z\right)$ is obtained following
the treatment in Sec. 4.3

$P\left(\nu|M,z\right)=\iint\mathrm{d}\hat{M_{\kappa}}\mathrm{d}\theta_{\mathrm{s}}P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})P(\hat{M_{\kappa}},\theta_{\mathrm{s}}|M,z,\boldsymbol{\mathbf{p}}),$
(39)

while $P(M,z_{\mathrm{cl}})$ is derived from the halo mass function
$\mathrm{d}N/\mathrm{d}M/\mathrm{d}z$ (Bocquet et al., 2016). To be consistent
with the richness–mass scaling relation, $P\left(\nu|M,z\right)$ and $P(M,z)$
are also derived under the Planck Collaboration et al. (2020) cosmology.

1 dNdRichness.pdf

Figure 9: Same as Fig. 8, number counts of weak lensing peaks’ optical
counterparts as a function of their optical richness given the scaling
relation $P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$ derived with different
allowed range of scattering $\Delta_{s}$. These results are derived with
$P(\nu|\hat{M_{\kappa}},\theta_{\mathrm{s}})$ for HSC-Y3 after we unblind the
cosmological analysis in Chiu et al. (2024). The green squares show the number
counts of the closest CAMIRA counterpart for weak-lensing peaks in the entire
HSC-Y3 footprint with a signal-to-noise ratio larger than $4.7$.