<figure><img src="image_1.png" /></figure>

HAL
open science

# 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, Vittorio Ghirardini, Ang Liu, Sebastian Grandis, Esra Bulbul, Y.

Emre Bahar, Johan Comparat, Sebastian Bocquet, Nicolas Clerc, Matthias

Klein, et al.

► To cite this version:

I-Non Chiu, Vittorio Ghirardini, Ang Liu, Sebastian Grandis, Esra Bulbul, et al. 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. Astron.Astrophys., 2022, 661, pp.A11. 10.1051/0004-6361/202141755. [hal-03312618](https://hal.science/hal-03312618)

HAL Id: hal-03312618

https://hal.science/hal-03312618

Submitted on 11 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# 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 (邱奕儻)¹,²,³✉, Vittorio Ghirardini⁴, Ang Liu⁴, Sebastian Grandis⁵, Esra Bulbul⁴, Y. Emre Bahar⁴, Johan Comparat⁴, Sebastian Bocquet⁵, Nicolas Clerc⁶, Matthias Klein⁵, Teng Liu⁴, Xiangchong Li⁷,⁸, Hironao Miyatake⁹,¹⁰,¹¹,⁸, Joseph Mohr⁵,⁴, Surhud More¹²,⁸, Masamune Oguri¹³,⁸,⁷, Nobuhiro Okabe¹⁴,¹⁵,¹⁶, Florian Pacaud¹⁷, Miriam E. Ramos-Ceja⁴, Thomas H. Reiprich¹⁷, Tim Schrabback¹⁷, and Keiichi Umetsu³

*(Affiliations can be found after the references)*

Received 9 July 2021 / Accepted 10 September 2021

## ABSTRACT

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 ≈140 deg², 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 ≲ z ≲ 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 η and the shear profile g+, of individual clusters through the count-rate-to-mass-and-redshift (η–M$_{500}$–z) relation and the weak-lensing-mass-to-mass-and-redshift (M$_{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¹³h⁻¹M☉ ≲ M$_{500}$ ≲ 10¹⁵h⁻¹M☉ with a median of ≈10¹⁴h⁻¹M☉ 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$_X$ and L$_b$), the emission-weighted temperature T$_X$, the mass of intra-cluster medium M$_g$, and the mass proxy Y$_X$, which is the product of T$_X$ and M$_g$. Except for L$_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 ≲1.7σ. 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 – large-scale structure of Universe – cosmology: observations – 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 XXIV 2016; Bocquet et al. 2015, 2019; de Haan et al. 2016) 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 Ω$_M$, the density parameter of matter in the Universe, and σ$_8$, the dispersion of linear density fluctuations on a comoving scale of 8h⁻¹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 ≈ 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 ≳ 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 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 to 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. 2014a,b; 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_M$, $\sigma_8$, and the dark energy equation of state $w$ at levels of a few percent.

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 \ge 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 \ge 1.1$ based on the HSC data sets was performed for a large sample of ($\ge 1700$) optically selected clusters over a footprint with an area of $\approx 140 \text{ deg}^2$ (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 \text{ deg}^2$ 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. 2022), 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$CDM cosmology with $\Omega_M = 0.3$, a mean baryon density $\Omega_b = 0.05$, the Hubble parameter $H_0 = h \times 100 \text{ km s}^{-1} \text{ Mpc}^{-1}$ with $h = 0.7$, $\sigma_8 = 0.8$, and a spectral index of the primordial power spectrum $n_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_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 $N(x, y^2)$ ($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.
<figure><img src="image_2.png" /><figcaption>**Fig. 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_{\text{cont}} < 0.2$ (see Sect. 2.3). These clusters are represented by the circles, color-coded by their observed redshifts $z_{\text{el}}$, with sizes proportional to their observed count rate $\eta$. The other clusters with $f_{\text{cont}} \ge 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.</figcaption></figure>

## 2.1. The eROSITA final equatorial depth survey

The eFEDS is a small survey covering an area of ≈140 deg². The depth of the eFEDS reaches an exposure time of ≈2.2 ksec, 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 ≈1 deg². 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 ≈18'' at 1.49 keV, with an average of ≈26'' 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'' with a typical positional uncertainty at the level of ≈4.6'' (1σ). The survey footprint consists of four scanning-mode units spanning a range of ≈126° to ≈146° (≈−3° to ≈+6°) 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 \ge 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 deg² 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σ limiting magnitudes of a 2'' 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'' at i-band (Mandelbaum et al. 2018a). 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. (2022). 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. (2018a). The shape measurement is rigorously calibrated against the image simulations (Mandelbaum et al. 2018b), such that the systematic uncertainty of the multiplicative bias meets the requirement of ≤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 ≈22 per arcmin². 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 ≈430 deg², in which the footprint of the GAMA09 field significantly overlaps with the eFEDS. As a result, the majority (≈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 ≈50% of those
currently covered by the S19A. This implies an improvement of
√2−1 ≈ 40% in terms of the signal-to-noise ratio (S/N) of weak-
lensing measurements from S16A to S19A.

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. (2022), to
which we refer readers for a complete description. In what fol-
lows, 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 keV. The source detection is performed with the eSSAUsers_201009 software (Brunner et al. 2022). 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 iter- atively performed to obtain a reliable estimation of background noise. A list of X-ray sources is created after the iterations. Then, the β-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 like- lihood $\mathcal{L}_{\text{det}}$, the extent likelihood $\mathcal{L}_{\text{ext}}$, and the observed count rate $\eta$ at 0.2–2.3 keV. 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 bias<fnref n="1" /> of the observed count rate $\eta$ in Sect. 4.1. Finally, the thresholds of $\mathcal{L}_{\text{det}} > 5$ and $\mathcal{L}_{\text{ext}} > 6$ are applied to create the catalog of cluster candidates, resulting in 542 sys- tems in total in the eFEDS. These 542 clusters and their X-ray properties are presented in Liu et al. (2022a).

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. (2022). Briefly, the combination of the MCMF run on
these two data sets yields, for each cluster, the photometric red-
shift estimate $z_{\text{cl}}$, the optical richness $\lambda$, and the estimator $f_{\text{cont}}$
of the probability that the cluster appears as a superposition of
galaxies along the line of sight by chance. The quantity of $f_{\text{cont}}$
indicates the level of contamination in the optical confirmation.
For instance, a value of $f_{\text{cont}} = 0.2$ means that there is a 20%
chance of detecting a candidate with the derived richness $\lambda$ at
the redshift $z_{\text{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. 2022). A cut on $f_{\text{cont}}$ is
equivalent to a reduction in the cluster contamination in the ini-
tial 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 ≈20% using simulations (Liu et al. 2022a). Applying

<figure><img src="image_3.png" /><figcaption>**Fig. 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_{\text{cont}} < 0.2$ (see Sect. 2.3) are shown as the circles, while the rest are marked by red triangles.</figcaption></figure>

the cut of $f_{\text{cont}} < 0.2$ results in a deduction of 18% (or 99 clus-
ters) in the sample size, which is consistent with the expectation
from the simulation. The residual contamination from nonex-
tended sources is suggested to be at a level of ≈4% with the
cut of $f_{\text{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_{\text{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$) rela-
tion 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 observ-
ables, which are the soft-band luminosity $L_X$ at the energy
band of 0.5–2 keV, the bolometric luminosity $L_b$, the emission-
weighted temperature $T_X$, the gas mass $M_g$, and the mass proxy
$Y_X$ of the ICM for each cluster. These measurements will be fully
presented in a forthcoming paper (Bahar et al. 2022), and the
details of the X-ray analysis are given in Ghirardini et al. (2021)
and Liu et al. (2022a). We therefore only outline the important
steps here.

We stress that these follow-up X-ray measurements ($L_X$, $L_b$, $T_X$, $M_g$, $Y_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. 2022). 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

<footnote n="1">¹ We 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}$.</footnote>
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_X$, the X-ray follow-up measurements ($L_X$, $L_b$, and $M_g$) are derived based on the fitting to the X-ray imaging extracted in the soft-energy band of 0.5–2 keV. 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_g$ is obtained by integrating the best-fit model to the cluster radius. Meanwhile, the luminosities ($L_X$ and $L_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 keV and 30 keV.

The temperature $T_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 ≈5 to ≈1000, with the first quartile, median, and the third quartile of 43, 75, and 147, respectively. Meanwhile, the 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_X$ (thus $Y_X$) for all clusters. In this work, we therefore restrict the modeling of the $T_X$ and $Y_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}_{\text{det}} > 50$ and an extent likelihood of $\mathcal{L}_{\text{ext}} > 50$, resulting in 64 clusters. With these cuts on $\mathcal{L}_{\text{det}}$ and $\mathcal{L}_{\text{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 kpc, within which we can measure the temperature with the highest S/N for a maximal number of clusters (Liu et al. 2022a). Given such a low number of X-ray photons, we do not expect significant differences between $T_X(R < R_{500})$ and $T_X(R < 500 \text{ kpc})$. As a result, we find that the difference between $T_X(R < R_{500})$ and $T_X(R < 500 \text{ kpc})$ for the subsample with $\mathcal{L}_{\text{det}} > 50$ and $\mathcal{L}_{\text{ext}} > 50$ has an inverse-variance-weighted mean of $\langle T_X(R < R_{500}) - T_X(R < 500 \text{ kpc}) \rangle = (0.06 \pm 0.36) \text{ keV}$, 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 kpc 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_X(R < R_{500})$, instead of $T_X(R < 500 \text{ kpc})$.

In summary, the modeling of the $T_X$ and $Y_X$ scaling relations is performed based on the subsample of 64 clusters with $\mathcal{L}_{\text{det}} > 50$ and $\mathcal{L}_{\text{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_X$, $L_b$, and $M_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–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_s$ arising from a galaxy cluster at redshift $z_{\text{cl}}$ depends on the distances between the cluster, the observer, and the source. This is characterized by the critical surface density $\Sigma_{\text{crit}}$ defined as

$$ \Sigma_{\text{crit}}(z_{\text{cl}}, z_s) = \frac{c^2}{4\pi G} \frac{1}{D_l(z_{\text{cl}})\beta(z_{\text{cl}}, z_s)}, \quad (1) $$

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

$$ \beta(z_{\text{cl}}, z_s) = \begin{cases} \frac{D_{ls}(z_{\text{cl}}, z_s)}{D_s(z_s)} & \text{if } z_s > z_{\text{cl}} \\ 0 & \text{if } z_s \le z_{\text{cl}} \end{cases}, \quad (2) $$

where $D_{ls}(z_{\text{cl}}, z_s)$ and $D_s(z_s)$ 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}, \quad (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_+(R)$ at the projected radius
R describes the differential surface mass density $\Delta\Sigma_m$ of the cluster at R with respect to the critical surface density inferred from the source, namely,

$$ \gamma_+(R) = \frac{\Delta\Sigma_m(R)}{\Sigma_{\text{crit}}}. \qquad (4) $$

Meanwhile, the convergence $\kappa(R)$ is the ratio of the projected surface mass density $\Sigma_m(R)$ to the critical surface density, as

$$ \kappa(R) = \frac{\Sigma_m(R)}{\Sigma_{\text{crit}}}. \qquad (5) $$

We note that $\Delta\Sigma_m(R) = \Sigma_m(<R) - \Sigma_m(R)$.

The shear profile, $g_+(R)$, 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_{\text{crit}}$ to interpret $\gamma_+$ and $\kappa$ (see Sect. 3.9 for the modeling of $\Sigma_{\text{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_{\text{cat},i} = m_{\text{true}} + dm_{1,i} + dm_{2,i}, \qquad (6) $$

where $m_{\text{cat},i}$ is the blinded multiplicative bias in the i-th shape catalog, $m_{\text{true}}$ is the true estimate of the multiplicative bias, and $dm_{1,i}$ and $dm_{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 $dm_{1,i}$ and $dm_{2,i}$ are encrypted and are different among these three blinded catalogs; the HSC analysis teams are also different. The term $dm_{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 $dm_{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 $dm_{2,i}$ until unblinding. The value of $dm_{2,i}$ varies from -0.1 to 0.1 randomly. Only one blinded catalog among the three has $dm_{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 $dm_{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 $dm_{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_X$, $L_b$, $M_g$, $T_X$, and $Y_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_X$, $L_b$, $T_X$, $M_g$, and $Y_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_{cl}$, the source selection is performed using the observed probability distribution of redshift, $P(z)$, such that a galaxy is considered as a source galaxy if

$$0.98 < \int_{z_{\text{min}}}^{\infty} P(z) dz \text{ and } z_{\text{MC}} < 2.5, \qquad (7)$$

where we set $z_{min} = z_{cl} + 0.2$, and $z_{MC}$ is the point estimator randomly sampled from $P(z)$. The contamination to lensing signals due to cluster members is suggested to be at the level of a few percent using our $P(z)$-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_{cl} = 0.05$, $z_{cl} = 0.25$, $z_{cl} = 0.50$, $z_{cl} = 0.78$, and $z_{cl} = 1.1$, respectively. At the median redshift of eFEDS clusters, $z_{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_m$. The latter require prior knowledge of redshift and the distance-to-redshift relation with assumed cosmological parameters to convert $g_+$ to $\Delta\Sigma_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_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_{cl}$, we derive the tangential shear profile $g_+$ in angular bins of clustercentric radius $\theta_i$,

$$g_+(\theta_i) = \frac{1}{2\Re(\theta_i)(1 + \Re(\theta_i))} \frac{\sum_{s \in \theta_i} w_s e_{+,s}}{\sum_{s \in \theta_i} w_s}, \qquad (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; $\Re(\theta_i)$ is the shear response (see also Mandelbaum et al. 2005); the factor $(1 + \Re(\theta_i))$ is the correction accounting for multiplicative shear bias, which is calibrated against the image simulations (Mandelbaum et al. 2018a,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, \qquad (9)$$

where $(e_1, e_2)$ 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_e^2 + e_{\text{rms}}^2}, \qquad (10)$$

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

The shear response $\Re(\theta_i)$ in the angular bin $\theta_i$ is calculated as

$$\Re(\theta_i) = 1 - \frac{\sum_{s \in \theta_i} w_s e_{\text{rms},s}^2}{\sum_{s \in \theta_i} w_s}, \qquad (11)$$

where the subscript $s$ runs over the source sample in the angular bin $\theta_i$. The factor $(1 + \Re(\theta_i))$ is calculated as

$$1 + \Re(\theta_i) = \frac{\sum_{s \in \theta_i} w_s (1 + m_s)}{\sum_{s \in \theta_i} w_s}, \qquad (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 $h^{-1}$Mpc and 3.5 $h^{-1}$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_+(\theta)$ 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(\theta)$ 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. \qquad (13)$$

The azimuthally averaged cross-shear profile $e_\times(R)$ 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(\theta)$ 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 $\mathcal{C}$ comprises a component characterizing the noise of shape measurements, $\mathcal{C}_{\text{shape}}$, and another component $\mathcal{C}_{\text{uLSS}}$ accounting for the uncertainty arising from uncorrelated large-scale structures around the cluster. That is,

$$ \mathcal{C} = \mathcal{C}_{\text{shape}} + \mathcal{C}_{\text{uLSS}}, \qquad (14) $$

in which

$$ \mathcal{C}_{\text{shape},i,j} = \sigma_{\text{shape}}^2(\theta_i)\delta_{i,j}, $$

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

$$ \sigma_{\text{shape}}^2(\theta_i) = \frac{1}{4\Re^2(\theta_i)(1 + \Re(\theta_i))^2 \sum_{s \in \theta_i} w_s}. \qquad (15) $$

We compute $\mathcal{C}_{\text{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_m^{\text{NL}}$ at the cluster redshift $z_{\text{cl}}$ using the code CAMB (Lewis et al. 2000) in order to derive the lensing power spectrum $C_l^{\text{KK}}$. 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^{\text{KK}}$. Finally, we compute $\mathcal{C}_{\text{uLSS}}$ by integrating the product of $C_l^{\text{KK}}$ and two second-order Bessel functions associated with the radial binning of shear profiles. This calculation of $\mathcal{C}_{\text{uLSS}}$ is independently repeated for all clusters using the cosmological parameters fixed to the fiducial cosmology.

We note that the term $\mathcal{C}_{\text{uLSS}}$ depends on the underlying cosmology. Ideally, the variation in $\mathcal{C}_{\text{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 $\mathcal{C}_{\text{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 $\mathcal{C}_{\text{uLSS}}$ so as to improve the speed of our calculation. We leave the improvement of a cosmology-dependent $\mathcal{C}_{\text{uLSS}}$ to future work.

We also note that the intrinsic variation of shear profiles at fixed cluster mass<fnref n="2" /> – 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).

### 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(z)$, 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(z)$ 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(z)$.

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(z; R)$ of selected source galaxies at the radial bin $R$ can be written as

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

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

$$ \int dz P(z; R) = (1 - f_{\text{cl}}(R)) \int dz P_{\text{bkg}}(z) + f_{\text{cl}}(R) \int dz P_{\text{cl}}(z; R), $$

<footnote n="2">² This is characterized by the term $C^{\text{int}}$ in Eq. (16) of Umetsu et al. (2020).</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 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_{\text{lolim}}$ and $z_{\text{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 $\le 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.</figcaption></figure>

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

$$f_{\text{cl}}(R) = 1 - \frac{\int_{z \in \mathfrak{Z}} dz P(z; R)}{\int_{z \in \mathfrak{Z}} dz P_{\text{bkg}}(z)} \quad (17)$$

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

In practice, in the interest of more precise measurements, we estimate $f_{\text{cl}}(R)$ 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 \le 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(z; R_i)$ of redshift observed at the radial bin $R_i$ as

$$P_j(z; R_i) = \frac{\sum_{s \in j, R_i} w_s P_s(z)}{\sum_{s \in j, R_i} w_s}, \quad (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

$$\langle P \rangle (z; R_i) = \frac{\sum_{j=1}^{N_{\text{cl}}} P_j(z; R_i) N_{\text{gal},j}(R_i)}{\sum_{j=1}^{N_{\text{cl}}} N_{\text{gal},j}(R_i)}, \quad (19)$$

where $j$ runs over the subsample of $N_{\text{cl}}$ clusters, and $N_{\text{gal},j}(R_i)$ is the number of observed galaxies at the radius bin $R_i$ around the cluster $j$. We then define the redshift interval $\mathfrak{Z}$

as $\{z|z < z_{\text{lolim}}\} \cup \{z|z > z_{\text{uplim}}\}$, where $z_{\text{lolim}} = \min\{z_j\}_{j=1}^{N_{\text{cl}}} - \Delta z$, $z_{\text{uplim}} = \max\{z_j\}_{j=1}^{N_{\text{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 $\langle P \rangle (z; R_i)$ over the interval $\mathfrak{Z}$.

Concomitently, we estimate $P_{\text{bkg}}(z)$ by repeating Eqs. (18) and (19) with the same scheme of radial binning and cluster redshift on a number of $N_{\text{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_{\text{cl}}(R)$ 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 $\le 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. 2022), 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) = (1 - f_{\text{mis}}) \times P_{\text{cen}}(x|\sigma_{\text{cen}}, \alpha) + f_{\text{mis}} \times P_{\text{mis}}(x|\sigma_{\text{mis}}). \quad (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_{\text{cen}}(x|\sigma_{\text{cen}}, \alpha)$,

which reads

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

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

$$P_{\text{mis}}(x|\sigma_{\text{mis}}) = \frac{1}{\sigma_{\text{mis}}} \left(\frac{x}{\sigma_{\text{mis}}}\right) \exp\left(-\frac{1}{2} \left(\frac{x}{\sigma_{\text{mis}}}\right)^2\right). \quad (22)$$

These two components are weighted by the miscentering fraction $f_{\text{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_{\text{cen}}(x)$ occurs, and the steeper fall-off toward $x=0$. Moreover, $P_{\text{cen}}(x)$ reduces to a half-normal distribution and zero, as $\alpha$ approaches zero and $+\infty$, respectively. Also, $P_{\text{cen}}$ is identical to $P_{\text{mis}}$ if $\alpha = 1$.

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_{\text{off}}/r_s$, where $R_{\text{off}}$ is the offset between the optical and X-ray centers, and $r_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_{\text{piv}} = 1.4 \times 10^{14} h^{-1} M_{\odot}$ at the pivotal redshift $z_{\text{piv}} = 0.35$. For each cluster, the offset $R_{\text{off}}$ is calculated first in the physical unit using the redshift-distance relation inferred in the fiducial cosmology. We then convert $R_{\text{off}}$ to $x$ using the fixed scale radius $r_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

$$(f_{\text{mis}}, \sigma_{\text{cen}}, \alpha, \sigma_{\text{mis}}) = (0.54 \pm 0.02, 0.17 \pm 0.01, 0.26 \pm 0.05, 0.61 \pm 0.03). \quad (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
<figure><img src="image_5.png" /><figcaption>Fig. 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.</figcaption></figure>

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$ kpc) 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}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).

### 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_{crit}$ in converting $g_+$ to $\Delta\Sigma_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_{cl})$ in the inferred mass of a cluster at the redshift $z_{cl}$ due to the bias in photo-z can be expressed in terms of the lensing efficiency $\beta$, as

$$ \frac{\beta_{\text{true}}}{\beta_{\text{obs}}} (z_{\text{cl}}) \equiv 1 + b_z(z_{\text{cl}}), \qquad (24) $$

where $\beta_{true}$ and $\beta_{obs}$ are the lensing efficiency inferred by the true and observed redshifts, respectively. In an ideal case, $\beta_{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_{cl})$. In what follows, we summarize the steps. First, the photometric redshifts $z_{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 catalog<fnref n="3" />. 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_{COSMOS30}$ as the true redshift, and the photometric redshift distribution $P_{HSC}(z)$ 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_{SOM}$, such that the weighted distributions of the properties match those of the

<footnote n="3">³ We 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/</footnote>
HSC sources (see also Medezinski et al. 2018b). Finally, we include the new weights $w_{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_{cl}$, we select the source sample, labeled S, in the COSMOS field according to the description in Sect. 3.3, and compute

$$ \beta_{\text{true/obs}}(z_{\text{cl}}) = \int \beta(z_{\text{cl}}, z) P_{\text{true/obs}}(z) dz, \qquad (25) $$

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

$$ P_{\text{obs}}(z) = \frac{\sum_{s \in S} w_{\text{SOM},s} w_s P_{\text{HSC},s}(z)}{\sum_{s \in S} w_{\text{SOM},s} w_s}. \qquad (26) $$

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

In practice, there is one additional factor that needs to be considered in this work: the weighting factor $w_{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_{\text{cl}}) = \frac{\langle \beta_{\text{S16A}}(z_{\text{cl}}) \rangle}{\langle \beta_{\text{S19A}}(z_{\text{cl}}) \rangle}, $$

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

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

Figure 5 shows the resulting $b_z(z_{cl})$ as a function of cluster redshift $z_{cl}$. As seen, the photo-z bias on the cluster mass is estimated to be $\le 2\%$ for clusters at redshift $z_{cl} \le 0.7$, and becomes $\le 6\%$ for $z_{cl} \ge 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,\text{S19A}}(z_{cl})$ is included in deriving the weak-lensing mass bias for eFEDS clusters in Sect. 4.2.

### 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_{cl}$, there are three input observables used in the modeling:

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

- the lensing covariance matrix $C$, derived as Eq. (14), and

<figure><img src="image_6.png" /><figcaption>Fig. 5. Photo-z bias $b_z(z_{cl})$ as a function of cluster redshift $z_{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_{cl})$ is applied to the calibration of weak-lensing mass bias for eFEDS clusters (see Sect. 4.2).</figcaption></figure>

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

$$ g_{+}^{\text{mod}}(\theta | M_{\text{WL}}, \mathbf{p}_{\text{mis}}, \mathbf{p}_{c}) = \frac{\gamma_{\text{mod}}(\theta | M_{\text{WL}}, \mathbf{p}_{\text{mis}}, \mathbf{p}_{c})}{1 - \kappa_{\text{mod}}(\theta | M_{\text{WL}}, \mathbf{p}_{\text{mis}}, \mathbf{p}_{c})} \times \left( 1 + \kappa_{\text{mod}}(\theta | M_{\text{WL}}, \mathbf{p}_{\text{mis}}, \mathbf{p}_{c}) \frac{\langle \beta^2(\mathbf{p}_{c}) \rangle}{\langle \beta(\mathbf{p}_{c}) \rangle^2} - 1 \right) \qquad (28) $$

in which

$$ \langle \beta(\mathbf{p}_c) \rangle = \int \beta(z_{\text{cl}}, z | \mathbf{p}_c) P(z) dz, \qquad (29) $$

$$ \langle \beta^2(\mathbf{p}_c) \rangle = \int \beta^2(z_{\text{cl}}, z | \mathbf{p}_c) P(z) dz, \qquad (30) $$

and the profiles $\kappa_{\text{mod}}$ and $\gamma_{\text{mod}}$ are calculated as

$$ \kappa_{\text{mod}}(\theta | M_{\text{WL}}, \mathbf{p}_{\text{mis}}, \mathbf{p}_{c}) = \frac{\sum_{m}^{\text{mod}} (\theta | M_{\text{WL}}, \mathbf{p}_{\text{mis}}, \mathbf{p}_{c})}{\sum_{\text{crit}} (\mathbf{p}_{c})}, \qquad (31) $$

and

$$ \gamma_{\text{mod}}(\theta | M_{\text{WL}}, \mathbf{p}_{\text{mis}}, \mathbf{p}_{c}) = \frac{\Delta \sum_{m}^{\text{mod}} (\theta | M_{\text{WL}}, \mathbf{p}_{\text{mis}}, \mathbf{p}_{c})}{\sum_{\text{crit}} (\mathbf{p}_{c})}, \qquad (32) $$

where

$$ \Sigma_{\text{crit}}(\mathbf{p}_c) = \frac{c^2}{4\pi G} \frac{1}{D_1(z_{\text{cl}}|\mathbf{p}_c)} \langle \beta \rangle. \qquad (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(z)$ used to calculate Eqs. (29) and (30) is estimated as

$$P(z) = \frac{\sum_{s \in j} w_s P_s(z)}{\sum_{s \in j} w_s}, \qquad (34)$$

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

Given the set of cosmological parameters $\mathbf{p}_c$, the calculation of $\Sigma_m^{\text{mod}}$ and $\Delta\Sigma_m^{\text{mod}}$ is performed in the radial bins in physical units by converting $\theta$ to $R$, as

$$R = D_1(z_{\text{cl}}|\mathbf{p}_c) \times \theta. \qquad (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_m^{\text{mod}}$ (thus in $\Delta\Sigma_m^{\text{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_{\text{mis}}$. That is,

$$\Delta\Sigma_m^{\text{mod}} = (1 - f_{\text{mis}}) \Delta\Sigma_m^{\text{cen}} + f_{\text{mis}} \Delta\Sigma_m^{\text{mis}}. \qquad (36)$$

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

$$\Sigma_m^{\text{mis}}(R) = \int dR_{\text{mis}} P(R_{\text{mis}}) \Sigma_m^{\text{cen}}(R|R_{\text{mis}}), \qquad (37)$$

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

$$\Sigma_m^{\text{cen}}(R|R_{\text{mis}}) = \int_0^{2\pi} \frac{d\phi}{2\pi} \Sigma_m^{\text{cen}} \left( \sqrt{R^2 + R_{\text{mis}}^2 + 2RR_{\text{mis}}\cos\phi} \right). \qquad (38)$$

When estimating $\Sigma_m^{\text{cen}}$ and $\Sigma_m^{\text{mis}}$ in Eqs. (36), we fix the concentration parameter, given a weak-lensing mass $M_{\text{WL}}$, according to the concentration-to-mass relation from Diemer & Kravtsov (2015). With the concentration parameter, we calculate the corresponding scale radius $r_s$ at the given $M_{\text{WL}}$ in order to convert the offset radius $R_{\text{mis}}$ in Eq. (37) to the dimensionless radius $x$, i.e., $x \equiv R_{\text{mis}}/r_s$. In this way, we evaluate $P(R_{\text{mis}})$ at $R_{\text{mis}}$ following the second component of the derived offset distribution $P(x)$, as in Eq. (20). That is, we ignore the centering component, $P_{\text{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_{\text{cen}} \approx 0.2$ in Eq. (21) is not significantly different from that of a perfectly centered model, $\Sigma_m^{\text{cen}}$, at the radial range of interest in this work. By doing so, we achieve a faster calculation of $\Delta\Sigma_m^{\text{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(R_{\text{mis}})$, 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_+(\theta)$ given a weak-lensing mass $M_{\text{WL}}$ and a set of parameters, $\mathbf{p} = \mathbf{p}_{\text{mis}} \cup \mathbf{p}_c$, while accounting for the covariance among radial bins, as

$$\ln P(g_+ | M_{\text{WL}}, z_{\text{cl}}, \mathbf{p}) = \frac{-1}{2} \times (g_+^{\text{mod}}(\theta | M_{\text{WL}}, z_{\text{cl}}, \mathbf{p}) - g_+(\theta))^{\text{T}} \cdot \mathbb{C}^{-1} \times (g_+^{\text{mod}}(\theta | M_{\text{WL}}, z_{\text{cl}}, \mathbf{p}) - g_+(\theta)) + \text{constant.} \qquad (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. (2022b) and in Sect. 3 of Liu et al. (2022a), 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 ≈40 for the sample studied in this work. Lowering this threshold would increase the false
<figure><img src="image_7.png" /><figcaption>Fig. 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(b_\eta)$. 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.</figcaption></figure>

positive rate of the resulting cluster catalogs from the simulations. A threshold of 30 counts would increase the overall count-rate bias by ≈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}(M_{500}, z) = \exp(A_{f}) \times \left(\frac{M_{500}}{M_{\text{piv}}}\right)^{B_{f} + \delta_{f} \ln\left(\frac{z}{z_{\text{piv}}}\right)} \left(\frac{z}{z_{\text{piv}}}\right)^{\gamma_{f}}, \quad (40)$$

with four parameters: {$A_{f}, B_{f}, \delta_{f}, \gamma_{f}$}. We use $M_{\text{piv}} = 1.4 \times 10^{14}h^{-1}M_{\odot}$ and $z_{\text{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 {$A_{f,i}, B_{f,i}, \delta_{f,i}, \gamma_{f,i}$}$_{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

$$ \begin{aligned} A_{f} &= 0.18 \pm 0.02, \\ B_{f} &= -0.16 \pm 0.03, \\ \delta_{f} &= -0.015 \pm 0.05, \\ \gamma_{f} &= 0.42 \pm 0.03. \end{aligned} \quad (41) $$

Finally, the constraints in Eq. (41) 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 Sect. 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. (2022a).

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.

## 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. 2014b). 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. (2022). 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_{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 Magneticum<fnref n="4" /> (Dolag et al., in prep.), 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 contamination<fnref n="5" />, (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_{WL}$. We note that the same radial binning is used in deriving the $M_{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

$$ \begin{aligned} \langle \ln b_{WL} (M_{500}, z) \rangle &= \left\langle \ln \left( \frac{M_{WL}}{M_{500}} \right) \right\rangle \\ &= \ln A_{WL} + B_{WL} \times \ln \left( \frac{M_{500}}{2 \times 10^{14} h^{-1} M_{\odot}} \right) + \gamma_{WL} \times \ln \left( \frac{1+z}{1+0.6} \right), \end{aligned} \quad (42) $$

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

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

<figure><img src="image_8.png" /><figcaption>Fig. 7. Best-fit weak-lensing mass bias $b_{WL}$ normalized to the pivotal mass $2 \times 10^{14} h^{-1} 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.2 h^{-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.</figcaption></figure>

The resulting constraints on the parameters are

$$ \begin{aligned} A_{WL} &= 0.896 \pm 0.076, \\ B_{WL} &= -0.055 \pm 0.021, \\ \gamma_{WL} &= -0.426 \pm 0.087, \\ \sigma_{WL} &= 0.296 \pm 0.047. \end{aligned} \quad (44) $$

We note that we do not observe significant redshift- or mass-dependent scatter $\sigma_{WL}$ in the weak-lensing mass bias; we therefore use the same $\sigma_{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_{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

$$ \begin{aligned} A_{WL} &= 0.903 \pm 0.057, \\ B_{WL} &= -0.057 \pm 0.022, \\ \gamma_{WL} &= -0.474 \pm 0.080, \\ \sigma_{WL} &= 0.238 \pm 0.032. \end{aligned} \quad (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.

<footnote n="4">⁴ http://www.magneticum.org/index.html</footnote>

<footnote n="5">⁵ Because 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.</footnote>
## 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_X$, $L_b$, $T_X$, $M_g$, and $Y_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, 2019; Liu et al. 2015; Chiu et al. 2016b, 2018a, 2020b,a; Bulbul et al. 2019; 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_X$ describing the normalization at the pivotal mass $M_{\text{piv}}$ and the pivotal redshift $z_{\text{piv}}$;

- $B_X$ describing the overall power-law index in mass;

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

- $\delta_X$ characterizing the redshift-dependent power-law index in the cluster mass;

- $\sigma_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,

$$ \begin{aligned} \langle \ln \mathcal{X} | M_{500} \rangle = \ln A_X &+ \left[ B_X + \delta_X \ln \left( \frac{1+z}{1+z_{\text{piv}}} \right) \right] \times \ln \left( \frac{M_{500}}{M_{\text{piv}}} \right) \\ &+ C_{\text{SS},X} \times \ln \left( \frac{E(z)}{E(z_{\text{piv}})} \right) + \gamma_X \times \ln \left( \frac{1+z}{1+z_{\text{piv}}} \right), \end{aligned} \quad (46) $$

with log-normal intrinsic scatter at fixed mass of

$$ \sigma_X = (\text{Var}[\ln \mathcal{X} | M_{500}])^{\frac{1}{2}}, \quad (47) $$

where $C_{\text{SS},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_X, L_b, M_g, T_X, Y_X\}$. We model the $\mathcal{X}-M_{500}-z$ relation separately for each X-ray observable. The self-similar prediction of $C_{\text{SS},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_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,

$$ \begin{aligned} \langle \ln \eta | M_{500} \rangle = \ln A_\eta &+ \left[ B_\eta + \delta_\eta \ln \left( \frac{1+z}{1+z_{\text{piv}}} \right) \right] \times \ln \left( \frac{M_{500}}{M_{\text{piv}}} \right) \\ &+ C_{\text{SS},\eta} \times \ln \left( \frac{E(z)}{E(z_{\text{piv}})} \right) + \gamma_\eta \times \ln \left( \frac{1+z}{1+z_{\text{piv}}} \right) \\ &- 2 \times \ln \left( \frac{D_L(z)}{D_L(z_{\text{piv}})} \right) + \ln (b_\eta (M_{500}, z)), \end{aligned} \quad (48) $$

with log-normal intrinsic scatter at fixed mass of

$$ \sigma_{\eta} = (\text{Var}[\ln \eta | M_{500}])^{\frac{1}{2}}, \quad (49) $$

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

### 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_{\text{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,

$$ \begin{aligned} P(g_+ | \eta, z_{\text{cl}}, \mathbf{p}) &= \frac{N(g_+, \eta | z_{\text{cl}}, \mathbf{p})}{N(\eta | z_{\text{cl}}, \mathbf{p})} \\ &= \frac{\int P(g_+, \eta | M_{500}, z_{\text{cl}}, \mathbf{p}) n(M_{500}, z_{\text{cl}}) dM_{500}}{\int P(\eta | M_{500}, z_{\text{cl}}, \mathbf{p}) n(M_{500}, z_{\text{cl}}) dM_{500}}, \end{aligned} \quad (50) $$
<table><caption>**Table 1.** Summary of the priors used in the modeling (see Sect. 5.4 for more details).</caption><thead><tr><th>Parameter</th><th>Prior</th></tr></thead><tbody><tr><td></td><td>The M<sub>WL</sub>-M<sub>500</sub>-z relation<br>Eq. (42)</td></tr><tr><td>A<sub>WL</sub></td><td>N(0.896, 0.076<sup>2</sup>) or N(0.903, 0.057<sup>2</sup>)</td></tr><tr><td>B<sub>WL</sub></td><td>N(-0.055, 0.021<sup>2</sup>) or N(-0.057, 0.022<sup>2</sup>)</td></tr><tr><td>γ<sub>WL</sub></td><td>N(-0.426, 0.087<sup>2</sup>) or N(-0.474, 0.08<sup>2</sup>)</td></tr><tr><td>σ<sub>WL</sub></td><td>N(0.293, 0.047<sup>2</sup>) or N(0.238, 0.032<sup>2</sup>)</td></tr><tr><td></td><td>The η-M<sub>500</sub>-z relation<br>Eq. (48) with C<sub>SS,η</sub> = 2</td></tr><tr><td>A<sub>η</sub></td><td>U(0, 0.5)</td></tr><tr><td>B<sub>η</sub></td><td>U(0, 5)</td></tr><tr><td>δ<sub>η</sub></td><td>U(-3, 3)</td></tr><tr><td>γ<sub>η</sub></td><td>U(-3, 3)</td></tr><tr><td>σ<sub>η</sub></td><td>N(0.3, 0.08<sup>2</sup>) and U(0.05, 0.8)</td></tr><tr><td>A<sub>f</sub></td><td>N(0.18, 0.02<sup>2</sup>)</td></tr><tr><td>B<sub>f</sub></td><td>N(-0.16, 0.03<sup>2</sup>)</td></tr><tr><td>δ<sub>f</sub></td><td>N(-0.015, 0.05<sup>2</sup>)</td></tr><tr><td>γ<sub>f</sub></td><td>N(0.42, 0.03<sup>2</sup>)</td></tr><tr><td></td><td>The X-M<sub>500</sub>-z relation<br>Eq. (46)</td></tr><tr><td>B<sub>X</sub></td><td>U(0, 5)</td></tr><tr><td>δ<sub>X</sub></td><td>U(-3, 3)</td></tr><tr><td>γ<sub>X</sub></td><td>U(-3, 3)</td></tr><tr><td>σ<sub>X</sub></td><td>U(0.05, 0.8)</td></tr><tr><td></td><td>X = L<sub>X</sub> with C<sub>SS,L<sub>X</sub></sub> = 2</td></tr><tr><td>A<sub>L<sub>X</sub></sub>/ergs</td><td>U(10<sup>41</sup>, 10<sup>46</sup>)</td></tr><tr><td></td><td>X = L<sub>b</sub> with C<sub>SS,L<sub>b</sub></sub> = 7/3</td></tr><tr><td>A<sub>L<sub>b</sub></sub>/ergs</td><td>U(10<sup>41</sup>, 10<sup>44</sup>)</td></tr><tr><td>A<sub>T<sub>X</sub></sub>/keV</td><td>X = T<sub>X</sub> with C<sub>SS,T<sub>X</sub></sub> = 2/3</td></tr><tr><td>A<sub>M<sub>g</sub></sub>/M<sub>☉</sub></td><td>U(0.5, 25)</td></tr><tr><td>A<sub>M<sub>g</sub></sub>/M<sub>☉</sub></td><td>X = M<sub>g</sub> with C<sub>SS,M<sub>g</sub></sub> = 0</td></tr><tr><td>A<sub>Y<sub>X</sub></sub>/keV·M<sub>☉</sub></td><td>U(5 × 10<sup>11</sup>, 5 × 10<sup>13</sup>)</td></tr><tr><td></td><td>X = Y<sub>X</sub> with C<sub>SS,Y<sub>X</sub></sub> = 2/3</td></tr><tr><td></td><td>U(10<sup>12</sup>, 10<sup>14</sup>)</td></tr><tr><td></td><td>Correlated scatter</td></tr><tr><td>ρ<sub>WL,X</sub>, ρ<sub>WL,η</sub>, ρ<sub>X,η</sub></td><td>U(-0.9, 0.9)</td></tr><tr><td></td><td>Miscentering</td></tr><tr><td>f<sub>mis</sub></td><td>N(0.54, 0.02<sup>2</sup>)</td></tr><tr><td>σ<sub>mis</sub></td><td>N(0.61, 0.03<sup>2</sup>)</td></tr><tr><td></td><td>Cosmological parameters</td></tr><tr><td>Ω<sub>M</sub></td><td>N(0.3, 0.016<sup>2</sup>)</td></tr><tr><td>σ<sub>8</sub></td><td>N(0.8, 0.014<sup>2</sup>)</td></tr><tr><td>H<sub>0</sub></td><td>N(70, 5.6<sup>2</sup>)</td></tr></tbody></table>

**Notes.** 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<sub>WL</sub>-M<sub>500</sub>-z relation describe the cases of modeling the lensing profile with R > 0.2h<sup>-1</sup>Mpc and R > 0.5h<sup>-1</sup>Mpc, respectively.

where N(η|z<sub>cl</sub>, p) and N(g<sub>+</sub>, η|z<sub>cl</sub>, p) denote the numbers of clusters observed with the count rate η and the observable set {g<sub>+</sub>, η}, respectively; n(M<sub>500</sub>, z<sub>cl</sub>) is the halo mass function evaluated at the true cluster mass M<sub>500</sub> at the cluster redshift z<sub>cl</sub> using the Bocquet et al. (2016) fitting formula with the inclusion of baryons; P(η|M<sub>500</sub>, z<sub>cl</sub>, p) and P(g<sub>+</sub>, η|M<sub>500</sub>, z<sub>cl</sub>, p) describe the probabilities of observing η and {g<sub>+</sub>, η}, respectively, for the cluster with the true mass M<sub>500</sub> at the redshift z<sub>cl</sub> given the parameter set p. In Eq. (50), 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(η|M<sub>500</sub>, z<sub>cl</sub>, p) in the denominator of Eq. (50) is calculated as

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

where the term P(η|M<sub>500</sub>, z<sub>cl</sub>, p) describes the log-normal distribution of the pristine count rate η̃ with the mean following Eq. (48) and the intrinsic scatter σ<sub>η</sub>. This characterizes the intrinsic distribution of the pristine count rate while including the observed bias b<sub>η</sub> calibrated by the simulations. The term P(η|η̃) accounts for the measurement uncertainty of η with respect to η̃.

The probability P(g<sub>+</sub>, η|M<sub>500</sub>, z<sub>cl</sub>, p) reads

$$P(g_+, \eta|M_{500}, z_{cl}, \mathbf{p}) = \iint d\tilde{\eta} dM_{WL} P(g_+|M_{WL}) P(\eta|\tilde{\eta}) P(M_{WL}, \tilde{\eta}|M_{500}, z_{cl}, \mathbf{p}), \quad (52)$$

in which the term P(g<sub>+</sub>|M<sub>WL</sub>) accounts for the measurement uncertainty of the observed g<sub>+</sub> given the weak-lensing mass M<sub>WL</sub>, which is evaluated using Eq. (39). The term P(M<sub>WL</sub>, η̃|M<sub>500</sub>, z<sub>cl</sub>, p) then describes a joint distribution of M<sub>WL</sub> and η̃ characterized by both the intrinsic scatter of these two quantities and the correlated scatter between them given the underlying true cluster mass M<sub>500</sub>. That is, P(M<sub>WL</sub>, η̃|M<sub>500</sub>, z<sub>cl</sub>, p) follows a multivariate log-normal distribution of (M<sub>WL</sub>, η̃) characterized by an intrinsic covariance matrix,

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

where ρ<sub>WL,η</sub> is the correlation coefficient between M<sub>WL</sub> and η̃.

In this way, the joint likelihood of a sample comprising N<sub>cl</sub> eFEDS clusters, given the set of observables {g<sub>+,i</sub>, η<sub>i</sub>}<sub>i=1</sub><sup>N<sub>cl</sub></sup>, is calculated as

$$\mathcal{L}(\mathbf{p}) = \prod_{i=1}^{N_{\mathrm{cl}}} P(g_{+,i}|\eta_i, \mathbf{p}). \quad (54)$$

In this framework, we have a parameter set p consisting of 19 free parameters, that is,

$$\mathbf{p} = \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}},$$

$$\mathbf{p}_{\eta} = \{A_{\eta}, B_{\eta}, \delta_{\eta}, \gamma_{\eta}, \sigma_{\eta}\},$$

$$\mathbf{p}_{\mathrm{F}} = \{A_{\mathrm{F}}, B_{\mathrm{F}}, \gamma_{\mathrm{F}}, \sigma_{\mathrm{F}}\},$$

$$\mathbf{p}_{\mathrm{WL}} = \{A_{\mathrm{WL}}, B_{\mathrm{WL}}, \gamma_{\mathrm{WL}}, \sigma_{\mathrm{WL}}\},$$

$$\mathbf{p}_{\mathrm{mis}} = \{f_{\mathrm{mis}}, \sigma_{\mathrm{mis}}\},$$

$$\mathbf{p}_{\mathrm{c}} = \{\Omega_{\mathrm{M}}, H_0, \sigma_{\mathrm{g}}\}, \quad (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.

### 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,

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

where $P(\mathcal{X}, g_+, \eta | M_{500}, z_{\mathrm{cl}}, \mathbf{p})$ 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(\mathcal{X}, g_+, \eta | M_{500}, z_{\mathrm{cl}}, \mathbf{p}) = \int \int \int d\tilde{\eta} dM_{\mathrm{WL}} d\tilde{\mathcal{X}} \\ P(\mathcal{X}|\tilde{\mathcal{X}}) P(g_+ | M_{\mathrm{WL}}) P(\eta | \tilde{\eta}) P(\tilde{\mathcal{X}}, M_{\mathrm{WL}}, \tilde{\eta} | M_{500}, z_{\mathrm{cl}}, \mathbf{p}). \quad (57) $$

In an analogy to Eq. (52), the term $P(\tilde{\mathcal{X}}, M_{\mathrm{WL}}, \tilde{\eta} | M_{500}, z_{\mathrm{cl}}, \mathbf{p})$ 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_{\mathrm{X}, \eta} \\ \sigma_{\mathrm{WL}} \sigma_{\eta} \rho_{\mathrm{WL}, \eta} & \sigma_{\mathcal{X}} \sigma_{\eta} \rho_{\mathrm{X}, \eta} & \sigma_{\eta}^2 \end{pmatrix}, \quad (58) $$

where $\rho_{\mathcal{X}, \mathcal{Y}}$ denotes the correlation coefficient between the observables X and Y. The term $P(\mathcal{X}|\tilde{\mathcal{X}})$ 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,

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

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

calculated as

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

where $P(\tilde{\mathcal{X}}, \tilde{\eta} | M_{500}, z_{\mathrm{cl}}, \mathbf{p})$ 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}. \quad (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. (50).

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

$$ \mathcal{L}(\mathbf{p}) = \prod_{i=1}^{n} P(\mathcal{X}_i, g_{+,i} | \eta_i, \mathbf{p}) \\ \times \prod_{j=1}^{m} P(g_{+,j} | \eta_j, \mathbf{p}) \times \prod_{k=1}^{n} P(\mathcal{X}_k | \eta_k, \mathbf{p}), \quad (62) $$

where $i$ runs over the clusters with the available observables of $\{\mathcal{X}, g_+, \eta\}$, $j$ runs over the clusters with only observed shear profiles but $\mathcal{X}$, and $k$ runs over the clusters with only $\mathcal{X}$ but $g_+$. In total, we have a parameter set $\mathbf{p}$ consisting of 26 free parameters, i.e.,

$$ \begin{aligned} \mathbf{p} &= \mathbf{p}_{\eta} \cup \mathbf{p}_{\mathrm{F}} \cup \mathbf{p}_{\mathcal{X}} \cup \mathbf{p}_{\mathrm{WL}} \cup \mathbf{p}_{\mathrm{mis}} \cup \mathbf{p}_{\rho} \cup \mathbf{p}_{\mathrm{c}}, \\ \mathbf{p}_{\mathcal{X}} &= \{A_{\mathcal{X}}, B_{\mathcal{X}}, \delta_{\mathcal{X}}, \gamma_{\mathcal{X}}, \sigma_{\mathcal{X}}\}, \\ \mathbf{p}_{\rho} &= \{\rho_{\mathrm{WL},\mathcal{X}}, \rho_{\mathrm{WL},\eta}, \rho_{\mathcal{X},\eta}\}, \end{aligned} \quad (63) $$

where $\mathbf{p}_{\eta}$, $\mathbf{p}_{\mathrm{F}}$, $\mathbf{p}_{\mathrm{WL}}$, $\mathbf{p}_{\mathrm{mis}}$, and $\mathbf{p}_{\mathrm{c}}$ are defined in Eq. (55); $\mathbf{p}_{\mathcal{X}}$ describes the $\mathcal{X}-M_{500}-z$ relation; $\mathbf{p}_{\rho}$ denotes the intrinsic correlation coefficients among the three observables.

It is worth noting that we expect a degeneracy between the parameters of the scaling relations obtained in a forward modeling. This is because we are constraining the scaling relations between the true cluster mass $M_{500}$ and the multiple observables $\{\mathcal{X}, \eta, g_+\}$, where the latter is fixed in the modeling. In this way, the posteriors of $M_{500}$ and the scaling relation parameters are simultaneously obtained by exploring the likelihood space. As a result, a change in the posterior of the normalization of one relation, for example, will coherently change those of the other scaling relations given a set of fixed observables. This is especially true for the parameters between the $\eta-M_{500}-z$ and $\mathcal{X}-M_{500}-z$ relations, because the degeneracy with the weak-lensing mass will be broken by applying informative priors on the parameters of the $M_{\mathrm{WL}}-M_{500}-z$ relation (see Sect. 4.2). We see below that the strongest degeneracy occurs in the case of $\mathcal{X} = M_g$ (see Sect. 6.4). Despite the degeneracy, we stress that the inclusion of the $\eta-M_{500}-z$ relation in modeling the $\mathcal{X}-M_{500}-z$ relations is still needed, because the count rate is used as the X-ray mass proxy. By including the $\eta-M_{500}-z$ relation, we fully account for the parameter degeneracy in constraining the $\mathcal{X}-M_{500}-z$ relation.

It is also important to make one remark: As in Eq. (62), we use all available 313 shear profiles in modeling the $\mathcal{X}-M_{500}-z$ relations, such that the constraining power from the weak-lensing mass calibration essentially remains the same for all cases of $\mathcal{X} = \{L_X, L_b, T_X, M_g, Y_X\}$. Even in the modeling of $T_X$ and $Y_X$ scaling
<table><caption>**Table 2.** Parameter constraints on the X-ray scaling relations.</caption><thead><tr><th>Model</th><th>A<sub>η</sub></th><th>B<sub>η</sub></th><th>δ<sub>η</sub></th><th>γ<sub>η</sub></th><th>σ<sub>η</sub></th><th>A<sub>X</sub></th><th>B<sub>X</sub></th><th>δ<sub>X</sub></th><th>γ<sub>X</sub></th><th>σ<sub>X</sub></th><th>ρ<sub>WL,X</sub></th><th>ρ<sub>WL,η</sub></th><th>ρ<sub>X,η</sub></th><th>N<sub>η,g</sub></th><th>N<sub>X,η,g</sub></th><th>N<sub>X,η</sub></th></tr></thead><tbody><tr><td>WL Mass Calibration</td><td>0.124<sup>+0.022</sup><sub>-0.019</sub></td><td>1.58<sup>+0.17</sup><sub>-0.14</sub></td><td>1.0<sup>+1.0</sup><sub>-1.4</sub></td><td>-0.44<sup>+0.81</sup><sub>-0.85</sub></td><td>0.301<sup>+0.089</sup><sub>-0.078</sub></td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.01<sup>+0.39</sup><sub>-0.55</sub></td><td>-</td><td>313</td><td>-</td></tr><tr><td>WL Mass Calibration (R > 0.2 h<sup>-1</sup>Mpc)</td><td>0.124<sup>+0.023</sup><sub>-0.021</sub></td><td>1.54<sup>+0.17</sup><sub>-0.14</sub></td><td>0.9<sup>+1.0</sup><sub>-1.1</sub></td><td>-0.04<sup>+0.82</sup><sub>-0.84</sub></td><td>0.294<sup>+0.081</sup><sub>-0.083</sub></td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.22<sup>+0.36</sup><sub>-0.50</sub></td><td>-</td><td>313</td><td>-</td></tr><tr><td>The L<sub>X</sub>-M<sub>500c</sub>-z Relation + WL Mass Calibration</td><td>0.124<sup>+0.021</sup><sub>-0.019</sub></td><td>1.50<sup>+0.15</sup><sub>-0.14</sub></td><td>0.3<sup>+1.0</sup><sub>-1.3</sub></td><td>-0.50<sup>+0.74</sup><sub>-0.84</sub></td><td>0.337<sup>+0.054</sup><sub>-0.069</sub></td><td>3.36<sup>+0.53</sup><sub>-0.49</sub></td><td>1.44<sup>+0.14</sup><sub>-0.11</sub></td><td>-0.07<sup>+1.26</sup><sub>-0.79</sub></td><td>-0.51<sup>+0.93</sup><sub>-0.75</sub></td><td>0.120<sup>+0.138</sup><sub>-0.060</sub></td><td>0.24<sup>+0.38</sup><sub>-0.067</sub></td><td>-0.21<sup>+0.41</sup><sub>-0.31</sub></td><td>0.35<sup>+0.29</sup><sub>-0.45</sub></td><td>0</td><td>313</td><td>121</td></tr><tr><td>The L<sub>b</sub>-M<sub>500c</sub>-z Relation + WL Mass Calibration</td><td>0.127<sup>+0.022</sup><sub>-0.018</sub></td><td>1.55<sup>+0.16</sup><sub>-0.14</sub></td><td>0.2<sup>+1.3</sup><sub>-1.0</sub></td><td>-0.67<sup>+0.99</sup><sub>-0.80</sub></td><td>0.305<sup>+0.062</sup><sub>-0.059</sub></td><td>9.2<sup>+1.6</sup><sub>-1.3</sub></td><td>1.59 ± 0.14</td><td>0.2<sup>+1.3</sup><sub>-1.1</sub></td><td>-0.45<sup>+1.00</sup><sub>-0.86</sub></td><td>0.102<sup>+0.143</sup><sub>-0.043</sub></td><td>0.61<sup>+0.26</sup><sub>-0.58</sub></td><td>-0.31<sup>+0.40</sup><sub>-0.34</sub></td><td>0.38<sup>+0.28</sup><sub>-0.67</sub></td><td>0</td><td>313</td><td>121</td></tr><tr><td>The T<sub>X</sub>-M<sub>500c</sub>-z Relation + WL Mass Calibration</td><td>0.130<sup>+0.020</sup><sub>-0.021</sub></td><td>1.61<sup>+0.17</sup><sub>-0.14</sub></td><td>0.6<sup>+1.2</sup><sub>-1.1</sub></td><td>-0.61<sup>+0.88</sup><sub>-0.90</sub></td><td>0.288<sup>+0.067</sup><sub>-0.089</sub></td><td>3.27<sup>+0.26</sup><sub>-0.31</sub></td><td>0.65 ± 0.11</td><td>-0.02<sup>+0.66</sup><sub>-0.70</sub></td><td>-1.03<sup>+0.54</sup><sub>-0.75</sub></td><td>0.069<sup>+0.061</sup><sub>-0.014</sub></td><td>0.47<sup>+0.36</sup><sub>-0.59</sub></td><td>0.00<sup>+0.27</sup><sub>-0.57</sub></td><td>0.65<sup>+0.22</sup><sub>-0.59</sub></td><td>266</td><td>47</td><td>17</td></tr><tr><td>The M<sub>b</sub>-M<sub>500c</sub>-z Relation + WL Mass Calibration</td><td>0.114<sup>+0.019</sup><sub>-0.017</sub></td><td>1.49<sup>+0.15</sup><sub>-0.16</sub></td><td>1.56<sup>+0.85</sup><sub>-1.12</sub></td><td>0.16<sup>+0.80</sup><sub>-0.85</sub></td><td>0.361<sup>+0.044</sup><sub>-0.061</sub></td><td>1.08 ± 0.13</td><td>1.190<sup>+0.099</sup><sub>-0.118</sub></td><td>0.40<sup>+0.78</sup><sub>-0.70</sub></td><td>0.32<sup>+0.59</sup><sub>-0.61</sub></td><td>0.074<sup>+0.063</sup><sub>-0.019</sub></td><td>0.08<sup>+0.60</sup><sub>-0.43</sub></td><td>-0.30<sup>+0.32</sup><sub>-0.38</sub></td><td>0.63<sup>+0.23</sup><sub>-0.49</sub></td><td>0</td><td>313</td><td>121</td></tr><tr><td>The Y<sub>X</sub>-M<sub>500c</sub>-z Relation + WL Mass Calibration</td><td>0.128<sup>+0.020</sup><sub>-0.015</sub></td><td>1.60<sup>+0.18</sup><sub>-0.14</sub></td><td>0.73<sup>+1.39</sup><sub>-0.79</sub></td><td>-0.36<sup>+0.75</sup><sub>-0.81</sub></td><td>0.224<sup>+0.070</sup><sub>-0.068</sub></td><td>6.9<sup>+1.5</sup><sub>-1.2</sub></td><td>1.78<sup>+0.24</sup><sub>-0.20</sub></td><td>-1.1<sup>+1.9</sup><sub>-1.2</sub></td><td>-1.43<sup>+1.20</sup><sub>-0.93</sub></td><td>0.106<sup>+0.171</sup><sub>-0.047</sub></td><td>0.12<sup>+0.53</sup><sub>-0.50</sub></td><td>-0.44<sup>+0.50</sup><sub>-0.36</sub></td><td>0.794<sup>+0.079</sup><sub>-0.344</sub></td><td>266</td><td>47</td><td>17</td></tr></tbody></table>

**Notes.** The first (second) row contains the result of the weak-lensing mass calibration, which is to simultaneously model the count rate $\eta$ and the shear profile $g_+$ in a radial range of $R > 0.5h^{-1}$Mpc ($R > 0.2h^{-1}$Mpc). The following rows present the joint modeling of the follow-up X-ray scaling relation and the weak-lensing mass calibration. The first column records the type of modeling. The second to sixth columns contain the best-fit parameters describing the $\eta-M_{500}-z$ relation. The seventh to eleventh show the parameter constraints of the target $X-M_{500}-z$ relation. The twelfth to fourteenth columns record the correlation coefficient in the intrinsic scatter among the underlying mass proxies. In each modeling, the number $N_{\eta, WL}$ ($N_{X, \eta, WL}$, $N_{X, \eta}$) of clusters used with the observable set of $\{\eta, g_+\}$ ($\{X, \eta, g_+\}$, $\{X, \eta\}$) is shown in the fifteenth, sixteenth, and seventeenth columns.

relations, where only 64 clusters are present with temperature measurements, the relative errors on the normalization ($A_{T_X}$ and $A_{Y_X}$) are much better than what would be obtained with the first term in Eq. (62) alone, which is due to the fact that we have fully used the available data sets in two surveys. The numbers of clusters used in each modeling with the combination of the data sets are presented in the last three columns of Table 2.

## 5.4. Statistical inference

Our statistical inference can be written in a generic form as follows. Letting **p** be a set of parameters describing the modeling and $\mathcal{D}$ a data vector, the posterior $P(\mathbf{p}|\mathcal{D})$ of **p** given the observed data can be written, through Bayes' theorem, as

$$P(\mathbf{p}|\mathcal{D}) \propto \mathcal{L}(\mathcal{D}|\mathbf{p}) \cdot \mathcal{P}(\mathbf{p}), \quad (64)$$

where $\mathcal{L}(\mathcal{D}|\mathbf{p})$ is the likelihood of observing $\mathcal{D}$ given the parameter **p**, and $\mathcal{P}(\mathbf{p})$ is the prior on **p**. The parameter space is explored by the Affine Invariant MCMC algorithm, which is carried out by emcee (Foreman-Mackey et al. 2013, 2019).

The following priors are adopted in the modeling: The constraints on the parameters describing the bias in the observed count rate and the weak-lensing mass, as in Eqs. (40) and (42), respectively, are used as the priors on these parameters. This effectively accounts for the bias in the observed quantities using the simulation calibration (see Sect. 4). When constraining the $\eta-M_{500}-z$ relation, we adopt a Gaussian prior on the log-normal scatter $\sigma_\eta$ derived from the simulation. Specifically, this includes two factors: First, the intrinsic scatter of the true count rate within each mass and redshift bin in the middle panel of Fig. 6 is quantified at a level of ≈22%. Second, the scatter of the residual in $\ln b_\eta$ between the model and the simulated data is at a level of ≈8%, as evaluated using all mass and redshift bins in the right panel of Fig. 6. These two lead to the intrinsic scatter of $\eta$ at a level of ≈30% in total. As a result, we use a Gaussian prior, $N(0.3, 0.08^2)$ with a dispersion of 0.08 on $\sigma_\eta$. We note that the scatter of the modeling residual at a level of 8% is used as the dispersion of the Gaussian prior as a conservative estimate.

When modeling other X-ray observable-to-mass-and-redshift relations, we adopt uniform priors on the parameters. The Gaussian priors on the parameters, $f_{\text{mis}}$ and $\sigma_{\text{mis}}$, are adopted following the constraints of the miscentering, as quantified in Sect. 3.7. We adopt the Gaussian priors on the cosmological parameters, i.e., $\Omega_M$, $\sigma_8$, and $H_0$, based on the

forecast from Grandis et al. (2019) if including the ground-based weak-lensing mass calibration. The adopted priors are summarized in Table 1.

# 6. Results and discussion

We first present and discuss the mass calibration in Sect. 6.1; we then turn to the constraint of the follow-up X-ray scaling relations in Sects. 6.2–6.5.

## 6.1. The $\eta-M_{500}-z$ relation and the cluster mass $M_{500}$

The count rate-to-mass-and-redshift ($\eta-M_{500}-z$) relation is obtained by simultaneously modeling the observed count rate $\eta$ and shear profile $g_+$, as detailed in Sect. 5.2. The parameter constraints are shown in Fig. 8 and listed in Table 2. We only present the parameters of the $\eta-M_{500}-z$ relation and correlated scatter in Table 2, as others are largely following the adopted Gaussian priors.

As clearly seen in Fig. 8, there is no statistical difference between the modeling of shear profiles at $R > 0.2h^{-1}$Mpc and $R > 0.5h^{-1}$Mpc. This consistency is not trivial: as we adopt the different simulation-calibrated priors on the parameters of $\{A_{\text{WL}}, B_{\text{WL}}, \gamma_{\text{WL}}, \sigma_{\text{WL}}\}$ for the cases of $R > 0.2h^{-1}$Mpc and $R > 0.5h^{-1}$Mpc (see Sect. 4.2), the consistent results suggest that the simulation-calibrated $M_{\text{WL}}-M_{500}-z$ relation accurately accounts for the weak-lensing mass bias, especially at the cluster core. Hereafter, we focus on the results of the modeling with $R > 0.5h^{-1}$Mpc, as the default analysis of the mass calibration in this work.

The resulting count-rate-to-mass-and-redshift relation is obtained as

$$ \left\langle \ln \left( \frac{\eta}{\text{counts/sec}} \middle| M_{500} \right) \right\rangle = \ln \left( 0.124^{+0.022}_{-0.019} + \left[ (1.58^{+0.17}_{-0.14}) + (1.0^{+1.0}_{-1.4}) \ln \left( \frac{1+z}{1+z_{\text{piv}}} \right) \right] \times \ln \left( \frac{M_{500}}{M_{\text{piv}}} \right) + 2 \times \ln \left( \frac{E(z)}{E(z_{\text{piv}})} \right) - 2 \times \ln \left( \frac{D_{\text{L}}(z)}{D_{\text{L}}(z_{\text{piv}})} \right) + (-0.44^{+0.81}_{-0.85}) \times \ln \left( \frac{1+z}{1+z_{\text{piv}}} \right) + \ln (b_{\eta}(M_{500}, z)) \right\rangle, \quad (65) $$
<figure><img src="image_9.png" /><figcaption>**Fig. 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}$.</figcaption></figure>

with the log-normal scatter of $\sigma_{\eta} = 0.301^{+0.089}_{-0.078}$. The last term, $\ln(b_{\eta}(M_{500}, z))$, characterizes the bias in the observed count rate with respect to the true count rate through the simulation calibration, as described in Sect. 4.1. By subtracting $\ln(b_{\eta}(M_{500}, z))$ from Eq. (65), the remaining terms describe the true count-rate-to-mass-and-redshift relation.

With the derived $\eta-M_{500}-z$ relation, we estimate the lensing-calibrated posterior of the cluster mass $M_{500}$ for each cluster at the redshift $z$ given the observed count rate $\eta$, as

$$P(M_{500}|\eta, z, \mathbf{p}) \propto P(\eta|M_{500}, z, \mathbf{p}) P(M_{500}|z, \mathbf{p}), \quad (66)$$

where $\mathbf{p}$ is the best-fit parameter of the modeling, and

$P(M_{500}|z, \mathbf{p})$ is the normalized mass function evaluated at the cluster redshift. We then randomly sample a mass estimate from the resulting posterior as the ensemble mass for individual clusters (see more applications in Bocquet et al. 2019; Chiu et al. 2020a). In this way, the individual cluster mass is statistically inferred from the ensemble population modeling. The same method but with a different terminology, the “mass forecasting”, has been used to estimate the mass of clusters based on the X-ray temperature in the XXL survey (Umetsu et al. 2020). We show the ensemble mass of each cluster in Fig. 9. The eFEDS sample spans a mass range between $\approx 10^{13}h^{-1}M_{\odot}$ and $\approx 10^{14.8}h^{-1}M_{\odot}$ in terms of $M_{500}$, of which low-mass systems ($M_{500} \lesssim 10^{14}h^{-1}M_{\odot}$) are mostly at low redshift ($z \lesssim 0.4$). The typical uncertainty in
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Cluster mass $M_{500}$ and redshift $z$ of each cluster color coded by the quantity $f_{\text{cont}}$. The mass is randomly sampled from the posterior of $M_{500}$ as an ensemble mass (see Eq. (66)).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 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_{\text{piv}} = 0.35$ (pivotal mass $M_{\text{piv}} = 1.4 \times 10^{14} h^{-1} M_{\odot}$).</figcaption></figure>

the cluster mass is at a level of ≈27% (≈32%, ≈24%) for a cluster with $M_{500} \approx 10^{14} h^{-1} M_{\odot}$ (≈10$^{13.6} h^{-1} M_{\odot}$, ≈10$^{14.4} h^{-1} M_{\odot}$). The mass $M_{500}$ of the eFEDS clusters is presented in Appendix C.

We show the observed count rate $\eta$ as a function of the cluster mass $M_{500}$ and redshift $z$, together with this best-fit model, in Fig. 10. When showing the mass-trend in the left panel, we divide the count rate by the redshift-dependent quantities at the cluster redshift. In this way, the count rate $\eta$ is re-normalized to the pivotal redshift. Similarly, we re-normalize the count rate to the pivotal mass when showing the redshift trend in the right panel of Fig. 10. As seen, the best-fit model provides a good description of the eFEDS clusters. Based on the result, the true count rate observed by eROSITA is $0.124^{+0.022}_{-0.019}$ counts per second for a cluster with the pivotal mass $M_{\text{piv}} = 1.4 \times 10^{14} h^{-1} M_{\odot}$ at the

pivotal redshift $z_{\text{piv}} = 0.35$. Moreover, the true count rate scales with mass and redshift as $\propto M_{500}^{1.58^{+0.17}_{-0.14}}$ and $\propto (1+z)^{-0.44^{+0.81}_{-0.85}}$, respectively, without significant cross-scaling ($\delta_{\eta} = 1.0^{+1.0}_{-1.4}$). We stress that there is no significant deviation from the expected redshift trend of the true count rate ($\gamma_{\eta} = -0.44^{+0.81}_{-0.85}$), which is consistent with what we observe in other X-ray scaling relations for the eFEDS sample. We continue our discussion of this perspective in Sects. 6.2–6.5.

For both the lensing modeling with and without the cluster core, the correlated scatter $\rho_{\text{WL},\eta}$ between the count rate and weak-lensing mass is constrained as statistically consistent with zero. This suggests that the true count rate and weak-lensing-inferred mass are nearly independent of each other. This is in line with the picture that the source of the intrinsic scatter in X-rays
<figure><img src="image_12.png" /><figcaption>Fig. 11. Stacked weak-lensing shear profile (black points), and the best-fit lensing models with (blue-dashed line) and without (red-dotted line) the miscentering component. The open circles indicate the inner regime of $R < 0.5 h^{-1} Mpc$, which are not used in the fiducial lensing model in this work.</figcaption></figure>

and in the optical is not strongly correlated. The former is related to various activities in the ICM caused by either the output of energetic sources or the thermal emission of clusters, while the latter is largely attributed to the halo orientation, concentration, and the projection of line-of-sight structures.

Finally, we show the stacked shear profile and the best-fit models in Fig. 11. To derive the stacked profile, we properly weight the lensing shear profiles $g_+$ of individual clusters at each radial bin. We use the diagonal term<fnref n="6" /> of the lensing covariance matrix (see Eq. (14)) as the radius-dependent weight for each cluster. We repeat the same stacking procedure for both observed shear profiles and those predicted by the best-fit parameters, as shown by the black points and blue-dashed line, respectively. We additionally derive the best-fit model without the miscentering by setting $f_{\text{mis}} = 0$, as shown by the red dotted line. It is seen that the best-fit model with the miscentering provides a good description of the stacked shear profile, and that the miscentered model comes into agreement with that of perfectly centered sources at $R \ge 0.5 h^{-1} Mpc$. This picture is also supported by the consistent constraints between the lensing modeling with and without the cluster core. However, we note that the best-fit miscentered model (blue-dashed) is slightly lower than the stacked shear profile (open circles) in the inner radial range $R < 0.5 h^{-1} Mpc$ at a level of $\le 1\sigma$, indicating that the miscentering of the X-ray clusters is over-corrected. This is not surprising, as the modeling of the miscentering is based on the offset distribution between the X-ray and optical centers, instead of the former and the true center, in which case the offset distribution is prone to the miscentering of the optical centers (see more discussions

in Sect. 3.7). We stress that the modeling of the miscentering is subdominant to our results, and that the use of simulations to calibrate the X-ray miscentering is clearly needed in a future work.

## 6.2. The $L_X-M_{500}-z$ relation

The luminosity-to-mass-and-redshift relation has been intensively studied (e.g., Reiprich & Böhringer 2002; Maughan 2007; Vikhlinin et al. 2009; Arnaud et al. 2010) and has been shown to have larger scatter at fixed mass than other X-ray mass proxies. This scatter is suggested to be related to complex mergers (Ritchie & Thomas 2002; Poole et al. 2008) or the nature of cluster cores, for example, the presence of a cool core (Pratt et al. 2009). Despite the large scatter, the X-ray luminosity is relatively easy to obtain with a high S/N, and thus becomes a convenient mass proxy for clusters with only a few tens of photon counts. Moreover, the luminosity has a direct link to the observed count rate $\eta$, which is used as the mass proxy in this work. Therefore, the luminosity provides straightforward access to the physical properties of eFEDS clusters.

The luminosity-to-mass-and-redshift relations predicted from the self-similar model (Kaiser 1986) are

$$L_X \propto M_{500} \times E(z)^2,$$

and

$$L_b \propto M_{500}^{4/3} \times E(z)^{7/3},$$

for the soft-band and bolometric luminosities, respectively. We refer readers to Böhringer et al. (2012) and Giodini et al. (2013) for a more recent and complete discussion of the self-similar scaling among X-ray observables. Here, we assume that the emissivity in the soft band does not depend on the temperature (see Fig. 1 in Mohr et al. 1999), such that the mass and redshift-dependence of the temperature does not enter into the soft-band luminosity-to-mass-and-redshift relation.

The parameter constraints from the joint modeling of the $L_X-M_{500}-z$ ($L_b-M_{500}-z$) relation and the weak-lensing mass calibration are shown in Fig. 12 (Fig. 13). The results are listed in Table 2. In this work, the $L_X-M_{500}-z$ relation of the eFEDS sample is derived as

$$ \left( \ln \left( \frac{L_X}{\text{ergs}} \middle| M_{500} \right) \right) = \ln(3.36^{+0.53}_{-0.49}) + \ln(10^{43}) \\ + \left[ (1.44^{+0.14}_{-0.13}) + (-0.07^{+1.26}_{-0.79}) \ln \left( \frac{1+z}{1+z_{\text{piv}}} \right) \right] \times \ln \left( \frac{M_{500}}{M_{\text{piv}}} \right) \\ + 2 \times \ln \left( \frac{E(z)}{E(z_{\text{piv}})} \right) + (-0.51^{+0.93}_{-0.75}) \times \ln \left( \frac{1+z}{1+z_{\text{piv}}} \right), \quad (67) $$

with the log-normal scatter of $\sigma_{L_X} = 0.120^{+0.138}_{-0.060}$ at fixed mass. Meanwhile, the $L_b-M_{500}-z$ relation is derived as

$$ \left( \ln \left( \frac{L_b}{\text{ergs}} \middle| M_{500} \right) \right) = \ln(9.2^{+1.6}_{-1.3}) + \ln(10^{43}) \\ + \left[ (1.59 \pm 0.14) + (0.2^{+1.3}_{-1.1}) \ln \left( \frac{1+z}{1+z_{\text{piv}}} \right) \right] \times \ln \left( \frac{M_{500}}{M_{\text{piv}}} \right) \\ + \frac{7}{3} \times \ln \left( \frac{E(z)}{E(z_{\text{piv}})} \right) + (-0.45^{+1.00}_{-0.86}) \times \ln \left( \frac{1+z}{1+z_{\text{piv}}} \right), \quad (68) $$

with the log-normal scatter of $\sigma_{L_b} = 0.102^{+0.143}_{-0.043}$. The soft-band (bolometric) luminosity is constrained as $3.36^{+0.53}_{-0.49} \times$

<footnote n="6">⁶ In the case that off-diagonal terms are large, the exact form including the correlation among radial bins must be used, as Eq. (136) in Umetsu (2020). In this work, the off-diagonal terms in the lensing covariance matrix are small, and we stack the measurements only for the purpose of visualizations. Thus, we use the diagonal approximation when producing Fig. 11.</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Parameter constraints in the joint modeling of the $L_X-M_{500}-z$ relation and the weak-lensing mass calibration. The weak-lensing mass calibration consists of the modeling of X-ray count rates and shear profiles at $R > 0.5 \ h^{-1}\mathrm{Mpc}$. The result of the weak-lensing mass calibration is shown in blue, while that of the joint modeling with the $L_X-M_{500}-z$ relation is shown in red. The marginalized posteriors of the parameters and the correlation between them are presented in the same way as in Fig. 8. The adopted priors are indicated by the dotted lines in the diagonal subplots.</figcaption></figure>

$10^{43}$ (9.2$^{+1.6}_{-1.3}$ × $10^{43}$) ergs at the pivotal mass $M_{\text{piv}} = 1.4 \times 10^{14}h^{-1}M_{\odot}$ at the pivotal redshift $z_{\text{piv}} = 0.35$. The intrinsic scatter of the soft-band ($\sigma_{L_X} = 0.120^{+0.138}_{-0.060}$) and bolometric ($\sigma_{L_b} = 0.102^{+0.143}_{-0.043}$) luminosity are statistically consistent with each other. We note that the asymmetric posterior of the scatter is mostly due to the lower bound of the imposed prior.

We discuss our results as follows. First, a mass-trend slope that is steeper than (statistically consistent with) the self-similar prediction is obtained for the $L_X-M_{500}-z$ ($L_b-M_{500}-z$) relation with $B_{L_X} = 1.44^{+0.14}_{-0.13}$ ($B_{L_b} = 1.59 \pm 0.14$), corresponding to a

significance level of ≈3.4σ (≈1.7σ). Second, there is no significant deviation in the redshift trend ($\gamma_{L_X} = -0.51^{+0.93}_{-0.75}$ and $\gamma_{L_b} = -0.45^{+1.00}_{-0.86}$) from the self-similar model ($\gamma = 0$). Finally, no sign of a redshift-dependent mass trend ($\delta_{L_X} = -0.07^{+1.26}_{-0.79}$ and $\delta_{L_b} = 0.2^{+1.3}_{-1.1}$) is revealed for the eFEDS sample. This picture can be better visualized in Figs. 14 and 15 for the $L_X-M_{500}-z$ and $L_b-M_{500}-z$ relations, respectively. In Figs. 14 and 15, we also re-normalize the luminosity to the pivotal mass and redshift, respectively, as done in a similar way in Fig. 10. As seen in Fig. 15, the mass and redshift trends of $L_b$ are
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Parameter constraints in the joint modeling of the $L_b-M_{500}-z$ relation and the weak-lensing mass calibration presented in an identical way to those shown in Fig. 12.</figcaption></figure>

statistically consistent with the self-similar predictions. For $L_X$ in Fig. 14, the mass-trend slope of the best-fit model (gray shaded region) is steeper than that predicted by the self-similar model (red dashed line) in the left panel, while the redshift trends of them show good consistency in the right panel. These are all in broad agreement with the results of Bulbul et al. (2019), despite the fact that their sample of clusters is selected by the SZE signature observed by the SPT. We stress that the SPT sample of Bulbul et al. (2019) only probes massive clusters with $M_{500} \gtrsim 10^{14.3}h^{-1}M_{\odot}$, while the eFEDS sample includes low-mass systems with $M_{500} \lesssim 10^{14}h^{-1}M_{\odot}$. That is, this picture still holds down to a group scale for luminosity-mass-redshift relations. Our result for $L_X$ suggests that (1) the mass-trend slope of eFEDS clusters is steeper than the self-similar prediction down to a mass scale of $M_{500} \approx 10^{13}h^{-1}M_{\odot}$, while the redshift trend is

in good agreement with the self-similarity out to redshift $z \approx 1.3$, and that (2) there is no significant difference in the luminosity-mass-redshift relation between the eROSITA-selected and the SZE-selected clusters.

In terms of correlated scatter, we find no clear correlation between the count rate and the luminosity ($\rho_{L_X,\eta} = 0.35^{+0.29}_{-0.45}$ and $\rho_{L_b,\eta} = 0.38^{+0.28}_{-0.67}$, see also Figs. 12 and 13). Meanwhile, no clear correlation is seen between the luminosity and the weak-lensing mass ($\rho_{WL,L_X} = 0.24^{+0.38}_{-0.67}$ and $\rho_{WL,L_b} = 0.61^{+0.26}_{-0.58}$). No correlation is shown between the count rate and weak-lensing mass either, as consistent with the result in Sect. 5.2. We note that we do expect positive $\rho_{L_X,\eta}$ and $\rho_{L_b,\eta}$, given a direct relation between the count rate and luminosity. However, the constraining power based on the eFEDS sample is generally not strong enough to precisely constrain the intrinsic correlated scatter. A larger
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Observed soft-band luminosity $L_X$ as a function of the cluster mass $M_{500}$ (left panel) and redshift $z$ (right panel). When showing the mass (redshift) trend in the left (right) panel, we re-normalize the observed quantity to the pivotal redshift $z_{piv} = 0.35$ (mass $M_{piv} = 1.4 \times 10^{14} h^{-1} M_{\odot}$) as done in a similar manner in order to produce Fig. 10. The best-fit model and observed clusters are shown in the same way as in Fig. 10. The pink circles indicate the weighted mean of the eFEDS sample.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Observed bolometric luminosity $L_b$ as a function of the cluster mass $M_{500}$ (left panel) and redshift $z$ (right panel). This plot is produced and shown in the same way as in Fig. 14.</figcaption></figure>

sample based on the first-year eROSITA survey will significantly improve these constraints in a future work.

## Comparisons with the literature

The $L_X-M_{500}-z$ relation of the eFEDS sample has a mass trend ($B_{L_X} = 1.44_{-0.13}^{+0.14}$) that is broadly consistent with the literature, including Pratt et al. (2009) who found a slope of $1.62 \pm 0.12$ using an X-ray-selected sample of clusters (REXCESS) at redshift $z \approx 0.1$, Vikhlinin et al. (2009) with a slope of $1.61 \pm 0.14$ from a cluster sample out to $z \approx 0.5$ using high-quality Chandra data; Andersson et al. (2011) who derived a slope of $1.16 \pm 0.20$ using a sample of 15 SPT-selected clusters with the redshift scaling fixed to $E(z)^{1.85}$; Lovisari et al. (2015) in which a slope of $1.39 \pm 0.05$ was obtained from a sample of local clusters and groups down to a mass scale of $\approx 10^{13} M_{\odot}$; Mantz et al. (2016) where a slope of $1.35 \pm 0.06$ is obtained for a sample of $\approx 200$

clusters selected from the ROSAT survey; and Bulbul et al. (2019) who derived a slope of $1.93_{-0.18}^{+0.19}$ using 59 SPT-selected clusters out to redshift $z \approx 1.3$ observed with the XMM-Newton telescope. Schellenberger & Reiprich (2017) obtained a slope of $1.35 \pm 0.07$ for an X-ray-flux-limited sample based on a joint data set from Chandra and ROSAT. Except for Bulbul et al. (2019), which has a steeper slope than ours at a level of $\approx 2\sigma$, the resulting mass trend of the soft-band luminosity of the eFEDS sample is in agreement with the literature at a level of $\lesssim 1.5\sigma$. We show the comparison of $L_X$ with the literature in the upper-left panel of Fig. 16, in which we additionally include the result of the simulated clusters (golden stars) from the C-Eagle cosmological hydrodynamical simulation (Barnes et al. 2017) and the self-similar model (blue line) with the normalization anchored to the best-fit value of the eFEDS sample. As seen, the mass scaling of the eFEDS sample is consistent with that found by these latter authors, including the simulated clusters, and shows
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Comparison between the eFEDS sample and the results from the literature for $L_X$ (upper-left), $T_X$ (upper-right), $M_g$ (lower-left), and $Y_X$ (lower-right). The eFEDS clusters are shown by the circles color coded by the redshift. The 68% confidence levels of the best-fit models at the pivotal redshift $z_{\text{piv}} = 0.35$ are plotted as the gray shaded regions. The self-similar prediction with the normalization anchored to the best-fit value of the eFEDS sample is shown by the blue line, while the results from the literature are indicated by different colors. The simulated clusters (Barnes et al. 2017) are marked as the golden stars. We apply a correction factor to the cluster mass from the literature to account for the systematic difference in $M_{500}$, following Appendix D.</figcaption></figure>

a steeper mass trend than the self-similar model ($B = 1$) at a level of $\approx 3\sigma$.

In terms of the bolometric luminosity $L_b$, we find that the resulting slope of the eFEDS sample ($B_{L_b} = 1.59 \pm 0.14$) is shallower than that from Pratt et al. (2009) and Bulbul et al. (2019), with $1.90 \pm 0.11$ and $2.12^{+0.23}_{-0.18}$, at a level of $\approx 1.6\sigma$ and $\approx 2.2\sigma$, respectively. It is worth mentioning that a mass scaling of $L_b \propto M_{500}^{1.52\pm0.10}$ was obtained by combining heterogeneous samples in Reichert et al. (2011), which is in good agreement with our results.

We constrain the deviation of $L_X$ and $L_b$ from the self-similar redshift trend as $\gamma_{L_X} = -0.51^{+0.93}_{-0.75}$ and $\gamma_{L_b} = -0.45^{+1.00}_{-0.86}$, respectively. This implies that the redshift trend of $L_X$ and $L_b$ closely follows the self-similar model, which is in excellent agreement with the result from the SPT-selected massive clusters out to a redshift of $z \approx 1.3$ (Bulbul et al. 2019). A consistent result is also reported in Mantz et al. (2016) with $\gamma_{L_X}^{0.65\pm0.38}$ for the soft-band luminosity. Moreover, Vikhlinin et al. (2009) found a redshift scaling of $L_X \propto E(z)^{1.85\pm0.42}$, which is also consistent with the self-similarity ($C_{SS,L_X} = 2$). By homogenizing the different
samples from the literature, Reichert et al. (2011) derived a redshift scaling<fnref n="7" /> of $L_b \propto E(z)^{1.24^{+0.16}_{-0.82}}$ out to a high redshift $z \approx 1.4$, which is significantly shallower than the self-similar prediction and the eFEDS results.

The log-normal scatter of the soft-band (bolometric luminosity is constrained as $\sigma_{L_x} = 0.120^{+0.138}_{-0.060}$ ($\sigma_{L_b} = 0.102^{+0.143}_{-0.043}$) at fixed mass. This is in broad agreement with, but mildly smaller than, the scatter of 0.3 in the soft-band luminosity derived from the simulated clusters in Barnes et al. (2017). Lovisari et al. (2015) and Vikhlinin et al. (2009) derived the scatter of 0.245 and 0.396, respectively, which are both mildly higher than ours in the soft-band luminosity. Meanwhile, a high value of $0.43 \pm 0.03$ was derived in Mantz et al. (2016) for the scatter of the soft-band luminosity. An interesting result from the FABLE simulation (Henden et al. 2019) suggests that the scatter of bolometric luminosity decreases with redshift from $\approx 0.25$ at redshift $z \approx 0$ to $\approx 0.13$ at redshift $z \approx 1$. Based on these comparisons, we conclude that the scatter of the X-ray luminosity in the eFEDS sample is in broad agreement with the literature.

### 6.3. The $T_x-M_{500}-z$ relation

The self-similar prediction of the $T_x-M_{500}-z$ relation reads

$$T_x \propto M_{500}^{3/2} E(z)^{3/2},$$

which is obtained by assuming that galaxy clusters reach the virial condition ($M/R \sim T$) with the conversion between the cluster mass and radius ($M_{500} \propto R_{500}^3 \rho_c$). The X-ray temperature was considered as a reliable mass proxy in early studies of galaxy clusters (Smith et al. 1979; Mohr & Evrard 1997; Mohr et al. 1999; Ikebe et al. 2002); therefore the temperature-to-mass relation was intensively investigated using the X-ray hydrostatic mass (e.g., Finoguenov et al. 2001; Vikhlinin et al. 2006) before the lensing mass became available (e.g., Lieu et al. 2016). In practice, galaxy clusters are not fully virialized and contain nonthermal emissions (O'Hara et al. 2006), resulting in scatter at fixed mass and deviation from self-similarity. Thus, studying the scaling relation of the temperature sheds light on this perspective.

For the eFEDS sample, the $T_x-M_{500}-z$ relation is constrained as

$$ \left( \ln \left( \frac{T_x}{\text{keV}} \middle| M_{500} \right) \right) = \ln \left( 3.27^{+0.26}_{-0.31} \right) + \left[ (0.65 \pm 0.11) + (-0.02^{+0.66}_{-0.70}) \ln \left( \frac{1+z}{1+z_{\text{piv}}} \right) \right] \times \ln \left( \frac{M_{500}}{M_{\text{piv}}} \right) + \frac{2}{3} \times \ln \left( \frac{E(z)}{E(z_{\text{piv}})} \right) + (-1.03^{+0.54}_{-0.75}) \times \ln \left( \frac{1+z}{1+z_{\text{piv}}} \right), \quad (69) $$

with log-normal scatter of $\sigma_{T_x} = 0.069^{+0.061}_{-0.014}$ at fixed mass. The parameter constraints of the $T_x-M_{500}-z$ relation are shown in Fig. 17 and are listed in Table 2. At the pivotal mass $M_{\text{piv}} = 2.5 \times 10^{14} h^{-1} M_{\odot}$ and the pivotal redshift $z_{\text{piv}} = 0.35$, the temperature of the eFEDS sample is constrained as $A_{T_x} = 3.27^{+0.26}_{-0.31}$ keV.

In terms of mass and redshift trends, the resulting $T_x-M_{500}-z$ relation shows (1) a mass trend ($B_{T_x} = 0.65 \pm 0.11$) that is in good agreement with the self-similarity ($B = 2/3$), (2) a redshift

trend ($\gamma_{T_x} = -1.03^{+0.54}_{-0.75}$) that is shallower than, but still statistically consistent with, the self-similar behavior at a significance of $1.7\sigma$, and (3) no clear redshift-dependent mass slope with $\delta_{T_x} = -0.02^{+0.66}_{-0.70}$. The results of a self-similar mass trend and a redshift scaling shallower than the self-similar prediction are not consistent with the picture we see in the scaling relation of the luminosity in Sect. 6.2, where we find a steeper mass trend than the self-similar slope with a nearly zero deviation from the self-similar redshift scaling. This picture can be visualized in Fig. 18, where the mass and redshift trends of the temperature are presented in the left and right panels, respectively.

A possible explanation to the shallower redshift trend could be as follows. In the right panel of Fig. 18, it is seen that the $T_x$ measurements above $\approx 4$ keV contain larger error bars and higher scatter than the low-temperature regime. This is not surprising, given that the capability of the eROSITA telescope to measure high temperatures ($\ge 5$ keV) is limited by the reduction of the sensitivity to the high-energy band at $> 3$ keV (Liu et al. 2022a). Consequently, the noisy measurements decrease the constraining power of the high-temperature eFEDS clusters, which are mainly the high-mass systems at high redshift. As a result, the modeling of the $T_x-M_{500}-z$ relation is weighted toward the low-mass eFEDS clusters or groups at low redshift ($z \le 0.3$). The shallow redshift trend could also be attributed to either the fact that the redshift-dependence of $T_x$ is intrinsically shallower than the self-similar prediction at a group scale, or that the temperature of the eFEDS clusters with the low mass ($M_{500} \le 10^{14} h^{-1} M_{\odot}$) at low redshift ($z \le 0.3$) is higher than the full population on average. Currently, the redshift trend is only marginally shallower than the self-similar prediction at a level of $\approx 1.7\sigma$. However, in the future, deep follow-up observations from either XMM-Newton or Chandra telescopes are clearly warranted for a thorough investigation of the $T_x-M_{500}-z$ relation for eROSITA-detected clusters.

The log-normal intrinsic scatter of $T_x$ for the eFEDS sample is constrained as $\sigma_{T_x} = 0.069^{+0.061}_{-0.014}$, with a correlation coefficient of $\rho_{T_x,\eta} = 0.65^{+0.22}_{-0.59}$ and $\rho_{\text{WL},T_x} = 0.47^{+0.36}_{-0.59}$ with the count rate and weak-lensing mass, respectively. Again, no clear correlated scatter between the count rate and weak-lensing mass is seen ($\rho_{\text{WL},\eta} = 0.00^{+0.27}_{-0.57}$) in the joint modeling of the $T_x-M_{500}-z$ relation and the mass calibration.

### Comparisons with the literature

We compare the $T_x-M_{500}-z$ relation of the eFEDS sample with the literature as follows. In Arnaud et al. (2007), the core-excised temperature is constrained to have a mass-dependent scaling with a slope of $0.58 \pm 0.03$ using a sample of ten nearby clusters; this slope becomes self-similar as $0.67 \pm 0.07$ if restricting the sample to clusters with temperature $T_x \ge 3$ keV (Arnaud et al. 2005). The mass trend from Vikhlinin et al. (2009) was derived with a self-similar slope of $0.65 \pm 0.04$ for the core-excised temperature-to-mass relation using a sample of clusters with $10^{14} \lesssim \frac{M_{500}}{h^{-1}M_{\odot}} \lesssim 10^{15}$ out to redshift $z \approx 0.5$. Based on a combined and local sample from the literature, Reichert et al. (2011) obtained a scaling of $T_x \propto M_{500}^{0.57 \pm 0.03}$, which is statistically consistent with the eFEDS result (within $1\sigma$). In Mahdavi et al. (2013), the core-excised temperature-to-mass relation was derived as $0.51^{+0.42}_{-0.16}$ ($0.70^{+0.11}_{-0.08}$) with respect to the weak-lensing (X-ray hydrostatic) mass using a sample of 50 massive clusters in the redshift range of $0.15 < z < 0.55$. The core-excised temperature-to-mass relation was found to be $0.58 \pm 0.01$ by Lovisari et al. (2015), which is shallower than the self-similar

<footnote n="7">⁷ These authors obtained $M_{500} \propto E(z)^{-0.93^{+0.62}_{-0.12}}$. For the redshift scaling of $L_b$ given a cluster mass, the slope is obtained by multiplying $-\frac{4}{3}$ to $-0.93^{+0.62}_{-0.12}$.</footnote>
<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Parameter constraints in the joint modeling of the $T_X-M_{500}-z$ relation and the weak-lensing mass calibration, presented in the same way as in Fig. 12.</figcaption></figure>

prediction at high significance. In Mantz et al. (2016), the mass trend was found to be $0.63 \pm 0.03$ for ROSAT-selected clusters in the range $10^{14} \lesssim \frac{M_{500}}{h^{-1}M_{\odot}} \lesssim 10^{15}$. Using a sample of 59 SPT clusters with $M_{500} \ge 3 \times 10^{14}h^{-1}M_{\odot}$ in the redshift range $0.2 \lesssim z \lesssim 1.3$, Bulbul et al. (2019) determined the mass trend as $0.79^{+0.09}_{-0.10}$ and $0.81^{+0.09}_{-0.08}$ for core-included and core-excised temperature, respectively. The mass-trend slope of the temperature measured within an aperture with a radius of 300 kpc was derived as $0.56^{+0.12}_{-0.10}$ for a sample of galaxy clusters and groups selected in the XXL survey (Lieu et al. 2016; Pierre et al. 2016); this slope became $0.60^{+0.04}_{-0.05}$ when combining the sample of galaxy groups from the COSMOS field (Kettula et al. 2013) and the sample of massive clusters from Mahdavi et al. (2013). A highly self-similar redshift trend of the temperature, $T_X \propto E(z)^{0.69 \pm 0.05}$, was

also derived in Reichert et al. (2011), which is consistent with the eFEDS result. The latest result from the XXL survey using the S16A HSC lensing data delivered a constraint on the mass trend of $T_X(< 300 \text{ kpc})$ as $0.85 \pm 0.31$ and $0.75 \pm 0.27$ when the redshift scaling is a free parameter and fixed to the self-similar prediction, respectively (Umetsu et al. 2020). In summary, the resulting mass trend of the eFEDS sample ($B_{T_X} = 0.65 \pm 0.11$) is statistically consistent with the literature at a level of $\le 1.5\sigma$. The upper-right panel of Fig. 16 contains the comparison, showing good consistency between the eFEDS sample and the previous work. In addition, we also show the simulated clusters from the C-Eagle cosmological simulation (Barnes et al. 2017) as the golden stars, which nicely follow the trend of the eFEDS sample.
<figure><img src="image_19.png" /><figcaption>Fig. 18. Observed X-ray temperature $T_X$ as a function of the cluster mass $M_{500}$ (left panel) and redshift $z$ (right panel). This plot is produced and shown in the same way as in Fig. 14.</figcaption></figure>

In terms of the redshift trend, the result of the eFEDS sample shows a mild deviation ($\gamma_{T_X} = -1.03^{+0.54}_{-0.75}$) from the self-similar behavior ($T_X \propto E(z)^{2/3}$) at a level of $\approx 1.7\sigma$. This redshift trend is consistent with the results of Bulbul et al. (2019), who obtained a scaling as $\propto E(z)^{2/3}(1+z)^{-0.22^{+0.29}_{-0.35}}$ for core-included temperature out to redshift $z \approx 1.3$. A similar picture is also suggested for the core-excised temperature in Bulbul et al. (2019), in which the redshift scaling is constrained as $E(z)^{4/3}(1+z)^{-0.30^{+0.27}_{-0.28}}$. The latest result from the XXL survey (Umetsu et al. 2020) obtained a redshift scaling of $\propto E(z)^{0.18\pm0.66}$ for the temperature estimated in an aperture with a radius of 300 kpc, which is much shallower than but still statistically consistent with the self-similar behavior, given the large error bar. In summary, the eFEDS sample shows a redshift scaling that is mildly shallower than the self-similar prediction but is still statistically consistent with the previous work.

The log-normal intrinsic scatter of $T_X$ is constrained as $\sigma_{T_X} = 0.069^{+0.061}_{-0.014}$ for the eFEDS sample, which is in broad agreement with the literature: In Bulbul et al. (2019), the scatter of the core-excised and core-included temperature was constrained as $0.13^{+0.05}_{-0.04}$ and $0.18 \pm 0.04$, respectively. Mantz et al. (2016) derived an intrinsic scatter of $0.161 \pm 0.019$ for the ROSAT sample. The sample of simulated clusters from Barnes et al. (2017) suggests a scatter of 0.14 for the core-included temperature. The intrinsic scatter of the core-excised temperature at fixed mass<fnref n="8" /> was suggested to lie in the range between $\approx 0.10$ and $\approx 0.16$ by Vikhlinin et al. (2009). A value of $0.23^{+0.17}_{-0.10}$ for the intrinsic scatter of $T_X$ was suggested by Mahdavi et al. (2013). Similarly, the intrinsic scatter was obtained as $0.30^{+0.12}_{-0.10}$ for the XXL clusters, and was reduced to $0.18^{+0.13}_{-0.11}$ if excluding disturbed clusters (Lieu et al. 2016). The simulation work in Henden et al. (2019) suggested the intrinsic scatter to be at a level of $\approx 0.23$, slightly decreasing toward high redshift.

## 6.4. The $M_g-M_{500}-z$ relation

The self-similar model predicts the $M_g-M_{500}-z$ relation as $M_g \propto M_{500}$,

meaning that the ICM mass fraction, $M_g/M_{500}$, is invariant among clusters at any redshift. This has an important implication: Galaxy clusters are considered as a “closed box”, inside which baryons cannot escape from the potential, meaning that the baryon fraction of clusters is equivalent to that of the Universe. In practice, the deviation from the universal baryon fraction implies an effect from baryonic physics (Lin et al. 2003) or the assembly history of galaxy clusters (Chiu et al. 2016b, 2018a).

The $M_g-M_{500}-z$ relation of the eFEDS sample is derived as

$$ \left\langle \ln \left( \frac{M_g}{M_{\odot}} \middle| M_{500} \right) \right\rangle = \ln(1.08 \pm 0.13) + \ln(10^{13}) + \left[ (1.190^{+0.099}_{-0.118}) + (0.40^{+0.78}_{-0.70}) \ln \left( \frac{1+z}{1+z_{\text{piv}}} \right) \right] \times \ln \left( \frac{M_{500}}{M_{\text{piv}}} \right) + (0.32^{+0.59}_{-0.61}) \times \ln \left( \frac{1+z}{1+z_{\text{piv}}} \right), \quad (70) $$

with a log-normal scatter of $\sigma_{M_{ICM}} = 0.074^{+0.063}_{-0.019}$ at fixed mass. The parameter constraints are shown in Fig. 19 and are listed in Table 2. The ICM mass of the eFEDS clusters is constrained to be $1.08 \pm 0.13 \times 10^{13} M_{\odot}$, corresponding to an ICM mass fraction of $\approx 0.053 \pm 0.007$ (with $h = 0.7$) at the pivotal mass $M_{\text{piv}} = 1.4 \times 10^{14} h^{-1} M_{\odot}$ and redshift $z_{\text{piv}} = 0.35$.

Based on these latter results, the scaling trends in the mass ($B_{M_{ICM}} = 1.190^{+0.099}_{-0.118}$) and redshift ($\gamma_{M_{ICM}} = 0.32^{+0.59}_{-0.61}$) are both steeper than, but still statistically consistent with, the self-similar predictions of $B = 1$ and $\gamma = 0$ at significance levels of $\lesssim 1\sigma$. Moreover, no redshift-dependent mass trend at a significant level is seen with $\delta_{M_{ICM}} = 0.40^{+0.78}_{-0.70}$. That is, our results for the eFEDS clusters suggest that the ICM mass $M_g$ largely follows the self-similar prediction in both the mass and redshift trends with only mild deviations.

The log-normal intrinsic scatter of $M_g$ at fixed cluster mass is constrained as $\sigma_{M_{ICM}} = 0.074^{+0.063}_{-0.019}$ with a mildly positive correlation coefficient $\rho_{M_g,\eta} = 0.63^{+0.23}_{-0.49}$ with that from the count rate $\eta$ at a level of $5 \approx 1\sigma$. This is in line with expectations, because the ICM mass and count rate are both derived mainly based on an integral of the X-ray surface brightness profile, as a pair of highly correlated quantities. Meanwhile, there is no significant correlated scatter between the ICM mass and the weak-lensing mass ($\rho_{WL,M_g} = 0.08^{+0.60}_{-0.43}$) and between the count

<footnote n="8">⁸ We make use of Eq. (4) in Evrard et al. (2014) to convert the intrinsic scatter of the mass at fixed temperature to that of the temperature at fixed mass.</footnote>
<figure><img src="image_20.png" /><figcaption>**Fig. 19.** Parameter constraints in the joint modeling of the $M_g–M_{500}–z$ relation and the weak-lensing mass calibration, presented in the same way as in Fig. 12.</figcaption></figure>

rate and weak-lensing mass ($\rho_{WL,\eta} = -0.30^{+0.32}_{-0.38}$), as consistent with other follow-up X-ray scaling relations.

We make an additional remark. Although the ICM mass has low intrinsic scatter at fixed mass and can be directly calculated based on the best-fit model of X-ray imaging, it is a less ideal mass proxy than the luminosity. This is because prior knowledge of the underlying cluster mass is needed to properly calculate the integrated ICM mass enclosed by $R_{500}$, given that $M_{500}(< r) \sim r$ for a large radius $r$. On the other hand, the luminosity is dominated by the cluster core and hence does not strongly rely on a prior on the cluster radius, despite the large scatter.

## Comparisons with the literature

We first compare the mass trend of $M_g$ between the eFEDS sample and the literature: The mass trend of the $M_g–M_{500}–z$ relation

was constrained as $1.25 \pm 0.06$ by Arnaud et al. (2007) based on a sample of ten nearby clusters. In Vikhlinin et al. (2009), the scaling of the ICM mass fraction is constrained as $M_g / M_{500} \propto 0.037 \pm 0.006 \ln M_{500}$, corresponding to a mass trend<fnref n="9" /> as $M_g \propto M_{500}^{1.15 \pm 0.02}$. Based on a representative sample of local clusters selected by the X-ray luminosity, Pratt et al. (2009) obtained a mass-trend slope of $1.21 \pm 0.03$ in a mass range of $10^{14}h^{-1}M_{\odot} \lesssim M_{500} \lesssim 10^{15}h^{-1}M_{\odot}$. Zhang et al. (2012) derived a mass scaling of $M_g \propto M_{500}^{1.38 \pm 0.36}$ using a sample of 19 nearby clusters with masses of $3 \times 10^{14}h^{-1}M_{\odot} \lesssim M_{500} \lesssim 10^{15}h^{-1}M_{\odot}$. Mahdavi et al. (2013) constrained the mass trend with a self-similar slope of $1.04 \pm 0.1$ for massive clusters with a temperature above 3 keV. A slope of $1.22 \pm 0.04$ was derived in Lovisari et al. (2015) across

<footnote n="9">⁹ This slope is obtained by re-fitting the data from Vikhlinin et al. (2009) in the functional form of $M_g / M_{500} = B(M_{500})$. See also Sect. 5.2 in Chiu et al. (2018a).</footnote>
<figure><img src="image_21.png" /><figcaption>**Fig. 20.** Observed ICM mass $M_g$ as a function of the cluster mass $M_{500}$ (left panel) and redshift $z$ (right panel). This plot is produced and shown in the same way as in Fig. 14.</figcaption></figure>

the mass range of $10^{13}h^{-1}M_{\odot} \lesssim M_{500} \lesssim 10^{15}h^{-1}M_{\odot}$. A mass-trend slope of $1.007 \pm 0.012$, which is highly consistent with self-similarity, was also suggested by Mantz et al. (2016) using the ROSAT-selected clusters. The 100 brightest galaxy clusters in the XXL survey with $M_{500} \ge 2 \times 10^{13}h^{-1}M_{\odot}$ resulted in a slope of $1.21^{+0.11}_{-0.10}$ (Eckert et al. 2016). The ICM mass of SPT-selected clusters at a mass range of $M_{500} \ge 3 \times 10^{14}h^{-1}M_{\odot}$ was constrained as $M_g \propto M_{500}^{1.33 \pm 0.09}$ and $M_g \propto M_{500}^{1.26^{+0.12}_{-0.09}}$ in Chiu et al. (2018a) and Bulbul et al. (2019), respectively. In summary, the mass trend of $M_g$ obtained from the eFEDS sample shows broad consistency with the previous work and no clear tendency to deviate from the self-similar behavior ($B=1$). We plot these comparisons in the lower-left panel of Fig. 16, where the simulated clusters (golden stars) from Barnes et al. (2017) nicely show an extrapolation from the eFEDS sample toward the high-mass end. As there is no significant redshift-dependent mass scaling ($\delta_{M_{ICM}} = 0.40^{+0.78}_{-0.70}$) of the eFEDS sample, our result suggests that this picture still holds out to redshift $z \approx 1.3$.

With the capability of eROSITA in discovering high-redshift clusters, we are able to constrain the redshift trend of the $M_g$-$M_{500}-z$ relation and make comparisons with previous work. The redshift trend of $M_g$ in the eFEDS sample is constrained to have a slope of $\gamma_{M_{ICM}} = 0.32^{+0.59}_{-0.61}$, showing no clear deviation from self-similarity ($\gamma=0$) out to redshift $z \approx 1.3$ (see Fig. 20). Based on the sample from the SPT, the redshift scaling of $M_g$ with respect to the pivotal mass was constrained as $M_g \propto (1+z)^{-0.15 \pm 0.22}$ and $M_g \propto (1+z)^{0.18^{+0.30}_{-0.31}}$ in Chiu et al. (2018a) and Bulbul et al. (2019), respectively. A steep mass trend as $M_g \propto E(z)^{1.76 \pm 1.22}$ with a large error bar was derived in Sereno et al. (2020) based on the XXL sample. Briefly, the ICM mass of the eFEDS sample shows a consistent picture with previous work, namely no significant deviation from self-similarity in the redshift trend is suggested.

For the eFEDS sample, the intrinsic scatter of the ICM mass at fixed cluster mass is constrained as $\sigma_{M_{ICM}} = 0.074^{+0.063}_{-0.019}$, where the asymmetric error bars are attributed to the lower bound of the adopted prior. This result is consistent with the SPT sample of massive clusters ($M_{500} \ge 3 \times 10^{14}h^{-1}M_{\odot}$) out to high redshift $z \approx 1.3$ with the scatter of $0.11 \pm 0.02$ (Chiu et al. 2018a) and $0.10^{+0.05}_{-0.07}$ (Bulbul et al. 2019) using the Chandra and XMM-Newton telescopes, respectively. Pratt et al. (2009) derived the scatter of the ICM mass fraction at fixed cluster mass as

0.12 ± 0.03. The simulation results from Henden et al. (2019) suggest the scatter of $0.12^{+0.03}_{-0.02}$ at redshift $z=0$, decreasing to ≈0.05 at redshift $z \ge 1$. A value of 0.13 was obtained for the intrinsic scatter of $M_g$ in Arnaud et al. (2007) using the conversion formula in Evrard et al. (2014). Similarly, the intrinsic scatter of $M_g$ at a level of $0.14 \pm 0.06$ was obtained in Mahdavi et al. (2013). The simulations in Henden et al. (2019) suggest a decreasing intrinsic scatter at the level of ≈0.27 at redshift $z \approx 0$ toward a value of ≈0.12 at redshift $z \ge 1$. Based on these comparisons, we conclude that the eFEDS sample has intrinsic scatter of $M_g$ that is broadly consistent with that found in previous studies, suggesting that the ICM mass can be used as a mass proxy that shows little scatter.

## 6.5. The $Y_X-M_{500}-z$ relation

The mass proxy, $Y_X$, is defined as a product of the temperature $T_X$ and the ICM mass $M_g$, namely, $Y_X \equiv T_X M_g$. The mass proxy $Y_X$ was first proposed in Kravtsov et al. (2006) and is demonstrated to be insensitive to baryonic feedback (see also Nagai et al. 2007). It is found to show remarkably low scatter and is widely used to infer the mass of clusters via X-rays (Vikhlinin et al. 2009; Arnaud et al. 2010). In addition, this quantity can be directly connected to the SZE observable, providing a window to study galaxy clusters in synergy with X-ray and SZE (Bonamente et al. 2008; Andersson et al. 2011; Benson et al. 2013). The $Y_X-M_{500}-z$ relation has a self-similar scaling as

$$Y_X \propto M_{500}^{5/3} E(z)^{2/3},$$

which is directly derived from the self-similar behavior of $T_X$ and $M_g$ by definition.

In this work, the $Y_X-M_{500}-z$ relation is derived as

$$ \begin{aligned} \left\langle \ln \left( \frac{Y_X}{M_{\odot} \times \text{keV}} \middle| M_{500} \right) \right\rangle &= \ln(6.9^{+1.5}_{-1.2}) + \ln(10^{13}) \\ &+ \left[ (1.78^{+0.24}_{-0.20}) + (-1.1^{+1.9}_{-1.2}) \ln \left( \frac{1+z}{1+z_{\text{piv}}} \right) \right] \times \ln \left( \frac{M_{500}}{M_{\text{piv}}} \right) \\ &+ \frac{2}{3} \times \ln \left( \frac{E(z)}{E(z_{\text{piv}})} \right) + (-1.43^{+1.20}_{-0.93}) \times \ln \left( \frac{1+z}{1+z_{\text{piv}}} \right), \end{aligned} \quad (71) $$
<figure><img src="image_22.png" /><figcaption>**Fig. 21.** Parameter constraints in the joint modeling of the $Y_X-M_{500}-z$ relation and the weak-lensing mass calibration, presented in the same way as in Fig. 12.</figcaption></figure>

with log-normal scatter of $\sigma_{Y_X} = 0.106_{-0.047}^{+0.171}$. The mass proxy $Y_X$ of the eFEDS sample is constrained as $6.9_{-1.2}^{+1.5} \times 10^{13} M_{\odot}$ keV at the pivotal mass $M_{\text{piv}} = 2.5 \times 10^{14} h^{-1} M_{\odot}$ and the pivotal redshift $z_{\text{piv}} = 0.35$. We show the parameter constraints in Fig. 21 with the results listed in Table 2. The mass and redshift trends of $Y_X$ for the eFEDS sample are shown in Fig. 22, following the normalization scheme in producing Fig. 14.

Based on our results, we find a self-similar mass trend ($B_{Y_X} = 1.78_{-0.20}^{+0.24}$) and no significant deviation ($\gamma_{Y_X} = -1.43_{-0.93}^{+1.20}$) from self-similarity ($Y_X \propto E(z)^{2/3}$) in the redshift scaling. Both the mass and redshift scaling are better visualized in Fig. 22. In addition, no significant redshift-dependent mass trend is seen ($\delta_{Y_X} = -1.1_{-1.2}^{+1.9}$), suggesting that the mass scaling holds out

to high redshift $z \approx 1.3$. The self-similar slopes in both the mass and redshift scaling are expected, as a combined result of the $T_X-M_{500}-z$ and $M_g-M_{500}-z$ relations. We note that the $Y_X$ measurements are more noisy and show greater scatter above the best-fit model (the gray area in Fig. 22), which is clearly due to the higher uncertainty at the high-temperature end (see Sect. 6.3).

The intrinsic scatter of $Y_X$ is constrained as $\sigma_{Y_X} = 0.106_{-0.047}^{+0.171}$ with a correlation coefficient of $\rho_{WL,Y_X} = 0.12_{-0.50}^{+0.53}$ which is consistent with zero. Interestingly, the correlation between the intrinsic scatter of $Y_X$ and $\eta$ is constrained to be positive ($\rho_{Y_X,\eta} = 0.794_{-0.344}^{+0.079}$) at a significance level that is higher than those of $\rho_{T_X,\eta}$ and $\rho_{M_g,\eta}$. This is expected, because both $\rho_{T_X,\eta}$ and $\rho_{M_g,\eta}$ are also constrained to be mildly positive.
<figure><img src="image_23.png" /><figcaption>**Fig. 22.** Observed $Y_X$ as a function of the cluster mass $M_{500}$ (left panel) and redshift $z$ (right panel). This plot is produced and shown in the same way as in Fig. 14.</figcaption></figure>

## Comparisons with the literature

Similarly to the X-ray scaling relations discussed in Sects. 6.2 to 6.4, we first compare the mass trend ($B_{Y_X} = 1.78^{+0.24}_{-0.20}$) of the $Y_X-M_{500}-z$ relation of the eFEDS sample with the results of other studies: In a pilot study, Arnaud et al. (2007) derived a slope of $1.82 \pm 0.09$ based on a sample of ten nearby clusters while fixing the redshift scaling to the self-similarity ($Y_X \propto E(z)^{2/3}$). Vikhlinin et al. (2009) constrained the slope of the mass trend as $1.75^{+0.10}_{-0.09}$ for ROSAT-selected clusters out to $z \approx 0.5$, assuming $Y_X \propto E(z)^{2/3}$. A mass-trend slope of $1.79^{+0.26}_{-0.20}$ using a sample of 50 massive clusters was derived in Mahdavi et al. (2013) with the fixed redshift scaling of $E(z)^{2/3}$. Later, Lovisari et al. (2015) derived a slope of $1.75 \pm 0.03$ using a joint sample of galaxy groups and massive clusters with temperatures above 3 keV. The ROSAT-selected sample was used in Mantz et al. (2016) to constrain the $Y_X-M_{500}-z$ relation, resulting in a slope of $1.63 \pm 0.04$. Using a sample of 59 massive ($M_{500} \ge 3 \times 10^{14} h^{-1} M_\odot$) clusters selected in the SPT, the mass scaling of $Y_X$ was constrained to be $2.02^{+0.16}_{-0.17}$ out to redshift $z \approx 1.3$. We compare these results with that of the eFEDS sample in the lower-right panel of Fig. 16, where we additionally show the simulated clusters (golden stars) from Barnes et al. (2017). As seen, the eFEDS sample has the mass scaling of $Y_X$ that is consistent with previous work.

As an interesting comparison, our result of $B_{Y_X} = 1.78^{+0.24}_{-0.20}$ is also in good agreement with the mass trend of the integrated pressure $Y_{SZ}$ inferred from the SZE, which was derived with a slope of $1.67 \pm 0.29$ by Andersson et al. (2011) using 15 SPT clusters with the redshift scaling fixed to $E(z)^{2/3}$.

In terms of the redshift trend, the eFEDS sample shows a scaling as $Y_X \propto E(z)^{2/3}(1+z)^{-1.43^{+1.20}_{-0.93}}$. Our result is consistent with that of Bulbul et al. (2019), in which $Y_X$ was derived with the scaling as $Y_X \propto E(z)^{2/3}(1+z)^{-0.17^{+0.47}_{-0.50}}$, showing no significant deviation from self-similarity. A cosmological analysis based on the abundance of SPT-selected clusters resulted in a slope of $0.310^{+0.209}_{-0.140}$ (Bocquet et al. 2019), which is statistically consistent with our result. Based on these comparisons, the mass proxy $Y_X$ of the eFEDS sample shows a self-similar redshift trend, as consistent with previous work.

The log-normal intrinsic scatter of $Y_X$ at fixed cluster mass is constrained as $\sigma_{Y_X} = 0.106^{+0.171}_{-0.047}$, which is in good agreement with the previous literature: The ROSAT-sample from Mantz et al. (2016) resulted in the scatter of $0.185 \pm 0.016$. A joint analysis of cluster abundance and the integrated pressure from both X-rays and SZE resulted in the intrinsic scatter of $0.184^{+0.087}_{-0.089}$ for the SPT-selected clusters (Bocquet et al. 2019). The sample of simulated clusters from Barnes et al. (2017) suggested the scatter of 0.23 for the $Y_X-M_{500}-z$ relation. A similar value of $\approx 0.23$ was also suggested by the simulations from Henden et al. (2019) at redshift $z \approx 0$, decreasing to $\approx 0.15$ for clusters at redshift $z \gtrsim 1.0$. The intrinsic scatter in the core-included and core-excised $Y_X$ of the SPT clusters at fixed mass was constrained as $0.16^{+0.07}_{-0.10}$ and $0.13^{+0.07}_{-0.08}$, respectively (Bulbul et al. 2019). These results, together with that from the eFEDS, imply that $Y_X$ serves as a low scattering proxy.

## 7. Systematic errors

The key and novel ingredients in this work are the simulation-based calibrations of the observed count rate $\eta$ and the weak-lensing mass $M_{WL}$, as described in Sect. 4. The calibrations are obtained by performing the end-to-end analysis on the simulations with properties that are as realistic as those observed, resulting in an empirical relation between the observed quantity and the underlying true mass. As a result, a systematic change in the simulation calibrations results in a systematic uncertainty of the scaling relation (and hence the cluster mass). Therefore, in this work we effectively marginalize over the systematic uncertainty by applying the informative priors, which are obtained from the simulation-based calibrations, on the parameters of Eqs. (40) and (42). That is, the systematic errors of the simulation calibration are marginalized over in our results.

One of the remaining systematic errors is from whether or not the cluster core is included in the modeling of weak-lensing shear profiles. As seen in Table 2, including the cluster core results in a negligible difference in the parameters of the $\eta-M_{500}-z$ relation compared to the current uncertainty. This, again, is supported by our end-to-end calibration of the weak-lensing mass bias presented in Sect. 4.2. Quantitatively, the cluster mass
would decrease by ≈0.01 dex if the cluster core were include ;
this is a subdominant systematic uncertainty in this work.

A source of potential systematic error could be the photo-
z bias in the 30-band COSMOS sample, if indeed there is one.
In this work, the photo-z from the 30-band COSMOS catalog
(Laigle et al. 2016) is used as the true redshift. Hence, a photo-
z bias in the COSMOS catalog would result in a biased lensing
signal; this is especially true for the faint-end regime (see also
Schrabback et al. 2010). This issue, the reliability of the photo-
z in the COSMOS catalog, was also thoroughly discussed in a
recent cosmic shear analysis (see Sect. 8 in Hildebrandt et al.
2020), in which the authors also estimated how the outlier frac-
tion of the photo-z may affect the mean redshift of the lensing
sources. Although the accuracy and the outlier fraction of the
COSMOS photo-z are estimated to be less than ≈1% for rel-
atively bright galaxies (i ≲ 21 mag; Laigle et al. 2016), the
outlier fraction is quantified to be at a level of ≈6% for galax-
ies with an i-band magnitude of between 23 and 24. However, it
is unclear and hard to quantify how the outlier fraction of the
COSMOS catalog affects the cluster mass estimate without a
representative spectroscopic sample out to high redshift. More-
over, our source selection effectively removes source galaxies
with low photo-z, which further affects the systematic impact
arising from catastrophic photo-z outliers. Here, we make an
assumption that the galaxies with available spectroscopic red-
shifts<fnref n="10" /> in the COSMOS field provide a representative sample of
those at the faint end and/or the high redshift, and try to estimate
the bias arising from the inaccuracy in the COSMOS photo-z
(see also Schrabback et al. 2018a,b; Raihan et al. 2020). Specif-
ically, we re-calculate Eq. (24) based on the COSMOS galaxy
sample with available spectroscopic redshifts, and replace the
observed and true redshifts with the COSMOS 30-band photo-z
and the spectroscopic redshift when evaluating Eq. (25), respec-
tively. The resulting photo-z bias $b_z$ ($z_{cl}$), which indicates the bias
arising from the inaccuracy of the COSMOS 30-band photo-z,
shows a behavior that is generally consistent with zero at $z \lesssim 0.6$
with a weak trend that slightly increases to a level of ≈1% at
$z \approx 0.8$ and then decreases toward a level of ≈ −3% at $z \approx 1.2$,
despite with large error bars. This result suggests that the inac-
curacy of the COSMOS photo-z may introduce systematic bias
to the cluster mass (primarily at high redshift, $z \ge 0.8$) at a level
of $\le 3\%$, assuming that the available spectroscopic sample in the
COSMOS is complete. For a study with a significantly larger
sample size from the eROSITA All-Sky Survey, a more detailed
investigation on the photo-z bias using an alternative method
(e.g., the clustering-based calibration; Gatti et al. 2022) is clearly
recommended.

Another source of systematic error comes from the fact that
we ignore the correlation between the uncertainty in the follow-
up X-ray quantities and that of the cluster mass. Because we
perform the early unblinding (see Sect. 3.2) to self-consistently
estimate the X-ray properties extracted at the “fixed” cluster
radius R<sub>500</sub>, a change in the cluster mass M<sub>500</sub> will lead to a
change in R<sub>500</sub> and the follow-up X-ray observables. We gauge
this correlated uncertainty as follows. In this work, the statisti-
cal uncertainty of each cluster mass, which is estimated as the
dispersion of the resulting posterior of M<sub>500</sub> (see Sect. 6.1), is

at the level of ≈32%, ≈27%, and ≈23% for a cluster with a
mass of ≈10¹³.⁶h⁻¹M☉, ≈10¹⁴h⁻¹M☉, and ≈10¹⁴.⁴h⁻¹M☉, corre-
sponding to a change in the radius at a level of ≈10%, ≈9%,
and ≈8%, respectively. The ICM mass M_g is more sensitive to
the cluster radius among the X-ray observables, because it is a
direct integration of the surface brightness profile up to R₅₀₀.
Assuming that the ICM density distribution n_e follows a β-model
(Cavaliere & Fusco-Femiano 1976) with a typical value of β =
2/3, the integrated profile of the ICM is ∝ ∫ n_e r² dr ∝ r. That is,
the change in R₅₀₀ leads to roughly the same amount of uncer-
tainty in M_g. Meanwhile, a typical measurement uncertainty of
the ICM mass is on the order of ≈50%, ≈42%, and ≈28% for an
eFEDS cluster with a mass of ≈10¹³.⁶h⁻¹M☉, ≈10¹⁴h⁻¹M☉, and
≈10¹⁴.⁴h⁻¹M☉, respectively. In comparison, the radius-induced
statistical uncertainty is much smaller than the current mea-
surement uncertainty and is thus negligible. In terms of the
systematic uncertainty in the cluster mass, which is expected to
be at the level of ≈7.6% as inferred from the weak-lensing mass
bias calibration (see Eq. (44) in Sect. 4.2), the corresponding
change in R₅₀₀ is at the level of ≈7.6%/3 ≈ 2.5%. On the other
hand, the uncertainties on the mean M_g at the mass bins cen-
tering at log((M₅₀₀/h⁻¹M☉) = 13.6, 14, and 14.4 with a bin width of
0.2 dex are approximately 50%/√30 ≈ 9%, 42%/√90 ≈ 4%,
and 28%/√64 ≈ 3.5%, respectively, where the numbers in the
square root are the numbers of the clusters in the bins. That is,
the change due to a systematic uncertainty in the cluster mass is
still subdominant compared to the current measurement uncer-
tainty. Therefore, we conclude that the uncertainty due to the
change in M₅₀₀ is expected to be subdominant for the modeling
of the X-ray scaling relations. We note that the luminosity, which
is also estimated based on the surface brightness profile, is dom-
inated by the cluster core<fnref n="11" /> and is thus less sensitive to the cluster
radius. An improvement overcoming these systematic errors is to
use an iterative way to estimate the X-ray observable from a pro-
file in each step of the likelihood exploration (e.g., Benson et al.
2013; Bocquet et al. 2019), but we defer this to a future work.

8. Conclusions

We present the first weak-lensing mass calibration of eROSITA-detected clusters selected in the eFEDS survey using the three-year weak-lensing shape catalog (S19A) from the HSC survey. We use a sample of 434 optically confirmed clusters at redshift 0.01 $\le z \le 1.3$ (median 0.35) with contamination from point sources at the level of $\approx 2\%$. The uniform coverage of the HSC survey enables us to extract the lensing shear profile for the majority ($\approx 70\%$, or 313 clusters) of the sample.

We first perform the weak-lensing mass calibration by simul-
taneously modeling the observed count rate η, as the X-ray
mass proxy, and the shear profile g_+ on the basis of individual
clusters through a Bayesian population modeling of the count-
rate-to-mass-and-redshift (η–M_500–z) and weak-lensing-mass-
to-mass-and-redshift (M_WL–M_500–z) relations, respectively. The
weak-lensing mass M_WL, which is obtained by fitting a mass
model to g+, is carefully calibrated against cosmological hydro-
dynamic simulations that are exclusively constructed according
to the observed lensing properties, including the multiplicative
bias of shape measurements as a function of source redshift,
the shape noise of g+, the uncertainty arising from uncorre-
lated large-scale structures, the cluster member contamination,

<footnote n="10">¹⁰ Specifically, we make use of the spectroscopic redshifts estimated from 3DHST (Skelton et al. 2014; Momcheva et al. 2016), C3R2 (Masters et al. 2019), DEIMOS (Hasinger et al. 2018), FMOS-COSMOS (Silverman et al. 2015), LEGA-C (Straatman et al. 2018), PRIMUS (Coil et al. 2011; Cool et al. 2013), SDSS (Alam et al. 2015; Pâris et al. 2018), and zCOSMOS (Lilly et al. 2009).</footnote>

<footnote n="11">¹¹ $L_X \propto \int n_e^2 \sqrt{\bar{T}r^2} dr \propto C^{-1/r}$, with C up to a constant. Here, we assume that the temperature is insensitive to the cluster radius.</footnote>