11institutetext:  Tsung-Dao Lee Institute, and Key Laboratory for Particle
Physics, Astrophysics and Cosmology, Ministry of Education, Shanghai Jiao Tong
University, Shanghai 200240, China 22institutetext: Department of Astronomy,
School of Physics and Astronomy, and Shanghai Key Laboratory for Particle
Physics and Cosmology, Shanghai Jiao Tong University, Shanghai 200240, China
33institutetext: Academia Sinica Institute of Astronomy and Astrophysics
(ASIAA), 11F of AS/NTU Astronomy-Mathematics Building, No.1, Sec. 4, Roosevelt
Rd, Taipei 10617, Taiwan 44institutetext: Max Planck Institute for
Extraterrestrial Physics, Giessenbachstrasse 1, 85748 Garching, Germany
55institutetext: Faculty of Physics, Ludwig-Maximilians-Universität,
Scheinerstr. 1, 81679, Munich, Germany 66institutetext: IRAP, Universite de
Toulouse, CNRS, UPS, CNES, Toulouse, France 77institutetext: Department of
Physics, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033 Japan
88institutetext: Kavli Institute for the Physics and Mathematics of the
Universe (WPI), The University of Tokyo Institutes for Advanced Study (UTIAS),
The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa-shi, Chiba, 277-8583,
Japan 99institutetext: Kobayashi-Maskawa Institute for the Origin of Particles
and the Universe (KMI), Nagoya University, Nagoya, 464-8602, Japan
1010institutetext: Institute for Advanced Research, Nagoya University, Nagoya,
464-8601, Japan 1111institutetext: Division of Physics and Astrophysical
Science, Graduate School of Science, Nagoya University, Nagoya 464-8602, Japan
1212institutetext: The Inter University Centre for Astronomy and Astrophysics,
Ganeshkhind, Pune 411007, India 1313institutetext: Research Center for the
Early Universe, University of Tokyo, Tokyo 113-0033, Japan 1414institutetext:
Core Research for Energetic Universe, Hiroshima University, 1-3-1, Kagamiyama,
Higashi-Hiroshima, Hiroshima 739-8526, Japan 1515institutetext: Physics
Program, Graduate School of Advanced Science and Engineering, Hiroshima
University, 1-3-1 Kagamiyama, HigashiHiroshima, Hiroshima 739-8526, Japan
1616institutetext: Hiroshima Astrophysical Science Center, Hiroshima
University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan
1717institutetext: Argelander-Institut für Astronomie (AIfA), Universität
Bonn, Auf dem Hügel 71, 53121 Bonn, Germany

#  The eROSITA Final Equatorial-Depth Survey (eFEDS):

X-ray observable-to-mass-and-redshift relations of galaxy clusters and groups
with weak-lensing mass calibration from
the Hyper Suprime-Cam Subaru Strategic Program survey

I-Non Chiu （邱奕儂）,,, e-mail<EMAIL_ADDRESS>Vittorio Ghirardini 44
Ang Liu 44 Sebastian Grandis 55 Esra Bulbul 44 Y. Emre Bahar 44
Johan Comparat 44 Sebastian Bocquet 55 Nicolas Clerc 66 Matthias Klein 55 Teng
Liu 44 Xiangchong Li, 7788 Hironao Miyatake,,, 991010111188 Joseph Mohr, 5544
Surhud More, 121288 Masamune Oguri,, 13138877 Nobuhiro Okabe,, 141415151616
Florian Pacaud 1717 Miriam E. Ramos-Ceja 44 Thomas H. Reiprich 1717 Tim
Schrabback 1717 Keiichi Umetsu 33

We present the first weak-lensing mass calibration and X-ray scaling relations
of galaxy clusters and groups selected in the _eROSITA_ Final Equatorial Depth
Survey (eFEDS) observed by Spectrum Roentgen Gamma/_eROSITA_ over a contiguous
footprint with an area of $\approx 140$ deg2, using the three-year (S19A)
weak-lensing data from the Hyper Suprime-Cam (HSC) Subaru Strategic Program
survey. In this work, we study a sample of $434$ optically confirmed galaxy
clusters (and groups) at redshift $0.01\lesssim z\lesssim 1.3$ with a median
of $0.35$, of which $313$ systems are uniformly covered by the HSC survey to
enable the extraction of the weak-lensing shear observable. In a Bayesian
population modeling, we perform a blind analysis for the weak-lensing mass
calibration by simultaneously modeling the observed count rate $\eta$ and the
shear profile $g_{+}$ of individual clusters through the count-rate-to-mass-
and-redshift ($\eta$–$M_{500}$–$z$) relation and the weak-lensing-mass-to-
mass-and-redshift ($M_{\mathrm{WL}}$–$M_{500}$–$z$) relation, respectively,
while accounting for the bias in these observables using simulation-based
calibrations. As a result, the count-rate-inferred and lensing-calibrated
cluster mass is obtained from the joint modeling of the scaling relations, as
the ensemble mass spanning a range of $10^{13}h^{-1}\mathrm{M}_{\odot}\lesssim
M_{500}\lesssim 10^{15}h^{-1}\mathrm{M}_{\odot}$ with a median of $\approx
10^{14}h^{-1}\mathrm{M}_{\odot}$ for the eFEDS sample. With the mass
calibration, we further model the X-ray observable-to-mass-and-redshift
relations, including the rest-frame soft-band and bolometric luminosity
($L_{\mathrm{X}}$ and $L_{\mathrm{b}}$), the emission-weighted temperature
$T_{\mathrm{X}}$, the mass of intra-cluster medium $M_{\mathrm{g}}$, and the
mass proxy $Y_{\mathrm{X}}$, which is the product of $T_{\mathrm{X}}$ and
$M_{\mathrm{g}}$. Except for $L_{\mathrm{X}}$ with a steeper dependence on the
cluster mass at a statistically significant level, we find that the other
X-ray scaling relations all show a mass trend that is statistically consistent
with the self-similar prediction at a level of $\lesssim 1.7\sigma$.
Meanwhile, all these scaling relations show no significant deviation from the
self-similarity in their redshift scaling. Moreover, no significant redshift-
dependent mass trend is present. This work demonstrates the synergy between
the _eROSITA_ and HSC surveys in preparation for the forthcoming first-year
_eROSITA_ cluster cosmology.

###### Key Words.:

Galaxies: clusters: general – Galaxies: clusters: intracluster medium –
Gravitational lensing: weak – Cosmology: large-scale structure of Universe –
Cosmology: observations – Cosmology: dark energy

##  1 Introduction

The probe based on the abundance of galaxy clusters plays a crucial role in
constraining cosmology by utilizing the sample of clusters selected in the
optical (Costanzi et al., 2019b, 2021; To et al., 2021), X-rays (Mantz et al.,
2015; Schellenberger & Reiprich, 2017; Pacaud et al., 2018), and millimeter-
wavelength bands (Planck Collaboration et al., 2016; Bocquet et al., 2015; de
Haan et al., 2016; Bocquet et al., 2019) through the Sunyaev-Zel’dovich (SZ;
Sunyaev & Zel’dovich, 1972) effect. The number density of galaxy clusters as a
function of mass, the so-called halo mass function, is sensitive to
$\Omega_{\mathrm{M}}$, the density parameter of matter in the Universe, and
$\sigma_{8}$, the dispersion of linear density fluctuations on a comoving
scale of $8h^{-1}\mathrm{Mpc}$. Therefore, measurements of halo mass functions
over a wide range of mass and redshift essentially constrain the history of
the cosmic expansion as well as the growth of large-scale structures, enabling
a cosmological tool that is independent of and is as competitive as other
methods.

A necessary ingredient in cluster cosmology is a large sample of galaxy
clusters with a well-understood selection function. With recent wide and deep
surveys in the optical, such as the Dark Energy Survey (The Dark Energy Survey
Collaboration, 2005), a large sample of galaxy clusters and groups has been
constructed out to redshift $z\approx 0.6$ and over a footprint with an area
of more than a thousand square degrees. A sample of optically selected
clusters is now even extended to a much higher redshift at $z\gtrsim 1$ (e.g.,
Oguri et al., 2018a). Although a sizable sample of optically selected clusters
has largely improved cosmological constraints, modeling the selection function
of the clusters is challenging, mainly because of the projection effect (Zu et
al., 2017; Costanzi et al., 2019a; Sunayama et al., 2020), and could result in
a systematic bias in cosmological parameters (DES Collaboration et al., 2020).

On the other hand, searching for clusters based on the intra-cluster medium
(ICM), either in X-rays or through the SZ effect, provides a clean sample for
cosmological studies. The large SZ survey carried out by the South Pole
Telescope (Bleem et al., 2015) has enabled the construction of a highly pure
and complete sample of clusters out to redshift $z\approx 1.8$, based on which
unbiased cosmological constraints have been obtained (Bocquet et al., 2019).
An even larger sample of SZ-selected clusters has recently been released by
the Atacama Cosmology Telescope survey (Hilton et al., 2021). However, current
SZ surveys are limited to massive clusters because of the sensitivity of
detectors, or only scan a small fraction of the sky, resulting in a sample of
relatively small size. To extend the success of the SZ-based cosmology to the
low-mass regime, it is therefore necessary to combine the SZ sample with those
from other surveys (e.g., Costanzi et al., 2021). In terms of medium-size
X-ray surveys with advanced telescopes, such as the _XXL_ survey (Pierre et
al., 2016), a large number of clusters down to a low-mass regime has be
obtained (Adami et al., 2018), although mainly at low redshift because of the
cosmological dimming of bremsstrahlung emissions from clusters. Meanwhile, the
first all-sky survey in X-rays carried out by the _ROSAT_ mission (Truemper,
1982) had delivered the largest X-ray cluster catalog ever across the whole
sky (Böhringer et al., 2001, 2004; Piffaretti et al., 2011; Boller et al.,
2016), but the poor resolution and shallow depth have limited cosmological
studies in practice. However, this situation is no longer true with the
ongoing all-sky X-ray survey carried out by the state-of-the-art telescope,
_eROSITA_.

The _eROSITA_ (extended ROentgen Survey with an Imaging Telescope Array;
Merloni et al., 2012; Predehl et al., 2021) is an X-ray space telescope on
board the Russian-German “Spectrum-Roentgen-Gamma” (SRG) satellite (Sunyaev et
al., 2021), which was successfully launched on July 13, 2019. The main goal of
the _eROSITA_ mission is to unveil the nature of dark energy by carrying out
the _eROSITA_ All-Sky Survey (eRASS), which is to map the whole sky in X-rays
for four years. By doing so, the eRASS will deliver the largest sample of ICM-
selected galaxy clusters to date and discover more than 100 000 clusters at
the end of the four-year all-sky survey (Borm et al., 2014). With this
revolutionary sample, cluster-based cosmological analyses carried out as part
of the _eROSITA_ mission will provide unprecedented power with which to
constrain cosmology (Pillepich et al., 2012, 2018).

The key to the success of _eROSITA_ , and cluster cosmology in general, is an
accurate mass calibration on an observable-to-mass-and-redshift scaling
relation, which links the observed mass proxy to the true cluster mass at the
cluster redshift in order to construct unbiased halo mass functions (Pratt et
al., 2019). The technique of weak gravitational lensing (hereafter weak
lensing) has been considered as the most optimal way to calibrate cluster
mass, with huge success in practice (Umetsu et al., 2014; von der Linden et
al., 2014; von der Linden et al., 2014; Hoekstra et al., 2015; Okabe & Smith,
2016; Schrabback et al., 2018a; Dietrich et al., 2019; Okabe et al., 2019;
McClintock et al., 2019). A forecast from Grandis et al. (2019) shows that the
inclusion of weak-lensing mass calibrations using advanced optical surveys
could significantly tighten the constraints obtained from the _eROSITA_
sample, with uncertainties on $\Omega_{\mathrm{M}}$, $\sigma_{8}$, and the
dark energy equation of state $w$ at levels of a few percent.

Figure 1:  Angular distributions of eFEDS clusters and the available weak-
lensing data set from the HSC survey. Out of 542 eFEDS clusters, there are
$434$ clusters as a secure sample, i.e., $f_{\mathrm{cont}}<0.2$ (see Sect.
2.3). These clusters are represented by the circles, color-coded by their
observed redshifts $z_{\mathrm{cl}}$, with sizes proportional to their
observed count rate $\eta$. The other clusters with $f_{\mathrm{cont}}\geq
0.2$ are represented by the triangles. The underlying gray points represent a
subset of sources that are randomly drawn from the HSC three-year (S19A; see
Sect. 2.2) weak-lensing catalog.

However, a weak-lensing analysis of high-redshift clusters is challenging,
requiring substantially deep imaging with good resolution to resolve
background sources at even higher redshift. It is made feasible by using
pointing observations from space (Jee et al., 2011; Schrabback et al., 2018a,
2021) or deep high-resolution near-infrared imaging (Schrabback et al., 2018b)
for small samples, but not for a large number of clusters. It is also
important to stress that a uniform weak-lensing mass calibration across the
whole redshift range is needed; otherwise systematic errors due to
inhomogeneity would be introduced (Chiu et al., 2016a). To achieve a high-
quality weak-lensing mass calibration for high-redshift clusters ($z\gtrsim
0.7$) that is as consistent as the calibration available for low redshift, the
ongoing Hyper Suprime-Cam (HSC) Subaru Strategic Program survey (Aihara et
al., 2018a) is currently the only available resource for a large sample over
the _eROSITA_ footprint. The weak-lensing mass calibration out to redshift
$z\gtrsim 1.1$ based on the HSC data sets was performed for a large sample of
($\gtrsim 1700$) optically selected clusters over a footprint with an area of
$\approx 140$ deg2 (Murata et al., 2019). More recently, a weak-lensing
analysis using the HSC data sets was carried out for an X-ray-selected sample
of galaxy groups and clusters detected in the _XXL_ survey (Umetsu et al.,
2020), clearly demonstrating the potential of the HSC weak-lensing data for
the mass calibration of _eROSITA_ clusters.

To fully explore the survey capability of _eROSITA_ and its synergy with the
HSC survey, the _eROSITA_ team began the _eROSITA_ Final Equatorial Depth
Survey (eFEDS) in a footprint with an area of $\approx 140$ deg2 significantly
overlapping the HSC survey. By design, the depth of the eFEDS reaches the
average full depth of the eRASS in the equatorial area, thereby serving as a
performance verification phase before the start of the main survey. In this
work, we make use of the latest HSC imaging (S20A) for the optical
confirmation and redshift measurement of clusters selected in the eFEDS field,
followed by the weak-lensing mass calibration using the three-year HSC data
(S19A Li et al., 2021), which are proprietary and will be released as an
incremental data set of the third HSC Public Data Release. The goal of this
paper is to study various X-ray observable-to-mass-and-redshift scaling
relations with the mass calibration using the HSC weak-lensing data, setting
the stage for extending this analysis to the main eRASS work in the future.

The structure of this paper is as follows. In Sect. 2, we describe the cluster
sample and the data used in the analysis. The full description of the weak-
lensing analysis is given in Sect. 3. The simulation-based calibration of the
weak-lensing mass and the observed count rate in X-rays is given in Sect. 4,
followed by the modeling of scaling relations presented in Sect. 5. The
results and discussions are given in Sect. 6. We discuss the potential
systematic uncertainty in this work in Sect. 7. Finally, conclusions are made
in Sect. 8. Throughout this paper, we assume a fiducial cosmology, namely a
flat $\Lambda\textrm{CDM}$ cosmology with $\Omega_{\mathrm{M}}=0.3$, a mean
baryon density $\Omega_{\mathrm{b}}=0.05$, the Hubble parameter $H_{0}=h\times
100$ km s -1 Mpc-1 with $h=0.7$, $\sigma_{8}=0.8$, and a spectral index of the
primordial power spectrum $n_{\mathrm{s}}=0.95$. While modeling the scaling
relations, we vary the cosmological parameters around those of the fiducial
cosmology, with the aim being to maintain the flexibility needed in a future
cosmological analysis. The cluster true mass $M_{500}$ is defined by a sphere
with a radius of $R_{500}$, such that the enclosed mass density is $500$ times
the critical density $\rho_{\mathrm{c}}(z)$ of the Universe at the cluster
redshift. All quoted errors represent the $68\%$ confidence level (i.e.,
$1\sigma$) throughout this work, unless otherwise stated. The notation
$\mathcal{N}(x,y^{2})$ ($\mathcal{U}(x,y)$) stands for a normal distribution
with the mean $x$ and the standard deviation $y$ (a uniform distribution
between $x$ and $y$).

##  2 The cluster sample and data sets

In Sect. 2.1, a brief description of the eFEDS is given. We then briefly
introduce the Subaru HSC survey and the weak-lensing data in Sect. 2.2. The
cluster sample selected in the eFEDS is described in Sect. 2.3, while their
X-ray measurements are presented in Sect. 2.4.

###  2.1 The eROSITA Final Equatorial Depth Survey

The eFEDS is a small survey covering an area of $\approx 140$ deg2. The depth
of the eFEDS reaches an exposure time of $\approx 2.2$ $k$sec, roughly the
average full-depth of the eRASS, and therefore serves as a performance
verification phase for _eROSITA_ science, by design.

