© ESO 2022

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, Bocquet et al. 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 σ₈, 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 cosmologi-cal 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 cos-mological 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 ≈ 1.8, based on which unbiased cosmological constraints have been obtained (Bocquet et al. 2019). An even larger sample of SZ-selected clusters has recently been released by the Atacama Cosmology Telescope survey (Hilton et al. 2021). However, current SZ surveys are limited to massive clusters because of the sensitivity of detectors, or only scan a small fraction of the sky, resulting in a sample of relatively small size. To extend the success of the SZ-based cosmology to the low-mass regime, it is therefore necessary to combine the SZ sample with those from other surveys (e.g., Costanzi et al. 2021). In terms of medium-size X-ray surveys with advanced telescopes, such as the XXL survey (Pierre et al. 2016), a large number of clusters down to a low-mass regime has be obtained (Adami et al. 2018), although mainly at low redshift because of the cosmological dimming of bremsstrahlung emissions from clusters. Meanwhile, the first all-sky survey in X-rays carried out by the ROSAT mission (Truemper 1982) had delivered the largest X-ray cluster catalog ever across the whole sky (Böhringer et al. 2001, 2004; Piffaretti et al. 2011; Boller et al. 2016), but the poor resolution and shallow depth have limited cosmological studies in practice. However, this situation is no longer true with the ongoing all-sky X-ray survey carried out by the state-of-the-art telescope, eROSITA.