The survey took place at an orbit around the L2 point in a scanning mode. The
field of view of _eROSITA_ is $\approx 1$ deg2. The imaging is collected by
seven separate CCDs, each with its own mirror, with good detector uniformity
(no chip gaps). The on-axis half-energy-width is $\approx 18\arcsec$ at $1.49$
$k$eV, with an average of $\approx 26\arcsec$ for the whole field of view
(Predehl et al., 2021). The imaging quality is excellent; the accuracy of a
source location is better than $10\arcsec$ with a typical positional
uncertainty at the level of $\approx 4.6\arcsec$ ($1\sigma$). The survey
footprint consists of four scanning-mode units spanning a range of $\approx
126^{\circ}$ to $\approx 146^{\circ}$ ($\approx-3^{\circ}$ to
$\approx+6^{\circ}$) in Right Ascension (Declination) in total. The footprint
of the eFEDS is visualized in Fig. 1.

The uniqueness of the eFEDS survey is that the footprint largely overlaps the
HSC survey, which provides weak-lensing data products of the highest quality
to date for a large sample of clusters. This is especially true for galaxy
clusters at high redshift ($z\gtrsim 0.7$), where a weak-lensing study of a
large sample from the ground is only feasible using the HSC survey. This
provides a timely opportunity to statistically calibrate the mass of eROSITA-
detected clusters out to redshifts beyond unity before the era of the Euclid
mission (Laureijs et al., 2010) or the Legacy Survey of Space and Time (LSST;
Ivezic et al., 2008) carried out by the Vera C. Rubin Observatory. Therefore,
the synergy between the eFEDS and the HSC survey presented in this paper lays
the foundation for future work on this topic.

###  2.2 The HSC survey and the weak-lensing data

The HSC survey (Aihara et al., 2018a) is an imaging survey carried out as part
of a Subaru Strategic Program with the goal being to map a sky area of 1100
deg2 through five broadband filters ($grizy$). This is done using the wide-
field camera Hyper Suprime-Cam (Miyazaki, 2015; Miyazaki et al., 2018)
installed on the 8.2 m Subaru Telescope. The main goal of the HSC survey is to
perform state-of-the-art weak-lensing studies, paving a way forward for
science with the upcoming LSST.

The HSC survey comprises a three-layer imaging scheme: WIDE, DEEP, and
UltraDEEP. The $5\sigma$ limiting magnitudes of a $2\arcsec$ aperture in the
WIDE layer are $26.5$ mag, $26.1$ mag, $25.9$ mag, $25.1$ mag, and $24.4$ mag
for $g$-, $r$-, $i$-, $z$-, and $y$-band, respectively. This represents the
deepest imaging survey at the achieved area to date. Moreover, the imaging
quality is remarkable with a mean seeing of $0.58\arcsec$ at $i$-band
(Mandelbaum et al., 2018c). With this high-quality imaging, together with a
unique combination of depth and area, the data from the HSC survey are
excellent for weak-lensing studies over a large sample of clusters.

In this work, we use the three-year shape catalog constructed from the second
HSC Public Data Release (S19A) for weak-lensing studies. This shape catalog is
not released yet, and the full details are described in Li et al. (2021). In
what follows, we provide a brief summary of the three-year shape catalog. The
shape measurement is carried out in the $i$-band imaging following the same
methodology in constructing the first-year shape catalog (S16A) with a
comprehensive description given in Mandelbaum et al. (2018c). The shape
measurement is rigorously calibrated against the image simulations (Mandelbaum
et al., 2018a), such that the systematic uncertainty of the multiplicative
bias meets the requirement of $\lesssim 1.7\%$. Moreover, various null tests,
as well as those on the map level (Oguri et al., 2018b), are all statistically
consistent with zero. Only galaxies that satisfy the weak-lensing flags with
an $i$-band magnitude of smaller than $24.5$ mag are contained in the shape
catalog. With the high-quality of imaging, the source density in the HSC
survey reaches $\approx 22$ per arcmin2. It is worth mentioning that the same
method of the shape measurement was used to construct the first-year weak-
lensing catalog, which has been used to derive unbiased cosmic shears and
cosmological constraints (Hikage et al., 2019; Hamana et al., 2020). The
three-year catalog covers an area of $\approx 430$ deg2, in which the
footprint of the GAMA09 field significantly overlaps with the eFEDS. As a
result, the majority ($\approx 72\%$) of the eFEDS clusters are covered by the
HSC shape catalog, from which the shear profile is derived.

It is worth mentioning that the biggest benefit of using the latest weak-
lensing catalog instead of the first-year one is the increase in the area
coverage. The number of eFEDS clusters covered by the S16A weak-lensing data
is only $\approx 50\%$ of those currently covered by the S19A. This implies an
improvement of $\sqrt{2}-1\approx 40\%$ in terms of the signal-to-noise ratio
of weak-lensing measurements from S16A to S19A.

Figure 2:  Observed count rate $\eta$ (in the unit of counts per second) and
the redshift $z$ of the 542 eFEDS clusters. The $434$ secure clusters with
$f_{\mathrm{cont}}<0.2$ (see Sect. 2.3) are shown as the circles, while the
rest are marked by red triangles.

###  2.3 The eFEDS cluster sample

The galaxy clusters in the eFEDS are identified by the source-detection
algorithm in the eSASS pipeline. The details of the X-ray source detection are
given in Brunner et al. (2021), to which we refer readers for a complete
description. In what follows, we give a brief summary of the cluster finding
algorithm in the eFEDS.

The source detection of galaxy clusters is performed on the X-ray imaging at
the energy band of $0.2-2.3$ $k$eV. The source detection is performed with the
eSSASusers_201009 software (Brunner et al., 2021). The tasks are based on the
sliding-box algorithm, which searches for sources that are brighter than the
calculated background in the X-ray images. This process is iteratively
performed to obtain a reliable estimation of background noise. A list of X-ray
sources is created after the iterations. Then, the $\beta$-model (Cavaliere &
Fusco-Femiano, 1976) convolved with the derived point spread function is
fitted to each source to derive the X-ray properties, including for example
the detection likelihood $\mathcal{L}_{\mathrm{det}}$, the extent likelihood
$\mathcal{L}_{\mathrm{ext}}$, and the observed count rate $\eta$ at $0.2-2.3$
$k$eV. We note that the observed count rate $\eta$ is estimated using the
total photon counts evaluated by integrating the best-fit model that depends
on the detection radius, which is a characteristic scale to facilitate the
detection of an extended source in the eSASS pipeline. Thus, the observed
count rate $\eta$ is different from that estimated within the cluster radius
$R_{500}$. We estimate the bias111We note that this bias refers to the
difference between the observed count rate (returned by the eSASS pipeline)
and the underlying true count rate enclosed by the cluster radius $R_{500}$.
of the observed count rate $\eta$ in Sect. 4.1. Finally, the thresholds of
$\mathcal{L}_{\mathrm{det}}>5$ and $\mathcal{L}_{\mathrm{ext}}>6$ are applied
to create the catalog of cluster candidates, resulting in 542 systems in total
in the eFEDS. These 542 clusters and their X-ray properties are presented in
Liu et al. (2021a).

The sample of 542 eFEDS clusters are then processed through the optical
confirmation by running the Multi-Component Matched Filter (MCMF) algorithm
(Klein et al., 2018). The MCMF is run on the imaging from both the HSC S20A
data set and the DESI Legacy Imaging Survey (Dey et al., 2019). Full details
of the optical confirmation are given in Klein et al. (2021). Briefly, the
combination of the MCMF run on these two data sets yields, for each cluster,
the photometric redshift estimate $z_{\mathrm{cl}}$, the optical richness
$\lambda$, and the estimator $f_{\mathrm{cont}}$ of the probability that the
cluster appears as a superposition of galaxies along the line of sight by
chance. The quantity of $f_{\mathrm{cont}}$ indicates the level of
contamination in the optical confirmation. For instance, a value of
$f_{\mathrm{cont}}=0.2$ means that there is a $20\%$ chance of detecting a
candidate with the derived richness $\lambda$ at the redshift
$z_{\mathrm{cl}}$ in a random field (see also Klein et al., 2019). A value
higher than this represents an even higher probability of a random
superposition (Klein et al., 2021). A cut on $f_{\mathrm{cont}}$ is equivalent
to a reduction in the cluster contamination in the initial catalog, under the
assumption that the contamination from X-rays and the optical is uncorrelated
(see also the discussion at length in Grandis et al., 2020). The contamination
from point sources in the initial catalog of eFEDS clusters is suggested to be
at a level of $\approx 20\%$ using simulations (Liu et al., 2021a). Applying
the cut of $f_{\mathrm{cont}}<0.2$ results in a deduction of $18\%$ (or $99$
clusters) in the sample size, which is consistent with the expectation from
the simulation. The residual contamination from nonextended sources is
suggested to be at a level of $\approx 4\%$ with the cut of
$f_{\mathrm{cont}}<0.2$ (Klein et al., 2019), as suggested by running the MCMF
algorithm on the sample of ROSAT clusters (Boller et al., 2016). In this work,
we apply the cut of $f_{\mathrm{cont}}<0.2$ to the parent sample to construct
a secure sample of $434$ clusters, which is used to study their scaling
relations. The sample in the plane of the observed count rate $\eta$ and
redshift $z$ is presented in Fig. 2.

In addition to the X-ray scaling relations, we also perform a modeling of the
richness-to-mass-and-redshift ($\lambda$–$M_{500}$–$z$) relation for the eFEDS
clusters following the exact methodology described in Sect. 5.3. The results
are shown in Appendix A.

###  2.4 The X-ray measurements from the eFEDS

In this work, we use five X-ray quantities as follow-up observables, which are
the soft-band luminosity $L_{\mathrm{X}}$ at the energy band of $0.5-2k$eV,
the bolometric luminosity $L_{\mathrm{b}}$, the emission-weighted temperature
$T_{\mathrm{X}}$, the gas mass $M_{\mathrm{g}}$, and the mass proxy
$Y_{\mathrm{X}}$ of the ICM for each cluster. These measurements will be fully
presented in a forthcoming paper (Bahar et al., in preparation), and the
details of the X-ray analysis are given in Ghirardini et al. (2021) and Liu et
al. (2021a). We therefore only outline the important steps here.

We stress that these follow-up X-ray measurements ($L_{\mathrm{X}}$,
$L_{\mathrm{b}}$, $T_{\mathrm{X}}$, $M_{\mathrm{g}}$, $Y_{\mathrm{X}}$) are
all estimated within the cluster radius ($<R_{500}$) including the cluster
core, because we do not have enough X-ray photons to derive these quantities
with the core excised (Bahar et al., in preparation). To derive the cluster
radius $R_{500}$, the cluster mass needs to be known before the follow-up
X-ray analysis. In this work, the cluster mass $M_{500}$ is inferred from a
joint modeling of the observed count rate and the weak-lensing shear profile
(see Sect. 5.2). As the cluster mass based on the lensing modeling is still
blinded at the time these X-ray observables are estimated, an early unblinding
exclusively to the _eROSITA_ team is carried out to estimate the “correct”
cluster radius $R_{500}$ (see more details in Sect. 3.2), which is then used
as the aperture size to derive the X-ray observables. After deriving the
follow-up X-ray observables, the _eROSITA_ team delivers them to the HSC team
for the modeling of the X-ray scaling relations. It is important to note that
the modeling of the follow-up X-ray scaling relations is solely led by the HSC
team, which has no exchange of any mass-related information with the _eROSITA_
team after the early unblinding. Therefore, the analysis of the X-ray scaling
relations in this work is still performed blindly to minimize confirmation
bias.

With the exception of $T_{\mathrm{X}}$, the X-ray follow-up measurements
($L_{\mathrm{X}}$, $L_{\mathrm{b}}$, and $M_{\mathrm{g}}$) are derived based
on the fitting to the X-ray imaging extracted in the soft-energy band of
$0.5-2k$eV. The fitting model is centered on the X-ray center of each cluster,
accounting for the sky and instrumental backgrounds. The ICM density profile
is parameterized by the Vikhlinin et al. (2006) model, which is then projected
onto the sky and fitted to the X-ray surface brightness profile. In this way,
the ICM mass $M_{\mathrm{g}}$ is obtained by integrating the best-fit model to
the cluster radius. Meanwhile, the luminosities ($L_{\mathrm{X}}$ and
$L_{\mathrm{b}}$) are obtained from a conversion factor depending on the
temperature and the cluster redshift. The conversion factor is derived at each
step of the Markov Chain Monte Carlo (MCMC) chain of the spectral analysis,
and the full distribution of the conversion factor is used to compute the
luminosity. This effectively marginalizes the derived luminosity over the
uncertainty arising from the spectral analysis, especially from the
temperature measurements. Moreover, this also allows us to estimate the
luminosity even for clusters whose temperatures are poorly constrained. We
note that this conversion factor weakly depends on the temperature; there is
only a factor of two difference between temperatures of 2 $k$eV and 30 $k$eV.

The temperature $T_{\mathrm{X}}$ is obtained by the X-ray spectral analysis
using the XSPEC (Arnaud, 1996). The spectral analysis is challenging in the
eFEDS observation because of the low number of X-ray photons. The net number
of X-ray counts obtained per cluster varies from $\approx 5$ to $\approx
1000$, with the first quartile, median, and the third quartile of $43$, $75$,
and $147$, respectively. Meanwhile, the signal-to-noise ratio (S/N) of such
X-ray data sets extracted at the cluster radius $R_{500}$ spans a range of
$0.02$ and $38.7$ with a median (mean) of $5.2$ ($6.5$).

Given the low photon counts, we are not able to robustly measure the
temperature $T_{\mathrm{X}}$ (thus $Y_{\mathrm{X}}$) for all clusters. In this
work, we therefore restrict the modeling of the $T_{\mathrm{X}}$ and
$Y_{\mathrm{X}}$ scaling relations to a subsample of eFEDS clusters, for which
the temperature can be reliably measured. Specifically, the subsample is
constructed by selecting clusters with a detection likelihood of
$\mathcal{L}_{\mathrm{det}}>50$ and an extent likelihood of
$\mathcal{L}_{\mathrm{ext}}>50$, resulting in $64$ clusters. With these cuts
on $\mathcal{L}_{\mathrm{det}}$ and $\mathcal{L}_{\mathrm{ext}}$, we construct
a highly pure and complete subsample of eFEDS clusters with a median (mean)
S/N of $13$ ($15$). We assess the temperature measurement by comparing the
temperatures obtained within $R_{500}$ to that within a fixed aperture with a
physical radius of $500$ $k$pc, within which we can measure the temperature
with the highest S/N for a maximal number of clusters (Liu et al., 2021a).
Given such a low number of X-ray photons, we do not expect significant
differences between $T_{\mathrm{X}}(R<R_{500})$ and
$T_{\mathrm{X}}(R<500\leavevmode\nobreak\ k\mathrm{pc})$. As a result, we find
that the difference between $T_{\mathrm{X}}(R<R_{500})$ and
$T_{\mathrm{X}}(R<500\leavevmode\nobreak\ k\mathrm{pc})$ for the subsample
with $\mathcal{L}_{\mathrm{det}}>50$ and $\mathcal{L}_{\mathrm{ext}}>50$ has
an inverse-variance-weighted mean of $\left\langle
T_{\mathrm{X}}(R<R_{500})-T_{\mathrm{X}}(R<500\leavevmode\nobreak\
k\mathrm{pc})\right\rangle=\left(0.06\pm 0.36\right)\leavevmode\nobreak\
k\mathrm{eV}$, showing no significant trend as a function of S/N. This is
consistent with Giles et al. (2016), where they found no clear difference
between the temperatures extracted within $R_{500}$ and $300$ $k$pc using the
spectra with low counts. We stress that the modeling of the X-ray scaling
relations throughout this work is based on the temperature measurements
$T_{\mathrm{X}}(R<R_{500})$, instead of
$T_{\mathrm{X}}(R<500\leavevmode\nobreak\ k\mathrm{pc})$.

In summary, the modeling of the $T_{\mathrm{X}}$ and $Y_{\mathrm{X}}$ scaling
relations is performed based on the subsample of $64$ clusters with
$\mathcal{L}_{\mathrm{det}}>50$ and $\mathcal{L}_{\mathrm{ext}}>50$, of which
there are $47$ systems with weak-lensing measurements. Meanwhile, we perform
the modeling of other X-ray follow-up scaling relations ($L_{\mathrm{X}}$,
$L_{\mathrm{b}}$, and $M_{\mathrm{g}}$) based on the full secure sample of
$434$ eFEDS clusters, of which there are $313$ systems with weak-lensing
measurements.

##  3 Weak-lensing analysis

In this section, we first provide a brief overview of the theory behind the
weak-lensing analysis in Sect. 3.1. We then describe the blinding strategy in
our weak-lensing analysis in Sect. 3.2. Finally, we describe the weak-lensing
analysis of the eFEDS clusters using the HSC data sets in Sects. 3.3 to 3.9.
Our goal is to extract the observed shear profiles for individual clusters and
to quantify the systematic errors that will be accounted for in our forward-
modeling approach. We note that we largely follow the weak-lensing analysis as
presented in Umetsu et al. (2020) and Miyatake et al. (2019), to which we
refer readers for more details.

###  3.1 Weak-lensing basic theory

A brief overview of weak gravitational lensing with emphasis on applications
to galaxy clusters is given in this section. We refer readers to Bartelmann &
Schneider (2001), Hoekstra et al. (2013), and Umetsu (2020) for more details.

Cosmic structures deflect light rays, resulting in effects of gravitational
lensing. Galaxy clusters in the limit of the thin lens approximation act as a
single lens embedded in a homogeneous universe, such that all sources behind
clusters are lensed. The lensing strength of a source at redshift
$z_{\mathrm{s}}$ arising from a galaxy cluster at redshift $z_{\mathrm{cl}}$
depends on the distances between the cluster, the observer, and the source.
This is characterized by the critical surface density $\Sigma_{\mathrm{crit}}$
defined as

$\Sigma_{\mathrm{crit}}\left(z_{\mathrm{cl}},z_{\mathrm{s}}\right)=\frac{c^{2}}{4\pi
G}\frac{1}{D_{\mathrm{l}}\left(z_{\mathrm{cl}}\right)\beta\left(z_{\mathrm{cl}},z_{\mathrm{s}}\right)}\,,$
(1)

in which $G$ is Newton’s constant, $c$ is the speed of light,
$D_{\mathrm{l}}\left(z_{\mathrm{cl}}\right)$ is the angular diameter distance
of the cluster, and $\beta$ is the lensing efficiency of the source,

$\beta\left(z_{\mathrm{cl}},z_{\mathrm{s}}\right)=\begin{cases}\frac{D_{\mathrm{l},\mathrm{s}}\left(z_{\mathrm{cl}},z_{\mathrm{s}}\right)}{D_{\mathrm{s}}\left(z_{\mathrm{s}}\right)}&\leavevmode\nobreak\
\mathrm{if}\leavevmode\nobreak\ z_{\mathrm{s}}>z_{\mathrm{cl}}\\\
0&\leavevmode\nobreak\ \mathrm{if}\leavevmode\nobreak\ z_{\mathrm{s}}\leq
z_{\mathrm{cl}}\end{cases}\,,$ (2)

where $D_{\mathrm{l},\mathrm{s}}\left(z_{\mathrm{cl}},z_{\mathrm{s}}\right)$
and $D_{\mathrm{s}}\left(z_{\mathrm{s}}\right)$ are the angular diameter
distances of the cluster–source and observer–source pairs, respectively.

In terms of weak-shear effects, gravitational lensing arising from clusters
distorts the imaging of background sources, resulting in a coherent distortion
in the tangential direction around the centers of the clusters. By measuring
the shape of background sources, this effect can be quantified statistically
by the quantity of the “reduced shear” $g_{+}$, defined as

$g_{+}=\frac{\gamma_{+}}{1-\kappa}\,,$ (3)

where $\gamma_{+}$ and $\kappa$ are the tangential shear and convergence of
the cluster, respectively. In an azimuthal average with respect to a cluster
center, the tangential shear $\gamma_{+}\left(R\right)$ at the projected
radius $R$ describes the differential surface mass density
$\Delta{\Sigma}_{\mathrm{m}}$ of the cluster at $R$ with respect to the
critical surface density inferred from the source, namely,

$\gamma_{+}\left(R\right)=\frac{\Delta{\Sigma}_{\mathrm{m}}\left(R\right)}{\Sigma_{\mathrm{crit}}}\,.$
(4)

Meanwhile, the convergence $\kappa\left(R\right)$ is the ratio of the
projected surface mass density ${\Sigma}_{\mathrm{m}}\left(R\right)$ to the
critical surface density, as

$\kappa\left(R\right)=\frac{{\Sigma}_{\mathrm{m}}\left(R\right)}{\Sigma_{\mathrm{crit}}}\,.$
(5)

We note that
$\Delta{\Sigma}_{\mathrm{m}}\left(R\right)={\Sigma}_{\mathrm{m}}\left(<R\right)-{\Sigma}_{\mathrm{m}}\left(R\right)$.

The shear profile, $g_{+}\left(R\right)$, can then be statistically obtained
by measuring the azimuthal average of tangential ellipticities over a large
sample of background sources around the center of clusters (see Sect. 3.4). To
infer the mass of a cluster from the observed shear profile, the redshift of
sources must be known to compute the critical surface density
$\Sigma_{\mathrm{crit}}$ to interpret $\gamma_{+}$ and $\kappa$ (see Sect. 3.9
for the modeling of $\Sigma_{\mathrm{crit}}$).

###  3.2 Blinding analysis

We carry out the weak-lensing analysis in a blind fashion to avoid
confirmation bias. The blinding strategy is decided by the HSC weak-lensing
team and has been widely used in cosmological analyses, such as those of
Hikage et al. (2019) and Hamana et al. (2020). A full description of the
blinding strategy is presented in Hikage et al. (2019), to which we refer the
reader for details. We summarize the blinding procedure below.

A two-level blinding strategy is used. The first level of blinding is at the
level of the catalog, while the second is at the level of the analysis itself.
In terms of the catalog level, we blind the true estimate of source shears by
perturbing the multiplicative bias, as

$m_{\mathrm{cat},i}=m_{\mathrm{true}}+\mathrm{d}m_{1,i}+\mathrm{d}m_{2,i}\,,$
(6)

where $m_{\mathrm{cat},i}$ is the blinded multiplicative bias in the $i$-th
shape catalog, $m_{\mathrm{true}}$ is the true estimate of the multiplicative
bias, and $\mathrm{d}m_{1,i}$ and $\mathrm{d}m_{2,i}$ are two random variables
for blinding. While many HSC analyses are being carried out in parallel to
this work, we prepare three blinded catalogs ($i\in\\{0,1,2\\}$) independently
from other HSC analyses. The random variables of $\mathrm{d}m_{1,i}$ and
$\mathrm{d}m_{2,i}$ are encrypted and are different among these three blinded
catalogs; the HSC analysis teams are also different. The term
$\mathrm{d}m_{1,i}$ can only be decrypted by the leader of the analysis, and
is removed before performing the analysis. This term is needed to avoid an
accidental unblinding by comparing the blinded catalogs of different analysis
teams. The term $\mathrm{d}m_{2,i}$ can only be decrypted by a “blinder-in-
chief”, once we are ready to unblind the analysis. The blinder-in-chief is not
involved in the analysis and is not aware of the values of $\mathrm{d}m_{2,i}$
until unblinding. The value of $\mathrm{d}m_{2,i}$ varies from $-0.1$ to $0.1$
randomly. Only one blinded catalog among the three has $\mathrm{d}m_{2,i}=0$,
as the true catalog. The analysis team has to run the identical analysis on
these three blinded catalogs, which is a computationally expensive element of
our blinding strategy. However, such a strategy comes with the advantage that
an end-to-end rerun is not needed, once these catalogs are unblinded. This
completes the catalog-level of blinding.

For the analysis-level blinding, we never compare the results among three
blinded catalogs until unblinding. This avoids an automatic unblinding arising
from an accidental comparison of the blinded catalogs with values of
$\mathrm{d}m_{2,i}$ all close to each other. Moreover, for eFEDS clusters with
available optical counterparts in the HSC survey, we never compare the mass
estimates or shear profiles of common clusters between these two surveys.
Also, we never compare the mass distribution of eFEDS clusters with those from
other cluster surveys until unblinding. These strategies keep our weak-lensing
analysis blinded.

It is worth mentioning that the first-year shape catalog of the HSC survey was
already publicly available by the time we initialized this work. However, we
stress that we neither run the analysis using the public catalog, nor make any
kind of comparisons between the blinded and public catalogs.

We unblind our analysis once the following criteria are met:

  * •

The analysis codes pass validation tests using mock catalogs of at least ten
times the size of the eFEDS sample, ensuring that the codes can recover the
input parameters within statistical uncertainties.

  * •

The posteriors of parameters from the sampled chains are converged.

  * •

The best-fit models provide a good description of observed shear profiles.

  * •

The systematic errors arising from the miscentering on the final results are
quantified.

  * •

The value $\mathrm{d}m_{1,i}$ has been correctly subtracted from the blinded
catalogs.

  * •

The selection bias due to resolution cuts and magnitude cuts has been
appropriately applied (see Hikage et al., 2019).

It is important to stress that the measurement of the follow-up X-ray
observables, i.e., $L_{\mathrm{X}}$, $L_{\mathrm{b}}$, $M_{\mathrm{g}}$,
$T_{\mathrm{X}}$, and $Y_{\mathrm{X}}$, is not performed in a blind analysis.
This is because the correct $M_{500}$ must be known to extract the X-ray
spectra enclosed by the corresponding cluster radius $R_{500}$. Ideally, the
X-ray spectral analysis should be performed in terms of profiles with fine
radial bins around the cluster center, and we calculate the X-ray measurement
at the corresponding cluster radius given a cluster mass at each iteration of
the likelihood exploration. However, the X-ray analysis, especially the
extraction of spectra observed by _eROSITA_ , is extremely time-consuming and
beyond what we can afford given the timescale. We therefore perform an early
unblinding on the X-ray analysis in this work. Specifically, the weak-lensing
mass calibration (in Sect. 5.2), which only involves the X-ray observable,
namely the count rate $\eta$, is performed first in a completely blinded way
as described above. After the weak-lensing mass calibration passes the
unblinding criteria, three blinded sets of cluster masses are produced based
on the three blinded shape catalogs. We then ask the blinder-in-chief to
unblind the cluster mass privately and only to the _eROSITA_ team, such that
they can extract the X-ray spectra within the “correct” $R_{500}$. A set of
the X-ray observables extracted within the “correct” radius is then delivered
to the HSC lensing team to perform the blind modeling of X-ray observable-to-
mass-and-redshift relations ($L_{\mathrm{X}}$, $L_{\mathrm{b}}$,
$T_{\mathrm{X}}$, $M_{\mathrm{g}}$, and $Y_{\mathrm{X}}$, see Sect. 5.3) for
each blinded catalog. After the early unblinding, we strictly forbid the
exchange of information between the HSC and _eROSITA_ teams by any means until
the official unblinding takes place. As the modeling of the scaling relations
is performed solely by the lensing team, this ensures that the analysis in
this work is still carried out blindly with an adjustment to minimize
excessive measurements in X-rays.

##### Post-unblinding analysis

Despite the careful treatment in the blinding analysis, we made changes to
correct the discovered errors in a post-unblinding analysis. We refer readers
to Appendix B for more details. We find no significant changes after
correcting the errors, and therefore the interpretations of this work are not
affected. We note that the post-unblinding analysis only affected the modeling
of the follow-up X-ray scaling relations but not the weak-lensing mass
calibration.

###  3.3 Source selection

We use the same source selection as in Umetsu et al. (2020) to select source
galaxies behind each cluster. Specifically, for a cluster at the redshift
$z_{\mathrm{cl}}$, the source selection is performed using the observed
probability distribution of redshift, $P\left(z\right)$, such that a galaxy is
considered as a source galaxy if

$0.98<\int_{{z}_{\mathrm{min}}}^{\infty}P\left(z\right)\mathrm{d}z\leavevmode\nobreak\
\leavevmode\nobreak\ \mathrm{and}\leavevmode\nobreak\ \leavevmode\nobreak\
{z}_{\mathrm{MC}}<2.5\,,$ (7)

where we set ${z}_{\mathrm{min}}=z_{\mathrm{cl}}+0.2$, and ${z}_{\mathrm{MC}}$
is the point estimator randomly sampled from $P\left(z\right)$. The
contamination to lensing signals due to cluster members is suggested to be at
the level of a few percent using our $P\left(z\right)$-based source selection,
as fully explored and quantified in Medezinski et al. (2018b). An independent
examination of the cluster contamination in this work verifies that the level
of contamination is at the level of $\lesssim 3\%$, on average (for more
details, see Sect. 3.6).

In this work, we use the photometric redshift estimated from the machine-
learning-based code, DEmP (Hsieh & Yee, 2014), which has been widely used in
obtaining not only redshifts but also other physical quantities, such as
stellar mass (Lin et al., 2017). Our source selection leads to source
densities of $\approx 15$, $\approx 11$, $\approx 6$, $\approx 2$, and
$\approx 0.3$ galaxies per square arcmin for a cluster at a redshift of
$z_{\mathrm{cl}}=0.05$, $z_{\mathrm{cl}}=0.25$, $z_{\mathrm{cl}}=0.50$,
$z_{\mathrm{cl}}=0.78$, and $z_{\mathrm{cl}}=1.1$, respectively. At the median
redshift of eFEDS clusters, $z_{\mathrm{piv}}=0.35$, the source density
reaches $\approx 9.6$ galaxies per square arcmin.

###  3.4 Tangential shear profiles

In what follows, we detail the procedure used to extract the observed
tangential shear profile, as the lensing observable, for each cluster. We
largely follow the analysis presented in Umetsu et al. (2020) with one
distinct difference: we measure the dimensionless tangential shear profile
$g_{+}$, as the quantity that can be directly observed in lensing observations
without knowing the redshift of sources, while Umetsu et al. (2020) measured
the differential surface mass density $\Delta{\Sigma}_{\mathrm{m}}$. The
latter require prior knowledge of redshift and the distance-to-redshift
relation with assumed cosmological parameters to convert $g_{+}$ to
$\Delta{\Sigma}_{\mathrm{m}}$. This results in a difficulty in the lensing
modeling in a cosmological analysis, where the cosmological parameters are
varied and hence change the observable $\Delta{\Sigma}_{\mathrm{m}}$.
Conversely, the tangential shear profile $g_{+}$ is directly observed and is
invariant among different cosmology. Therefore, the use of $g_{+}$ as a
lensing observable enables a clean forward-modeling approach, which will be
needed in a cosmological analysis in the future.

For each cluster at redshift $z_{\mathrm{cl}}$, we derive the tangential shear
profile $g_{+}$ in angular bins of clustercentric radius $\theta_{i}$,

$g_{+}\left(\theta_{i}\right)=\frac{1}{2\mathfrak{R(\theta_{i})}\left(1+\mathfrak{K(\theta_{i})}\right)}\frac{\sum\limits_{s\in\theta_{i}}{w}_{s}{e_{+}}_{,s}}{\sum\limits_{s\in\theta_{i}}{w}_{s}}\,,$
(8)

where the subscript $s$ runs over the source sample that is selected according
to the description in Sect. 3.3; ${w}_{s}$ and ${e_{+}}_{,s}$ are the lensing
weight and tangential ellipticity of the $s$-th source galaxy, respectively;
$\mathfrak{R}(\theta_{i})$ is the shear response (see also Mandelbaum et al.,
2005); the factor $\left(1+\mathfrak{K(\theta_{i})}\right)$ is the correction
accounting for multiplicative shear bias, which is calibrated against the
image simulations (Mandelbaum et al., 2018c, b).

The tangential ellipticity $e_{+}$ of a source galaxy around the cluster
center is calculated as

$e_{+}=-\cos(2\phi)e_{1}-\sin(2\phi)e_{2}\,,$ (9)

where $\left(e_{1},e_{2}\right)$ is a two-component estimate of ellipticity
measured in the cartesian coordinate of the sky defined in the HSC survey, and
$\phi$ is the positional angle measured from the Right Ascension direction to
the line connecting the cluster center and the source galaxy. The lensing
weight $w$ is evaluated as

$w=\frac{1}{{\sigma_{\mathrm{e}}}^{2}+{e_{\mathrm{rms}}}^{2}}\,,$ (10)

where $\sigma_{\mathrm{e}}$ and $e_{\mathrm{rms}}$ are the measurement
uncertainty and the root-mean-square estimate of the ellipticity per
component, respectively.

The shear response $\mathfrak{R}(\theta_{i})$ in the angular bin $\theta_{i}$
is calculated as

$\mathfrak{R}(\theta_{i})=1-\frac{\sum\limits_{s\in\theta_{i}}{w}_{s}{{e_{\mathrm{rms}}}_{,s}}^{2}}{\sum\limits_{s\in\theta_{i}}{w}_{s}}\,,$
(11)

where the subscript $s$ runs over the source sample in the angular bin
$\theta_{i}$. The factor $\left(1+\mathfrak{K(\theta_{i})}\right)$ is
calculated as

$1+\mathfrak{K}(\theta_{i})=\frac{\sum\limits_{s\in\theta_{i}}{w}_{s}\left(1+m_{s}\right)}{\sum\limits_{s\in\theta_{i}}{w}_{s}}\,,$
(12)

where $m_{s}$ is the multiplicative bias of the $s$-th source galaxy.

The angular binning in clustercentric radii is defined as follows. For each
cluster, we first use a logarithmic binning between $0.2\leavevmode\nobreak\
h^{-1}\mathrm{Mpc}$ and $3.5\leavevmode\nobreak\ h^{-1}\mathrm{Mpc}$ with ten
steps in physical units. Then, the angular binning is obtained by dividing the
physical radius by the angular diameter distance evaluated at the cluster
redshift in the fiducial cosmology. In this way, we ensure that roughly the
same portion of the cluster profiles at different redshifts is probed. The use
of the angular-binned shear profile $g_{+}\left(\theta\right)$ enables us to
self-consistently calculate the cluster mass profile at any redshift-inferred
distance in a forward-modeling approach where the cosmological parameters vary
(see Sect. 3.9).

For each cluster, we additionally compute the B-mode component of the shear
profile $g_{\times}\left(\theta\right)$ according to Eq. (8), after replacing
$e_{+}$ with $e_{\times}$, which is defined as

$e_{\times}=\sin(2\phi)e_{1}-\cos(2\phi)e_{2}\,.$ (13)