The eROSITA (extended ROentgen Survey with an Imaging Telescope Array; Merloni et al. 2012; Predehl et al. 2021) is an X-ray space telescope on board the Russian-German “Spectrum-Roentgen-Gamma” (SRG) satellite (Sunyaev et al. 2021), which was successfully launched on July 13, 2019. The main goal of the eROSITA mission is to unveil the nature of dark energy by carrying out the eROSITA All-Sky Survey (eRASS), which is to map the whole sky in X-rays for four years. By doing so, the eRASS will deliver the largest sample of ICM-selected galaxy clusters to date and discover more than 100 000 clusters at the end of the four-year all-sky survey (Borm et al. 2014). With this revolutionary sample, cluster-based cosmological analyses carried out as part of the eROSITA mission will provide unprecedented power with which to constrain cosmology (Pillepich et al. 2012, Pillepich et al. 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 Ω_M, σ₈, 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, Schrabback et al. 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 ≳ 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 ≳ 1.1 based on the HSC data sets was performed for a large sample of (≳1700) optically selected clusters over a footprint with an area of ≈140 deg² (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 ≈ l40 deg² 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 ΛCDM cosmology with Ω_M = 0.3, a mean baryon density Ω_b = 0.05, the Hubble parameter H₀ = h × 100 km s⁻¹ Mpc⁻¹ with h = 0.7, σ₈ = 0.8, and a spectral index of the primordial power spectrum n_s = 0.95. While modeling the scaling relations, we vary the cosmologi-cal parameters around those of the fiducial cosmology, with the aim being to maintain the flexibility needed in a future cos-mological analysis. The cluster true mass M₅₀₀ is defined by a sphere with a radius of R₅₀₀, such that the enclosed mass density is 500 times the critical density ρ_c(z) of the Universe at the cluster redshift. All quoted errors represent the 68% confidence level (i.e., 1σ) throughout this work, unless otherwise stated. The notation $\mathcal{N}\left(x,y^2\right)\left(\mathcal{U}\left(x,y\right)\right)$ 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.

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., fcont < 0.2 (see Sect. 2.3). These clusters are represented by the circles, color-coded by their observed redshifts Zcl, with sizes proportional to their observed count rate η. The other clusters with fcont ≥ 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.

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 ≈l 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 red-shift (z ≳ 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.2m 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 $\sqrt{2-1}\approx 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 follows, we give a brief summary of the cluster finding algorithm in the eFEDS.

The source detection of galaxy clusters is performed on the X-ray imaging at the energy band of 0.2-2.3 keV. The source detection is performed with the eSSASusers_20l009 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 likelihood ${\mathcal{L}}_{\text{det}}$, the extent likelihood ${\mathcal{L}}_{\text{ext}}$, and the observed count rate η at 0.2–2.3 keV. We note that the observed count rate η 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 η is different from that estimated within the cluster radius R₅₀₀. We estimate the bias¹ of the observed count rate η 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 systems 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_cı, the optical richness λ, and the estimator f_cont of the probability that the cluster appears as a superposition of galaxies along the line of sight by chance. The quantity of f_C0nt indicates the level of contamination in the optical confirmation. For instance, a value of f_C0nt = 0.2 means that there is a 20% chance of detecting a candidate with the derived richness λ at the redshift z_c1 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_cont is equivalent to a reduction in the cluster contamination in the initial catalog, under the assumption that the contamination from X-rays and the optical is uncorrelated (see also the discussion at length in Grandis et al. 2020). The contamination from point sources in the initial catalog of eFEDS clusters is suggested to be at a level of ≈20% using simulations (Liu et al. 2022a). Applying the cut of f_cont < 0.2 results in a deduction of 18% (or 99 clusters) in the sample size, which is consistent with the expectation from the simulation. The residual contamination from nonex-tended sources is suggested to be at a level of ≈4% with the cut of f_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_C0nt < 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 η 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 (λ-M₅₀₀−Z) relation for the eFEDS clusters following the exact methodology described in Sect. 5.3. The results are shown in Appendix A.

Fig. 2 Observed count rate η (in the unit of counts per second) and the redshift z of the 542 eFEDS clusters. The 434 secure clusters with fcont < 0.2 (see Sect. 2.3) are shown as the circles, while the rest are marked by red triangles.

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₅₀₀) 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₅₀₀, the cluster mass needs to be known before the follow-up X-ray analysis. In this work, the cluster mass M₅₀₀ 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₅₀₀ (see more details in Sect. 3.2), which is then used as the aperture size to derive the X-ray observables. After deriving the follow-up X-ray observables, the eROSITA team delivers them to the HSC team for the modeling of the X-ray scaling relations. It is important to note that the modeling of the follow-up X-ray scaling relations is solely led by the HSC team, which has no exchange of any mass-related information with the eROSITA team after the early unblinding. Therefore, the analysis of the X-ray scaling relations in this work is still performed blindly to minimize confirmation bias.

With the exception of T_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₅₀₀ 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₅₀₀ 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₅₀₀) and T_X(R < 500 kpc). As a result, we find that the difference between T_X(R < R₅₀₀) and T_X(R < 500 kpc) for the subsample with ${\mathcal{L}}_{\text{det}}>50$ and ${\mathcal{L}}_{\text{ext}}>50$ has an inverse-variance-weighted mean of 〈T_X(R < R₅₀₀) − T_X(R < 500 kpc)〉 = (0.06 ± 0.36) 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₅₀₀ 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₅₀₀), instead of T_X(R < 500 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_cl depends on the distances between the cluster, the observer, and the source. This is characterized by the critical surface density Σ_crit defined as $${\displaystyle {\sum }_{\text{crit}}\left(z_{\text{cl}},z_s\right)}=\frac{c^2}{4\pi G}\frac{1}{D_1\left(z_{\text{cl}}\right)\beta \left(z_{\text{cl}},z_s\right)},$$(1)

in which G is Newton’s constant, c is the speed of light, D₁ (z_cl) is the angular diameter distance of the cluster, and β is the lensing efficiency of the source, $$\beta \left(z_{\text{cl}},z_s\right)\{\begin{array}{cc}\frac{D_{\text{l,s}}\left(z_{\text{cl}},z_s\right)}{D_{\text{s}}\left(z_s\right)}& {\text{ifz}}_s>z_{\text{cl}}\\ 0& {\text{ifz}}_s\le z_{\text{cl}}\end{array},$$(2)

where D_l,s (z_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-k},$$(3)

where γ₊ and κ are the tangential shear and convergence of the cluster, respectively. In an azimuthal average with respect to a cluster center, the tangential shear γ₊ (R) at the projected radius R describes the differential surface mass density ΔΣ_m of the cluster at R with respect to the critical surface density inferred from the source, namely, $${\gamma }_{+}\left(R\right)=\frac{\Delta {\displaystyle {\sum }_{\text{m}}\left(R\right)}}{{\displaystyle {\sum }_{\text{crit}}}}.$$(4)

Meanwhile, the convergence κ (R) is the ratio of the projected surface mass density Σ_m (R) to the critical surface density, as $$\kappa \left(R\right)=\frac{{\displaystyle {\sum }_{\text{m}}\left(R\right)}}{{\displaystyle {\sum }_{\text{crit}}}}.$$(5)

We note that ΔΣ_m (R) = Σ_m (< R) − Σ_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 Σ_crit to interpret γ₊ and κ (see Sect. 3.9 for the modeling of Σ_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}}+\text{d}{m}_{1,i}+\text{d}{m}_{2,i},$$(6)

where m_cat,i is the blinded multiplicative bias in the i-th shape catalog, m_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 ∈ {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₅₀₀ must be known to extract the X-ray spectra enclosed by the corresponding cluster radius R₅₀₀. 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 η, 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₅₀₀. 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_c1, 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.68<{\displaystyle {\int }_{z_{\min }}^{\infty }P\left(z\right)\text{d}z}\text{and}{\text{z}}_{\text{MC}}<2.5,$$(7)

where we set z_min = z_c1 + 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 ≲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 ≈15, ≈11, ≈6, ≈2, and ≈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 ≈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 dimen-sionless 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 ΔΣ_m. The latter require prior knowledge of redshift and the distance-to-redshift relation with assumed cosmological parameters to convert g₊ to ΔΣ_m. This results in a difficulty in the lensing modeling in a cosmologi-cal analysis, where the cosmological parameters are varied and hence change the observable ΔΣ_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_ci, we derive the tangential shear profile g₊ in angular bins of clustercentric radius θ_i, $$g_{+}\left({\theta }_i\right)=\frac{1}{2\Re \left({\theta }_{\text{i}}\right)\left(1+ℌ\left({\theta }_{\text{i}}\right)\right)}\frac{{\displaystyle \sum _{s\in {\theta }_i}{w}_s{e}_{+,s}}}{{\displaystyle \sum _{s\in {\theta }_i}{w}_s}},$$(8)

where the subscript s runs over the source sample that is selected according to the description in Sect. 3.3; w_s and e_+,s are the lens-ing weight and tangential ellipticity of the s-th source galaxy, respectively; ℜ(θ_i) is the shear response (see also Mandelbaum et al. 2005); the factor (1 + ℜ(θ_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 \left(2\varphi \right)e_1-\sin \left(2\varphi \right)e_2,$$(9)

where (e₁, e₂) is a two-component estimate of ellipticity measured in the cartesian coordinate of the sky defined in the HSC survey, and ϕ 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},$$(10)

where σ_e and e_rms are the measurement uncertainty and the root-mean-square estimate of the ellipticity per component, respectively.

The shear response ℜ(θ_i) in the angular bin θ_i is calculated as $$\Re \left({\theta }_i\right)=1-\frac{{\displaystyle \sum _{s\in {\theta }_i}{w}_s{e}_{\text{rms,}s}{}^2}}{{\displaystyle \sum _{s\in {\theta }_i}{w}_s}},$$(11)

where the subscript s runs over the source sample in the angular bin θ_i. The factor (1 + ℜ(θ_i)) is calculated as $$0+ℌ\left({\theta }_i\right)=\frac{{\displaystyle \sum _{s\in {\theta }_i}{w}_s\left(1+m_s\right)}}{{\displaystyle \sum _{s\in {\theta }_i}{w}_s}},$$(12)

where m_s is the multiplicative bias of the s-th source galaxy.

The angular binning in clustercentric radii is defined as follows. For each cluster, we first use a logarithmic binning between 0.2 h⁻¹Mpc and 3.5 h⁻¹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₊ (θ) 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_× (θ) according to Eq. (8), after replacing e₊ with e_×, which is defined as $$e_{\times }=\sin \left(2\varphi \right)e_1-\cos \left(2\varphi \right)e_2.$$(13)

The azimuthally averaged cross-shear profile e_× (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_× (θ) 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 ℂ comprises a component characterizing the noise of shape measurements, ℂ_shape, and another component ℂ_uLss accounting for the uncertainty arising from uncorrelated large-scale structures around the cluster. That is, $$ℂ={ℂ}_{\text{shape}}+{ℂ}_{\text{uLSS}},$$(14)

in which $${ℂ}_{{\text{shape}}_{i,j}}={\sigma }_{\text{shape}}^2\left({\theta }_i\right){\delta }_{i,j},$$

where δ_i,j is a Kronecker delta function, and $\Re \left({\theta }_i\right)$ in the angular bin θ_i is evaluated as $${\sigma }_{\text{shape}}^2\left({\theta }_i\right)=\frac{1}{4\Re \left({\theta }_i\right){\left(1+ℌ\left({\theta }_i\right)\right)}^2{\displaystyle \sum _{s\in {\theta }_i}{w}_s}}.$$(15)

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

We note that the term ℂ_uLSS depends on the underlying cosmology. Ideally, the variation in ℂ_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 ℂ_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 ℂ_uLSS so as to improve the speed of our calculation. We leave the improvement of a cosmology-dependent ℂ_uLSS to future work.

We also note that the intrinsic variation of shear profiles at fixed cluster mass² − 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 red-shift 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 red-shift 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 red-shift 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\left(z;R\right)=\left(1-f_{\text{cl}}\left(R\right)\right)P_{\text{bkg}}\left(z\right)+f_{\text{cl}}\left(R\right)P_{\text{cl}}\left(z;R\right),$$(16)

where f_cl (R) is the amount of cluster contamination at the radial bin R, P_bkg is the redshift distribution of galaxies at backgrounds, and P_cl is the redshift distribution of cluster members that lead to contamination. By construction, we have 0 ≤ f_cl (R) ≤ 1. We can write Eq. (16) in the integral of redshift, $${\displaystyle \int \text{d}zP\left(z;R\right)=\left(1-f_{\text{cl}}\left(R\right)\right){\displaystyle \int \text{d}z{P}_{\text{bkg}}}\left(z\right)+f_{\text{cl}}\left(R\right){\displaystyle \int \text{d}z{P}_{\text{cl}}\left(z;R\right)}},$$

and then integrate over the redshift interval ℑ, where the cluster component vanishes, and finally arrive at $$f_{\text{cl}}\left(R\right)=1\frac{{\displaystyle \underset{z\in \mathfrak{J}}{\int }\text{d}zP\left(z;R\right)}}{{\displaystyle \underset{z\in \mathfrak{J}}{\int }\text{d}z{P}_{\text{bkg}}\left(z\right)}}.$$(17)

In this way, the cluster contamination at the radial bin R can be derived for a given set of observed P (z; R) and P_bkg (z).

In practice, in the interest of more precise measurements, we estimate f_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 ≤ 0.35) and high-redshift (z > 0.35), and two count-rate bins, low-rate (η < 0.1) and high-rate (η > 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\left(z;R_i\right)=\frac{{\displaystyle \sum _{s\in j,R_i}{w}_s{P}_s}\left(z\right)}{{\displaystyle \sum _{s\in j,R_i}{w}_s}},$$(18)

where s runs over the source galaxies which are located at the radius bin R_i around the cluster j. In each subsample of clusters, we then derive the number-weighted average of redshift distributions observed at the radius bin R_i as $$\langle P\rangle \left(z;R_i\right)=\frac{{\displaystyle \sum _{j=1}^{N_{\text{d}}}{P}_j\left(z;R_i\right)N_{gal,j}\left(R_i\right)}}{{\displaystyle \sum _{j=1}^{N_{\text{d}}}{N}_{\text{gal,j}}\left(R_i\right)}},$$(19)

where j runs over the subsample of N_cl clusters, and N_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 ℑ as $\left\{z|z<z_{\text{lolim}}\right\}\cup \left\{z|z>z_{\text{uplim}}\right\}$, where $z_{\text{lolim}}=\text{min}{\left\{z_j\right\}}_{j=1}^{N_{\text{cl}}}-\text{Δ}z,z_{\text{uplim}}=\text{max}{\left\{z_j\right\}}_{j=1}^{N_{\text{cl}}}+\text{Δ}z$, and Z_j is the redshift of the j-th cluster. We set Δ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 Δz = 0.2; we verified that the results with Δz = 0.4 are statistically consistent with those based on Δz = 0.2, although with much larger error bars. In this way, we can estimate the numerator in Eq. (17) by integrating 〈P〉 (z;R_i) over the interval ℑ.

Concomitently, we estimate P_bkg (z) by repeating Eqs. (18) and (19) with the same scheme of radial binning and cluster red-shift on a number of N_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_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-η 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 η, 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 ≲3% (lσ), 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).

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 (zlolim and zuplim) 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 ≤3% (1σ), which is used to calibrate the weak-lensing mass bias (see Sect. 4.2). In these three panels, the numbers in the parentheses indicate the total numbers of the clusters used in the binning.

3.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\left(x\right)=\left(1-f_{\text{mis}}\right)\times P_{\text{cen}}\left(x|{\sigma }_{\text{cen}},\alpha \right)+f_{\min }\times P_{\text{mis}}\left(x|{\sigma }_{\text{mis}}\right).$$(20)

In Eq. (20), we assume that the distribution P(x) can be decomposed into two components: one is the centering component characterized by a modified Rayleigh distribution P_cen(x|σ_cen, α),

which reads $$P_{\text{cen}}\left(x|{\sigma }_{\text{cen}},\alpha \right)=\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)},$$(21)

where Γ is the gamma function; and the other is the miscentering component as modeled by a standard Rayleigh distribution, $$P_{\text{mis}}\left(x|{\sigma }_{\text{mis}}\right)=\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).$$(22)