The azimuthally averaged cross-shear profile $e_{\times}\left(R\right)$ is
expected to vanish if the signal is due to weak lensing, therefore enabling a
null-test for our lensing signals. We verified that the stacked profile of
$g_{\times}\left(\theta\right)$ in this work is indeed statistically
consistent with zero.

###  3.5 Covariance matrices

For each cluster, we calculate the covariance matrix that will be used to
model the observed $g_{+}$, by following the prescription in Umetsu et al.
(2020). Specifically, the covariance matrix $\mathbb{C}$ comprises a component
characterizing the noise of shape measurements, $\mathbb{C}_{\mathrm{shape}}$,
and another component $\mathbb{C}_{\mathrm{uLSS}}$ accounting for the
uncertainty arising from uncorrelated large-scale structures around the
cluster. That is,

$\mathbb{C}=\mathbb{C}_{\mathrm{shape}}+\mathbb{C}_{\mathrm{uLSS}}\,,$ (14)

in which

${\mathbb{C}_{\mathrm{shape}}}_{i,j}=\sigma_{\mathrm{shape}}^{2}(\theta_{i})\delta_{i,j}\,,$

where $\delta_{i,j}$ is a Kronecker delta function, and
${\sigma_{\mathrm{shape}}}^{2}(\theta_{i})$ in the angular bin $\theta_{i}$ is
evaluated as

$\sigma_{\mathrm{shape}}^{2}(\theta_{i})=\frac{1}{4\mathfrak{R}^{2}(\theta_{i})\left(1+\mathfrak{K}(\theta_{i})\right)^{2}\sum\limits_{s\in\theta_{i}}{w}_{s}}\,.$
(15)

We compute $\mathbb{C}_{\mathrm{uLSS}}$ by following the prescription in
Appendix A of Miyatake et al. (2019) (see also Hoekstra, 2003). Specifically,
for each cluster, we compute the nonlinear matter power spectrum
$P^{\mathrm{NL}}_{\mathrm{m}}$ at the cluster redshift $z_{\mathrm{cl}}$ using
the code CAMB (Lewis et al., 2000) in order to derive the lensing power
spectrum $C_{l}^{\kappa\kappa}$. We use the full redshift distribution of the
sources by stacking the observed redshift distributions of the source sample
in order to infer the lensing weight function in calculating
$C_{l}^{\kappa\kappa}$. Finally, we compute $\mathbb{C}_{\mathrm{uLSS}}$ by
integrating the product of $C_{l}^{\kappa\kappa}$ and two second-order Bessel
functions associated with the radial binning of shear profiles. This
calculation of $\mathbb{C}_{\mathrm{uLSS}}$ is independently repeated for all
clusters using the cosmological parameters fixed to the fiducial cosmology.

We note that the term $\mathbb{C}_{\mathrm{uLSS}}$ depends on the underlying
cosmology. Ideally, the variation in $\mathbb{C}_{\mathrm{uLSS}}$ due to a
change of cosmological parameters needs to be taken into account in each
iteration step of a forward modeling. However, in this work, the covariance
matrix is dominated by the shape noise in the fitting range of radii of
interest. Moreover, we adopt Gaussian priors on the cosmological parameters
with the mean values used in the fiducial cosmology. Hence, the cosmology-
dependence of $\mathbb{C}_{\mathrm{uLSS}}$ is not expected to be a dominant
factor in the final results of our work (see Kodwani et al., 2019). Therefore,
in this work we ignore the cosmological dependence of
$\mathbb{C}_{\mathrm{uLSS}}$ so as to improve the speed of our calculation. We
leave the improvement of a cosmology-dependent $\mathbb{C}_{\mathrm{uLSS}}$ to
future work.

We also note that the intrinsic variation of shear profiles at fixed cluster
mass222This is characterized by the term $C^{\mathrm{int}}$ in Eq. (16) of
Umetsu et al. (2020). —due to for example the triaxiality or concentration of
halos— is accounted for by modeling the intrinsic scatter of the weak-lensing
mass bias (see Sect. 4.2). Therefore, we do not include this term of intrinsic
variation in Eq. (14).

Figure 3:  Comparison of the redshift distributions between cluster and random
fields at low (left panel) and high (middle panel) redshifts. In the left and
middle panels, the color scheme ranging from dark red to blue indicates the
scale of the clustercentric radius from the core to large radii, while the
black dashed lines are the redshift distributions of random fields. Meanwhile,
the vertical gray regions indicate the redshift range of the cluster samples,
while the upper and lower limits ($z_{\mathrm{lolim}}$ and
$z_{\mathrm{uplim}}$) of the redshift interval in Eq. (17) are shown by green
dotted lines. In the right panel, the profiles of cluster contamination in
each subsample are shown with the color scheme according to the sample
binning. The resulting profile of the total sample is in black. The horizontal
gray region indicates the conservative estimate of cluster contamination at a
level of $\leq 3\%$ ($1\sigma$), which is used to calibrate the weak-lensing
mass bias (see Sect. 4.2). In these three panels, the numbers in the
parentheses indicate the total numbers of the clusters used in the binning.

###  3.6 Cluster member contamination

Because of the uncertainty in determining redshifts, a sample of
photometrically selected sources might contain galaxies that are not behind
clusters, are therefore not lensed, and thus do not contribute to lensing
signals. The inclusion of these nonbackground galaxies dilutes lensing
signals, and therefore introduces bias in inferred cluster masses. For the
contamination due to foreground galaxies, this bias can be statistically
accounted for by properly weighting the lensing signal according to the
observed redshift distribution $P\left(z\right)$, with a calibration against
random fields where no clusters are present. However, cluster environments are
extremely biased fields, of which the observed redshift distribution cannot be
directly interpreted by those of random fields. Moreover, the cluster
contamination is expected to be more significant around the center of
clusters, resulting in a radial dependence of bias in lensing signals. This
needs to be accounted for; otherwise significant bias will be introduced in
the inferred mass (Medezinski et al., 2018b).

In this section, we examine the contamination of lensing signals due to
cluster members by following the prescription in McClintock et al. (2019) (see
also Gruen et al., 2014; Melchior et al., 2017; Varga et al., 2019; Chiu et
al., 2020b). Specifically, this method is based on a decomposition of the
observed redshift distributions into a cluster component and the other from
fore- and background (hereafter background, for simplicity). By comparing the
observed $P\left(z\right)$ between the cluster and random fields, the
contamination arising from cluster members can be quantified. We stress that
this methodology has been validated by simulations (Varga et al., 2019),
showing that the underlying true cluster contamination can be recovered based
on the observed $P\left(z\right)$.

It is important to note that the source selection is different between
McClintock et al. (2019) and this work: these latter authors made use of all
galaxies around clusters with available shape measurements, and properly
assigned lensing weights to galaxies according to their observed redshift
distributions, such that foreground galaxies and cluster members did not
contribute to the lensing signals. Then, they accounted for the residual
contamination due to cluster members by quantifying the excess observed
redshift distributions at the cluster redshift with respect to random fields.
In this work, on the other hand, we first minimize the cluster contamination
by employing a stringent source selection based on the observed redshift
distribution (see Sect. 3.3); we then measure the lensing signal of these
secure sources. Consequently, our strategy leads to a clean sample with the
cluster member contamination at the level of only a few percent (Medezinski et
al., 2018b; Umetsu et al., 2020). Importantly, our source selection still
results in high densities of source galaxies thanks to the deep imaging of the
HSC survey.

In this work, we conduct a nonparametric approach to quantify the cluster
contamination by comparing the observed redshift distributions between the
cluster and random fields. In what follows, we express the clustercentric
radius in the physical unit $R$ through the redshift–distance relation in the
fiducial cosmology. We assume that the observed redshift distribution
$P\left(z;R\right)$ of selected source galaxies at the radial bin $R$ can be
written as

$P\left(z;R\right)=\left(1-f_{\mathrm{cl}}\left(R\right)\right)P_{\mathrm{bkg}}\left(z\right)+f_{\mathrm{cl}}\left(R\right)P_{\mathrm{cl}}\left(z;R\right)\,,$
(16)

where $f_{\mathrm{cl}}\left(R\right)$ is the amount of cluster contamination
at the radial bin $R$, $P_{\mathrm{bkg}}$ is the redshift distribution of
galaxies at backgrounds, and $P_{\mathrm{cl}}$ is the redshift distribution of
cluster members that lead to contamination. By construction, we have $0\leq
f_{\mathrm{cl}}\left(R\right)\leq 1$. We can write Eq. (16) in the integral of
redshift,

$\int\mathrm{d}zP\left(z;R\right)=\left(1-f_{\mathrm{cl}}\left(R\right)\right)\int\mathrm{d}zP_{\mathrm{bkg}}\left(z\right)+f_{\mathrm{cl}}\left(R\right)\int\mathrm{d}zP_{\mathrm{cl}}\left(z;R\right)\,,$

and then integrate over the redshift interval $\mathfrak{I}$, where the
cluster component vanishes, and finally arrive at

$f_{\mathrm{cl}}\left(R\right)=1-\frac{\int\limits_{z\in\mathfrak{I}}\mathrm{d}zP\left(z;R\right)}{\int\limits_{z\in\mathfrak{I}}\mathrm{d}zP_{\mathrm{bkg}}\left(z\right)}\,.$
(17)

In this way, the cluster contamination at the radial bin $R$ can be derived
for a given set of observed $P\left(z;R\right)$ and
$P_{\mathrm{bkg}}\left(z\right)$.

In practice, in the interest of more precise measurements, we estimate
$f_{\mathrm{cl}}\left(R\right)$ according to Eq. (17) by stacking clusters in
different redshift and count-rate bins. Namely, the clusters are stacked in
two redshift bins, low-redshift ($z\leq 0.35$) and high-redshift ($z>0.35$),
and two count-rate bins, low-rate ($\eta<0.1$) and high-rate ($\eta>0.1$).
This leads to four subsamples in total based on the observed redshifts and
count rates. We note that we have tried a finer binning in redshift and find
consistent results. For each individual cluster, labeled $j$, we first derive
the weighted probability distribution $P_{j}\left(z;R_{i}\right)$ of redshift
observed at the radial bin $R_{i}$ as

$P_{j}\left(z;R_{i}\right)=\frac{\sum\limits_{s\in
j,R_{i}}{w}_{s}P_{s}(z)}{\sum\limits_{s\in j,R_{i}}{w}_{s}}\,,$ (18)

where $s$ runs over the source galaxies which are located at the radius bin
$R_{i}$ around the cluster $j$. In each subsample of clusters, we then derive
the number-weighted average of redshift distributions observed at the radius
bin $R_{i}$ as

$\left\langle
P\right\rangle\left(z;R_{i}\right)=\frac{\sum\limits_{j=1}^{N_{\mathrm{cl}}}P_{j}\left(z;R_{i}\right)N_{\mathrm{gal},j}\left(R_{i}\right)}{\sum\limits_{j=1}^{N_{\mathrm{cl}}}N_{\mathrm{gal},j}\left(R_{i}\right)}\,,$
(19)

where $j$ runs over the subsample of $N_{\mathrm{cl}}$ clusters, and
$N_{\mathrm{gal},j}\left(R_{i}\right)$ is the number of observed galaxies at
the radius bin $R_{i}$ around the cluster $j$. We then define the redshift
interval $\mathfrak{I}$ as
$\left\\{z|z<{z}_{\mathrm{lolim}}\right\\}\cup\left\\{z|z>{z}_{\mathrm{uplim}}\right\\}$,
where
${z}_{\mathrm{lolim}}=\mathtt{min}\left\\{{z}_{j}\right\\}_{j=1}^{N_{\mathrm{cl}}}-\Delta
z$ ,
${z}_{\mathrm{uplim}}=\mathtt{max}\left\\{{z}_{j}\right\\}_{j=1}^{N_{\mathrm{cl}}}+\Delta
z$, and ${z}_{j}$ is the redshift of the $j$-th cluster. We set $\Delta
z=0.2$, which is consistent with our source selection (see Sect. 3.3). We note
that our results are not sensitive to the current choice of $\Delta z=0.2$; we
verified that the results with $\Delta z=0.4$ are statistically consistent
with those based on $\Delta z=0.2$, although with much larger error bars. In
this way, we can estimate the numerator in Eq. (17) by integrating
$\left\langle P\right\rangle\left(z;R_{i}\right)$ over the interval
$\mathfrak{I}$.

Concomitently, we estimate $P_{\mathrm{bkg}}\left(z\right)$ by repeating Eqs.
(18) and (19) with the same scheme of radial binning and cluster redshift on a
number of $N_{\mathrm{cl}}$ pointing on random fields. This allows us to
calculate the contribution arising from pure background sources in a fully
consistent way as in cluster fields. Finally, we compute
$f_{\mathrm{cl}}\left(R\right)$ in each subsample by following Eq. (17).

The results of cluster contamination are shown in the right panel of Fig. 3,
where the contamination profiles of four subsamples are in blue, orange,
green, and red. These are all statistically consistent with zero
contamination, with a mild exception that the subsample of low-$\eta$ clusters
at high redshift shows large variation but also with large error bars. This is
mainly due to noisy measurements of photo-$z$ for high-redshift sources. We
further repeat the procedure but without binning in $\eta$, resulting in two
subsamples at low ($0<z<0.35$) and high ($0.35<z<1.20$) redshifts; their
contamination profiles are shown in purple and brown, respectively. No clear
contamination is seen. We additionally show the number-weighted average of
redshift distributions (i.e., Eq. (19)) of these two subsamples in the left
and middle panels, where the redder colors represent inner radii, together
with those estimated from random fields in dashed lines. As seen in these
figures, they all show no signs of cluster contamination. This is in perfect
agreement with previous HSC work (Medezinski et al., 2018b, a; Miyatake et
al., 2019; Umetsu et al., 2020), demonstrating a highly pure source sample for
weak-lensing studies.

Based on these results, we conclude that the level of cluster contamination in
this work is at the level of $\lesssim 3\%$ ($1\sigma$), which is a
conservative estimate, shown by the gray shaded region in the right panel of
Fig. 3. This estimation of cluster contamination will be taken into account
when deriving the weak-lensing mass bias (see Sect. 4.2).

###  3.7 Cluster miscentering

In this work, the cluster center is defined as the center of X-ray emissions.
However, X-ray centers do not necessarily represent the true center of the
total mass distribution of galaxy clusters. The miscentering between the
“observed” center and the “true” center causes bias in shear profiles with
respect to perfectly centered clusters (Johnston et al., 2007a, b), especially
at small clustercentric radii. Consequently, the lensing modeling needs to
account for this miscentering effect; otherwise bias in the inferred mass will
be introduced.

It has been shown that miscentering has significant effects on the observed
properties of optically selected clusters (Biesiadzinski et al., 2012; Sehgal
et al., 2013; Rozo & Rykoff, 2014). Therefore, extensive efforts have been
made to quantify the offset between the optical center, usually defined by the
BCG, and other center proxies, such as the peak of projected mass
distributions inferred from lensing (Oguri et al., 2010; Zitrin et al., 2012),
the second brightest cluster galaxy (Hoshino et al., 2015), and the ICM-based
center defined in X-rays (Lin & Mohr, 2004; Mahdavi et al., 2013; Lauer et
al., 2014; Zhang et al., 2019) or at millimeter wavelength (Song et al., 2012;
Saro et al., 2015; Bleem et al., 2020). It is worth mentioning that the center
inferred from weak-lensing mass maps is largely affected by the noise arising
from shape measurements (Dietrich et al., 2012), and therefore it is extremely
difficult to determine the true center of overall cluster mass distributions
from observations alone. We note that it is possible to account for the
modeling bias when centering on the weak-lensing mass peak for clusters whose
lensing signals are detected at high significance (Sommer et al., 2021), which
is unfortunately not the case in our study.

In this work, we gauge the miscentering in the eFEDS clusters using the offset
between the X-ray centers and the BCGs, which are used as the optical centers.
The BCGs are identified by the MCMF as the brightest galaxies that are
consistent with the red-sequence prediction at the cluster redshift. We note
that the complex nature of BCGs may result in incorrect identification of the
cluster center in the optical. For example, the BCGs are suggested to scatter
more from the red sequence than the typical bright galaxies (van der Burg et
al., 2014; Kravtsov et al., 2018), and some BCGs may not be red in the group
scale (Liu et al., 2012). To assess the impact of incorrectly identified
optical centers, we re-derive the offset distribution (see below) by replacing
the optical center with the peak of the galaxy density map for the clusters
whose BCGs are significantly different from both the peak of galaxy density
and X-ray maps. We find that the resulting difference is negligible.

The offset distribution $P(x)$ between the X-ray and optical centers as a
function of the dimensionless radius $x$ is characterized by a composite model
(see also Saro et al., 2015; Zhang et al., 2019):

$P(x)=\left(1-f_{\mathrm{mis}}\right)\times
P_{\mathrm{cen}}\left(x|{\sigma}_{\mathrm{cen}},\alpha\right)+f_{\mathrm{mis}}\times
P_{\mathrm{mis}}\left(x|{\sigma}_{\mathrm{mis}}\right)\,.$ (20)

In Eq. (20), we assume that the distribution $P(x)$ can be decomposed into two
components: one is the centering component characterized by a modified
Rayleigh distribution $P_{\mathrm{cen}}(x|{\sigma}_{\mathrm{cen}},\alpha)$,
which reads

$P_{\mathrm{cen}}(x|{\sigma}_{\mathrm{cen}},\alpha)=\frac{\frac{1}{{\sigma}_{\mathrm{cen}}}\left(\frac{x}{{\sigma}_{\mathrm{cen}}}\right)^{\alpha}\exp\left(-\frac{1}{2}\left(\frac{x}{{\sigma}_{\mathrm{cen}}}\right)^{2}\right)}{2^{\frac{\alpha-1}{2}}\Gamma\left(\frac{\alpha+1}{2}\right)}\,,$
(21)

where $\Gamma$ is the gamma function; and the other is the miscentering
component as modeled by a standard Rayleigh distribution,

$P_{\mathrm{mis}}(x|{\sigma}_{\mathrm{mis}})=\frac{1}{{\sigma}_{\mathrm{mis}}}\left(\frac{x}{{\sigma}_{\mathrm{mis}}}\right)\exp\left(-\frac{1}{2}\left(\frac{x}{{\sigma}_{\mathrm{mis}}}\right)^{2}\right)\,.$
(22)

These two components are weighted by the miscentering fraction
$f_{\mathrm{mis}}$ in Eq. (20). We note that the parameter $\alpha$ is
required to be positive and characterizes the shape of the centering
component: The smaller $\alpha$ is the smaller radius where the peak of
$P_{\mathrm{cen}}(x)$ occurs, and the steeper fall-off toward $x=0$. Moreover,
$P_{\mathrm{cen}}(x)$ reduces to a half-normal distribution and zero, as
$\alpha$ approaches zero and $+\infty$, respectively. Also, $P_{\mathrm{cen}}$
is identical to $P_{\mathrm{mis}}$ if $\alpha=1$.

Figure 4:  Distribution of offsets between the X-ray and optically defined
centers of eFEDS clusters. The blue histogram is the observed offset
distribution in the unit of the scale radius. The black solid line represents
the best-fit model, i.e., Eq. (20), which consists of centering and
miscentering components, as indicated by the red and green dashed lines,
respectively.

Our goal is to derive a universal offset distribution, based on which we model
the miscentering of the shear profile for individual clusters. We define the
dimensionless radius as $x\equiv R_{\mathrm{off}}/r_{\mathrm{s}}$, where
$R_{\mathrm{off}}$ is the offset between the optical and X-ray centers, and
$r_{\mathrm{s}}$ is the scale radius of a Navarro-Frenk-White (hereafter NFW;
Navarro et al., 1997) profile with a concentration parameter evaluated using
the Diemer & Kravtsov (2015) fitting formula for a halo with the pivotal mass
$M_{\mathrm{piv}}=1.4\times 10^{14}h^{-1}\mathrm{M}_{\odot}$ at the pivotal
redshift $z_{\mathrm{piv}}=0.35$. For each cluster, the offset
$R_{\mathrm{off}}$ is calculated first in the physical unit using the
redshift–distance relation inferred in the fiducial cosmology. We then convert
$R_{\mathrm{off}}$ to $x$ using the fixed scale radius $r_{\mathrm{s}}$. In
this way, we effectively assume a universal mass and redshift for each cluster
in deriving the offset distribution $P(x)$, as an ensemble behavior. Later,
when modeling the miscentering in Sect. 3.9, we do vary the scale radius for a
given weak-lensing mass according to the concentration-to-mass relation at
each step of the likelihood exploration. Ideally, different core radii based
on individual cluster masses should be used in deriving the universal offset
distribution, which could be achieved by an iterative process after obtaining
the cluster mass. However, we show below that the resulting posteriors of the
cluster mass $M_{500}$ obtained with and without the modeling of the cluster
core are consistent with each other (see Sect. 6.1), suggesting that our
modeling of the miscentering is sufficiently adequate in this work.

Figure 4 shows the offset distribution (blue histogram) and the best-fit
miscentering model (black solid line) as a function of $x$. The best-fit
parameters for Eq. (20) are

$\left(f_{\mathrm{mis}},{\sigma}_{\mathrm{cen}},\alpha,{\sigma}_{\mathrm{mis}}\right)=\\\
\left(0.54\pm 0.02,0.17\pm 0.01,0.26\pm 0.05,0.61\pm 0.03\right)\,.$ (23)

We stress the following particularities in modeling the miscentering of eFEDS
clusters. First, in this work we statistically correct for the miscentering
for individual clusters assuming a universal distribution, $P(x)$, as directly
determined from the data using Eq. (20). For each cluster, this is done by
taking an average of a data array representing a set of shear profiles with
different offset radii $x$ weighted by the distribution $P(x)$ (see more
details in Sect. 3.9). This approach could not be optimal on the basis of
individual clusters; this is because a shear profile with average miscentering
is not necessarily the best description of each cluster, of which the
miscentering is just one realization from $P(x)$. However, this approach is
statistically correct, because we simultaneously model all observed shear
profiles in an ensemble manner, in which the miscentering distribution follows
the derived $P(x)$. Moreover, we follow the identical modeling of the
miscentering, including the distribution $P(x)$, in calibrating the weak-
lensing mass bias using simulations (see Sect. 4.2), which ensures that we
correctly infer the underlying true cluster mass in this framework.

Second, we determine the distribution $P(x)$ as the offset between the X-ray
center and the BCG, while the required knowledge for modeling the miscentering
is the offset between the X-ray and “true” centers of a cluster potential.
That is, the distribution $P(x)$ contains the intrinsic miscentering of both
the X-ray center and the BCG with respect to the true center, as well as their
respective measurement uncertainties at both wavelengths. Among these, the
dominant component is the intrinsic miscentering of the BCG with respect to
the true center, because (1) the X-ray center is considered as a much better
tracer of the total potential than other proxies (Lin & Mohr, 2004; Mahdavi et
al., 2013; Lauer et al., 2014; Zhang et al., 2019) and (2) the uncertainty on
the measurement of the optical center is much smaller than the characteristic
miscentering radius. Therefore, the miscentering effect estimated in this work
merely serves as an upper limit on X-ray miscentering. It was suggested that
X-ray centers could bias the mass estimate even when the modeling of the
cluster core ($R<500$ $k$pc) is excluded (Schrabback et al., 2021). To
accurately determine the intrinsic miscentering of eFEDS clusters, the most
promising way is to utilize an end-to-end simulation suite, in which we
compare the true center of halos with the X-ray center identified by the same
cluster-finding algorithm. A dedicated effort using simulations to determine
the miscentering effect for eFEDS clusters is warranted, but we leave this for
a future analysis.

To quantify the systematic errors arising from the modeling of miscentering,
we compare the final results with and without the modeling of shear profiles
in the cluster core ($R<0.5h^{-1}\mathrm{Mpc}$), that is including or
excluding the three innermost radial bins, respectively. We find that the
miscentering effect is a subdominant factor in this work (see Sect. 7).

Figure 5:  Photo-$z$ bias $b_{z}(z_{\mathrm{cl}})$ as a function of cluster
redshift $z_{\mathrm{cl}}$. The result derived from the first-year HSC data is
shown by the red crosses, which are re-scaled to the three-year HSC data shown
by the black points based on the difference in the photo-$z$ between the data
sets (see Sect. 3.8). We also show a reference at the level of $1\%$ by the
gray shaded region. The resulting $b_{z}(z_{\mathrm{cl}})$ is applied to the
calibration of weak-lensing mass bias for eFEDS clusters (see Sect. 4.2).

###  3.8 Photometric redshift bias

The cluster mass is estimated from lensing signals with photo-$z$, which is
used to calculate the critical surface density $\Sigma_{\mathrm{crit}}$ in
converting $g_{+}$ to $\Delta{\Sigma}_{\mathrm{m}}$ (see Sect. 3.1).
Therefore, a bias in photo-$z$ would ultimately bias the inferred cluster
mass.

Given an observed shear profile, the bias $b_{z}(z_{\mathrm{cl}})$ in the
inferred mass of a cluster at the redshift $z_{\mathrm{cl}}$ due to the bias
in photo-$z$ can be expressed in terms of the lensing efficiency $\beta$, as

$\frac{{\beta}_{\mathrm{true}}}{{\beta}_{\mathrm{obs}}}\left(z_{\mathrm{cl}}\right)\equiv
1+b_{z}(z_{\mathrm{cl}})\,,$ (24)

where ${\beta}_{\mathrm{true}}$ and ${\beta}_{\mathrm{obs}}$ are the lensing
efficiency inferred by the true and observed redshifts, respectively. In an
ideal case, ${\beta}_{\mathrm{true}}$ should be estimated from a sample of
galaxies with redshifts secured by spectroscopic observations. Moreover, the
distributions of the properties of this sample (e.g., of galaxy magnitudes or
colors) should be consistent with those of the source sample, and should also
be independent of the photo-$z$ calibration. However, in practice this sample
is difficult to obtain, especially for our source sample with the limiting
depth of $i<24.5$ mag in the HSC survey. Therefore, we make use of the galaxy
sample observed in the COSMOS field (Ilbert et al., 2009; Laigle et al., 2016)
with high-quality photometric redshifts to assess the photo-$z$ bias in this
work.

We follow the procedure detailed in Miyatake et al. (2019) to quantify the
photo-$z$ bias $b_{z}(z_{\mathrm{cl}})$. In what follows, we summarize the
steps. First, the photometric redshifts ${z}_{\mathrm{COSMOS30}}$ estimated by
the 30-band COSMOS photometry (Laigle et al., 2016) are used to represent the
true redshifts. Second, as the COSMOS field is also observed by the HSC
survey, we additionally prepare a catalog of the COSMOS field that is
processed in the same configuration as the HSC WIDE layer, resulting in a HSC
WIDE-depth version of the COSMOS catalog333We note that we only use a subset
of galaxies in the COSMOS field that were not used to train the photo-$z$ of
the source sample in the HSC survey. This HSC WIDE-depth version of the COSMOS
catalog can be acquired at https://hsc-
release.mtk.nao.ac.jp/doc/index.php/s17a-wide-cosmos/. We then estimate the
photometric redshift of the galaxies in this HSC WIDE-depth COSMOS catalog in
the same way as in the HSC survey. This ensures the homogeneity of the
photo-$z$ between the COSMOS reference catalog and the HSC data. In this way,
for each galaxy in the COSMOS field, we can assess the 30-band COSMOS
photo-$z$ ${z}_{\mathrm{COSMOS30}}$ as the true redshift, and the photometric
redshift distribution $P_{\mathrm{HSC}}\left(z\right)$ estimated by the HSC
survey. Moreover, we can select the galaxy sample in the COSMOS field in the
same way as our source selection using the photometric redshift (see Sect.
3.3).

Next, we employ a re-weighting technique (for more details and some caveats of
this method, see Bonnett et al., 2016; Gruen & Brimioulle, 2017; Hikage et
al., 2019), such that the observed properties of the galaxy population in the
COSMOS field match those of the HSC source sample. Briefly, the HSC source
sample is classified into cells of a self-organizing map (SOM; Masters et al.,
2015) based on the $i$-band magnitude and four colors. We then classify the
COSMOS galaxies according to this SOM, and compute their new weights
$w_{\mathrm{SOM}}$, such that the weighted distributions of the properties
match those of the HSC sources (see also Medezinski et al., 2018b). Finally,
we include the new weights $w_{\mathrm{SOM}}$ into the lensing weights $w$ in
calculating the weighted distribution of redshift to infer the lensing
efficiency. In practice, for a cluster at redshift $z_{\mathrm{cl}}$, we
select the source sample, labeled $\mathsf{S}$, in the COSMOS field according
to the description in Sect. 3.3, and compute

${\beta}_{\mathrm{true}/\mathrm{obs}}(z_{\mathrm{cl}})=\int\beta\left(z_{\mathrm{cl}},z\right)P_{\mathrm{true}/\mathrm{obs}}\left({z}\right)\mathrm{d}{z}\,,$
(25)

where $P_{\mathrm{true}}(z)$ is the distribution of true redshift
${z}_{\mathrm{COSMOS30}}$ in the source sample $\mathsf{S}$, and

$P_{\mathrm{obs}}\left(z\right)=\frac{\sum\limits_{s\in\mathsf{S}}{w_{\mathrm{SOM}}}_{,s}{w}_{s}{P_{\mathrm{HSC}}}_{,s}\left(z\right)}{\sum\limits_{s\in\mathsf{S}}{w_{\mathrm{SOM}}}_{,s}{w}_{s}}\,.$
(26)

Consequently, we compute the photo-$z$ bias $b_{z}(z_{\mathrm{cl}})$ following
Eq. (24).

In practice, there is one additional factor that needs to be considered in
this work: the weighting factor $w_{\mathrm{SOM}}$ in the COSMOS catalog is
specifically tuned to match the property of the sources observed in the first-
year weak-lensing data (S16A), while the three-year (S19A) data are used in
this work. Therefore, the difference in the photo-$z$ between the first- and
three-year data needs to be accounted for when calculating the photo-$z$ bias.
In this work, the difference in the photo-$z$ is taken into account by
multiplying the ratio of the lensing efficiency between the S19 and S16 data
sets. Specifically, the ratio is calculated as

$R_{\beta}(z_{\mathrm{cl}})=\frac{\left\langle\beta_{\mathrm{S16A}}(z_{\mathrm{cl}})\right\rangle}{\left\langle\beta_{\mathrm{S19A}}(z_{\mathrm{cl}})\right\rangle}\,,$

where the $\beta_{\mathrm{S16A}}$ and $\beta_{\mathrm{S19A}}$ are the lensing
efficiency calculated using the photo-$z$ from the S16A and S19A data sets,
respectively, and the bracket $\left\langle\cdot\right\rangle$ stands for the
mean for a selected sample of sources $\mathsf{S}$ given a cluster redshift
$z_{\mathrm{cl}}$. In this way, the final photo-$z$ bias $b_{z,\mathrm{S19A}}$
used in this work is that from the S16A data set multiplying a re-scaling
factor, namely

$1+b_{z,\mathrm{S19A}}(z_{\mathrm{cl}})=R_{\beta}(z_{\mathrm{cl}})\times\left(1+b_{z,\mathrm{S16A}}(z_{\mathrm{cl}})\right)\,.$
(27)

Figure 5 shows the resulting $b_{z}(z_{\mathrm{cl}})$ as a function of cluster
redshift $z_{\mathrm{cl}}$. As seen, the photo-$z$ bias on the cluster mass is
estimated to be $\lesssim 2\%$ for clusters at redshift
$z_{\mathrm{cl}}\lesssim 0.7$, and becomes $\lesssim 6\%$ for
$z_{\mathrm{cl}}\gtrsim 0.8$. There is no clear difference in the photo-$z$
bias between the S16A and S19A data sets, suggesting that their photo-$z$
performances are consistent with each other. This is in close agreement with
previous HSC studies, in which the mean photo-$z$ bias was estimated to be
$\approx 2\%$ and $\approx 0.9\%$ for clusters selected in the ACTpol
(Miyazaki et al., 2018) and XXL (Umetsu et al., 2020) surveys, respectively.
We note that this resulting photo-$z$ bias
$b_{z,\mathrm{S19A}}(z_{\mathrm{cl}})$ is included in deriving the weak-
lensing mass bias for eFEDS clusters in Sect. 4.2.

Figure 6:  Bias $b_{\eta}$, which is defined as the ratio of the observed
count rate $\eta$ to the true count rate $\hat{}\eta$, from the simulation and
the best-fit model (middle and left panels, respectively). These are presented
as a function of true cluster mass $M_{500}$ at each redshift, color coded
according to $\ln\left(b_{\eta}\right)$. A clear gradient depending on both
mass and redshift is seen. The best-fit model, as described in Eq. (40),
provides a good description of the simulation, as suggested by their
difference as shown in the right panel.

###  3.9 Modeling of shear profiles

In this section, we describe the modeling of the observed shear profile of
individual clusters, which is jointly fitted in the likelihood described in
Sect. 5.

For an eFEDS cluster at redshift $z_{\mathrm{cl}}$, there are three input
observables used in the modeling:

  * •

the observed shear profile ${g_{+}}\left(\theta\right)$ as a function of
angular radius $\theta$, estimated in Eq. (8),

  * •

the lensing covariance matrix $\mathbb{C}$, derived as Eq. (14), and

  * •

the observed redshift distribution $P\left(z\right)$, estimated in Eq. (34).

Given a weak-lensing mass $M_{\mathrm{WL}}$ and a set of parameters
$\mathbf{p}_{\mathrm{mis}}=\\{f_{\mathrm{mis}},{\sigma}_{\mathrm{mis}}\\}$
describing the miscentering distribution, together with a set of cosmological
parameters $\mathbf{p}_{\mathrm{c}}$ that are used to calculate the redshift-
inferred distance, the observed shear profile of an eFEDS cluster is modeled
as (Seitz & Schneider, 1997)