These two components are weighted by the miscentering fraction f_mis in Eq. (20). We note that the parameter α is required to be positive and characterizes the shape of the centering component: The smaller α is the smaller radius where the peak of P_cen(x) occurs, and the steeper fall-off toward x = 0. Moreover, P_cen(x) reduces to a half-normal distribution and zero, as α approaches zero and +∞, respectively. Also, P_cen is identical to P_mis if α = 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 = R_off/r_s, where R_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_piv = 1.4 × 10¹⁴h⁻¹ M_⊙ at the pivotal redshift z_piv = 0.35. For each cluster, the offset R_off is calculated first in the physical unit using the redshift−distance relation inferred in the fiducial cosmology. We then convert R_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₅₀₀ 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

$$\begin{array}{c}\left(F_{\min },{\sigma }_{\text{con}},\alpha ,{\sigma }_{\text{mis}}\right)\\ =\left(0.54\pm 0.02,0.17\pm 0.01,0.26\pm 0.05,0.61\pm 0.03\right).\end{array}$$(23)

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

Second, we determine the distribution P(x) as the offset between the X-ray center and the BCG, while the required knowledge for modeling the miscentering is the offset between the X-ray and “true” centers of a cluster potential. That is, the distribution P(x) contains the intrinsic miscentering of both the X-ray center and the BCG with respect to the true center, as well as their respective measurement uncertainties at both wavelengths. Among these, the dominant component is the intrinsic miscentering of the BCG with respect to the true center, because (1) the X-ray center is considered as a much better tracer of the total potential than other proxies (Lin & Mohr 2004; Mahdavi et al. 2013; Lauer et al. 2014; Zhang et al. 2019) and (2) the uncertainty on the measurement of the optical center is much smaller than the characteristic miscentering radius. Therefore, the miscentering effect estimated in this work merely serves as an upper limit on X-ray miscentering. It was suggested that X-ray centers could bias the mass estimate even when the modeling of the cluster core (R < 500 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⁻¹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).

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 mis-centering components, as indicated by the red and green dashed lines, respectively.

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 Σ_crit in converting g₊ to ΔΣ_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_cı) 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 ß, as $$\frac{{\beta }_{\text{tiue}}}{{\beta }_{\text{obs}}}\left(z_{\text{cl}}\right)\equiv 1+b_z\left(z_{\text{cl}}\right),$$(24)