$\begin{aligned}
g_{+}^{\mathrm{mod}}\left(\theta|M_{\mathrm{WL}},\mathbf{p}_{\mathrm{mis}},\mathbf{p}_{\mathrm{c}}\right)=\frac{{\gamma}_{\mathrm{mod}}\left(\theta|M_{\mathrm{WL}},\mathbf{p}_{\mathrm{mis}},\mathbf{p}_{\mathrm{c}}\right)}{1-{\kappa}_{\mathrm{mod}}\left(\theta|M_{\mathrm{WL}},\mathbf{p}_{\mathrm{mis}},\mathbf{p}_{\mathrm{c}}\right)}\times\\\
\left(1+{\kappa}_{\mathrm{mod}}\left(\theta|M_{\mathrm{WL}},\mathbf{p}_{\mathrm{mis}},\mathbf{p}_{\mathrm{c}}\right)\left(\frac{\left\langle{\beta}^{2}\left(\mathbf{p}_{\mathrm{c}}\right)\right\rangle}{\left\langle\beta\left(\mathbf{p}_{\mathrm{c}}\right)\right\rangle^{2}}-1\right)\right)\end{aligned}\,,$
(28)

in which

$\left\langle\beta\left(\mathbf{p}_{\mathrm{c}}\right)\right\rangle=\int\beta\left(z_{\mathrm{cl}},z|\mathbf{p}_{\mathrm{c}}\right)P\left(z\right)\mathrm{d}z\,,$
(29)
$\left\langle{\beta}^{2}\left(\mathbf{p}_{\mathrm{c}}\right)\right\rangle=\int{\beta}^{2}\left(z_{\mathrm{cl}},z|\mathbf{p}_{\mathrm{c}}\right)P\left(z\right)\mathrm{d}z\,,$
(30)

and the profiles ${\kappa}_{\mathrm{mod}}$ and ${\gamma}_{\mathrm{mod}}$ are
calculated as

${\kappa}_{\mathrm{mod}}\left(\theta|M_{\mathrm{WL}},\mathbf{p}_{\mathrm{mis}},\mathbf{p}_{\mathrm{c}}\right)=\frac{{\Sigma}_{\mathrm{m}}^{\mathrm{mod}}\left(\theta|M_{\mathrm{WL}},\mathbf{p}_{\mathrm{mis}},\mathbf{p}_{\mathrm{c}}\right)}{{\Sigma_{\mathrm{crit}}}\left(\mathbf{p}_{\mathrm{c}}\right)}\,,$
(31)

and

${\gamma}_{\mathrm{mod}}\left(\theta|M_{\mathrm{WL}},\mathbf{p}_{\mathrm{mis}},\mathbf{p}_{\mathrm{c}}\right)=\frac{\Delta{\Sigma}_{\mathrm{m}}^{\mathrm{mod}}\left(\theta|M_{\mathrm{WL}},\mathbf{p}_{\mathrm{mis}},\mathbf{p}_{\mathrm{c}}\right)}{{\Sigma_{\mathrm{crit}}}\left(\mathbf{p}_{\mathrm{c}}\right)}\,,$
(32)

where

$\Sigma_{\mathrm{crit}}\left(\mathbf{p}_{\mathrm{c}}\right)=\frac{c^{2}}{4\pi
G}\frac{1}{D_{\mathrm{l}}\left(z_{\mathrm{cl}}|\mathbf{p}_{\mathrm{c}}\right)\left\langle\beta\right\rangle}\,.$
(33)

Given that our source selection leads to a highly pure sample with cluster
contamination that is consistent with zero without a clear radial trend (see
Fig. 3), we use a radially independent redshift distribution to calculate the
critical surface density for all radii. The observed redshift distribution
$P\left(z\right)$ used to calculate Eqs. (29) and (30) is estimated as

$P\left(z\right)=\frac{\sum\limits_{s\in j}{w}_{s}P_{s}(z)}{\sum\limits_{s\in
j}{w}_{s}}\,,$ (34)

where the index $s$ runs over the sources with $R<3.5\leavevmode\nobreak\
h^{-1}\mathrm{Mpc}$ in the fiducial cosmology.

Given the set of cosmological parameters $\mathbf{p}_{\mathrm{c}}$, the
calculation of ${\Sigma}_{\mathrm{m}}^{\mathrm{mod}}$ and
$\Delta{\Sigma}_{\mathrm{m}}^{\mathrm{mod}}$ is performed in the radial bins
in physical units by converting $\theta$ to $R$, as

$R=D_{\mathrm{l}}\left(z_{\mathrm{cl}}|\mathbf{p}_{\mathrm{c}}\right)\times\theta\,.$
(35)

This allows us to self-consistently compute the lensing profiles in every step
of the likelihood exploration, given the variation of cosmological parameters.

The miscentering is accounted for in ${\Sigma}_{\mathrm{m}}^{\mathrm{mod}}$
(thus in $\Delta{\Sigma}_{\mathrm{m}}^{\mathrm{mod}}$) as a composite model
comprising a profile of a perfectly centered component and the profile of a
miscentered component, weighted by a relative normalization associated with
$f_{\mathrm{mis}}$. That is,

$\displaystyle\Delta{\Sigma}_{\mathrm{m}}^{\mathrm{mod}}$
$\displaystyle=\left(1-f_{\mathrm{mis}}\right)\Delta{\Sigma}_{\mathrm{m}}^{\mathrm{cen}}+f_{\mathrm{mis}}\Delta{\Sigma}_{\mathrm{m}}^{\mathrm{mis}}\,.$
(36)

The profile of ${\Sigma}_{\mathrm{m}}^{\mathrm{cen}}$ represents the surface
mass density of a perfectly centered halo, which is evaluated using a NFW
model. Meanwhile, the profile of ${\Sigma}_{\mathrm{m}}^{\mathrm{mis}}$ is
modeled as an average of a set of NFW models with miscentering, weighted by a
miscentering distribution $P\left(R_{\mathrm{mis}}\right)$, as

${\Sigma}_{\mathrm{m}}^{\mathrm{mis}}\left(R\right)=\int\mathrm{d}R_{\mathrm{mis}}P\left(R_{\mathrm{mis}}\right){\Sigma}_{\mathrm{m}}^{\mathrm{cen}}\left(R|R_{\mathrm{mis}}\right)\,,$
(37)

where ${\Sigma}_{\mathrm{m}}^{\mathrm{cen}}\left(R|R_{\mathrm{mis}}\right)$ is
the surface mass density of a halo with a central offset of $R_{\mathrm{mis}}$
azimuthally averaged over the positional angle $\phi$, expressed as (Yang et
al., 2006; Johnston et al., 2007a)

${\Sigma}_{\mathrm{m}}^{\mathrm{cen}}\left(R|R_{\mathrm{mis}}\right)=\int\limits_{0}^{2\pi}\frac{\mathrm{d}\phi}{2\pi}{\Sigma}_{\mathrm{m}}^{\mathrm{cen}}\left(\sqrt{R^{2}+R_{\mathrm{mis}}^{2}+2RR_{\mathrm{mis}}\cos\phi}\right)\,.$
(38)

When estimating ${\Sigma}_{\mathrm{m}}^{\mathrm{cen}}$ and
${\Sigma}_{\mathrm{m}}^{\mathrm{mis}}$ in Eqs. (36), we fix the concentration
parameter, given a weak-lensing mass $M_{\mathrm{WL}}$, according to the
concentration-to-mass relation from Diemer & Kravtsov (2015). With the
concentration parameter, we calculate the corresponding scale radius
$r_{\mathrm{s}}$ at the given $M_{\mathrm{WL}}$ in order to convert the offset
radius $R_{\mathrm{mis}}$ in Eq. (37) to the dimensionless radius $x$, i.e.,
$x\equiv R_{\mathrm{mis}}/r_{\mathrm{s}}$. In this way, we evaluate
$P\left(R_{\mathrm{mis}}\right)$ at $R_{\mathrm{mis}}$ following the second
component of the derived offset distribution $P\left(x\right)$, as in Eq.
(20). That is, we ignore the centering component, $P_{\mathrm{cen}}$, in the
offset distribution when calculating the profiles. This is a reasonable
approximation, because a miscentered shear profile with a miscentering
distribution at the level of ${\sigma}_{\mathrm{cen}}\approx 0.2$ in Eq. (21)
is not significantly different from that of a perfectly centered model,
${\Sigma}_{\mathrm{m}}^{\mathrm{cen}}$, at the radial range of interest in
this work. By doing so, we achieve a faster calculation of
$\Delta{\Sigma}_{\mathrm{m}}^{\mathrm{mod}}$ in Eq. (36).

We note that we statistically account for the effect of miscentering in the
modeling of observed shear profiles by assuming a universal distribution
$P\left(R_{\mathrm{mis}}\right)$, as an average of all clusters. This approach
might not be optimal for individual clusters. However, this is statistically
correct in terms of accounting for the miscentering in an average behavior,
which is the main focus of this work. We refer readers to Sect. 3.7 for
further discussion of this topic.

Finally, we compute the log probability of observing the shear profile
$g_{+}\left(\theta\right)$ given a weak-lensing mass $M_{\mathrm{WL}}$ and a
set of parameters,
$\mathbf{p}=\mathbf{p}_{\mathrm{mis}}\cup\mathbf{p}_{\mathrm{c}}$, while
accounting for the covariance among radial bins, as

$\ln
P\left(g_{+}|M_{\mathrm{WL}},z_{\mathrm{cl}},\mathbf{p}\right)=\frac{-1}{2}\times\\\
\left(g_{+}^{\mathrm{mod}}\left(\theta|M_{\mathrm{WL}},z_{\mathrm{cl}},\mathbf{p}\right)-g_{+}\left(\theta\right)\right)^{\mathrm{T}}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\scriptscriptstyle\bullet$}}}}}\leavevmode\nobreak\
\mathbb{C}^{-1}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\scriptscriptstyle\bullet$}}}}}\leavevmode\nobreak\
\\\
\left(g_{+}^{\mathrm{mod}}\left(\theta|M_{\mathrm{WL}},z_{\mathrm{cl}},\mathbf{p}\right)-g_{+}\left(\theta\right)\right)+\mathrm{constant}\,.$
(39)

##  4 Simulation calibrations

In this work, we use simulations to calibrate the relation between the
observed and true underlying quantities. These quantities are the X-ray count
rate and the weak-lensing inferred mass, which are described in Sects. 4.1 and
4.2, respectively.

###  4.1 The X-ray count-rate bias

The observed count rate estimated by the eSASS pipeline is biased with respect
to that enclosed by the cluster radius $R_{500}$, which is referred to as the
true count rate $\hat{\eta}$. The main reason for this is that the observed
count rate is estimated based on the best-fit model, which depends on the
scale of extendedness EXT, which in turn is the best-fit core size that best
facilitates the cluster detection and does not reflect the true underlying
cluster core radius. When constraining the true count-rate-to-mass-and-
redshift relation, this bias needs to be accounted for.

In this work, we use large simulations to empirically calibrate the bias
between the observed and true count rates at a given true cluster mass and
redshift. The simulation setup of point sources and clusters is fully
described in Comparat et al. (2019) and Comparat et al. (2020), respectively,
to which we refer readers for more details. In particular, the detail of the
simulation scheme applied to eFEDS is described in Liu et al. (2021b) and in
Sect. 3 of Liu et al. (2021a), which is specifically for clusters. For each
simulated cluster, we first evaluate the true count rate $\hat{\eta}$ enclosed
by $R_{500}$ given the cluster true mass $M_{500}$. We note that the
$K$-correction is included in calculating $\hat{\eta}$. Meanwhile, we run the
identical eSASS pipeline on the locations of simulated clusters to measure
their observed count rate $\eta$. In this way, we can derive the bias,

$b_{\eta}\equiv\frac{\eta}{\hat{\eta}}\,,$

as a function of cluster true mass and redshift. We only use clusters with
true counts larger than $40$ when calculating $b_{\eta}$, because the purity
and completeness of the cluster sample drops significantly below this
threshold. The threshold at $40$ counts is motivated by the count distribution
of eFEDS clusters: the minimum number of observed counts is $\approx 40$ for
the sample studied in this work. Lowering this threshold would increase the
false positive rate of the resulting cluster catalogs from the simulations. A
threshold of $30$ counts would increase the overall count-rate bias by
$\approx 4\%$, which is still smaller than the modeling residual (see below)
and hence is not a dominant factor in our analysis. The results are shown in
the middle panel of Fig. 6.

We empirically calibrate the bias by a power-law function in mass and
redshift, as

$b_{\eta}\left(M_{500},z\right)=\exp\left(A_{\mathrm{f}}\right)\times\left(\frac{M_{500}}{M_{\mathrm{piv}}}\right)^{B_{\mathrm{f}}+\delta_{\mathrm{f}}\ln\left(\frac{z}{z_{\mathrm{piv}}}\right)}\left(\frac{z}{z_{\mathrm{piv}}}\right)^{\gamma_{\mathrm{f}}}\,,$
(40)

with four parameters:
$\left\\{A_{\mathrm{f}},B_{\mathrm{f}},\delta_{\mathrm{f}},\gamma_{\mathrm{f}}\right\\}$.
We use $M_{\mathrm{piv}}=1.4\times 10^{14}h^{-1}\mathrm{M}_{\odot}$ and
$z_{\mathrm{piv}}=0.35$ in Eq. (40). The best-fit model describing the bias is
shown in the left panel of Fig. 6. This model provides a good description of
the simulated data, as suggested by their difference in the right panel of
Fig. 6. This difference, which is the modeling residual of the count-rate
bias, can be well approximated by a Gaussian distribution with a standard
deviation of $0.08$. This amount of modeling inaccuracy is absorbed into the
Gaussian prior applied to the intrinsic scatter of observed count rates (see
Sect. 5.4). We repeat the end-to-end procedure above for three different sets
of simulations, which results in three sets of best-fit parameters. These
three sets of simulations differ from one another in terms of the underlying
N-body dark-matter-only simulations, the area coverage, and the abundance of
point-source populations (at a level of $25\%$). Exhaustive runs of various
simulations spanning a complete list of configurations are too time-consuming
and are beyond the scope of this work. The goal of employing these different
configurations is to gauge the inaccuracy of the derived count-rate bias as a
systematic uncertainty that will be marginalized over in the analysis.
Therefore, we take the mean and standard error of the three sets of the best-
fit parameters $\left\\{{A_{\mathrm{f}}}_{,i}\leavevmode\nobreak\
,{B_{\mathrm{f}}}_{,i}\leavevmode\nobreak\
,{\delta_{\mathrm{f}}}_{,i}\leavevmode\nobreak\
,{\gamma_{\mathrm{f}}}_{,i}\right\\}_{i=1,2,3}$ separately inferred from the
three simulations as the priors that will be marginalized over. The mean and
standard error of these four parameters among the three simulations are

$\displaystyle A_{\mathrm{f}}$ $\displaystyle=0.18\pm 0.02\,,$ $\displaystyle
B_{\mathrm{f}}$ $\displaystyle=-0.16\pm 0.03\,,$
$\displaystyle\delta_{\mathrm{f}}$ $\displaystyle=-0.015\pm 0.05\,,$
$\displaystyle\gamma_{\mathrm{f}}$ $\displaystyle=0.42\pm 0.03\,.$ (41)

Finally, the constraints in Eq. (4.1) are used as the priors on these
parameters in the modeling of scaling relations, effectively accounting for
the bias in $\eta$.

It is worth mentioning that the synthetic clusters in the simulations are
drawn from a library that is constructed based on real clusters, and that the
resulting X-ray properties of the synthetic clusters have been shown to be in
excellent agreement with the observed relations among luminosities,
temperatures, and masses (Comparat et al., 2020). Moreover, the halo shape,
gas profile, and dynamical state of synthetic clusters are assigned to each
simulated halo depending on the cluster mass and redshift, with dedicated
treatments to capture the signature of cool-core clusters in terms of
emissivity at cores (see Section 3.2 in Comparat et al., 2020). The resulting
catalog of synthetic clusters has been shown to be in good agreement with the
eFEDS clusters, and has been used to construct the selection function in Liu
et al. (2021a).

It is important to stress that the modeling of the “observed” count rate
$\eta$, as the main quantity of interest, is not affected by the resulting
count-rate bias $b_{\eta}$. This is because the parameters used to
characterize the $\eta$–$M_{500}$–$z$ relation are let free (see Eq. (48) in
Sect. 5.1), such that a change in the resulting count-rate bias will be
absorbed into the posteriors of the parameters given a fixed $\eta$. As a
result, the resulting count-rate bias will not change the overall modeling of
the “observed” count-rate-to-mass-and-redshift relation, but only affects the
inferred “true” count rate $\hat{\eta}$ given a true cluster mass $M_{500}$ at
a redshift $z$. However, including the modeling of the count-rate bias is
still necessary, because it gives a flexible and complete functional form with
which to empirically calibrate the $\eta$–$M_{500}$–$z$ relation.

Figure 7:  Best-fit weak-lensing mass bias $b_{\mathrm{WL}}$ normalized to the
pivotal mass $2\times 10^{14}h^{-1}\mathrm{M}_{\odot}$ as a function of
redshift. This plot shows the results of the weak-lensing mass calibration
including the modeling of the shear profiles at $R>0.2h^{-1}\mathrm{Mpc}$. The
results of the four snapshots in the hydrodynamics simulation are shown by the
black points, while the best-fit redshift trend is indicated by the blue line.
We note that we account for the correlated uncertainties on the calibrated
bias among the four snapshots of redshift when obtaining the best-fit
relation.

###  4.2 Weak-lensing mass bias

Weak lensing is expected to be the most direct way to probe the total mass of
clusters, given that it does not rely on any assumption on the dynamical state
of halos. However, the weak-lensing inferred cluster mass could still be
biased due to inaccurate assumptions (e.g., Pratt et al., 2019; Umetsu et al.,
2020). For example, clusters are triaxial (Clowe et al., 2004; Corless & King,
2007; Chiu et al., 2018b), such that the assumption of the spherical symmetry
could introduce bias and additional scatter in the inferred mass. Another
example is that the density profile of clusters could deviate from the assumed
functional form used in modeling due to the presence of correlated large-scale
structures or substructures (Clowe et al., 2004; King & Corless, 2007).
Moreover, the weak-lensing inferred mass depends on the fitting range of radii
(Applegate et al., 2014), the accuracy of source photo-$z$ (Kelly et al.,
2014), and the effect of miscentering (Rozo et al., 2011; von der Linden et
al., 2014). Therefore, the mass inferred from weak lensing is not free from
bias and needs to be calibrated.

In this work, we calibrate the weak-lensing inferred mass by simulations,
similarly to the methodology in Schrabback et al. (2018a), Dietrich et al.
(2019), and Sommer et al. (2021). We note that a similar approach was also
adopted for shear-selected clusters (Chen et al., 2020). The goal is to
calibrate both the weak-lensing mass bias $b_{\mathrm{WL}}$ as a function of
true cluster mass and redshift, and the intrinsic scatter at fixed mass. To
properly account for the various systematic and statistical effects impacting
weak-lensing mass measurements, we follow the method presented by Grandis et
al. (2021). We use the Magneticum444http://www.magneticum.org/index.html
(Dolag et al., in preparation), a suite of hydrodynamical structure formation
simulations, to create a library of halo density contrasts at different halo
masses and redshifts. We then use the systematic and statistical properties of
our lensing survey described in Sect. 3 to transform the density contrasts
into synthetic shear profiles. These properties include (1) the multiplicative
bias as a function of the source redshift, (2) the redshift distribution of
selected sources at a given cluster redshift, (3) the cluster member
contamination555Because there is no clear cluster contamination in this work
(see Sect. 3.6), we set the upper bound of the cluster contamination to $6\%$
at a level of $2\sigma$ in the inner most bin, following a projected NFW model
toward large radii., (4) the miscentering distribution of the eFEDS clusters,
and (5) the photo-$z$ bias of lensing sources observed in the three-year weak-
lensing data. The shear profiles are then drawn including the scatter due to
uncorrelated large-scale structures, as described in Eq. (14). Finally, these
shear profiles are fitted with our shear model following the identical
procedure described in Sect. 3.9 to extract the weak-lensing mass
$M_{\mathrm{WL}}$. We note that the same radial binning is used in deriving
the $M_{\mathrm{WL}}$ from both synthetic shear profiles and those from the
HSC data. Given the true halo mass $M_{500}$ from the simulation, we then fit
a weak-lensing mass to true halo mass relation at a given cluster redshift by
a log-normal distribution, resulting in a set of best-fit mean and variance.
The posteriors on these parameters are derived by drawing independent
realizations of the shear library while varying the weak-lensing parameters,
as described in Grandis et al. (2021).

We perform the same procedure above at four different cluster redshifts:
$0.25$, $0.47$, $0.78$, and $1.18$. With the results at these four snapshots
of redshift, we perform a joint fit of the weak-lensing mass bias as a
function of the cluster true mass and redshift, which is described as

$\left\langle\ln
b_{\mathrm{WL}}\left(M_{500},z\right)\right\rangle\equiv\left\langle\ln\left(\frac{M_{\mathrm{WL}}}{M_{500}}\right)\right\rangle=\\\
\ln A_{\mathrm{WL}}+B_{\mathrm{WL}}\times\ln\left(\frac{M_{500}}{2\times
10^{14}h^{-1}\mathrm{M}_{\odot}}\right)+\gamma_{\mathrm{WL}}\times\ln\left(\frac{1+z}{1+0.6}\right)\,,$
(42)

with log-normal intrinsic scatter $\sigma_{\mathrm{WL}}$ at fixed mass,

$\sigma_{\mathrm{WL}}\equiv\left(\mathrm{Var}\left[\ln\left(b_{\mathrm{WL}}|M_{500}\right)\right]\right)^{\frac{1}{2}}\,.$
(43)

The resulting constraints on the parameters are

$\displaystyle A_{\mathrm{WL}}$ $\displaystyle=0.896\pm 0.076\,,$
$\displaystyle B_{\mathrm{WL}}$ $\displaystyle=-0.055\pm 0.021\,,$
$\displaystyle\gamma_{\mathrm{WL}}$ $\displaystyle=-0.426\pm 0.087\,,$
$\displaystyle\sigma_{\mathrm{WL}}$ $\displaystyle=0.296\pm 0.047\,.$ (44)

We note that we do not observe significant redshift- or mass-dependent scatter
$\sigma_{\mathrm{WL}}$ in the weak-lensing mass bias; we therefore use the
same $\sigma_{\mathrm{WL}}$ for all clusters in this work. The weak-lensing
mass bias obtained from the simulations and the best-fit relation as a
function of redshift are shown in Fig. 7, in which we can see that the model
provides a good description of the resulting $b_{\mathrm{WL}}$.

We also repeat the calibration procedure for the modeling without the cluster
core (i.e., excluding the three inner most bins), for which the resulting
constraints are

$\displaystyle A_{\mathrm{WL}}$ $\displaystyle=0.903\pm 0.057\,,$
$\displaystyle B_{\mathrm{WL}}$ $\displaystyle=-0.057\pm 0.022\,,$
$\displaystyle\gamma_{\mathrm{WL}}$ $\displaystyle=-0.474\pm 0.080\,,$
$\displaystyle\sigma_{\mathrm{WL}}$ $\displaystyle=0.238\pm 0.032\,.$ (45)

These constraints are used as the priors on the parameters of the weak-
lensing-mass-bias-to-mass relation, which accounts for the bias of weak-
lensing masses in the likelihood modeling.

Table 1:  Summary of the priors used in the modeling (see Sect. 5.4 for more
details). The first and second columns represent the name of and the prior on
the parameters, respectively. The former and latter priors on the parameters
of the $M_{\mathrm{WL}}$–$M_{500}$–$z$ relation describe the cases of modeling
the lensing profile with $R>0.2h^{-1}\mathrm{Mpc}$ and
$R>0.5h^{-1}\mathrm{Mpc}$, respectively.

Parameter | Prior
---|---
The $M_{\mathrm{WL}}$–$M_{500}$–$z$ relation
Eq. (42)
$A_{\mathrm{WL}}$ |  $\mathcal{N}(0.896,0.076^{2})$ or $\mathcal{N}(0.903,0.057^{2})$
$B_{\mathrm{WL}}$ |  $\mathcal{N}(-0.055,0.021^{2})$ or $\mathcal{N}(-0.057,0.022^{2})$
$\gamma_{\mathrm{WL}}$ |  $\mathcal{N}(-0.426,0.087^{2})$ or $\mathcal{N}(-0.474,0.08^{2})$
$\sigma_{\mathrm{WL}}$ |  $\mathcal{N}(0.293,0.047^{2})$ or $\mathcal{N}(0.238,0.032^{2})$
The $\eta$–$M_{500}$–$z$ relation
Eq. (48) with $C_{\mathrm{SS},\eta}=2$
$A_{\eta}$ | $\mathcal{U}(0,0.5)$
$B_{\eta}$ | $\mathcal{U}(0,5)$
$\delta_{\eta}$ | $\mathcal{U}(-3,3)$
$\gamma_{\eta}$ | $\mathcal{U}(-3,3)$
$\sigma_{\eta}$ |  $\mathcal{N}(0.3,0.08^{2})$ and $\mathcal{U}(0.05,0.8)$
$A_{\mathrm{f}}$ | $\mathcal{N}(0.18,0.02^{2})$
$B_{\mathrm{f}}$ | $\mathcal{N}(-0.16,0.03^{2})$
$\delta_{\mathrm{f}}$ | $\mathcal{N}(-0.015,0.05^{2})$
$\gamma_{\mathrm{f}}$ | $\mathcal{N}(0.42,0.03^{2})$
The $\mathcal{X}$–$M_{500}$–$z$ relation
Eq. (46)
$B_{\mathcal{X}}$ | $\mathcal{U}(0,5)$
$\delta_{\mathcal{X}}$ | $\mathcal{U}(-3,3)$
$\gamma_{\mathcal{X}}$ | $\mathcal{U}(-3,3)$
$\sigma_{\mathcal{X}}$ | $\mathcal{U}(0.05,0.8)$
$\mathcal{X}=L_{\mathrm{X}}$ with $C_{\mathrm{SS},L_{\mathrm{X}}}=2$
$\frac{A_{L_{\mathrm{X}}}}{\mathrm{ergs}}$ | $\mathcal{U}(10^{41},10^{46})$
$\mathcal{X}=L_{\mathrm{b}}$ with $C_{\mathrm{SS},L_{\mathrm{b}}}=7/3$
$\frac{A_{L_{\mathrm{b}}}}{\mathrm{ergs}}$ | $\mathcal{U}(10^{41},10^{44})$
$\mathcal{X}=T_{\mathrm{X}}$ with $C_{\mathrm{SS},T_{\mathrm{X}}}=2/3$
$\frac{A_{T_{\mathrm{X}}}}{k\mathrm{eV}}$ | $\mathcal{U}(0.5,25)$
$\mathcal{X}=M_{\mathrm{g}}$ with $C_{\mathrm{SS},M_{\mathrm{g}}}=0$
$\frac{A_{M_{\mathrm{g}}}}{\mathrm{M}_{\odot}}$ | $\mathcal{U}(5\times 10^{11},5\times 10^{13})$
$\mathcal{X}=Y_{\mathrm{X}}$ with $C_{\mathrm{SS},Y_{\mathrm{X}}}=2/3$
$\frac{A_{Y_{\mathrm{X}}}}{k\mathrm{eV}\cdot\mathrm{M}_{\odot}}$ | $\mathcal{U}(10^{12},10^{14})$
Correlated scatter
$\rho_{\mathrm{WL},\mathcal{X}}$, $\rho_{\mathrm{WL},\eta}$,$\rho_{\mathcal{X},\eta}$ | $\mathcal{U}(-0.9,0.9)$
Miscentering
$f_{\mathrm{mis}}$ | $\mathcal{N}(0.54,0.02^{2})$
${\sigma}_{\mathrm{mis}}$ | $\mathcal{N}(0.61,0.03^{2})$
Cosmological parameters
$\Omega_{\mathrm{M}}$ | $\mathcal{N}(0.3,0.016^{2})$
$\sigma_{8}$ | $\mathcal{N}(0.8,0.014^{2})$
$H_{0}$ | $\mathcal{N}(70,5.6^{2})$

##  5 Modeling of X-ray-observable-to-mass-and-redshift relations

In this section, we describe the forward modeling of X-ray observables as
functions of the cluster mass and redshift. Given the X-ray observables and
weak-lensing shear profiles of a sample of eFEDS clusters, our goal is to
model the X-ray-observable-to-mass-and-redshift relations by simultaneously
fitting these observables on the basis of individual clusters.

We perform two types of modeling. First, we model the count rate $\eta$ as a
function of cluster mass and redshift together with the mass calibration,
using the observed shear profiles from the HSC data. Second, we model the
follow-up X-ray observables ($L_{\mathrm{X}}$, $L_{\mathrm{b}}$,
$T_{\mathrm{X}}$, $M_{\mathrm{g}}$, and $Y_{\mathrm{X}}$) as a function of
cluster mass and redshift, jointly with the modeling of the count rate and the
weak-lensing mass calibration. These two types of modeling are described in
Sects. 5.2 and 5.3, respectively.

It is important to stress that these fitting frameworks account for Malmquist
bias arising from the selection function of eFEDS clusters, and for Eddington
(1913) bias due to the steep gradient of the halo mass function. We have fully
verified and extensively used these modeling strategies in previous work
(Bocquet et al., 2015; Liu et al., 2015; Chiu et al., 2016b, 2018a; Bulbul et
al., 2019; Bocquet et al., 2019; Chiu et al., 2020b, a; Schrabback et al.,
2021), to which we refer readers for more details.

Although the eFEDS clusters are selected by both the detection likelihood and
extent likelihood, instead of the observed count rate $\eta$ as the main X-ray
proxy used in this work, we note that our strategy to model the scaling
relation still provides the unbiased result, as demonstrated in Grandis et al.
(2020).

It is important to stress that we do not impose any additional cut on the
follow-up observables, including the weak-lensing shear profile, such that the
selection of the cluster sample does not depend on them and is only determined
by the selection used in constructing the initial cluster catalog. Therefore,
the likelihoods presenting the weak-lensing mass calibration and the modeling
of the X-ray scaling relations (in Sects. 5.2 and 5.3, respectively) are
complete and free from the selection bias. Moreover, the inclusion of the
modeling of the correlated scatter among the X-ray follow-up observable, the
weak-lensing mass, and the count rate accounts for the possible selection
biases in the X-ray scaling relations.

By using mock catalogs that are at least ten times larger than the eFEDS
sample, we verified that the modeling codes can recover the input parameters
within statistical uncertainties that are $\sqrt{10}\approx 3$ times smaller
than those expected for the real sample. The mock catalogs are generated
following the procedure in Grandis et al. (2019), and we perform the modeling
on them in an identical way to that used on the real data. This ensures that
the results in this work are robust and unbiased.

###  5.1 Functional forms of X-ray scaling relations

We first summarize the functional forms used to describe the X-ray observable-
to-mass-and-redshift relations. The adopted form follows that in Bulbul et al.
(2019), in which there are five parameters describing the behavior of an
observable as a function of mass and redshift. Specifically, for an X-ray
observable labeled $\mathcal{X}$, we have

  * •

$A_{\mathcal{X}}$ describing the normalization at the pivotal mass
$M_{\mathrm{piv}}$ and the pivotal redshift $z_{\mathrm{piv}}$;

  * •

$B_{\mathcal{X}}$ describing the overall power-law index in mass;

  * •

$\gamma_{\mathcal{X}}$ describing the deviation from the self-similar redshift
evolution (Kaiser, 1986) of the observable $\mathcal{X}$;

  * •

$\delta_{\mathcal{X}}$ characterizing the redshift-dependent power-law index
in the cluster mass;

  * •

$\sigma_{\mathcal{X}}$ describing the log-normal scatter of $\mathcal{X}$ at
fixed mass.

In this way, the X-ray observable-to-mass-and-redshift
($\mathcal{X}$–$M_{500}$–$z$) relation reads,

$\left\langle\ln\mathcal{X}|M_{500}\right\rangle=\ln A_{\mathcal{X}}+\\\
\left[B_{\mathcal{X}}+\delta_{\mathcal{X}}\ln\left(\frac{1+z}{1+z_{\mathrm{piv}}}\right)\right]\times\ln\left(\frac{M_{500}}{M_{\mathrm{piv}}}\right)+\\\
C_{\mathrm{SS},\mathcal{X}}\times\ln\left(\frac{E(z)}{E(z_{\mathrm{piv}})}\right)+\gamma_{\mathcal{X}}\times\ln\left(\frac{1+z}{1+z_{\mathrm{piv}}}\right)\,,$
(46)

with log-normal intrinsic scatter at fixed mass of

$\sigma_{\mathcal{X}}\equiv\left(\mathrm{Var}\left[\ln\mathcal{X}|M_{500}\right]\right)^{\frac{1}{2}}\,,$
(47)

where $C_{\mathrm{SS},\mathcal{X}}$ is the power-law index of the observable
$\mathcal{X}$ in redshift predicted by the self-similar model. In this work,
we have five follow-up X-ray observables, namely
$\mathcal{X}=\\{L_{\mathrm{X}},L_{\mathrm{b}},M_{\mathrm{g}},T_{\mathrm{X}},Y_{\mathrm{X}}\\}$.
We model the $\mathcal{X}$–$M_{500}$–$z$ relation separately for each X-ray
observable. The self-similar prediction of $C_{\mathrm{SS},\mathcal{X}}$ for
each of them is listed in Table 1.

For the count-rate-to-mass-and-redshift ($\eta$–$M_{500}$–$z$) relation, we
additionally include two factors: the scaling depending on the luminosity
distance $D_{\mathrm{L}}$ to the redshift of clusters (Grandis et al., 2019),
and the count-rate bias $b_{\eta}$ calibrated by the simulations (see Sect.
4.1). The resulting $\eta$–$M_{500}$–$z$ relation reads,

$\left\langle\ln\eta|M_{500}\right\rangle=\ln A_{\eta}+\\\
\left[B_{\eta}+\delta_{\eta}\ln\left(\frac{1+z}{1+z_{\mathrm{piv}}}\right)\right]\times\ln\left(\frac{M_{500}}{M_{\mathrm{piv}}}\right)+\\\
C_{\mathrm{SS},\eta}\times\ln\left(\frac{E(z)}{E(z_{\mathrm{piv}})}\right)+\gamma_{\eta}\times\ln\left(\frac{1+z}{1+z_{\mathrm{piv}}}\right)\\\
-2\times\ln\left(\frac{D_{\mathrm{L}}\left(z\right)}{D_{\mathrm{L}}\left(z_{\mathrm{piv}}\right)}\right)+\ln\left(b_{\eta}\left(M_{500},z\right)\right)\,,$
(48)

with log-normal intrinsic scatter at fixed mass of