where ß_true and β_obs are the lensing efficiency inferred by the true and observed redshifts, respectively. In an ideal case, β_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³. 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 $${\beta }_{\text{true/obs}}\left(z_{\text{cl}}\right){\displaystyle \int \beta \left(z_{\text{cl}},z\right)}{P}_{\text{true/obs}}\left(z\right)\text{d}z,$$(25)

where P_trae(z) is the distribution of true redshift z_COSMOS30 in the source sample S, and $$P_{\text{obs}}\left(z\right)=\frac{{\displaystyle \sum _{s\in \text{S}}{w}_{{\text{SOM,}}_s{w}_s{P}_{\text{HSC,s}}\left(z\right)}}}{{\displaystyle \sum _{s\in \text{S}}{w}_{{\text{SOM,}}_s{w}_s}}}.$$(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 }\left(z_{\text{cl}}\right)=\frac{\langle {\beta }_{\text{S16A}}\left(z_{\text{cl}}\right)\rangle }{\langle {\beta }_{\text{S19A}}\left(z_{\text{cl}}\right)\rangle },$$

where the β_sl6A and β_s19A are the lensing efficiency calculated using the photo-z from the S16A and S19A data sets, respectively, and the bracket 〈·〉 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,SI9A used in this work is that from the S16A data set multiplying a re-scaling factor, namely $$1+b_{z,\text{S19A}}\left(z_{\text{cl}}\right)=R_{\beta }\left(z_{\text{cl}}\right)\times \left(1+b_{z,\text{S16A}}\left(z_{\text{cl}}\right)\right).$$(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 ≲2% for clusters at redshift z_cl ≲ 0.7, and becomes ≲6% for Z_cl ≳ 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 ≈2% and ≈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,_sl9a(z_cl) is included in deriving the weak-lensing mass bias for eFEDS clusters in Sect. 4.2.

Fig. 5 Photo-z bias bz (zcl) as a function of cluster redshift zcl. 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 bz(zcl) is applied to the calibration of weak-lensing mass bias for eFEDS clusters (see 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_cı, there are three input observables used in the modeling:

• the observed shear profile g+ (θ) as a function of angular radius θ, estimated in Eq. (8),

• the lensing covariance matrix ℂ, derived as Eq. (14), and

• the observed redshift distribution P (z), estimated in Eq. (34). Given a weak-lensing mass M_WL and a set of parameters P_miS = {f_mis,σ_mis} describing the miscentering distribution, together with a set of cosmological parameters 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)

$$\begin{array}{c}{g}_{+}^{\mathrm{mod}}\left(\theta |M_{\text{WL,Pmis,Pc}}\right)=\frac{{\gamma }_{\mathrm{mod}}\left(\theta |M_{\text{WL},\text{Pmis,Pc}}\right)}{1-{\kappa }_{\mathrm{mod}}\left(\theta |M_{\text{WL,Pmis},\text{Pc}}\right)}\\ \times {\left(1+{\kappa }_{\mathrm{mod}}\left(\theta |M_{\text{WL,Pmis,Pc}}\right)\left(\frac{\langle {\beta }^2\left(\text{Pc}\right)\rangle }{{\langle {\beta }^2\left(\text{Pc}\right)\rangle }^2}-1\right)\right)}^{,}\end{array}$$(28)

in which $$\langle \beta \left(\text{Pc}\right)\rangle ={\displaystyle \int \beta \left(z_{\text{cl,z}|\text{Pc}}\right)}P\left(z\right)dz,$$(29)

$$\langle \beta \left(\text{Pc}\right)\rangle ={\displaystyle \int {\beta }^2\left(z_{\text{cl,z}|\text{Pc}}\right)}P\left(z\right)dz,$$(30)

and the profiles κ_mod and γ_mod are calculated as $${\kappa }_{\mathrm{mod}}\left(\theta |M_{\text{WL,Pmin,Pc}}\right)=\frac{{\displaystyle {\sum }_{\text{m}}^{\text{mad}}\left(\theta |M_{\text{WL,Pmin,Pc}}\right)}}{{\displaystyle {\sum }_{\text{crit}}\left(\text{Pc}\right)}},$$(31)

and $${\gamma }_{\mathrm{mod}}\left(\theta |M_{\text{WL}},{\text{p}}_{\text{mis}},{\text{p}}_{\text{c}}\right)\frac{{\displaystyle {\sum }_{\text{m}}^{\text{mad}}\left(\theta |M_{\text{WL}},{\text{p}}_{\text{mis}},{\text{p}}_{\text{c}}\right)}}{{\displaystyle {\sum }_{\text{creit}}\left({\text{p}}_{\text{c}}\right)}},$$(32)

where $${\displaystyle {\sum }_{\text{crit}}\left({\text{p}}_{\text{c}}\right)}=\frac{c^2}{4\pi G}\frac{1}{D_1\left(z_{\text{cl}}|{\text{p}}_{\text{c}}\right)\langle \beta \rangle }.$$(33)

Given that our source selection leads to a highly pure sample with cluster contamination that is consistent with zero without a clear radial trend (see Fig. 3), we use a radially independent red-shift 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\left(z\right)=\frac{{\displaystyle \sum _{s\in j}{w}_s{P}_s\left(z\right)}}{{\displaystyle \sum _{s\in j}{w}_s}},$$(34)

where the index s runs over the sources with R < 3.5 h⁻¹Mpc in the fiducial cosmology.

Given the set of cosmological parameters p_c, the calculation of ${\text{Σ}}_{\text{m}}^{\text{mod}}$ and ${\text{ΔΣ}}_{\text{m}}^{\text{mod}}$ is performed in the radial bins in physical units by converting θ to R, as $$R=D_1\left(z_{\text{cl}}|{\text{p}}_{\text{c}}\right)\times \theta .$$(35)

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

The miscentering is accounted for in ${\text{Σ}}_{\text{m}}^{\text{mod}}$ (thus in ${\text{ΔΣ}}_{\text{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_mis. That is, $$\Delta {\displaystyle {\sum }_{\text{m}}^{\mathrm{mod}}=\left(1-f_{mis}\right)\Delta {\displaystyle {\sum }_{\text{m}}^{\text{cen}}+}{f}_{\min }}\Delta {\displaystyle {\sum }_{\text{m}}^{\text{mis}}.}$$(36)

The profile of ${\text{Σ}}_{\text{m}}^{\text{cen}}$ represents the surface mass density of a perfectly centered halo, which is evaluated using a NFW model. Meanwhile, the profile of ${\text{Σ}}_{\text{m}}^{\text{mis}}$ is modeled as an average of a set of NFW models with miscentering, weighted by a miscentering distribution P (R_mis), as $${\displaystyle {\sum }_{\text{m}}^{\text{mis}}\left(R\right)}={\displaystyle \int \text{d}{R}_{\text{mis}}P\left(R_{\text{mis}}\right)}{\displaystyle {\sum }_{\text{m}}^{\text{cen}}\left(R|R_{\text{mis}}\right)},$$(37)

where ${\text{Σ}}_{\text{m}}^{\text{cen}}\left(R|R_{\text{mis}}\right)$ is the surface mass density of a halo with a central offset of R_mis azimuthally averaged over the positional angle ϕ, expressed as (Yang et al. 2006; Johnston et al. 2007a) $${\displaystyle {\sum }_{\text{m}}^{\text{cen}}\left(R|R_{\min }\right)={\displaystyle \underset{0}{\overset{2\pi }{\int }}\frac{\text{d}\varphi }{2\pi }{\displaystyle {\sum }_{\text{m}}^{\text{cen}}\left(\sqrt{R^2+R_{\text{mis}}^2+2R{R}_{\mathrm{mi}\text{s}}\cos \varphi }\right).}}}$$(38)

When estimating ${\text{Σ}}_{\text{m}}^{\text{cen}}$ and ${\text{Σ}}_{\text{m}}^{\text{mis}}$ in Eqs. (36), we fix the concentration parameter, given a weak-lensing mass M_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_WL in order to convert the offset radius R_mis in Eq. (37) to the dimensionless radius x, i.e., x = R_mis/r_s. In this way, we evaluate P (R_mis) at R_mis following the second component of the derived offset distribution P (x), as in Eq. (20). That is, we ignore the centering component, P_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 σ_cen ≈ 0.2 in Eq. (21) is not significantly different from that of a perfectly centered model, ${\text{Σ}}_{\text{m}}^{\text{cen}}$, at the radial range of interest in this work. By doing so, we achieve a faster calculation of ${\text{ΔΣ}}_{\text{m}}^{\text{mod}}$ in Eq. (36).

We note that we statistically account for the effect of miscen-tering in the modeling of observed shear profiles by assuming a universal distribution P (R_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 miscen-tering 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₊ (θ) given a weak-lensing mass M_WL and a set of parameters, p = p_mis ∪ p_c, while accounting for the covariance among radial bins, as $$\begin{array}{c}\ln P\left(g_{+}|M_{\text{WL,}}{z}_{\text{cl}}\text{,p}\right)=\frac{-1}{2}\\ \times {\left(g_{+}^{\mathrm{mod}}\left(\theta |M_{\text{WL,}}{z}_{\text{cl}},\text{p}\right)-g_{+}\left(\theta \right)\right)}^{\text{T}}·{ℂ}^{-1}·\\ \left(g_{+}^{\mathrm{mod}}\left(\theta |M_{\text{WL}},z_{\text{cl}},\text{p}\right)-g_{+}\left(\theta \right)\right)+\text{constant}\text{.}\end{array}$$(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₅₀₀, which is referred to as the true count rate $\widehat{\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 $\widehat{\eta }$ enclosed by R₅₀₀ given the cluster true mass M₅₀₀. We note that the K-correction is included in calculating $\widehat{\eta }$. Meanwhile, we run the identical eSASS pipeline on the locations of simulated clusters to measure their observed count rate η. In this way, we can derive the bias, $$b_{\eta }\equiv \frac{\eta }{\widehat{\eta }},$$

as a function of cluster true mass and redshift. We only use clusters with true counts larger than 40 when calculating b_η, 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 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 }\left(M_{500},z\right)=\exp \left(A_{\text{f}}\right)\times {\left(\frac{M_{500}}{M_{\text{piv}}}\right)}^{B_{\text{f}}+{\delta }_{\text{f}}\ln {\left(\frac{z}{z_{\text{piv}}}\right)}^{{\gamma }_{{}_{\text{f}}}}}{\left(\frac{z}{z_{\text{piv}}}\right)}^{{\gamma }_{{}_{\text{f}}}},$$(40)

with four parameters: {A_f, B_f, δ_f, γ_f}. We use M_piv = 1.4 × 10¹⁴h⁻¹M_⊙ and z_piv = 0.35 in Eq. (40). The best-fit model describing the bias is shown in the left panel of Fig. 6. This model provides a good description of the simulated data, as suggested by their difference in the right panel of Fig. 6. This difference, which is the modeling residual of the count-rate bias, can be well approximated by a Gaussian distribution with a standard deviation of 0.08. This amount of modeling inaccuracy is absorbed into the Gaussian prior applied to the intrinsic scatter of observed count rates (see Sect. 5.4). We repeat the end-to-end procedure above for three different sets of simulations, which results in three sets of best-fit parameters. These three sets of simulations differ from one another in terms of the underlying N-body dark-matter-only simulations, the area coverage, and the abundance of point-source populations (at a level of 25%). Exhaustive runs of various simulations spanning a complete list of configurations are too time-consuming and are beyond the scope of this work. The goal of employing these different configurations is to gauge the inaccuracy of the derived count-rate bias as a systematic uncertainty that will be marginalized over in the analysis. Therefore, we take the mean and standard error of the three sets of the best-fit parameters ${\left\{A_{\text{f},i},B_{\text{f},i},{\delta }_{\text{f},i},{\gamma }_{\text{f},i}\right\}}_{i=1,2,3}$ separately inferred from the three simulations as the priors that will be marginalized over. The mean and standard error of these four parameters among the three simulations are $$\begin{array}{c}{A}_{\text{f}}=0.18\pm 0.02,\\ B_{\text{f}}=-0.16\pm 0.03,\\ {\delta }_{\text{f}}=-0.015\pm 0.05,\\ {\gamma }_{\text{f}}=0.42\pm \mathrm{0.03.}\end{array}$$(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 η.

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 η, as the main quantity of interest, is not affected by the resulting count-rate bias b_η. This is because the parameters used to characterize the η−M₅₀₀−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 η. 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 $\widehat{\eta }$ given a true cluster mass M₅₀₀ 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 η−M₅₀₀−z relation.

Fig. 6 Bias bηη, which is defined as the ratio of the observed count rate η to the true count rate $\widehat{\eta }$, from the simulation and the best-fit model (middle and left panels, respectively). These are presented as a function of true cluster mass M500 at each redshift, color coded according to In (bη). 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.

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 bWL 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⁴ (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⁵, (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 largescale 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 MWL. We note that the same radial binning is used in deriving the MWL from both synthetic shear profiles and those from the HSC data. Given the true halo mass M500 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{array}{c}\langle \ln b_{\text{WL}}\left(M_{500},z\right)\rangle \equiv \langle \ln \left(\frac{M_{\text{WL}}}{M_{500}}\right)\rangle \\ =\ln A_{\text{WL}}+B_{\text{WL}}\times \ln \left(\frac{M_{500}}{2\times {10}^{14}{h}^{-1}{M}_{\odot }}\right)+{\gamma }_{\text{WL}}\times \ln \left(\frac{1+z}{1+0.6}\right),\end{array}$$(42)

with log-normal intrinsic scatter σ_WL at fixed mass, $${\sigma }_{\text{WL}}\equiv {\left(\text{Var}\left[\ln \left(b_{\text{WL}}|M_{500}\right)\right]\right)}^{\frac{1}{2}}.$$(43)

The resulting constraints on the parameters are $$\begin{array}{c}{A}_{\text{WL}}=0.896\pm 0.076,\\ B_{\text{WL}}=-0.055\pm 0.021,\\ {\gamma }_{\text{WL}}=-0.426\pm 0.087,\\ {\sigma }_{\text{WL}}=0.296\pm \mathrm{0.047.}\end{array}$$(44)

We note that we do not observe significant redshift- or mass-dependent scatter σ_WL in the weak-lensing mass bias; we therefore use the same σ_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{array}{c}{A}_{\text{WL}}=0.903\pm 0.057,\\ B_{\text{WL}}=-0.057\pm 0.022,\\ {\gamma }_{\text{WL}}=-0.474\pm 0.080,\\ {\sigma }_{\text{WL}}=0.238\pm \mathrm{0.032.}\end{array}$$(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.

Fig. 7 Best-fit weak-lensing mass bias bWL normalized to the pivotal mass 2 × 1014h−1M⊙ as a function of redshift. This plot shows the results of the weak-lensing mass calibration including the modeling of the shear profiles at R > 0.2h−1Mpc. 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.

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 η 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 η 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 χ, we have

• A_χ describing the normalization at the pivotal mass M_piv and the pivotal redshift z_piv;

• B_χ describing the overall power-law index in mass;

• γ_χ describing the deviation from the self-similar redshift evolution (Kaiser 1986) of the observable χ;

• δ_χ characterizing the redshift-dependent power-law index in the cluster mass;

• σ_χ describing the log-normal scatter of χ at fixed mass.

In this way, the X-ray observable-to-mass-and-redshift (χ−M₅₀₀−z) relation reads, $$\begin{array}{c}\langle \text{ln}\chi |M_{500}\rangle =\ln A_{\chi }\\ +\left[B_{\chi }+{\delta }_{\chi }\ln \left(\frac{1+z}{1+z_{\text{piv}}}\right)\right]\times \ln \left(\frac{M_{500}}{M_{\text{piv}}}\right)\\ +C_{\text{SS,}\chi }\times \ln \left(\frac{E\left(z\right)}{E\left(z_{\text{Piv}}\right)}\right)+\gamma \chi \times \ln \left(\frac{1+z}{1+z_{\text{piv}}}\right),\end{array}$$(46)

with log-normal intrinsic scatter at fixed mass of $${\sigma }_{\chi }\equiv {\left(\text{Var}\left[\ln \chi |M_{500}\right]\right)}^{\frac{1}{2}},$$(47)

where C_SS,χ is the power-law index of the observable χ in redshift predicted by the self-similar model. In this work, we have five follow-up X-ray observables, namely χ = {L_X, L_b, M_g, T_X, Y_X}. We model the χ−M₅₀₀−z relation separately for each X-ray observable. The self-similar prediction of C_SS,χ for each of them is listed in Table 1.

For the count-rate-to-mass-and-redshift (η−M₅₀₀−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_η calibrated by the simulations (see Sect. 4.1). The resulting η−M₅₀₀−z relation reads, $$\begin{array}{c}\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_{SS,\eta }\times \ln \left(\frac{E\left(z\right)}{E\left(z_{\text{piv}}\right)}\right)+{\gamma }_{\eta }\times \ln \left(\frac{1+z}{1+z_{\text{piv}}}\right)\\ -2\times \ln \left(\frac{D_{\text{L}}\left(z\right)}{D_{\text{L}}\left(z_{\text{piv}}\right)}\right)+\ln \left(b_{\eta }\left(M_{500},z\right)\right),\end{array}$$(48)

with log-normal intrinsic scatter at fixed mass of $${\sigma }_{\eta }\equiv {\left(\text{Var}\left[\ln \eta |M_{500}\right]\right)}^{\frac{1}{2}},$$(49)

Throughout this work, except for T_X and Y_X, we use the pivotal mass M_piv = 1.4 × 10¹⁴h⁻¹ M_⊙ and the pivotal redshift z_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_Piv = 2.5 × 10¹⁴h⁻¹ M_⊙ for χ ∈ {T_X, Y_X}.

5.2 Weak-lensing mass calibration

Our goal is to model the observed X-ray count rate η 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_cl with an observed X-ray count rate η and a shear profile g₊, we compute the probability of observing g₊ given η and a set of the parameter p, namely, $$\begin{array}{c}P\left(g_{+}|\eta ,z_{\text{cl}},\text{p}\right)=\frac{N\left(g_{+},\eta |z_{\text{cl}},\text{p}\right)}{N\left(\eta |z_{\text{cl}},\text{p}\right)}\\ =\frac{{\displaystyle \int P\left(g_{+},\eta |M_{500},z_{\text{cl}},\text{p}\right)n\left(M_{500},z_{\text{cl}}\right)\text{d}{M}_{500}}}{{\displaystyle \int P\left(\eta |M_{500},z_{\text{cl}},\text{p}\right)n\left(M_{500},z_{\text{cl}}\right)\text{d}{M}_{500}}},\end{array}$$(50)

where N (η|z_cl, p) and N η|z_cl, p) denote the numbers of clusters observed with the count rate η and the observable set {g₊, η}, respectively; n (M₅₀₀, z_cl) is the halo mass function evaluated at the true cluster mass M₅₀₀ at the cluster redshift z_cl using the Bocquet et al. (2016) fitting formula with the inclusion of baryons; P (η|M₅₀₀, z_cl, p) and P (g₊, η|M₅₀₀, z_cl, p) describe the probabilities of observing η and{g₊, η}, respectively, for the cluster with the true mass M₅₀₀ at the redshift zcl 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₅₀₀, z_cl, p) in the denominator of Eq. (50) is calculated as $$P\left(\eta |M_{500},z_{\text{cl}},\text{p}\right)={\displaystyle \int \text{d}\tilde{\eta }P\left(\eta |\tilde{\eta }\right)P\left(\tilde{\eta }|M_{500},z_{\text{cl}},\text{p}\right),}$$(51)

where the term $P\left(\tilde{\eta }|M_{500,z_{c\text{l}}},\text{p}\right)$ describes the log-normal distribution of the pristine count rate $\tilde{\eta }$ with the mean following Eq. (48) and the intrinsic scatter σ_η. This characterizes the intrinsic distribution of the pristine count rate while including the observed bias b_η calibrated by the simulations. The term $P\left(\eta |\tilde{\eta }\right)$ accounts for the measurement uncertainty of η with respect to $\tilde{\eta }$. The probability P(g₊, η|M₅₀₀, z_cl, p) reads $$\begin{array}{c}P\left(g_{+},\eta |M_{500},z_{\text{cl}},\text{p}\right)\\ ={\displaystyle \iint \text{d}\tilde{\eta }\text{d}{M}_{\text{WL}}P\left(g_{+}|M_{\text{WL}}\right)P\left(\eta |\tilde{\eta }\right)P\left(M_{\text{WL}},\tilde{\eta }|M_{500},z_{\text{cl}},\text{p}\right),}\end{array}$$(52)

in which the term P (g₊ | M_WL) accounts for the measurement uncertainty of the observed g₊ given the weak-lensing mass M_WL, which is evaluated using Eq. (39). The term $P\left(M_{\text{WL}},\tilde{\eta }|M_{500}{}_{,z_{cl},}\text{p}\right)$ then describes a joint distribution of M_WL and $\tilde{\eta }$ characterized by both the intrinsic scatter of these two quantities and the correlated scatter between them given the underlying true cluster mass M₅₀₀. That is, $P\left(M_{\text{WL}},\tilde{\eta }|M_{500}{}_{,z_{c\text{l}},}\text{p}\right)$ follows a multivariate log-normal distribution of $M_{\text{WL}},\tilde{\eta }$ characterized by an intrinsic covariance matrix, $${\displaystyle {\sum }_{M_{\text{WL,}\eta }}=\left(\begin{array}{cc}{\sigma }_{\text{WL}}^2& {\sigma }_{\text{WL}}{\sigma }_{\eta }{\rho }_{\text{WL,}\eta }\\ {\sigma }_{\text{WL}}{\sigma }_{\eta }{\rho }_{\text{WL,}\eta }& {\sigma }_{\eta }^2\end{array}\right)},$$(53)

where ${\rho }_{{\text{WL,}}_{\eta }}$ is the correlation coefficient between M_WL and $\widehat{\eta }$.

In this way, the joint likelihood of a sample comprising N_cl eFEDS clusters, given the set of observables ${\left\{g_{+,i,}{\eta }_i\right\}}_{i=1}^{N_{c\text{l}}}$, is calculated as $$\mathcal{L}\left(\text{p}\right)={\displaystyle \prod _{1=1}^{N_{\text{cl}}}P\left(g_{+j}|{\eta }_i,\text{p}\right).}$$(54)

In this framework, we have a parameter set p consisting of 19 free parameters, that is, $$\begin{array}{l}\text{p}={\text{p}}_{\eta }\cup {\text{p}}_{\text{F}}\cup {\text{p}}_{\text{WL}}\cup {\text{p}}_{\text{mis}}\cup \left\{{\rho }_{\text{WL,}\eta }\right\}\cup {\text{p}}_{\text{c}},\hfill \\ {\text{p}}_{\eta }=\left\{A_{\eta },B_{\eta },{\delta }_{\eta },{\gamma }_{\eta },{\sigma }_{\eta }\right\},\hfill \\ {\text{p}}_{\text{F}}=\left\{A_{\text{f}},B_{\text{f}},{\gamma }_{\text{f}},{\sigma }_{\text{f}}\right\},\hfill \\ {\text{P}}_{\text{WL}}=\left\{A_{\text{WL}},B_{\text{WL}},{\gamma }_{\text{WL}},{\sigma }_{\text{WL}}\right\},\hfill \\ {\text{p}}_{\text{mis}}=\left\{f_{\text{mis}},{\sigma }_{\text{mis}}\right\},\hfill \\ {\text{p}}_{\text{c}}=\left\{{\Omega }_{\text{M}},H_{\text{0}},{\sigma }_{\text{8}}\right\},\hfill \end{array}$$(55)

where p_η describes the η−M₅₀₀−z relation as in Eqs. (48); p_F and p_WL describe the b_η−M₅₀₀−z and M_WL−M₅₀₀−z relations in Eqs. (40) and (42), respectively; and p_mis characterizes the miscentering distribution as in Eq. (20); p_c consists of the cos-mological parameters, i.e., Ω_m, H₀, and σ₈, which we vary in the modeling below.

5.3 Joint modeling of X-ray scaling relations

A joint modeling of an X-ray observable χ 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_cl with an available set of observables {χ, g₊, η}, we evaluate the probability of observing the cluster with χ and g₊, given η and the set of parameters p. That is, $$\begin{array}{ll}P\left(\chi ,g_{+}|\eta ,z_{\text{cl}},\text{p}\right)\hfill & =\frac{N\left(\chi ,g_{+},\eta |z_{\text{cl}},\text{p}\right)}{N\left(\eta |z_{\text{cl}},\text{p}\right)}\hfill \\ \hfill & =\frac{{\displaystyle \int P\left(\chi ,g_{+},\eta |M_{500},z_{\text{cl}},\text{p}\right)n\left(M_{500},z_{\text{cl}}\right)\text{d}{M}_{500}}}{{\displaystyle \int P\left(\eta |M_{500},z_{\text{cl}},\text{p}\right)n\left(M_{500},z_{\text{cl}}\right)\text{d}{M}_{500}}},\hfill \end{array}$$(56)

where P (χ, g₊, η|M₅₀₀, z_cl, p) is the joint probability of observing the cluster with the observables of {χ, g₊,η} at the true cluster mass M₅₀₀ and redshift z_cl, given the parameter p. This term is calculated as $$\begin{array}{l}P\left(\chi ,g_{+},\eta |M_{500,}{z}_{\text{cl}},\text{p}\right)={\displaystyle \int {\displaystyle \int {\displaystyle \int \text{d}\tilde{\eta }}\text{d}{M}_{\text{WL}}\text{d}\tilde{\chi }}}\hfill \\ P\left(\chi |\tilde{\chi }\right)P\left(g_{+}|M_{\text{WL}}\right)P\left(\eta |\tilde{\eta }\right)P\left(\tilde{\chi },M_{\text{WL}},\tilde{\eta }|M_{500},z_{\text{cl}},\text{p}\right).\hfill \end{array}$$(57)

In an analogy to Eq. (52), the term $P\left(\tilde{\mathcal{X}},M_{\text{WL}},\tilde{\eta }|M_{500}{}_{,z_{c\text{l}},}\text{p}\right)$ describes the intrinsic distribution of $\tilde{\mathcal{X}},M_{\text{WL}},\tilde{\eta }$ given the true cluster mass M₅₀₀, assuming a multivariate log-normal distribution characterized by an intrinsic covariance matrix, $${\displaystyle {\sum }_{M_{\text{WL}},\chi ,\eta }=\left(\begin{array}{ccc}{\sigma }_{\text{WL}}^2& {\sigma }_{\text{WL}}{\sigma }_{\chi }{\rho }_{\text{WL,}\chi }& {\sigma }_{\text{WL,}\eta }\\ {\sigma }_{\text{WL}}{\sigma }_{\chi }{\rho }_{\text{WL,}}{}_{\chi }& {\sigma }_{\chi }^2& {\sigma }_{\chi }{\sigma }_{\eta }{\rho }_{\chi ,\eta }\\ {\sigma }_{\text{WL}}{\sigma }_{\eta }{\rho }_{\text{WL,}}{}_{\eta }& {\sigma }_{\chi }{\sigma }_{\eta }{\rho }_{\chi ,\eta }& {\sigma }_{\eta }^2\end{array}\right),}$$(58)

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

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 χ given the count rate η, namely, $$\begin{array}{ll}P\left(\chi |\eta ,z_{\text{cl}},\text{p}\right)\hfill & =\frac{N\left(\chi ,\eta |z_{\text{cl}},\text{p}\right)}{N\left(\eta |z_{\text{cl}},\text{p}\right)}\hfill \\ \hfill & =\frac{{\displaystyle \int P\left(\chi ,\eta |M_{500},z_{\text{cl}},\text{p}\right)n\left(M_{500},z_{\text{cl}}\right)\text{d}{M}_{500}}}{{\displaystyle \int P\left(\eta |M_{500},z_{\text{cl}},\text{p}\right)n\left(M_{500},z_{\text{cl}}\right)\text{d}{M}_{500}}}.\hfill \end{array}$$(59)

We note that Eq. (59) is identical to Eq. (50) with a substitution of g₊ with χ. Accordingly, the probability P (χ, η|M₅₀₀, z_cl, p) is calculated as $$\begin{array}{l}P\left(\chi ,\eta |M_{500},z_{\text{cl}},\text{p}\right)\hfill \\ ={\displaystyle \int {\displaystyle \int \text{d}\tilde{\eta }\text{d}\tilde{\chi }P\left(\chi |\tilde{\chi }\right)P\left(\eta |\tilde{\eta }\right)P\left(\tilde{\chi },\tilde{\eta }|M_{500},z_{\text{cl}},\text{p}\right),}}\hfill \end{array}$$(60)

where $P\left(\tilde{\mathcal{X}},\tilde{\eta }|M_{500}{}_{,z_{c\text{l}},}\text{p}\right)$ describes the intrinsic distribution of $\left\{\tilde{\mathcal{X}},\tilde{\eta }\right\}$ characterized as a two-dimensional log-normal distribution with a covariance matrix of $${\displaystyle {\sum }_{\chi ,\eta }=\left(\begin{array}{cc}{\sigma }_{\chi }^2& {\sigma }_{\chi }{\sigma }_{\eta }{\rho }_{\chi ,}{}_{\eta }\\ {\sigma }_{\chi }{\sigma }_{\eta }{\rho }_{\chi ,\eta }& {\sigma }_{\eta }^2\end{array}\right).}$$(61)