$\sigma_{\eta}\equiv\left(\mathrm{Var}\left[\ln\eta|M_{500}\right]\right)^{\frac{1}{2}}\,,$
(49)

Throughout this work, except for $T_{\mathrm{X}}$ and $Y_{\mathrm{X}}$, we use
the pivotal mass $M_{\mathrm{piv}}=1.4\times 10^{14}h^{-1}\mathrm{M}_{\odot}$
and the pivotal redshift $z_{\mathrm{piv}}=0.35$ as the median values of our
eFEDS cluster sample. As the scaling relations of $T_{\mathrm{X}}$ and
$Y_{\mathrm{X}}$ are derived based on the subsample (see Sect. 2.4) with a
relatively high mass scale, we instead use the pivotal mass
$M_{\mathrm{piv}}=2.5\times 10^{14}h^{-1}\mathrm{M}_{\odot}$ for
$\mathcal{X}\in\left\\{T_{\mathrm{X}},Y_{\mathrm{X}}\right\\}$.

###  5.2 Weak-lensing mass calibration

Our goal is to model the observed X-ray count rate $\eta$ of eFEDS clusters as
a function of cluster mass and redshift with the weak-lensing calibration
using the observed shear profile of individual clusters.

Specifically, given a cluster at redshift $z_{\mathrm{cl}}$ with an observed
X-ray count rate $\eta$ and a shear profile $g_{+}$, we compute the
probability of observing $g_{+}$ given $\eta$ and a set of the parameter
$\mathbf{p}$, namely,

$\displaystyle P\left(g_{+}|\eta,z_{\mathrm{cl}},\mathbf{p}\right)$
$\displaystyle=\frac{N\left(g_{+},\eta|z_{\mathrm{cl}},\mathbf{p}\right)}{N\left(\eta|z_{\mathrm{cl}},\mathbf{p}\right)}$
$\displaystyle=\frac{\int
P\left(g_{+},\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)n\left(M_{500},z_{\mathrm{cl}}\right)\mathrm{d}M_{500}}{\int
P\left(\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)n\left(M_{500},z_{\mathrm{cl}}\right)\mathrm{d}M_{500}}\,,$
(50)

where $N\left(\eta|z_{\mathrm{cl}},\mathbf{p}\right)$ and
$N\left(g_{+},\eta|z_{\mathrm{cl}},\mathbf{p}\right)$ denote the numbers of
clusters observed with the count rate $\eta$ and the observable set
$\\{g_{+},\eta\\}$, respectively; $n\left(M_{500},z_{\mathrm{cl}}\right)$ is
the halo mass function evaluated at the true cluster mass $M_{500}$ at the
cluster redshift $z_{\mathrm{cl}}$ using the Bocquet et al. (2016) fitting
formula with the inclusion of baryons;
$P\left(\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)$ and
$P\left(g_{+},\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)$ describe the
probabilities of observing $\eta$ and $\\{g_{+},\eta\\}$, respectively, for
the cluster with the true mass $M_{500}$ at the redshift $z_{\mathrm{cl}}$
given the parameter set $\mathbf{p}$. In Eq. (5.2), the inclusion of the halo
mass function is needed to account for the Eddington (1913) bias. We note that
the difference in the final results obtained between the dark-matter-only and
hydro-simulation mass functions is expected to be negligible given the current
sample size (Bocquet et al., 2016; Castro et al., 2021).

The probability $P\left(\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)$ in the
denominator of Eq. (5.2) is calculated as

$P\left(\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)=\int\mathrm{d}\tilde{}\eta
P\left(\eta|\tilde{}\eta\right)P\left(\tilde{}\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)\,,$
(51)

where the term $P\left(\tilde{}\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)$
describes the log-normal distribution of the pristine count rate
$\tilde{}\eta$ with the mean following Eq. (48) and the intrinsic scatter
$\sigma_{\eta}$. This characterizes the intrinsic distribution of the pristine
count rate while including the observed bias $b_{\eta}$ calibrated by the
simulations. The term $P\left(\eta|\tilde{}\eta\right)$ accounts for the
measurement uncertainty of $\eta$ with respect to $\tilde{}\eta$.

The probability $P\left(g_{+},\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)$
reads

$P\left(g_{+},\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)=\\\
\int\int\mathrm{d}\tilde{}\eta\mathrm{d}M_{\mathrm{WL}}P\left(g_{+}|M_{\mathrm{WL}}\right)P\left(\eta|\tilde{}\eta\right)P\left(M_{\mathrm{WL}},\tilde{}\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)\,,$
(52)

in which the term $P\left(g_{+}|M_{\mathrm{WL}}\right)$ accounts for the
measurement uncertainty of the observed $g_{+}$ given the weak-lensing mass
$M_{\mathrm{WL}}$, which is evaluated using Eq. (39). The term
$P\left(M_{\mathrm{WL}},\tilde{}\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)$
then describes a joint distribution of $M_{\mathrm{WL}}$ and $\tilde{}\eta$
characterized by both the intrinsic scatter of these two quantities and the
correlated scatter between them given the underlying true cluster mass
$M_{500}$. That is,
$P\left(M_{\mathrm{WL}},\tilde{}\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)$
follows a multivariate log-normal distribution of
$\left(M_{\mathrm{WL}},\tilde{}\eta\right)$ characterized by an intrinsic
covariance matrix,

$\Sigma_{M_{\mathrm{WL}},\eta}=\begin{pmatrix}\sigma_{\mathrm{WL}}^{2}&\sigma_{\mathrm{WL}}\sigma_{\eta}\rho_{\mathrm{WL},\eta}\\\
\sigma_{\mathrm{WL}}\sigma_{\eta}\rho_{\mathrm{WL},\eta}&\sigma_{\eta}^{2}\end{pmatrix}\,,$
(53)

where $\rho_{\mathrm{WL},\eta}$ is the correlation coefficient between
$M_{\mathrm{WL}}$ and $\hat{\eta}$.

In this way, the joint likelihood of a sample comprising $N_{\mathrm{cl}}$
eFEDS clusters, given the set of observables
$\\{{g_{+}}_{,i},\eta_{i}\\}_{i=1}^{N_{\mathrm{cl}}}$, is calculated as

$\mathcal{L}\left(\mathbf{p}\right)=\prod\limits_{i=1}^{N_{\mathrm{cl}}}P\left({g_{+}}_{,i}|\eta_{i},\mathbf{p}\right)\,.$
(54)

In this framework, we have a parameter set $\mathbf{p}$ consisting of 19 free
parameters, that is,

$\displaystyle\mathbf{p}$
$\displaystyle=\mathbf{p}_{\eta}\cup\mathbf{p}_{\mathrm{F}}\cup\mathbf{p}_{\mathrm{WL}}\cup\mathbf{p}_{\mathrm{mis}}\cup\\{\rho_{\mathrm{WL},\eta}\\}\cup\mathbf{p}_{\mathrm{c}}\,,$
$\displaystyle\mathbf{p}_{\eta}$
$\displaystyle=\\{A_{\eta},B_{\eta},\delta_{\eta},\gamma_{\eta},\sigma_{\eta}\\}\,,$
$\displaystyle\mathbf{p}_{\mathrm{F}}$
$\displaystyle=\\{A_{\mathrm{f}},B_{\mathrm{f}},\gamma_{\mathrm{f}},\sigma_{\mathrm{f}}\\}\,,$
$\displaystyle\mathbf{p}_{\mathrm{WL}}$
$\displaystyle=\\{A_{\mathrm{WL}},B_{\mathrm{WL}},\gamma_{\mathrm{WL}},\sigma_{\mathrm{WL}}\\}\,,$
$\displaystyle\mathbf{p}_{\mathrm{mis}}$
$\displaystyle=\\{f_{\mathrm{mis}},{\sigma}_{\mathrm{mis}}\\}\,,$
$\displaystyle\mathbf{p}_{\mathrm{c}}$
$\displaystyle=\\{\Omega_{\mathrm{M}},H_{0},\sigma_{8}\\},\,$ (55)

where $\mathbf{p}_{\eta}$ describes the $\eta$–$M_{500}$–$z$ relation as in
Eqs. (48); $\mathbf{p}_{\mathrm{F}}$ and $\mathbf{p}_{\mathrm{WL}}$ describe
the $b_{\eta}$–$M_{500}$–$z$ and $M_{\mathrm{WL}}$–$M_{500}$–$z$ relations in
Eqs. (40) and (42), respectively; and $\mathbf{p}_{\mathrm{mis}}$
characterizes the miscentering distribution as in Eq. (20);
$\mathbf{p}_{\mathrm{c}}$ consists of the cosmological parameters, i.e.,
$\Omega_{\mathrm{M}}$, $H_{0}$, and $\sigma_{8}$, which we vary in the
modeling below.

Figure 8:  Parameter constraints in the joint modeling of X-ray count rates
and weak-lensing shear profiles. The results from the modeling of shear
profiles at $R>0.2$ $h^{-1}\mathrm{Mpc}$ and $R>0.5$ $h^{-1}\mathrm{Mpc}$ are
shown in green and blue, respectively. The marginalized posteriors of the
parameters are presented in the diagonal subplots, where the priors are
indicated by the dotted lines. The off-diagonal subplots contain the
correlations between the parameters. The contours display the $1\sigma$ and
$2\sigma$ confidence levels. We only show the priors in the case of fitting
shear profiles at $R>0.5$ $h^{-1}\mathrm{Mpc}$, for clarity. As seen, there is
no clear evidence of a difference in the results between the modeling of shear
profiles at $R>0.2$ $h^{-1}\mathrm{Mpc}$ and $R>0.5$ $h^{-1}\mathrm{Mpc}$.

Figure 9:  Cluster mass $M_{500}$ and redshift $z$ of each cluster color coded
by the quantity $f_{\mathrm{cont}}$. The mass is randomly sampled from the
posterior of $M_{500}$ as an ensemble mass (see Eq. (66)).

Figure 10:  Observed count rate $\eta$ as a function of the cluster mass
$M_{500}$ (left panel) and redshift $z$ (right panel). The circles in the left
and right panels represent the measurements of $\eta$ color-coded by the
cluster redshift and mass, respectively. The gray shaded regions are the best-
fit models. In the left (right) panel, we normalize the observed count rate to
that at the pivotal redshift $z_{\mathrm{piv}}=0.35$ (pivotal mass
$M_{\mathrm{piv}}=1.4\times 10^{14}h^{-1}\mathrm{M}_{\odot}$).

###  5.3 Joint modeling of X-ray scaling relations

A joint modeling of an X-ray observable $\mathcal{X}$ scaling relation and the
weak-lensing mass calibration (as described in Sect. 5.2) is described as
follows. Given an eFEDS cluster at redshift $z_{\mathrm{cl}}$ with an
available set of observables $\\{\mathcal{X},g_{+},\eta\\}$, we evaluate the
probability of observing the cluster with $\mathcal{X}$ and $g_{+}$, given
$\eta$ and the set of parameters $\mathbf{p}$. That is,

$\displaystyle
P\left(\mathcal{X},g_{+}|\eta,z_{\mathrm{cl}},\mathbf{p}\right)$
$\displaystyle=\frac{N\left(\mathcal{X},g_{+},\eta|z_{\mathrm{cl}},\mathbf{p}\right)}{N\left(\eta|z_{\mathrm{cl}},\mathbf{p}\right)}\,$
$\displaystyle=\frac{\int
P\left(\mathcal{X},g_{+},\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)n\left(M_{500},z_{\mathrm{cl}}\right)\mathrm{d}M_{500}}{\int
P\left(\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)n\left(M_{500},z_{\mathrm{cl}}\right)\mathrm{d}M_{500}}\,,$
(56)

where
$P\left(\mathcal{X},g_{+},\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)$ is
the joint probability of observing the cluster with the observables of
$\\{\mathcal{X},g_{+},\eta\\}$ at the true cluster mass $M_{500}$ and redshift
$z_{\mathrm{cl}}$, given the parameter $\mathbf{p}$. This term is calculated
as

$P\left(\mathcal{X},g_{+},\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)=\int\int\int\mathrm{d}\tilde{}\eta\mathrm{d}M_{\mathrm{WL}}\mathrm{d}\tilde{}\mathcal{X}\\\
P\left(\mathcal{X}|\tilde{}\mathcal{X}\right)P\left(g_{+}|M_{\mathrm{WL}}\right)P\left(\eta|\tilde{}\eta\right)P\left(\tilde{}\mathcal{X},M_{\mathrm{WL}},\tilde{}\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)\,.$
(57)

In an analogy to Eq. (52), the term
$P\left(\tilde{}\mathcal{X},M_{\mathrm{WL}},\tilde{}\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)$
describes the intrinsic distribution of
$\\{\tilde{}\mathcal{X},M_{\mathrm{WL}},\tilde{}\eta\\}$ given the true
cluster mass $M_{500}$, assuming a multivariate log-normal distribution
characterized by an intrinsic covariance matrix,

$\Sigma_{M_{\mathrm{WL}},\mathcal{X},\eta}=\begin{pmatrix}\sigma_{\mathrm{WL}}^{2}&\sigma_{\mathrm{WL}}\sigma_{\mathcal{X}}\rho_{\mathrm{WL},\mathcal{X}}&\sigma_{\mathrm{WL}}\sigma_{\eta}\rho_{\mathrm{WL},\eta}\\\
\sigma_{\mathrm{WL}}\sigma_{\mathcal{X}}\rho_{\mathrm{WL},\mathcal{X}}&\sigma_{\mathcal{X}}^{2}&\sigma_{\mathcal{X}}\sigma_{\eta}\rho_{\mathcal{X},\eta}\\\
\sigma_{\mathrm{WL}}\sigma_{\eta}\rho_{\mathrm{WL},\eta}&\sigma_{\mathcal{X}}\sigma_{\eta}\rho_{\mathcal{X},\eta}&\sigma_{\eta}^{2}\end{pmatrix}\,,$
(58)

where $\rho_{\mathrm{X},\mathrm{Y}}$ denotes the correlation coefficient
between the observables $\mathrm{X}$ and $\mathrm{Y}$. The term
$P\left(\mathcal{X}|\tilde{}\mathcal{X}\right)$ evaluates the probability of
observing $\mathcal{X}$ with the measurement uncertainty given the intrinsic
value of $\tilde{}\mathcal{X}$ predicted at the true cluster mass $M_{500}$.

If a cluster is not covered by the HSC survey (thus does not have an observed
weak-lensing shear profile), we evaluate the probability of observing this
cluster with the X-ray observable $\mathcal{X}$ given the count rate $\eta$,
namely,

$\displaystyle P\left(\mathcal{X}|\eta,z_{\mathrm{cl}},\mathbf{p}\right)$
$\displaystyle=\frac{N\left(\mathcal{X},\eta|z_{\mathrm{cl}},\mathbf{p}\right)}{N\left(\eta|z_{\mathrm{cl}},\mathbf{p}\right)}$
$\displaystyle=\frac{\int
P\left(\mathcal{X},\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)n\left(M_{500},z_{\mathrm{cl}}\right)\mathrm{d}M_{500}}{\int
P\left(\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)n\left(M_{500},z_{\mathrm{cl}}\right)\mathrm{d}M_{500}}\,.$
(59)

We note that Eq. (5.3) is identical to Eq. (5.2) with a substitution of
$g_{+}$ with $\mathcal{X}$. Accordingly, the probability
$P\left(\mathcal{X},\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)$ is
calculated as

$P\left(\mathcal{X},\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)=\\\
\int\int\mathrm{d}\tilde{}\eta\mathrm{d}\tilde{}\mathcal{X}P\left(\mathcal{X}|\tilde{}\mathcal{X}\right)P\left(\eta|\tilde{}\eta\right)P\left(\tilde{}\mathcal{X},\tilde{}\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)\,,$
(60)

where
$P\left(\tilde{}\mathcal{X},\tilde{}\eta|M_{500},z_{\mathrm{cl}},\mathbf{p}\right)$
describes the intrinsic distribution of
$\\{\tilde{}\mathcal{X},\tilde{}\eta\\}$ characterized as a two-dimensional
log-normal distribution with a covariance matrix of

$\Sigma_{\mathcal{X},\eta}=\begin{pmatrix}\sigma_{\mathcal{X}}^{2}&\sigma_{\mathcal{X}}\sigma_{\eta}\rho_{\mathcal{X},\eta}\\\
\sigma_{\mathcal{X}}\sigma_{\eta}\rho_{\mathcal{X},\eta}&\sigma_{\eta}^{2}\end{pmatrix}\,.$
(61)

If a cluster does not have the X-ray observable $\mathcal{X}$ but simply a
weak-lensing shear profile, the likelihood reduces to the pure mass
calibration, as in Eq. (5.2).

Finally, the likelihood of the joint modeling given a sample of eFEDS clusters
is calculated as

$\mathcal{L}\left(\mathbf{p}\right)=\prod\limits_{i=1}P\left(\mathcal{X}_{i},{g_{+}}_{,i}|\eta_{i},\mathbf{p}\right)\\\
\times\prod\limits_{j=1}P\left({g_{+}}_{,j}|\eta_{j},\mathbf{p}\right)\times\prod\limits_{k=1}P\left(\mathcal{X}_{k}|\eta_{k},\mathbf{p}\right)\,,$
(62)

where $i$ runs over the clusters with the available observables of
$\\{\mathcal{X},g_{+}\eta\\}$, $j$ runs over the clusters with only observed