© The Authors 2025

1. Introduction

The extended ROentgen Survey with an Imaging Telescope Array (eROSITA) X-ray telescope (Predehl et al. 2021), the primary soft-X-ray instrument of the Russian-German Spectrum-Roentgen-Gamma (SRG) space observatory (Sunyaev et al. 2021) launched in 2019, has revolutionized the study of galaxy clusters with its groundbreaking achievements in carrying out the eROSITA All-Sky Survey (hereafter eRASS; Merloni et al. 2012, 2024), which provides the deepest all-sky imaging in X-rays to date. The main goal of the eROSITA survey is to advance cosmological studies by constructing the largest sample of galaxy clusters and investigating the growth of their populations over cosmic time. The population growth of galaxy clusters is closely related to the cosmic structure formation (Bardeen et al. 1986; Allen et al. 2011; Kravtsov & Borgani 2012) and places stringent constraints on cosmological parameters, such as the mean matter fraction Ω_m, the degree of density-field inhomogeniety σ₈, and the equation of state of dark energy (Weinberg et al. 2013; Huterer et al. 2015).

Prior to the eROSITA era, cluster cosmology played a central role in cosmological analyses across multiple wavelengths, including X-rays (Reiprich & Böhringer 2002; Vikhlinin et al. 2009a; Mantz et al. 2015; Schellenberger & Reiprich 2017; Garrel et al. 2022), the optical (Costanzi et al. 2019a; To et al. 2021; Sunayama et al. 2024), and the millimetre wavelength (Salvati et al. 2022; Bocquet et al. 2024a) via the Sunyaev–Zeldovich effect (SZE; Sunyaev & Zel’dovich 1970; Sunyaev & Zel’dovich 1972). Moreover, samples selected by weak gravitational lensing (hereafter weak lensing or WL) have gained increasing attention in recent years as a promising probe (Chiu et al. 2024; Chen et al. 2025). Cluster cosmology heavily relies on the accurate determination of the cluster mass (Pratt et al. 2019), for which significant progress has been made in utilizing weak lensing as a reliable method (Okabe et al. 2010; Umetsu et al. 2014; von der Linden et al. 2014; Applegate et al. 2014; Okabe & Smith 2016; Schrabback et al. 2018; Dietrich et al. 2019). In the past decade, the deployments of Stage III weak lensing surveys, such as the Dark Energy Survey (DES; The Dark Energy Survey Collaboration 2005), the Kilo-Degree Survey (KiDS; Kuijken et al. 2015), and the Hyper Suprime-Cam (HSC) Subaru Strategic Program (Aihara et al. 2018a), have enabled the weak-lensing mass calibration of sizable samples using large and homogeneous datasets with a consistent methodology (Simet et al. 2017; McClintock et al. 2019; Murata et al. 2019; Bellagamba et al. 2019; Miyatake et al. 2019; Umetsu et al. 2020; Chiu et al. 2022; Shirasaki et al. 2024; Bocquet et al. 2024b). In particular, Chiu et al. (2023) conducted the first cosmological study using eROSITA clusters selected during the performance verification phase of the eRASS survey in synergy with the HSC weak-lensing mass calibration, demonstrating the potential of eROSITA-based cluster cosmology when it is combined with the Stage III weak lensing surveys.

In 2024, the first eROSITA Data Release (DR1; Merloni et al. 2024) from the German eROSITA Consortium (hereafter eROSITA-DE) delivered the largest sample of X-ray-selected clusters to date (Bulbul et al. 2024; Kluge et al. 2024) based on the first all-sky scan (eRASS1). This milestone also led to the tightest cosmological constraints obtained by far from cluster abundance (Ghirardini et al. 2024), enabled by joint weak-lensing mass calibrations from the surveys of HSC (this work), DES (Grandis et al. 2024), and KiDS (Kleinebreil et al. 2025). By leveraging the study of Ghirardini et al. (2024), constraints on various non-concordance cosmological models were also obtained (Artis et al. 2024, 2025). Importantly, the X-ray selection function of eRASS1 clusters has been precisely characterized (Clerc et al. 2024) using eRASS digital-twin simulations (Seppi et al. 2022). The combination of the largest cluster sample selected via the X-ray emission of intracluster medium (ICM), the accurate mass calibration with the state-of-the-art weak lensing datasets, and the precise characterization of the selection function has enabled unprecedented precision in cluster studies using the eROSITA sample, for example active galactic nucleus feedback in cluster cores (Bahar et al. 2024), the halo assembly bias of superclusters (Liu et al. 2024), the properties of halo clustering (Seppi et al. 2024), and the halo concentration and morphology of massive clusters (Okabe et al. 2025).

We note that eROSITA clusters are primarily selected in X-rays with minimal dependence on the characteristics of their galaxy populations, thus making them an ideal sample for studying the formation and evolution of galaxies hosted by massive halos. X-ray-selected samples also offer two advantages over other selection methods. First, they probe a mass regime generally much lower than that of SZE-selected samples, thanks to the deep and high-resolution imaging in X-rays. Second, by tracing the distribution of hot ICM, X-ray-selected samples are not affected by the projection effect, which currently poses severe challenges for the modelling of optically selected clusters (Costanzi et al. 2019b; Sunayama et al. 2020).

By taking full advantage of exquisite imaging from the HSC survey, this study aims to (1) perform the weak-lensing mass calibration, (2) measure the stellar mass M_⋆, BCG of their brightest cluster galaxies (BCGs), and (3) derive the stellar-mass-to-mass-and-redshift (M_⋆, BCG–M–z) relation for eROSITA clusters in the first-year (eRASS1) sample. Given the mass range of the eRASS1 clusters, the M_⋆, BCG–M–z relation corresponds to the high-mass end of the stellar-mass-to-halo-mass relation, a terminology widely used in the literature. Similarly to Leauthaud et al. (2012) and Coupon et al. (2015), a key feature of this work is that we simultaneously model the weak-lensing mass and the stellar-mass-to-halo-mass relation on a basis of individual clusters, allowing the direct determination of the mass and redshift trends of the M_⋆, BCG–M–z relation without relying on the abundance-matching method (Vale & Ostriker 2004; Kravtsov et al. 2004; Conroy et al. 2006), which depends on numerical simulations. The major strengths of this work are (1) the uniform X-ray selection of eRASS1 clusters extending to redshift z ≈ 1 without the selection bias toward galaxy properties and (2) the well-quantified selection function and the halo mass determination, ensuring that the resulting cosmological constraints are in acceptable agreement with those from other independent probes, thereby forming a cosmology-verified cluster sample suitable for astrophysical studies. During this study we identified 124 eRASS1 clusters within the common footprint of the eRASS and HSC surveys, covering an area of ≈500 deg².

This paper is organized as follows. In Sect. 2 we describe the eRASS1 clusters and the HSC datasets used in this work. The measurements of weak-lensing observables and the BCG stellar mass are presented in Sect. 3. The modelling is described in Sect. 4. We present and discuss the results in Sect. 5, and draw our conclusions in Sect. 6. Unless stated otherwise, the halo mass adopts the definition of M_500c, which is the mass enclosed by a sphere with a radius R_500c wherein the average mass density is 500 times the cosmic critical density ρ_c(z) at the cluster redshift z. In this work we also make use of the symbol M for the halo mass interchangeably with M_500c. We define the notation 𝒩(x,y²) (𝒰(a,b)) as a Gaussian distribution with the mean x and variance of y² (flat distribution between a and b). We adopt a concordance flat ΛCDM cosmological model with the standard cosmological parameters (Ω_m = 0.3, σ₈ = 0.8, and H₀ = 70 km/s/Mpc), which we allow to vary within reasonable ranges in the modelling.

2. The cluster sample and data

We describe the eRASS1 sample of galaxy clusters in Sect. 2.1. The HSC datasets used in this work are summarized in Sect. 2.2.

2.1. The cluster sample

We used the sample of galaxy clusters overlapping the footprint of the HSC survey from the cluster catalogue released in the eROSITA-DE DR1 (Bulbul et al. 2024). Specifically, the selection of the clusters is identical to that of the cosmological subsample used in the eRASS1 cosmological analysis (Ghirardini et al. 2024) with (1) the extent likelihood L_ext of L_ext > 6, (2) the galaxy richness λ of λ > 3, and (3) the redshift range of 0.1 < z < 0.8 determined from the optical confirmation (Kluge et al. 2024). The extent likelihood L_ext, a quantity returned by the eROSITA pipeline, describes the normalized logarithmic likelihood of a detected source being a spatially extended object rather than a point source. In the interest of uniformity, the photometric redshifts of the galaxy clusters are consistently adopted with the accuracy of δz/(1+z) < 0.005 calibrated using spectroscopic samples (Kluge et al. 2024). With such high accuracy, the uncertainty of the cluster redshift is negligible for the purpose of this work. As a result, there are 124 eRASS1 clusters in the common footprint between the eRASS and HSC surveys, as the final sample studied in this work. The sky distribution of the sample is shown in Fig. 1.

Fig. 1. Footprint of the HSC-Y3 weak lensing dataset (grey dots) and the sky distributions of the eRASS1 clusters together with those categorized based on the availability of weak-lensing shear profile g+ and/or the BCG stellar mass estimate M⋆, BCG. The eRASS1 clusters with available g+ are marked as blue circles (23 clusters), those with M⋆, BCG as green diamonds (28 clusters), and those with both measurements as red squares (73 clusters). The other eRASS1 clusters with neither g+ nor M⋆, BCG are shown as grey crosses, and are excluded from this study.

Leveraging the realistic mock simulations (Seppi et al. 2022, see also Comparat et al. 2020), the initial contamination of the cluster sample is estimated to be at a level of ≈6% with the X-ray selection of L_ext > 6 alone. After the optical confirmation, the resulting contamination is reduced to ≈5% (Kluge et al. 2024; Bulbul et al. 2024), which is subdominant compared to the Poissonian noise of our sample (at a level of ≈9%), and therefore negligible in this work. An independent modelling of the same weak-lensing shear profiles while accounting for the sample contamination in Kleinebreil et al. (2025) returned a result that is statistically consistent with ours (see Sect. 5.1), reinforcing the picture that the contamination in our sample is subdominant in this work.

The observed and corrected X-ray photon count rate (hereafter count rate C_R) in the 0.2–2.3 keV band is used as the mass proxy (Bulbul et al. 2024). The X-ray count rate C_R has been shown as a reliable mass proxy for eROSITA-detected clusters with well-characterized scaling with the cluster mass (Chiu et al. 2022; Grandis et al. 2024; Kleinebreil et al. 2025), although with large scatter (Ghirardini et al. 2024) that deserves being further studied in a future work. For each cluster, the selection function – the probability $\mathcal{I}({\widehat{C}}_{\mathrm{R}}|z_{\mathrm{cl}},\mathcal{H})$ of the cluster being detected – is modelled as a function of the intrinsic count rate ${\widehat{C}}_{\mathrm{R}}$ and evaluated at the cluster redshift z_cl and the sky location $\mathcal{H}$ (Clerc et al. 2024, see also Eq. (1) in Ghirardini et al. 2024). The intrinsic count rate ${\widehat{C}}_{\mathrm{R}}$ refers to the observed one C_R before it is affected by observational noise. The selection function is included in the modelling of the scaling relation in Sect. 4.2.

The distribution of the clusters in the observable space of the X-ray count rate C_R and redshift is shown in the upper panel of Fig. 2, where the colors represent the availability of the measurements described in the following subsection.

Fig. 2. Distributions of the observed count rate CR (top panel) and the BCG stellar mass M⋆, BCG (left panel) as functions of the cluster redshift for the eRASS1 clusters studied in this work. The clusters are colour-coded as in Fig. 1. We note that a clear dependence of the BCG stellar mass on the cluster redshift is revealed in the bottom panel, which arises from the X-ray selection primarily favoring massive clusters at high redshift.

2.2. The HSC datasets

We used the data from the HSC survey to (1) measure the weak-lensing shear profiles and hence the weak-lensing mass, and to (2) extract the stellar mass of the BCGs. In what follows, a brief description about the datasets is provided.

In terms of the weak-lensing measurements, we used the latest HSC Three-Year (Y3) shape catalogues (Li et al. 2022) that were constructed using the i-band imaging acquired during 2014 and 2019 with a mean seeing of 0.59 arcsec. The catalogue has a limiting i-band magnitude cut of i < 24.5 mag, resulting in an average density at a level of ≈20 galaxies/arcmin². The shape measurement was rigorously calibrated against image simulations (see e.g. Mandelbaum et al. 2018a), delivering the multiplicative bias m_γ at a level of |m_γ|< 9 × 10⁻³ independently of redshift up to z ≈ 3. Various null tests were examined to be consistent with zero or within the requirements for cosmic-shear analyses, ensuring sufficiently accurate weak-lensing measurements for cluster studies. The HSC-Y3 shape catalogue covers a footprint area of ≈500 deg² and is divided into six fields, namely GAMA09H, GAMA15H, XMM, VVDS, HECTOMAP, and WIDE12H. The clusters selected in eRASS1 are only located in the fields of GAMA09H, GAMA15H, and WIDE12H, of which the datasets are used in this work.

The photometric redshift (photo-z) of the weak-lensing source sample is needed to interpret the weak-lensing measurements. Specifically, we used the photo-z estimated by the machine-learning-based code DEmP (Direct Empirical Photometric method; Hsieh & Yee 2014), whose performances were fully quantified in Nishizawa et al. (2020). In short, the bias, scatter, and the outlier fraction of the photo-z estimates for the source galaxies are estimated to be at levels of |Δz|≈0.003(1+z), σ_Δz ≈ 0.019(1+z), and ≈5.4%, respectively, where Δz is defined as the difference between the photometric and spectroscopic redshifts. The full distribution P(z) of the photo-z estimates was used to model the weak-lensing signals.

In Kluge et al. (2024), the BCGs of eRASS1 clusters were identified as the brightest passive member galaxies within a characteristic radius. In addition, the authors estimated that approximately 85% of eRASS1 clusters’ BCGs coincide with the optical centers. To ensure a homogeneous identification of BCGs, we visually select each cluster’s brightest galaxy as the BCG using the HSC image. Specifically, we selected the BCG of each eRASS1 cluster as the brightest galaxy whose color is similar to that of the galaxy overdensity near the X-ray centre. As a result, we find that 25 (out of 101) eRASS1 clusters in our sample have their BCGs that differ from those automatically selected in Kluge et al. (2024). As a validation, we compared the photo-z estimates (photoz_best) of the BCGs from the HSC survey with the cluster redshifts (Z_LAMBDA) and calculated the BCG photo-z bias defined as Δz_BCG ≡ (photoz_best−Z_LAMBDA)/(1+Z_LAMBDA). We find that the bias Δz_BCG has a median value of 0.0028 and a standard deviation of 0.015¹, indicating that the BCG identification is robust in this work.

To estimate the stellar mass M_⋆, BCG of BCGs, we used the forced cmodel photometry at the grizY broadband queried from the HSC Public Data Release 3 (PDR3). No attempts to include the intracluster light (ICL) in M_⋆, BCG were made. We excluded the BCGs with pixels (in one of the grizY images) that are contaminated, saturated, or masked due to bright stars, or that have failures in the cmodel photometric fitting. This leads to a sample of BCGs with clean cmodel photometry. Aiming for the high signal-to-noise ratio and uniform HSC photometry, we only included the BCGs with at least two exposures in all grizY bands. Appendix A contains the flags used in querying the photometry of the BCGs. Appendix B shows the cutout images of these BCGs studied in this work. The identified sky locations of the BCGs are tabulated in Table B.1.

It is important to note that the cmodel algorithm adopts a composite model, where the bulge and disk components are described by de Vaucouleurs (1948) and exponential profiles, respectively. Consequently, the cmodel photometry does not fully capture the stellar distribution of the extended, diffuse ICL at large radii around BCGs. As a result, the cmodel photometry is expected to underestimate the “total” flux of a BCG+ICL system (see also Akino et al. 2022). In fact, the exact definition of the “total” stellar mass of BCGs becomes nuanced, as they are embedded within the surrounding ICL. A common approach to obtain the “total” stellar mass of a BCG+ICL system is to measure the flux within a large aperture, followed by a conversion from the observed light to the mass with an assumed mass-to-light ratio. In Huang et al. (2018a), they estimated the stellar mass profiles of BCGs using a series of well defined apertures and found that the stellar masses estimated from the cmodel photometry were systematically lower than those enclosed within 100 kpc, precisely because the cmodel photometry does not account for the extended ICL component at large radii. By stacking the light profiles of ≈3000 clusters in Chen et al. (2022), the authors found that the stellar mass of a BCG alone is well approximated by the stellar mass of the bulge/disk component enclosed within ≈50 kpc, beyond which the BCG component slowly transitions to the ICL-dominated regime at ≳100 kpc. In this study, where the cmodel photometry is used, we effectively estimate only the stellar mass of the bulge/disk components of a BCG+ICL system (within ≲50 kpc).

We note that the photometric catalogue of the BCGs is constructed independently of the HSC-Y3 shape catalogue, such that the availability of the optical photometry for the BCGs is not subject to the selection of the weak-lensing sources. Consequently, some clusters have available HSC grizY photometry to estimate the BCG stellar mass but no data to obtain the weak-lensing mass, and vice versa. The final sample of 124 eRASS1 clusters includes 23 and 28 systems with only weak-lensing data and BCG photometry, respectively; the remaining 73 clusters have both datasets.

The distributions of the X-ray count rate C_R, the BCG stellar mass M_⋆, BCG (which we derive in Sect. 3.2), and the cluster redshift z_cl are shown in Fig. 2 with colors indicating the availability of the HSC datasets.

3. Measurements

For individual clusters with the available datasets, we extracted two kinds of measurements, the weak-lensing shear profile (Sect. 3.1) and the BCG stellar mass (Sect. 3.2).

3.1. Weak-lensing measurements

To extract the weak-lensing measurements, we followed the identical procedure as detailed in Chiu et al. (2022), in which they used the same HSC-Y3 datasets to derive the tangential reduced shear profiles g₊(θ) as a function of the angular radius θ to calibrate the halo mass of clusters selected in the eROSITA Final Equatorial-Depth Survey (eFEDS; Liu et al. 2022). Later, the same weak-lensing analysis was used to derive cosmological constraints using the eFEDS cluster abundance (Chiu et al. 2023). Here, the weak-lensing data products produced in this work have been closely examined against other Stage III surveys (Kleinebreil et al. 2025) and used in the eRASS1 cosmological analyses (Ghirardini et al. 2024). In what follows, we provide a summary of the weak-lensing measurements.

On a basis of individual clusters, the weak-lensing measurements are composed of the tangential reduced shear profile g₊(θ) (with the measurement uncertainty) and the corresponding redshift distribution of the source galaxies. A galaxy was selected as a weak-lensing source of a cluster at redshift z_cl based on its full redshift distribution P(z), requiring it to be securely located at the background with ∫_{zcl + 0.2}P(z)dz > 0.98. With this “P-cut” selection (Oguri 2014; Medezinski et al. 2018), the resulting source sample is extremely clean without cluster member contamination out to redshift z ≈ 1.2 (see Fig. 3 in Chiu et al. 2022). Following Eq. (8) in Chiu et al. (2022), the shear profile g₊(θ) was extracted around the eROSITA-defined X-ray centre with ten logarithmic radial bins between 0.2 h⁻¹ Mpc and 3.5 h⁻¹ Mpc calculated in a fixed flat ΛCDM cosmology with Ω_m = 0.3 and h = 0.7. The three innermost bins were discarded in the weak-lensing mass modelling to avoid systematics, resulting in total seven logarithmic bins at a radial range of 0.5 h⁻¹ Mpc ≲ R < 3.5 h⁻¹ Mpc. We note that we converted the unit of the radial binning to the angular space, i.e. in arcmin, and self-consistently fitted the model of the cluster mass profile at a correct physical scale inferred from the redshift-distance relation while the cosmological parameters are varying in the forward modelling. The uncertainty of the tangential shear profile is expressed by the weak-lensing covariance matrix (see Eq. (14) in Chiu et al. 2022); however, we only include the shape noise, as opposed to that in Chiu et al. (2022) containing the scatter due to uncorrelated large-scale structures (Hoekstra 2003). The scatter raised from uncorrelated large-scale structures is accounted for in calibrating the weak-lensing mass bias (see Sect. 4.3 and Grandis et al. 2024). The stacked and lensing-weight-weighted P(z) of the source galaxies for each cluster is used to infer the lensing efficiency (see Eq. (34) in Chiu et al. 2022).

We describe the modelling of these shear profiles to calibrate the cluster halo mass in Sect. 4.3.

3.2. Measurements of BCG stellar masses

Similarly to the procedure in Chiu et al. (2018, see also Erfanianfar et al. 2019), we estimated the BCG stellar mass M_⋆, BCG using the SED (spectral energy distribution) fitting technique implemented in the χ² template-fitting code, Le Phare (Arnouts et al. 1999; Ilbert et al. 2006)². We describe the analysis steps below.

We used the five broadband photometry grizY from the PDR3 database of the HSC survey. For each BCG, the template fitting is fixed to the cluster redshift z_cl. The SED fitting has two passes. The first is to account for the zero-point offsets in the photometry. Specifically, we built the galaxy templates from the COSMOS field (Ilbert et al. 2009) with the Prevot et al. (1984) and Calzetti et al. (2000) dust-extinction laws³ with E(B−V) in a range between 0 and 0.5. Once the best-fit template was determined for each BCG, a systematic zero-point offset of each band was computed to minimize the differences between the models and observed magnitudes among all the BCGs. This process was iterated until a convergence is reached. The second pass of the SED fitting is to obtain the BCG stellar mass by fitting the galaxy templates predicted by the Stellar Population Synthesis (SPS). We built the SPS templates using the Bruzual & Charlot (2003) model with the Chabrier (2003) initial mass function (IMF). Three metallicities (Z = 0.004, 0.008, 0.02) are used and configured with the exponentially decaying star-formation rate characterized by an e-folding timescale τ. We used six values of τ up to 5 Gyr (i.e. $\frac{\tau }{\mathrm{Gyr}}=0.1,0.3,1,2,3,5$) because BCGs are expected to be dominated by old stellar populations without recent star-forming activities. This is supported by our data, as extending the timescale τ to 30 Gyr to include star-forming galaxy templates does not change the final results (see Sect. 5). Moreover, the galaxy templates are configured with the Calzetti et al. (2000) dust-extinction law of E(B−V) = 0, 0.1, 0.2, 0.3. The ages of these galaxy templates can have a maximum value of 14 Gyr. When fitting the Bruzual & Charlot (2003) templates, we apply the zero-point offsets obtained from the first pass to the individual bands and add the same amounts of systematic uncertainties to the measurement uncertainties in quadrature. We note that this effectively weakens the constraining power of the bands with a larger offset. Last, we increased the flux uncertainties (after adding the systematic uncertainties) by a factor of 6.5, which is suggested by our data, such that the median of the resulting reduced χ_red² of the BCGs is unity. This does not change the best-fit stellar mass but only increases the measurement uncertainty of the final mass estimate. We note that the statistical uncertainties of the photometric measurements of the BCGs in the grizY bands are at a level of 10⁻³ mag, and that the SED templates predicted by SPS are not perfect. As a result, even small discrepancies between the model and observation can lead to significantly large χ_red². By increasing the flux uncertainties while requiring the median χ_red² of the sample to be unity, we effectively account for the uncertainty in the derived stellar mass arising from the limitations of the SED templates. The measurement uncertainty of M_⋆, BCG is obtained by bootstrapping the Bruzual & Charlot (2003) template fitting by 1000 times, assuming that the flux uncertainties are Gaussian.

We estimated the systematic uncertainty in our stellar mass measurements as follows. We constructed a catalogue of galaxies observed in the COSMOS field from the HSC data base, estimate their stellar masses at known redshifts using the same method mentioned above, and compare the resulting stellar masses M_⋆, HSC to those in the latest COSMOS2020 catalogue (M_{⋆, COSMOS}) estimated using 30-band photometry (Weaver et al. 2022). In the stellar mass and redshift range occupied by the BCGs of the eRASS1 clusters ($M_{\star }\gtrsim {10}^{10.7}{\mathrm{M}}_{\odot }$ and z ≲ 1), we find a systematic offset in the stellar mass at a level of ≈0.1 dex, i.e. ⟨log(M_⋆, HSC/M_⊙)−log(M_{⋆, COSMOS}/M_⊙)⟩ ≈ 0.1 dex. Moreover, this systematic offset does not depend on the redshift. Therefore, we expect that (1) there exists a systematic uncertainty at a level of ≈0.1 dex in our M_⋆, BCG estimates with respect to those measured with the full configuration of SED templates on high-quality data such as in the COSMOS field, and that (2) the systematic uncertainty only affects the normalization but not the mass or redshift trends of the scaling relation.

The resulting BCG stellar mass estimates M_⋆, BCG are plotted as a function of the cluster redshift in Fig. 2. We present the results of the SED fitting in Appendix B. The measurements of the BCG setllar mass M_⋆, BCG of individual eRASS1 clusters are contained in Table B.1.

4. Modelling

Our goal is to constrain the scaling relations of the observables, including the X-ray count rate C_R and the BCG stellar mass M_⋆, BCG. We first provide an overall framework of the modelling in Sect. 4.1, and refer to Chiu et al. (2022, see also Bocquet et al. 2019 and Chiu et al. 2023) for more details. Next, we describe the modelling of the scaling relations in Sect. 4.2, followed by the description of the weak-lensing modelling in Sect. 4.3. Finally, we give the statistical inference to obtain the parameter posteriors in Sect. 4.4.

4.1. Modelling framework

The cluster sample is constructed based on the X-ray selection (L_ext > 6) with the aid of the optical confirmation to remove the contamination with unphysically low richness (λ ≤ 3), together with the selection of the cluster redshift at 0.1 < z < 0.8 (see Sect. 2.1 and Sect. 2.1 in Ghirardini et al. 2024). The extent likelihood L_ext is mainly subject to the parameters tuned to optimize in the cluster-finding algorithm and lacks a straightforward association with the physical properties of clusters. Therefore, it is not an ideal mass proxy. Instead, we use the observed X-ray count rate C_R as the mass proxy. Meanwhile, the intrinsic X-ray count rate ${\widehat{C}}_{\mathrm{R}}$ acts as an X-ray “selection” observable. The observed C_R and intrinsic ${\widehat{C}}_{\mathrm{R}}$ are related through the measurement uncertainty, and both can serve as selection observables.

In short, there are three selection observables $({\widehat{C}}_{\mathrm{R}},\lambda ,z)$, where the uncertainty of the cluster redshift z is negligible, and the richness selection mainly removes obvious contamination with an insignificant impact on the resulting sample (Ghirardini et al. 2024). That is, the cluster selection is primarily determined in X-rays.

The follow-up observations of individual clusters are then conducted to obtain the follow-up observables, which are the weak-lensing shear profile g₊ and the BCG stellar mass M_⋆, BCG in this study. It is extremely important to note that the follow-up observables are uniformly obtained for every cluster and that the lack of these follow-up observables for some clusters does not rely on any preferences for X-ray properties. In other words, no additional selections are made on the weak-lensing shear profile g₊ and the BCG stellar mass M_⋆, BCG.

We denote the selection and follow-up observables by the symbols of 𝒮 and 𝒪, where 𝒮 ≡ {C_R,λ,z} and 𝒪 ≡ {g₊,M_⋆, BCG}, respectively. For the scaling relations parametrized by a parameter vector p, we maximize the probability of observing the follow-up observable 𝒪 for a cluster selected with the selection observable 𝒮 as

$$\begin{array}{c}\hfill P(\mathcal{O}|\mathcal{S},\mathbf{p})=\frac{N(\mathcal{O},\mathcal{S}|\mathbf{p})}{N(\mathcal{S}|\mathbf{p})},\end{array}$$(1)

where the label N represents the number of clusters with the observables of 𝒮 and/or 𝒪 evaluated at the parameter vector p.

Equation (1) is referred to as the “mass calibration likelihood” widely used in the population modelling of galaxy clusters (e.g. Chiu et al. 2022, 2023; Ghirardini et al. 2024). Consider the selection function as a function of only the observed quantity 𝒮 rather than the intrinsic one $\widehat{\mathcal{S}}$. Because the selection of the clusters purely depends on 𝒮 instead of 𝒪, the selection-function factors of both the numerator and denominator in Eq. (1) are cancelled⁴ and have no effect (see Eq. (3) in Chiu et al. 2023). That is, the selection function of the cluster sample is not important for the mass calibration likelihood, i.e. Eq. (1), as long as no additional selection is made on the follow-up observable 𝒪, as in this study.

In this study and other eRASS1 work, however, the selection function $\mathcal{I}\left({\widehat{C}}_{\mathrm{R}}|z,\mathcal{H}\right)$ is calibrated at the sky location ℋ and redshift z of each cluster in terms of the intrinsic observable, namely the intrinsic count rate ${\widehat{C}}_{\mathrm{R}}$ (Clerc et al. 2024). Therefore, the selection function is implicitly used in calculating both the numerator and denominator of Eq. (1) and cannot be canceled. We note that the selection function $\mathcal{I}\left({\widehat{C}}_{\mathrm{R}}|z,\mathcal{H}\right)$ does not depend on the parameter vector p, which is distinct from the self-calibrated selection function used in Chiu et al. (2023).

Because (1) we ignore the uncertainty of the cluster redshift and (2) the richness selection of λ > 3 has no impact on the resulting X-ray selected sample (Ghirardini et al. 2024), Eq. (1) can be approximated with 𝒮 = {C_R,z} as (e.g. Eq. (56) in Chiu et al. 2022)

$$\begin{array}{cc}\hfill P(\mathcal{O}|\mathcal{S},\mathbf{p})& \approx P(M_{\star ,\mathrm{BCG}},g_{+}|C_{\mathrm{R}},z,\mathbf{p})\hfill \\ \hfill & =\frac{\int P(M_{\star ,\mathrm{BCG}},g_{+},C_{\mathrm{R}}|M,z,\mathbf{p})n(M|z,\mathbf{p})\mathrm{d}M}{\int P\left(C_{\mathrm{R}}|M,z,\mathbf{p}\right)n(M|z,\mathbf{p})\mathrm{d}M},\hfill \end{array}$$(2)

where n(M|z,p) is the halo mass function evaluated at the cluster redshift z and the parameter vector p.

In Eq. (2), the factor P(C_R|M,z,p) describes the probability of observing the X-ray count rate C_R given the halo mass M at the redshift z evaluated with the parameter vector p. It includes two components making the observed count rate C_R scatter around the mean value predicted at the fixed cluster mass and redshift. The first is the intrinsic scatter in the count rate, and the second is the measurement uncertainty due to the observational noise. Additionally, we include the selection function in terms of the intrinsic count rate ${\widehat{C}}_{\mathrm{R}}$ at the redshift z and the sky location ℋ, i.e. ${\mathcal{I}}_z\left({\widehat{C}}_{\mathrm{R}}\right)\equiv \mathcal{I}\left({\widehat{C}}_{\mathrm{R}}|z,\mathcal{H}\right)$. Therefore, we further deduce the probability P(C_R|M,z,p) as

$$\begin{array}{c}\hfill P\left(C_{\mathrm{R}}|M,z,\mathbf{p}\right)={\displaystyle \int P}\left(C_{\mathrm{R}}|{\widehat{C}}_{\mathrm{R}}\right){\mathcal{I}}_z\left({\widehat{C}}_{\mathrm{R}}\right)P\left({\widehat{C}}_{\mathrm{R}}|M,z,\mathbf{p}\right)\mathrm{d}{\widehat{C}}_{\mathrm{R}},\end{array}$$(3)

where $P\left({\widehat{C}}_{\mathrm{R}}|M,z,\mathbf{p}\right)$ describes the intrinsic count rate ${\widehat{C}}_{\mathrm{R}}$ scattering around the mean value predicted at the fixed cluster mass M and redshift z, and the $P\left(C_{\mathrm{R}}|{\widehat{C}}_{\mathrm{R}}\right)$ accounts for the scatter due to the measurement uncertainty. The count rate-to-mass-and-redshift (C_R–M–z) relation is modelled in Sect. 4.2. We assume a log-normal and constant intrinsic scatter σ_X in the X-ray count rate, as no significant variation of σ_X as a function of M and z is suggested (Chiu et al. 2022, 2023; Ghirardini et al. 2024).

Similarly, the factor P(M_⋆, BCG,g₊,C_R|M,z,p) in the numerator of Eq. (2) describes the probability of observing the observables {M_⋆, BCG,g₊,C_R} at the fixed cluster mass M and redshift z given the parameter vector p, which reads

$$\begin{array}{cc}& P(M_{\star ,\mathrm{BCG}},g_{+},C_{\mathrm{R}}|M,z,\mathbf{p})\hfill \\ \hfill & ={\displaystyle \int \int \int [}P\left(M_{\star ,\mathrm{BCG}}|{\widehat{M}}_{\star ,\mathrm{BCG}}\right)P(g_{+}|M_{\mathrm{WL}},z,\mathbf{p})P\left(C_{\mathrm{R}}|{\widehat{C}}_{\mathrm{R}}\right)\times \hfill \\ \hfill & {\mathcal{I}}_z\left({\widehat{C}}_{\mathrm{R}}\right)P({\widehat{M}}_{\star ,\mathrm{BCG}},M_{\mathrm{WL}},{\widehat{C}}_{\mathrm{R}}|M,z,\mathbf{p})]\mathrm{d}{\widehat{M}}_{\star ,\mathrm{BCG}}\mathrm{d}{M}_{\mathrm{WL}}\mathrm{d}{\widehat{C}}_{\mathrm{R}}.\hfill \end{array}$$(4)

In Eq. (4), the probability $P({\widehat{M}}_{\star ,\mathrm{BCG}},M_{\mathrm{WL}},{\widehat{C}}_{\mathrm{R}}|M,z,\mathbf{p})$ describes the joint intrinsic distribution of the observables $\}{\widehat{M}}_{\star ,\mathrm{BCG}},M_{\mathrm{WL}},{\widehat{C}}_{\mathrm{R}}\{$ at the fixed cluster mass M and redshift z, the factor $P\left(M_{\star ,\mathrm{BCG}}|{\widehat{M}}_{\star ,\mathrm{BCG}}\right)$ takes the measurement uncertainty of the BCG stellar mass into account, and P(g₊|M_WL,z,p) accounts for the measurement uncertainty of the shear profile, expressed by the weak-lensing covariance matrix (Sect. 3.1). The selection function ${\mathcal{I}}_z\left({\widehat{C}}_{\mathrm{R}}\right)$ is also included to calculate Eq. (4).

The correlations between the intrinsic scatters of the observable pairs of the X-ray count rate and the BCG stellar mass (ρ_X, BCG), the X-ray count rate and the weak-lensing mass (ρ_X, WL), and the BCG stellar mass and the weak-lensing mass (ρ_BCG, WL) are included in calculating the joint intrinsic distribution, $P({\widehat{M}}_{\star ,\mathrm{BCG}},M_{\mathrm{WL}},{\widehat{C}}_{\mathrm{R}}|M,z,\mathbf{p})$.

We provide the following remarks for Eq. (4). First, P(g₊|M_WL,z,p) models the observed shear profile g₊(θ) by a mass profile given the weak-lensing mass M_WL and the redshift z. This requires a redshift-distance relation, which has dependence on cosmological parameters, to compare the mass profile in the physical scale with the observed shear profile in the angular scale. Moreover, the mass model includes nuisance parameters to account for the miscentring effect, for example. We therefore keep the parameter vector p in P(g₊|M_WL,z,p) to account for the measurement uncertainty of the shear profile. Second, the weak-lensing mass M_WL, the mass estimate that best describes the observed shear profile, is not equal to the cluster halo mass M that is used to parametrize the halo mass function. This is because the halo model is not a perfect description of observed galaxy clusters; for example, the clusters have substructures or triaxiality that are not included in the model (see Sect. 4.2 in Chiu et al. 2022). This leads to a biased weak-lensing mass M_WL with respect to the cluster halo mass M. Such weak-lensing mass bias is accounted for by the weak-lensing mass-to-halo mass relation given in Sect. 4.3. Third, the measurement uncertainty of the BCG stellar mass is modelled as a Gaussian noise in terms of the ten-based logarithmic BCG stellar mass, i.e. log(M_⋆, BCG).

For the clusters with only the observed shear profile g₊, the modelling of the BCG stellar mass is omitted in Eq. (4), finally leading to

$$\begin{array}{cc}\hfill P(\mathcal{O}|\mathcal{S},\mathbf{p})& \approx P(g_{+}|C_{\mathrm{R}},z,\mathbf{p})\hfill \\ \hfill & =\frac{\int P(g_{+},C_{\mathrm{R}}|M,z,\mathbf{p})n(M|z,\mathbf{p})\mathrm{d}M}{\int P\left(C_{\mathrm{R}}|M,z,\mathbf{p}\right)n(M|z,\mathbf{p})\mathrm{d}M}.\hfill \end{array}$$(5)

For those clusters with only the BCG stellar mass estimates, we have

$$\begin{array}{cc}\hfill P(\mathcal{O}|\mathcal{S},\mathbf{p})& \approx P(M_{\star ,\mathrm{BCG}}|C_{\mathrm{R}},z,\mathbf{p})\hfill \\ \hfill & =\frac{\int P(M_{\star ,\mathrm{BCG}},C_{\mathrm{R}}|M,z,\mathbf{p})n(M|z,\mathbf{p})\mathrm{d}M}{\int P\left(C_{\mathrm{R}}|M,z,\mathbf{p}\right)n(M|z,\mathbf{p})\mathrm{d}M}.\hfill \end{array}$$(6)

4.2. Modelling of the scaling relations

In this work, we use three scaling relations to model the cluster observables as functions of the halo mass and redshift. They are the count rate-to-mass-and-redshift (C_R–M–z) relation, the BCG stellar-mass-to-mass-and-redshift (M_⋆, BCG–M–z) relation, and the weak-lensing mass bias-to-mass-and-redshift (b_WL–M–z or M_WL–M–z) relation. These scaling relations are constrained at the pivotal mass M_piv and redshift z_piv, where we use $M_{\mathrm{piv}}=1.4\times {10}^{14}{h}^{-1}{\mathrm{M}}_{\odot }$ and z_piv = 0.35 throughout this work. The correlated intrinsic scatters among these observables are also included in the modelling (see Sect. 4.1).

4.2.1. The count rate-to-mass-and-redshift relation

For the C_R–M–z relation, we use the same functional form as in the eRASS1 cosmological analysis (Ghirardini et al. 2024), namely

$$\begin{array}{cc}& \langle ln(\frac{{\widehat{C}}_{\mathrm{R}}}{\mathrm{counts}/\sec }|M,z)\rangle =ln\left(A_{\mathrm{X}}\right)\hfill \\ \hfill & +(B_{\mathrm{X}}+{\delta }_{\mathrm{X}}ln\left(\frac{1+z}{1+z_{\mathrm{piv}}}\right))ln\left(\frac{M}{M_{\mathrm{piv}}}\right)+2ln\left(\frac{E\left(z\right)}{E\left(z_{\mathrm{piv}}\right)}\right)\hfill \\ \hfill & -2ln\left(\frac{D_{\mathrm{L}}\left(z\right)}{D_{\mathrm{L}}\left(z_{\mathrm{piv}}\right)}\right)+{\gamma }_{\mathrm{X}}ln\left(\frac{1+z}{1+z_{\mathrm{piv}}}\right),\hfill \end{array}$$(7)

where A_X is the normalization with a unit of counts/sec at the pivotal mass M_piv and redshift z_piv, B_X is the power-law index of the mass trend with a redshift-dependent slope parametrized by δ_X, γ_X describes the power-law scaling of the redshift trend, D_L denotes the luminosity distance, and E(z) accounts for the self-similar redshift evolution of the X-ray luminosity (Kaiser & Silk 1986; Böhringer et al. 2012). We adopt a log-normal and constant intrinsic scatter, ${\sigma }_{\mathrm{X}}\equiv \sqrt{\mathrm{Var}(ln{\widehat{C}}_{\mathrm{R}}|M,z)}$, around the mean value predicted by Eq. (7) at a fixed mass M and redshift z.

4.2.2. The BCG stellar-mass-to-mass-and-redshift relation

For the M_⋆, BCG–M–z relation, or the stellar-mass-to-halo-mass relation, we use the power-law functional form as

$$\begin{array}{cc}& \langle ln(\frac{{\widehat{M}}_{\star ,\mathrm{BCG}}}{{\mathrm{M}}_{\odot }}|M,z)\rangle =A_{\mathrm{BCG}}ln\left(10\right)\hfill \\ \hfill & +(B_{\mathrm{BCG}}+{\delta }_{\mathrm{BCG}}(\frac{\frac{z}{z_{\mathrm{piv}}}}{\frac{1+z}{1+z_{\mathrm{piv}}}}-1))ln\left(\frac{M}{M_{\mathrm{piv}}}\right)+{\gamma }_{\mathrm{BCG}}ln\left(\frac{1+z}{1+z_{\mathrm{piv}}}\right),\hfill \end{array}$$(8)

where A_BCG is the normalization in a unit of dex, B_BCG is the power-law index of the mass trend with a redshift-dependent slope parametrized by δ_BCG, and γ_BCG describes the power-law slope of the redshift scaling. We assume a log-normal and constant intrinsic scatter σ_BCG around the mean value predicted by Eq. (8). Although A_BCG is in units of dex, we express the intrinsic scatter in terms of the natural log, i.e. ${\sigma }_{\mathrm{BCG}}=\sqrt{\mathrm{Var}(ln{\widehat{M}}_{\star ,\mathrm{BCG}}|M,z)}$, in the interest of consistency with the other scaling relations. We note that the parametrization of the redshift-dependent mass slope in Eq. (8), i.e. $(\left(\frac{z}{z_{\mathrm{piv}}}\right)/\left(\frac{1+z}{1+z_{\mathrm{piv}}}\right)-1)$, is a result of following that used in Moster et al. (2013, see their Eq. 14). We have tried a different functional form, ln(1+z)/(1+z_piv), for the redshift-dependent mass slope and find no significant impact on the final results and our conclusions.

4.2.3. The weak-lensing mass bias-to-mass-and-redshift

For the weak-lensing mass bias $b_{\mathrm{WL}}\equiv \frac{M_{\mathrm{WL}}}{M}$, we follow the identical functional form developed in Grandis et al. (2024) and also applied to the cosmological analysis in Ghirardini et al. (2024), as

$$\begin{array}{c}\hfill \langle ln(b_{\mathrm{WL}}|M,z)\rangle =ln{b}_{\mathrm{z}}(z|{\delta }_{\mathrm{WL}},{\gamma }_{\mathrm{WL}})+B_{\mathrm{WL}}ln\left(\frac{M}{M_{\mathrm{piv}}}\right),\end{array}$$(9)

where B_WL describes the mass trend of the weak-lensing mass bias, and lnb_z(z|δ_WL,γ_WL) is the redshift-dependent normalization of the mass bias parametrized by two parameters, δ_X and γ_X, as

$$\begin{array}{cc}& ln{b}_{\mathrm{z}}(z|{\delta }_{\mathrm{WL}},{\gamma }_{\mathrm{WL}})=\hfill \\ \hfill & \mathfrak{I}(z|{\mathbf{z}}_{\mathrm{sim}},{\mu }_0)+{\delta }_{\mathrm{WL}}\mathfrak{I}(z|{\mathbf{z}}_{\mathrm{sim}},{\mu }_1)+{\gamma }_{\mathrm{WL}}\mathfrak{I}(z|{\mathbf{z}}_{\mathrm{sim}},{\mu }_2),\hfill \end{array}$$(10)

in which the notation ℑ(x₀|x,y) stands for the linear interpolation at x₀ among the function that has the values y at x. In Grandis et al. (2024), they calibrated the weak-lensing mass bias using simulated clusters extracted at four snapshots of redshift, which we label them by z_sim in Eq. (10), and derived the mass bias μ₀ with the two principle components μ₁ and μ₂ characterizing the uncertainty. In this way, the distribution of the weak-lensing mass bias as a function of the cluster redshift can be evaluated by directly interpolating among the redshift snapshots z_sim with the parameters of δ_WL and γ_WL modelled by Gaussian distributions of 𝒩(0,1²). While the mean weak-lensing mass bias at the halo mass M and redshift z is calculated with Eqs. (9) and (10), we model the intrinsic scatter of the weak-lensing mass around the mean value by a log-normal and constant scatter as ${\sigma }_{\mathrm{WL}}\equiv \sqrt{\mathrm{Var}\left(b_{\mathrm{WL}}|M,z\right)}$. This modelling is different from those in Grandis et al. (2024) and Kleinebreil et al. (2025), where the intrinsic scatter of the weak-lensing mass was modelled as a function of the halo mass and redshift. This modification is consistent with the modelling of the cluster miscentring, for which we also adopt a fixed model globally for all eRASS1 clusters (see Sect. 4.3 for details). As seen in Sect. 5.1, the results of our modelling with a statistical approach fully agree with those of Kleinebreil et al. (2025), which used a method dependent on individual clusters.

4.3. Weak-lensing modelling

In this subsection, we describe (1) the modelling of the shear profile to obtain the weak-lensing mass M_WL, and (2) the calibration of the weak-lensing mass bias b_WL. For the former, we closely follow the fitting procedure detailed in Chiu et al. (2022) with minor modifications to be consistent with that adopted in Ghirardini et al. (2024). Meanwhile, the latter task has been quantified in Grandis et al. (2024, see also Grandis et al. 2021) and applied to the weak-lensing mass calibration of eRASS1 clusters (e.g. Kleinebreil et al. 2025) as well as the cosmological analyses (Ghirardini et al. 2024). We briefly summarize them, as follows.

Given a cluster at redshift z_cl and a background population with a stacked (and lensing-weight-weighted) redshift distribution P(z), the model $g_{+}^{\mathrm{mod}}$ of the tangential reduced shear profile as a function of the projected (physical) radius R is calculated with a parameter vector p⁵ as

$$\begin{array}{cc}\hfill g_{+}^{\mathrm{mod}}(R|\mathbf{p})& =\frac{{\gamma }_{+}^{\mathrm{mod}}(R|\mathbf{p})}{1-{\kappa }^{\mathrm{mod}}(R|\mathbf{p})}\hfill \\ \hfill & \times (1+{\kappa }^{\mathrm{mod}}(R|\mathbf{p})(\frac{\langle {\beta }^2\left(\mathbf{p}\right)\rangle }{{\langle \beta \left(\mathbf{p}\right)\rangle }^2}-1)),\hfill \end{array}$$(11)

in which β is the lensing efficiency depending on the angular diameter distance⁶,

$$\begin{array}{c}\hfill \beta \left(\mathbf{p}\right)=\{\begin{array}{cc}\frac{D_{\mathrm{A}}\left(z_{\mathrm{cl}},z|\mathbf{p}\right)}{D_{\mathrm{A}}(z|\mathbf{p})}\hfill & \text{if}z>z_{\mathrm{cl}}\hfill \\ 0\hfill & \text{else}\hfill \end{array},\end{array}$$(12)

and ${\gamma }_{+}^{\mathrm{mod}}$ and κ^mod are the models of the tangential shear and convergence, respectively, which read

$$\begin{array}{cc}\hfill {\gamma }_{+}^{\mathrm{mod}}(R|\mathbf{p})& =\frac{\mathrm{\Delta }{\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{mod}}(R|\mathbf{p})}{{\mathrm{\Sigma }}_{\mathrm{crit}}\left(z_{\mathrm{cl}}|\mathbf{p}\right)},\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill {\kappa }^{\mathrm{mod}}(R|\mathbf{p})& =\frac{{\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{mod}}(R|\mathbf{p})}{{\mathrm{\Sigma }}_{\mathrm{crit}}\left(z_{\mathrm{cl}}|\mathbf{p}\right)},\hfill \end{array}$$(14)

where ${\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{mod}}(R|\mathbf{p})$ is the model of the surface mass density at the projected radius R, Σ_crit is the critical surface mass density,

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{crit}}\left(z_{\mathrm{cl}}|\mathbf{p}\right)=\frac{c^2}{4\pi G}\frac{1}{D_{\mathrm{A}}\left(z_{\mathrm{cl}}|\mathbf{p}\right)\langle \beta \left(\mathbf{p}\right)\rangle },\end{array}$$(15)

and $\mathrm{\Delta }{\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{mod}}(R|\mathbf{p})$ is the excess surface mass density,

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{mod}}(R|\mathbf{p})\equiv \left[\frac{2}{R^2}{\displaystyle {\int }_0^R}{\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{mod}}(x|\mathbf{p})x\mathrm{d}x\right]-{\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{mod}}(R|\mathbf{p}).\end{array}$$(16)

In Eq. (11), the quantities of ⟨β(p)⟩ and ⟨β²(p)⟩ are obtained with the weights of the background redshift distribution P(z),

$$\begin{array}{cc}\hfill \langle \beta \left(\mathbf{p}\right)\rangle & ={\displaystyle \int \beta }\left(\mathbf{p}\right)P\left(z\right)\mathrm{d}z\hfill \end{array}$$(17)

$$\begin{array}{cc}\hfill \langle {\beta }^2\left(\mathbf{p}\right)\rangle & ={\displaystyle \int {\beta }^2}\left(\mathbf{p}\right)P\left(z\right)\mathrm{d}z.\hfill \end{array}$$(18)

The projected and physical radius R is determined as R = D_A(z_cl|p) × θ, where the θ is the clustercentric radius in a unit of radian. The comparison between the observed shear profile and the model is in the space of the angular radius θ, and therefore the redshift-distance relation must be used to correctly calculate the physical radius at a fixed θ and a given parameter vector p.

The model of the projected surface mass density ${\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{mod}}$ is composed of two components, the perfectly centered profile ${\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{cen}}$ and that with the miscentring ${\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{mis}}$. They are weighted by the fraction f_mis of the miscentered clusters, in a statistical approach, as

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{mod}}(R|\mathbf{p})=(1-f_{\mathrm{mis}})\times {\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{cen}}(R|\mathbf{p})+f_{\mathrm{mis}}\times {\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{mis}}(R|\mathbf{p}),\end{array}$$(19)

where ${\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{cen}}$ is modelled by a spherical Navarro–Frenk–White (hereafter NFW; Navarro et al. 1997) profile, and ${\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{mis}}$ is modelled by a miscentered NFW model (see e.g. Johnston et al. 2007) with a characteristic and dimensionless miscentring scale σ_mis (see Eq. (20) in Chiu et al. 2022), which is a varied parameter in p.

There are minor modifications from the previous study, Chiu et al. (2022). First, a different halo concentration-to-mass (c–M) relation is used when calculating the NFW profile. Here, we use the mean value predicted by the c–M relation from Ragagnin et al. (2021) at a given weak-lensing mass M_WL and the cluster redshift as the halo concentration, while they used the c–M relation from Diemer (2018). Second, we only include shape noises in the weak-lensing covariance matrix, as the measurement uncertainty of the observed shear profile. The scatter (denoted as ℂ_LSS) in the shear profile due to uncorrelated large-scale structures along the line of sight is accounted for when calibrating the weak-lensing mass bias. To be exact, the scatter ℂ_LSS calculated from the formula in Hoekstra (2003) is included in fitting the weak-lensing shear profile of simulated clusters when calibrating b_WL (see Sect. 2.1.6 in Grandis et al. 2021). Consequently, the resulting scatter of b_WL naturally includes the effect of uncorrelated line-of-sight structures. In Chiu et al. (2022), the scatter ℂ_LSS was calculated and initially included in the weak-lensing covariance matrix as a component of the measurement uncertainty. These two merely differ in the philosophy of the methodology, not the final results. Third, we update the miscentring model to that determined for the eRASS1 sample. Specifically, the miscentring of the eRASS1 clusters was calibrated using the synthetic X-ray eROSITA images of clusters in the eROSITA digital-twin simulations (Grandis et al. 2024, see also simulations produced in Comparat et al. 2019, 2020; Seppi et al. 2022), including both intrinsic and observational effects. As a result, the miscentring of clusters was expressed by the probability f_mis of being miscentered and a characteristic miscentring scale σ_mis that depends on the extent EXT and the detection likelihood L_det, two quantities returned by the eRASS detection pipeline.

The purpose of the first two aforementioned changes is to use the identical c–M relation and the measurement uncertainty as in Grandis et al. (2024), where the dedicated and updated calibration of the weak-lensing mass bias was carried out. Meanwhile, for the third modification we implement the miscentring model in a statistical approach to have a much speedier calculation. That is, we first calculate the distribution of the characteristic miscentring scale σ_mis of our eRASS1 sample based on the miscentring calibration (see Sect. 4.2 in Grandis et al. 2024), model it as a log-normal distribution, and apply the result as an informative Gaussian prior on lnσ_mis in Eq. (19) for all the clusters consistently. Similarly, we also do so for the parameter of the miscentring fraction f_mis for our eRASS1 sample. The resulting Gaussian priors on f_mis and lnσ_mis are 𝒩(0.32,0.043²) and 𝒩(−0.85,0.18²), respectively. This approach is equivalent to statistically correcting for the miscentring effect at a population level. Because we also model the scatter of the weak-lensing mass bias in the same way (see the last paragraph), we do not expect a significant impact on the final results from this modification. Indeed, our result is in excellent agreement with Grandis et al. (2024) and Kleinebreil et al. (2025), as seen in Sect. 5.1.

As fully described in Grandis et al. (2024, see also Grandis et al. 2021), the weak-lensing mass bias, b_WL ≡ M_WL/M, is calibrated by repeating the same modelling on the synthetic shear profiles of clusters extracted from the cosmological hydrodynamical TNG300 simulations (Pillepich et al. 2018), of which the underlying halo mass M used to parametrize the halo mass function is known. The bias and intrinsic scatter of the weak-lensing mass are derived, specifically tailored to include the observational effects of (1) the misfitting due to the halo triaxiality or shape, (2) the miscentring of the eROSITA X-ray centre, (3) the cluster member contamination⁷, (4) the photo-z bias of the background sources, (5) the calibration uncertainty of the galaxy shape measurement. Importantly, the systematic uncertainties of the resulting weak-lensing mass bias and scatter include baryonic effects. This is achieved by quantifying the systematic difference in the results of the b_WL calibrations between the same suite of the dark-matter-only and hydrodynamical simulations. In the end, the constraints on the parameters of Eqs. (9) and (10) are obtained and applied as informative priors in the final modelling. The resulting constraints are tabulated in Table A.3. of Kleinebreil et al. (2025).

In terms of the intrinsic scatter σ_WL of the weak-lensing mass bias, Grandis et al. (2024) parametrized σ_WL as a function of the cluster mass M and redshift z. In this work, we adopt a constant and average weak-lensing mass scatter σ_WL without modelling its mass and redshift dependence. This is a reasonable approach, given that we do not observe significant mass-dependent scatter and that the scatter σ_WL is statistically constant out to redshift z ≈ 0.8 (see Table A.3 in Kleinebreil et al. 2025). In practice, we calculate the distribution of the weak-lensing bias scatter of our eRASS1 sample based on the resulting model reported in Kleinebreil et al. (2025), and model it as a Gaussian distribution that we apply as an informative prior on σ_WL in the final modelling. For our sample, we find that σ_WL can be well described by 𝒩(0.23,0.03²).

4.4. Statistical inference

We adopt the Bayesian framework to explore the posterior P(p|𝒟) of the parameter vector p given the observed data vectors 𝒟,

$$\begin{array}{c}\hfill P(\mathbf{p}|\mathcal{D})\propto \mathcal{L}(\mathcal{D}|\mathbf{p})\mathcal{P}\left(\mathbf{p}\right),\end{array}$$(20)

where ℒ(𝒟|p) is the likelihood of observing 𝒟 at a given p, and 𝒫(p) is the prior on p. In this work, the data vector 𝒟 consists of both the selection and follow-up observables 𝒮 and 𝒪 of individual clusters, i.e. 𝒟 = {𝒮_i,𝒪_i|i∈eRASS1 clusters}. We use the importance nested sampler Multinest (Feroz et al. 2009) implemented in the CosmoSIS framework (Zuntz et al. 2015) to explore the likelihood.

With the technique of the population modelling, i.e. jointly fitting individual clusters at the same time, we calculate the log-likelihood with Eqs. (2), (5), and (6) as

$$\begin{array}{cc}& ln\mathcal{L}\left(\mathcal{D}|\mathbf{p}\right)={\displaystyle \sum _i}lnP({M_{\star ,\mathrm{BCG}}}_i,{g_{+}}_i|{C_{\mathrm{R}}}_i,{z_{\mathrm{cl}}}_i,\mathbf{p})\hfill \\ \hfill & +{\displaystyle \sum _j}lnP({g_{+}}_j|{C_{\mathrm{R}}}_j,{z_{\mathrm{cl}}}_j,\mathbf{p})+{\displaystyle \sum _k}lnP({M_{\star ,\mathrm{BCG}}}_k|{C_{\mathrm{R}}}_k,{z_{\mathrm{cl}}}_k,\mathbf{p}),\hfill \end{array}$$(21)

where z_cl is the cluster redshift (and fixed to its photometric redshift), i runs over the clusters with both the measurements of the weak-lensing shear profile g₊ and the BCG stellar mass M_⋆, BCG, and j and k run over those with only g₊ and M_⋆, BCG, respectively. The parameter vector p includes the parameters⁸ of

• {A_X,B_X,δ_X,γ_X,σ_X} of the C_R–M–z relation in Eq. (7),

• {A_BCG,B_BCG,δ_BCG,γ_BCG,σ_BCG} of the M_⋆, BCG–M–z relation in Eq. (8),

• {B_WL,δ_WL,γ_WL,σ_WL} of the M_WL–M–z relation in Eq. (9),

• {f_mis,lnσ_mis} of the weak-lensing miscentring modelling in Eq. (19), and

• {#x03A9;_m,h,σ₈,Ω_b,n_s} of the cosmological parameters in a flat ΛCDM model.

We adopt the following priors. For the cosmological parameters, we apply flat priors, namely 𝒰(0.10,0.99), 𝒰(0.03,0.07), 𝒰(0.50,0.90), 𝒰(0.94,1.0), and 𝒰(0.45,1.15) to Ω_m, Ω_b, h, n_s, and σ₈, respectively. Additionally, Gaussian informative priors of 𝒩(0.3,0.02²), 𝒩(0.8,0.02²), and 𝒩(0.7,0.05²) are applied to Ω_m, σ₈, and h, respectively. With these informative priors, the cosmological parameters of the flat ΛCDM model are effectively anchored to the widely accepted concordance values with the uncertainties of (Ω_m,σ₈,h) at levels of those obtained from the state-of-the-art constraints.

For the parameters of the weak-lensing mass bias, we apply a Gaussian prior 𝒩(0.047,0.021²) on the mass-dependent power-law index B_WL calibrated in Grandis et al. (2024) (and tabulated in Table A.3 of Kleinebreil et al. 2025). A unit Gaussian prior 𝒩(0,1²) is applied on both γ_WL and δ_WL, by construction. A Gaussian prior 𝒩(0.23,0.03²), informed by the eRASS1 sample we study in this work (see Sect. 4.3), is applied on the intrinsic scatter σ_WL of the weak-lensing mass bias.

We apply Gaussian priors of 𝒩(0.32,0.043²) and 𝒩(−0.85,0.18²) to the parameters of the weak-lensing miscentring model, f_mis and lnσ_mis, respectively (see Sect. 4.3). We required that the determinant of the correlated intrinsic scatter matrix be positive when jointly fitting the observables C_R, g₊, and M_⋆, BCG of a cluster.

For the other parameters, only flat priors were used. We summarize the adopted priors in Table 1. These priors are similar to those used in Chiu et al. (2022) with a major difference: We only applied the flat prior to the intrinsic scatter σ_X of the X-ray count rate in this work, while Chiu et al. (2022, and the cosmological analysis in Chiu et al. 2023) additionally applied a Gaussian prior 𝒩(0.3,0.08²) to σ_X.

5. Results and discussions

In a forward-modelling framework, we constrain both the C_R–M–z and M_⋆, BCG–M–z relations while accounting for the weak-lensing mass bias and scatter via the b_WL–M–z relation. In Sect. 5.1 we first present the weak-lensing calibrated C_R–M–z relation without the modelling of M_⋆, BCG–M–z relation, as an analysis of purely weak-lensing mass calibration. Finally, We present the M_⋆, BCG–M–z relation and discuss its implication in Sect. 5.2.

5.1. The C_R–M–z relation

Regardless of the BCG stellar mass estimates, there are 96 eRASS1 clusters with available weak-lensing shear profiles. Sampling the joint likelihood of those 96 clusters, i.e. Eq. (5), without the modelling of the M_⋆, BCG–M–z relation delivers the constraints on the parameters of the C_R–M–z relation in Eq. (7) as

$$\begin{array}{cc}& A_{\mathrm{X}}=0.{093}_{-0.036}^{+0.029}\hfill \\ \hfill & B_{\mathrm{X}}=1.{50}_{-0.30}^{+0.20}\hfill \\ \hfill & {\delta }_{\mathrm{X}}=0.7\pm 1.5\hfill \\ \hfill & {\gamma }_{\mathrm{X}}=-0.8_{-1.1}^{+1.5}\hfill \\ \hfill & {\sigma }_{\mathrm{X}}=0.{56}_{-0.23}^{+0.17}.\hfill \end{array}$$(22)

These numbers are listed in Table 2. The marginalized posteriors of and the two-dimensional joint posteriors among these parameters are shown as the blue contours in Fig. 3, while we present the full results of the weak-lensing-only modelling in Appendix C. The correlated intrinsic scatter between the X-ray count rate and the weak-lensing mass is constrained as ${\rho }_{\mathrm{WL},\mathrm{X}}=-0.{05}_{-0.63}^{+0.45}$, suggesting that the correlation between the X-ray and WL properties at a fixed halo mass is consistent with zero or ill-constrained in this work.

Fig. 3. Parameter marginalized posteriors (diagonal) and two-dimensional joint posteriors (off-diagonal) from different modelling are presented. The results of the weak-lensing mass calibration alone are shown as the blue contours, while the inclusions of the modelling of the M⋆, BCG–M–z relation are indicated by the red contours. We also indicate the results from Kleinebreil et al. (2025), based on the identical cluster sample and the weak-lensing data g+, as the green contours.

In Fig. 3 we find excellent agreement between this work and Kleinebreil et al. (2025, the green contours). We note that Kleinebreil et al. (2025) performed the weak-lensing mass calibration of the same sample of the 96 eROSITA clusters overlapping the HSC-Y3 footprint based on the identical weak-lensing shear profiles extracted in this work. The main difference is that they modelled (1) the miscentring on a per-cluster basis depending on the X-ray detection likelihood L_det and extent EXT and (2) the redshift-dependent intrinsic scatter of the weak-lensing mass. Meanwhile, we derive the Gaussian priors for the sample and applied them to the parameters at a population level in the modelling (see Sect. 4.3). As a result, these two results obtained from independent codes are fully consistent, again ensuring the robustness of the weak-lensing modelling.

In Fig. 4 we present the stacked shear profiles and best-fit models in four sub-samples defined by the observed count rates and redshift in the left panel. We use the median count rate (C_Rmed ≈ 0.487) and redshift (z_med ≈ 0.27) of the eRASS1 clusters to split the sample into four bins. The stacked shear profile for each sub-sample is calculated as the weighted average of the individual shear profiles g₊, where the weights are defined as the inverse square of the radially dependent shape noise. The best-fit models, with their uncertainties fully marginalized over all parameters, are shown as the grey shaded regions. Similarly, we produce the stacked profile of all the 96 eRASS1 clusters in the right panel. As seen, we find good agreement between the stacked profiles and our best-fit models, suggesting that the weak-lensing modelling provides excellent description of the data.

Fig. 4. Observed shear profiles g+ (black data points) and best-fit models with fully marginalized uncertainties (grey shaded regions). The four subplots in the left panel present the stacked profiles of the subsamples defined with respect to the median values of the observed count rate (CR = 0.487) and the redshift (z = 0.27). The stacked profile of all 96 eRASS1 clusters is presented in the right panel. The numbers of the clusters in the (sub)samples and the reduced χred2 are shown in the lower left corners. In both left and right panels, the radial values (the x-axis) of the measurements are obtained with the inverse-variance weights of the clusters in the (sub)samples.

In Fig. 5 we present the mass and redshift trend of the C_R–M–z relation re-scaled at the pivotal redshift z_piv = 0.35 and mass $M_{\mathrm{piv}}=1.4\times {10}^{14}{h}^{-1}{\mathrm{M}}_{\odot }$ in the left and right panels, respectively. The data points are eRASS1 clusters with their halo mass M randomly sampled from the mass posterior P(M|C_R,g₊,p) with the uncertainties fully marginalized over the posteriors of the parameters p (see Eq. (66) in Chiu et al. 2022). To show the mass trend (left panel), on the y-axis we remove the redshift dependence by dividing the observed count rate C_R by the redshift-dependent factors evaluated at the best-fit parameters; similarly, we remove the mass dependence of C_R to show the redshift trend in the right panel. As seen in the left panel, our results suggest a mass trend with a power-law index of $1.{50}_{-0.30}^{+0.20}$, which is significantly steeper than the self-similar prediction at the soft X-ray band (B_X = 1; the red dashed line). Meanwhile, no significant departure from the self-similar prediction of the redshift trend is revealed in the right panel (${\gamma }_{\mathrm{X}}=-0.8_{-1.1}^{+1.5}$), despite the large error bar. Additionally, we find no clear evidence for the cross-scaling between the halo mass and redshift (δ_X = 0.7 ± 1.5). The mass trend of the X-ray count rate steeper than the self-similar prediction implies that non-gravitational process (e.g. feedback from cluster cores) plays a key role in the X-ray emission. Moreover, this mass scaling this mass scaling neither depends on redshift nor deviates from the self-similar redshift evolution, implying that the non-gravitational mechanisms in massive clusters are already in place at high redshift (z ≈ 0.8 for our cluster sample). This picture is in excellent agreement with previous studies of eROSITA-selected samples (Chiu et al. 2022) and SZE-selected clusters (Bulbul et al. 2019) which probed an even higher redshift to z ≈ 1.3.

Fig. 5. CR–M–z relation of the eRASS1 clusters and the best-fit model in Eq. (22). Its fully marginalized uncertainties are shown as grey shaded regions. The left panel shows the mass trend normalized at the pivotal redshift zpiv = 0.35 after accounting for the redshift-dependent factors, while the right panel similarly presents the redshift scaling normalized at the pivotal mass $M_{\mathrm{piv}}=1.4\times {10}^{14}{h}^{-1}{\mathrm{M}}_{\odot }$. In both panels, the eRASS1 clusters are colour-coded based on the redshift (left panel) and posterior-sampled halo mass (right panel). For comparison, the results from Grandis et al. (2024, green dash-dotted lines) and Kleinebreil et al. (2025, purple lines), which are based on the weak-lensing mass calibration from the DES and KiDS surveys, respectively, are also plotted. The observed CR–M–z relation reveals a mass trend that is steeper than the self-similar prediction (red dashed line), while its redshift trend remains statistically consistent with the self-similar scaling. The normalization of the self-similar predictions (red dashed lines) is fixed to the best-fit AX.

In Fig. 5 the black dotted lines represent the result from Chiu et al. (2022), which were obtained based on ≈300 eROSITA-selected clusters and groups in the eFEDS survey using the same HSC-Y3 WL data. We find satisfactory agreement between Chiu et al. (2022) and this work with the following remarks. First, Chiu et al. (2022) used a slightly different functional form to model the C_R–M–z relation. Specifically, they included a simulation-calibrated correction factor for the observed count rate to probe the underlying true count rate relation in eFEDS (see Sect. 4.1 in Chiu et al. 2022). Consequently, their resulting mass and redshift trends are scaled as ∝M^{B_X − 0.16} and ∝z^0.42(1+z)^γ_X, respectively, where they contained $B_{\mathrm{X}}=1.{58}_{-0.14}^{+0.17}$ and ${\delta }_{\mathrm{X}}=-0.{44}_{-0.85}^{+0.81}$. In this work, we directly relate the observed count rate C_R to the halo mass without the correction factor. In Fig. 5 we plot the overall mass and redshift trends of Chiu et al. (2022) as $\propto M^{B_{\mathrm{X}}-0.16}\stackrel{\propto }{\sim }{M}^{1.42}$ and $\propto z^{0.42}{(1+z)}^{{\gamma }_{\mathrm{X}}}\stackrel{\propto }{\sim }{z}^{0.42}{(1+z)}^{-0.44}$, respectively. Second, as the major difference between Chiu et al. (2022) and the updated analysis of eRASS clusters in this study (also Grandis et al. 2024; Kleinebreil et al. 2025), the former applied a Gaussian prior 𝒩(0.3,0.08²) to the intrinsic scatter σ_X, as opposed to this work in which only an uninformative prior 𝒰(0.05,2.0) is adopted. Removing the Gaussian prior reveals a significantly large intrinsic scatter σ_X ≳ 0.55 as preferred by the count rate observed by eROSITA (see also Grandis et al. 2024; Kleinebreil et al. 2025). Meanwhile, there exists a strong degeneracy between the normalization A_X and σ_X (as seen in Fig. 3), leading to a lower A_X with increasing σ_X in this work. As seen in Fig. 5, we obtain a normalization A_X mildly lower than Chiu et al. (2022, $A_{\mathrm{X}}\approx 0.{148}_{-0.023}^{+0.026}$) at a level of ≈1.6σ.

It is worth mentioning that the large intrinsic scatter σ_X results in significant Eddington (1913) bias, as evident from the offset between the sampled masses P(M|C_R,g₊,p) and the best-fit C_R–M–z relation in Fig. 5. This is explained by the posterior P(M|C_R,g₊,p) evaluated as P(M|C_R,g₊,p) ∝ P(C_R,g₊|M,p)P(M|p), where P(M|p) is the normalized halo mass function. At a fixed halo mass, the large scatter σ_X of the observed count rate leads to a wide distribution P(C_R,g₊|M,p) so that the up-scatter contribution from the low-mass end of the halo mass function P(M|p) becomes significant. Consequently, the sampled masses (data points) are expected to be lower than those following the best-fit C_R–M–z relation (grey region) at a given fixed count rate in the left panel. The Eddington bias is accounted for in this work and other eROSITA analyses, thus the resulting scaling relations are unbiased. An even larger intrinsic scatter of σ_X ≈ 1 was found in Ghirardini et al. (2024), where the authors combined the weak-lensing mass calibration and the abundance of eRASS1 clusters to constrain cosmological parameters. It is surprising that the intrinsic scatter of the eROSITA observed count rate, obtained from either weak lensing alone or a combination with the cluster abundance, is significantly larger than that (≈0.3) of X-ray luminosity in previous studies (Pratt et al. 2009; Vikhlinin et al. 2009b; Lovisari et al. 2015; Mantz et al. 2016; Bulbul et al. 2019; Chiu et al. 2022; Bahar et al. 2022; Akino et al. 2022). This implies that significantly large scatter is introduced in measuring the X-ray count rate in the eROSITA data pipeline, increasing the intrinsic scatter at a fixed halo mass from a scale of 30% in the X-ray luminosity to a factor of exp(1) ≈ 2.7 in the count rate. A future study investigating this topic is warranted.

We show the comparisons with the weak-lensing-calibrated results of eRASS1 clusters using data from the DES (Grandis et al. 2024, the green lines) and KiDS (Kleinebreil et al. 2025, the purple lines) surveys in Fig. 5. We observe excellent agreement with the DES result as they constrained A_X = 0.088 ± 0.02, B_X = 1.62 ± 0.14, δ_X = −0.85 ± 0.93, γ_X = −0.32 ± 0.69, and σ_X = 0.61 ± 0.19. Meanwhile, mild discrepancy at ≲1.5σ is seen between this work the KiDS result with their error bars generally larger than ours, as they obtained A_X = 0.17 ± 0.042, B_X = 1.68 ± 0.27, δ_X = 1.57 ± 1.74, γ_X = −2.60 ± 1.45, and only the 1σ upper limit on the intrinsic scatter as σ_X < 0.85.

Fig. 6. As in Fig. 3, but for the marginalized posteriors and two-dimensional joint posteriors of the CR–M–z and M⋆, BCG–M–z relation.

It is worth noting that Kleinebreil et al. (2025) performed a direct comparison of the weak-lensing measurements between the KiDS and HSC surveys using 48 common clusters. They found excellent agreement in the observed shear profile g₊ but a significant discrepancy in the inferred excess surface mass density ΔΣ_m. Specifically, ΔΣ_m inferred from HSC is statistically higher than that from KiDS by ≈30% at a level of ≈4σ. Meanwhile, only marginal difference (at a level of ≈2.3σ) was seen in g₊ and ΔΣ_m between DES and KiDS. Because the source redshift is required to derive ΔΣ_m from the observed g₊, the agreement in g₊ but not in ΔΣ_m between the KiDS and HSC surveys implies a potential systematic uncertainty existing in the source photo-z calibration. Recall that g₊ ≈ γ₊ = ΔΣ_m/Σ_crit ∝ ΔΣ_m × β, where β is the lensing efficiency, which monotonically increases with source redshift. The higher ΔΣ_m in the HSC survey compared to KiDS suggests that the photo-z of the source sample is underestimated in HSC and/or overestimated in KiDS. However, the photo-z in the HSC survey is rigorously calibrated against spectroscopic samples and shows no bias at redshift z ≲ 1.2 (Rau et al. 2023). In addition, the agreement between the DES and KiDS surveys suggests no significant bias in the photo-z calibration of their source samples, which are primarily at relatively low redshift (z ≲ 1). A possible cause is that the photo-z of the HSC source sample at high redshift z ≳ 1.2 is biased low. This aligns well with the finding of the HSC-Y3 cosmic-shear analyses (Li et al. 2023; Dalal et al. 2023; Sugiyama et al. 2023), which suggest that the photo-z of the HSC sources at z ≳ 1.2 is biased low by ≈0.2 at a level of ≈2σ. However, the photo-z bias at z ≳ 1.2 at a level of ≈0.2 alone cannot explain the ≈30% difference in Σ_m for the eRASS1 clusters, which are primarily at redshift z ≈ 0.35. A detailed investigation of this systematic uncertainty is beyond the scope of this paper and is clearly warranty in future work, with a focus on the photo-z calibration at high redshift z ≳ 1.2.

We provide the mean of the mass posterior, $\langle log\left(\frac{M_{500\mathrm{c}}}{h^{-1}{\mathrm{M}}_{\odot }}\right)\rangle$, and its sampled mass in Table B.1. The mass posterior is calculated as P(M|C_R,𝒪,p), where 𝒪 = {g₊,M_⋆, BCG}, {M_⋆, BCG}, or {g₊}, depending on the data availability.

5.2. The M_⋆, BCG–M–z relation

To begin with, we verify that further including the modelling of the BCG stellar mass relation, as in Eq. (21), yields a consistent C_R–M–z relation. This is clearly evident from the excellent agreement between the blue (WL only) and red (WL+BCG) contours in Fig. 3, suggesting that the cluster mass estimates, informed by the X-ray count rate with an absolute weak-lensing calibration, are robust to constrain the BCG stellar mass relation.

As listed in Table 2, we obtain the constraints on the M_⋆, BCG–M–z relation, i.e. Eq. (8), as

$$\begin{array}{cc}\hfill A_{\mathrm{BCG}}& =11.{547}_{-0.074}^{+0.083}\hfill \\ \hfill [1.5pt]B_{\mathrm{BCG}}& =0.{184}_{-0.16}^{+0.068}\hfill \\ \hfill [1.5pt]{\delta }_{\mathrm{BCG}}& =0.22\pm 0.29\hfill \\ \hfill {\gamma }_{\mathrm{BCG}}& =1.{49}_{-0.75}^{+0.88}\hfill \\ \hfill [1.5pt]{\sigma }_{\mathrm{BCG}}& =0.{554}_{-0.049}^{+0.040}.\hfill \end{array}$$(23)

The parameter marginalized posteriors and two-dimensional joint posteriors are presented in Fig. 6 (see the full results in Appendix C). Meanwhile, we constrain the correlated intrinsic scatter as ρ_WL, X = −0.04 ± 0.43, ${\rho }_{\mathrm{WL},M_{\star ,\mathrm{BCG}}}=-0.{32}_{-0.49}^{+0.23}$, and ${\rho }_{\mathrm{X},M_{\star ,\mathrm{BCG}}}=-0.{10}_{-0.27}^{+0.24}$, suggesting no significant intrinsic correlation among the BCG stellar mass, X-ray count rate, and the weak-lensing mass at a fixed halo mass.

It is noticeable that there exists large intrinsic scatter in the BCG stellar mass M_⋆, BCG at a fixed halo mass, constrained as ${\sigma }_{\mathrm{BCG}}=0.{554}_{-0.049}^{+0.040}$ in this work. This corresponds to an intrinsic scatter of ≈0.24 ± 0.02 dex in a unit of log(M_⋆, BCG/M_⊙), which is consistent with previous studies (e.g. Leauthaud et al. 2012; Kravtsov et al. 2018; Erfanianfar et al. 2019; Akino et al. 2022). Unlike the C_R–M–z relation, the M_⋆, BCG–M–z relation has a shallow mass slope ($B_{\mathrm{BCG}}=0.{184}_{-0.16}^{+0.068}$), so that the scatter $\sqrt{\mathrm{Var}\left(lnM|M_{\star ,\mathrm{BCG}}\right)}$ in the logarithmic of the halo mass at a fixed BCG stellar mass is at a level of ≈σ_BCG/B_BCG ≈ 3, corresponding to nearly a factor of exp(3) ≈ 20 in the halo mass M. That is, it is not feasible to adopt the BCG stellar mass as a precise mass proxy on a per-cluster basis.

Interestingly, Golden-Marx et al. (2022) found that the intrinsic scatter σ_BCG of the BCG stellar mass at a fixed halo mass could be largely reduced by accounting for their local environments, as indicated by the magnitude gap between the BCG and the fourth brightest cluster galaxy (see also Golden-Marx & Miller 2018, 2019). Moreover, Golden-Marx et al. (2025) measured a non-zero correlation between the magnitude gap and the central stellar mass (BCG+ICL), supporting a scenario in which BCGs primarily assemble their masses at late times through ex situ processes, such as merging. A similar correlation between the central stellar mass and the halo concentration was also found in Zu et al. (2022). Therefore, it is possible to adopt the BCG stellar mass M_⋆, BCG as a reliable mass proxy for eROSITA clusters if the formation history of BCGs is properly taken into account (e.g. Huang et al. 2022).

In Fig. 7 we present the mass and redshift trends of the M_⋆, BCG–M–z relation in the left and right panels, respectively, with the re-normalization as similarly done in Fig. 5. Additionally, we plot the running inverse-variance-weighted means and the standard deviations of the measurements as the pink circles, using five sub-sample bins of equal cluster numbers. As seen, the best-fit model provides a good description of the data with the noticeably large intrinsic scatter at a fixed halo mass. Our results suggest no or mildly increasing redshift trend of the BCG stellar mass at a fixed halo mass (${\gamma }_{\mathrm{BCG}}=1.{49}_{-0.75}^{+0.88}$) at a level of ≈2σ. However, it is in contrast to the physical picture of the hierarchical structure formation, where galaxies at late time are expected to be more massive that those at earlier time, on average (see e.g. Moster et al. 2013). This is primarily caused by the strong degeneracy between the parameters B_BCG and γ_BCG (see Fig. 6), which arises from a combination of the large intrinsic scatter σ_BCG and the X-ray selection of the cluster sample. While the X-ray selection favors high-mass clusters at high redshift, the excessively large scatter σ_BCG (corresponding to an ≈80% variation in M_⋆, BCG at a fixed halo mass M) prevents us from distinguishing between two scenarios: (1) BCGs at a fixed halo mass are intrinsically increasing with redshift or (2) more massive BCGs at high redshift appear in the sample purely due to the selection bias (i.e. the Malmquist bias; Malmquist 1922). In the absence of low-mass clusters at high redshift, it is challenging to simultaneously constrain both the mass and redshift trends given the size of the current sample.

Fig. 7. As in Fig. 5, but for the M⋆, BCG–M–z relation of the eRASS1 clusters and the best-fit model in Eq. (23) with its fully marginalized uncertainties (grey shaded regions). We also plot the running inverse-variance-weighted means and standard deviations as the pink open circles, using five bins of an equal number of clusters.

In light of the strong degeneracy between B_BCG and γ_BCG, we additionally constrain the M_⋆, BCG–M–z relation with an informative prior applied to the parameter γ_BCG. The informative prior is derived as follows. Based on Moster et al. (2013), where the authors used an abundance-match method and accounted for the Eddington bias in the halo mass estimates, they constrained the M_⋆, BCG–M–z relation across wide ranges of halo mass and redshift (see their Table 1). We take an eRASS1 cluster with the mean of its mass posterior⁹ $\overline{M}$, place it at a series of equal-step redshift bins z_j between 0 and 1, and calculate its corresponding BCG stellar mass ${\langle M_{\star ,\mathrm{BCG}}|\overline{M},z_j\rangle }_{\mathrm{Moster}+13}$ (without the intrinsic scatter and measurement uncertainty) following the Moster et al. (2013) relation assuming that it is at each redshift z_j. We then re-fit the functional form $A{\left(\frac{1+z}{1+z_{\mathrm{piv}}}\right)}^{\gamma }$, which is the redshift scaling of the M_⋆, BCG–M–z relation, to the BCG stellar masses ${\langle M_{\star ,\mathrm{BCG}}|\overline{M},z_j\rangle }_{\mathrm{Moster}+13}$ at different redshifts z_j using a χ² minimization. This gives the best-fit power-law index γ of the redshift scaling for this eRASS1 cluster with the mass $\overline{M}$. We repeat the aforementioned procedure for all eRASS1 clusters, resulting a distribution of the best-fit γ which represents the distribution of the power-law indices for our eRASS1 sample at their mass scales. The distribution is obtained as a nearly Gaussian 𝒩(−0.4,0.14²). This result is free from the Malmquist bias because we do not include the intrinsic scatter and measurement uncertainty in calculating ${\langle M_{\star ,\mathrm{BCG}}|\overline{M},z_j\rangle }_{\mathrm{Moster}+13}$. In the end, we apply a Gaussian prior 𝒩(−0.4,0.14²) to the parameter γ_BCG and obtain an improved constraint on the M_⋆, BCG–M–z relation as

$$\begin{array}{c}\hfill \begin{array}{cc}{A}_{\mathrm{BCG}}& =11.{433}_{-0.055}^{+0.066}\\ B_{\mathrm{BCG}}& =0.38\pm 0.11\\ {\delta }_{\mathrm{BCG}}& =0.62\pm 0.25\\ {\gamma }_{\mathrm{BCG}}& =-0.36\pm 0.13\\ {\sigma }_{\mathrm{BCG}}& =0.{560}_{-0.055}^{+0.044}\end{array}.\end{array}$$(24)

These results are also listed in Table 2. The mass and redshift trends of the M_⋆, BCG–M–z relation with the Gaussian prior on γ_BCG are presented in Fig. C.2.

In Fig. 8 we present the parameter marginalized posteriors and two-dimensional joint posteriors in comparison with those without the Gaussian prior on γ_BCG. As seen, the constraints with (blue contours) and without (red contours) the Gaussian prior are statistically consistent with each other at a level of ≈1σ. However, the degeneracy between B_BCG and γ_BCG is broken by including the Gaussian prior, leading to a ≈1σ steeper mass trend (B_BCG = 0.38 ± 0.11) with smaller error bars. The steeper mass trend (B_BCG = 0.38 ± 0.11) is also in better agreement with previous studies (e.g. Gonzalez et al. 2013; van der Burg et al. 2014; Chiu et al. 2016; Kravtsov et al. 2018; Akino et al. 2022), suggesting that our constraints based on an X-ray-selected sample could be prone to the B_BCG–γ_BCG degeneracy. Interestingly, once the B_BCG–γ_BCG degeneracy is broken, there exists a cross-scaling between the mass and redshift (δ_BCG = 0.62 ± 0.25) at a level of ≈2.5σ. This parameter was constrained as −0.085 ± 0.045 in Moster et al. (2013)¹⁰, which has an opposite trend that is inconsistent with ours at a mildly significant level (≈2.8σ). A larger sample with more low-mass clusters at high redshift will shed light on this. The interpretations of the other parameters A_BCG and σ_BCG remain unaffected by the Gaussian prior.

Fig. 8. Comparison between the M⋆, BCG–M–z relations obtained with (blue contours) and without (red contours) the Gaussian prior applied on γBCG. The Gaussian prior is obtained based on the overall mass and redshift distributions of the eRASS1 clusters studied in this work, following the assumed relation from Moster et al. (2013) (see the text for details).

In what follows, we discuss potential systematic uncertainties arising from the SED fitting. In our fiducial analysis, we make an assumption that the BCGs are dominated by passively evolving galaxy populations, so we only fit the SPS model with the star-formation timescale τ up to 5 Gyr and low dust extinction E(B−V) = 0, 0.1, 0.2, 0.3. We examine the systematics regarding this assumption. Namely, we allow a new SPS model with τ up to 30 Gyr configured with the dust-extinction coefficient of E(B−V) = 0, 0.1, 0.2, 0.3, 0.4, 0.5 and repeat the whole analysis. We find the constraints of A_BCG = 11.498 ± 0.085, $B_{\mathrm{BCG}}=0.{192}_{-0.18}^{+0.057}$, ${\delta }_{\mathrm{BCG}}=0.{37}_{-0.32}^{+0.29}$, ${\gamma }_{\mathrm{BCG}}=0.{91}_{-0.77}^{+0.97}$, and ${\sigma }_{\mathrm{BCG}}=0.{566}_{-0.051}^{+0.043}$, which are in excellent agreement with the fiducial results. This reinforces the picture that the BCGs are dominated by passively evolving populations. On the other hand, we also examine the systematics due to the inclusion of the mid-infrared photometry. Specifically, we perform the SED fitting with the broadband photometry at 3.4 μm and 4.6 μm from WISE, and repeat the subsequent analysis. We obtain a result that is in agreement with our fiducial results at a level of ≲1σ, further ensuring the robustness of our SED fitting. We provide a more detailed description about the SED fitting including the WISE photometry in Appendix D.

We now turn into the comparisons of the M_⋆, BCG–M–z relation with previous studies. These comparisons are illustrated in Figs. 9 and 10, which represent the trends in mass and redshift, respectively, at the pivotal values of redshift (z_piv = 0.35) and mass ($M_{\mathrm{piv}}=1.4\times {10}^{14}{h}^{-1}{\mathrm{M}}_{\odot }$). For the re-normalization, we adopt the best-fit parameters in Eq. (24), which are obtained with the Gaussian prior on γ_BCG. We note that the best-fit parameters in Eq. (24) incorporates the external information of the Gaussian prior, which is derived based on the overall mass and redshift distributions of eRASS1 clusters informed by the previous work of Moster et al. (2013). As such, our focus in this subsection is on the consistency between the eRASS1 measurements and the literature.

Fig. 9. Comparisons of the mass trend of the BCG stellar mass between the eRASS1 clusters and previous observational (left panel) and simulation-based (right panel) studies. In both panels, the mass trend is re-normalized at the pivotal redshift (zpiv = 0.35) using the best-fit M⋆, BCG–M–z relation as in Eq. (24). The eRASS1 clusters are indicated by black stars, while previous studies are shown according to the lists to the right of each panel.

Fig. 10. As in Fig. 9, but for the comparisons of the redshift trend of the BCG stellar mass between the eRASS1 clusters and previous observational (left panel) and simulation-based (right panel) studies, and are re-normalized at the pivotal mass ($M_{\mathrm{piv}}=1.4\times {10}^{14}{h}^{-1}{\mathrm{M}}_{\odot }$).

In Figs. 9 and 10, the left panels compare the results with previous observational studies, while the right panels compare them with simulation-based results. The eRASS1 cluster measurements from this work, with posterior-sampled halo masses, are indicated by black stars.

For previous observational results, we include the following:

• Lidman et al. (2012, grey points) compiled a heterogeneous sample of clusters spanning at 0.1 ≲ z ≲ 1.6 with the halo mass estimated from either X-ray mass proxies or the line-of-sight velocity dispersion, and derived the stellar mass from the K-band luminosity.

• Gonzalez et al. (2013, purple triangles) studied 12 X-ray-selected clusters at local Universe (z ≈ 0.1) and estimated their BCG stellar masses including the ICL.

• van der Burg et al. (2014, red open circles) measured the BCG stellar mass using 11-band SED fitting and the halo mass from the line-of-sight velocity dispersion for a sample of high-redshift clusters (0.9 ≲ z ≲ 1.4) selected in the near-infrared.

• Chiu et al. (2016, green crosses) studied 14 SZE-selected massive clusters at redshift 0.6 ≲ z ≲ 1.3 and derived the BCG stellar mass using SED fitting to data with wavelength from the optical to mid-infrared.

• Kravtsov et al. (2018, brown pluses) derived the BCG stellar mass (using r-band luminosity) and the halo mass (using X-ray mass proxies) of 9 nearby clusters at z ≲ 0.1.

• DeMaio et al. (2018, pink pluses) studied the stellar mass of BCGs and ICL for 23 clusters at redshift 0.3 ≲ z ≲ 0.9 using the imaging from Hubble Space Telescope.

• Erfanianfar et al. (2019, yellow open circles) studied the M_⋆, BCG–M–z relation of 438 X-ray-selected clusters out to z ≈ 0.65 based on the BCG stellar mass estimated from SED fitting to SDSS and WISE photometry.

• Akino et al. (2022, red dash-dotted lines) measured the BCG stellar mass for 136 X-ray-selected clusters drawn from the XXL survey (Adami et al. 2018) using the 5-band grizY photometry from HSC.

As suggested by Chen et al. (2022, see also their detailed profile of the central stellar mass), we use the stellar mass measurements within < 50 kpc for both Kravtsov et al. (2018) and DeMaio et al. (2018) to roughly separate the BCG stellar mass M_⋆, BCG from the ICL. For these comparison samples, we consistently correct for the systematics in the halo and stellar mass estimates following the prescription in Chiu et al. (2018). Specifically, the halo mass M is multiplied by 0.96 (1.12) if it is estimated from the line-of-sight velocity dispersion (X-ray mass proxies), while the stellar mass is all corrected to that inferred by the Chabrier (2003) IMF.

In the comparisons, we consider the following previous simulation-based results, including those based on the abundance-matching technique which assigns mass M from N-body simulations to observed halos statistically:

• Both Moster et al. (2013, blue shaded regions) and Behroozi et al. (2013, red shaded regions) determined the stellar-mass-to-halo-mass relation out to high redshift z ≳ 4 using the observed stellar mass functions in the SDSS field (e.g. Baldry et al. 2008; Li & White 2009; Bouwens et al. 2012; Moustakas et al. 2013), each with a slightly different functional form.

• Girelli et al. (2020, the cyan dotted lines) studied the stellar-mass-to-halo-mass relation in the COSMOS field using a modified functional form from Moster et al. (2013) with a focus on its redshift evolution.

All the above-mentioned studies with the abundance-matching technique assumed a fixed intrinsic scatter of the central stellar mass given a halo mass at a level ranging from ≈0.15 dex to ≈0.20 dex. The 1σ confidence regions of Moster et al. (2013) and Behroozi et al. (2013) in Figs. 9 and 10 are obtained from a full marginalization of all parameters. Meanwhile, we directly take the data points from cosmological hydrodynamic simulations (at redshift z ≲ 0.1) of

• the IllustrisTNG simulation (Pillepich et al. 2018, green shaded regions), and

• the Three Hundred Project (Cui et al. 2022, pink points with error bars).

No other corrections are applied to the halo mass M and the BCG stellar mass M_⋆, BCG in the IllustrisTNG and Three Hundred projects in Figs. 9 and 10. We use the BCG stellar mass within 30 kpc for IllustrisTNG (Pillepich et al. 2018) due to the lack of that within 50 kpc.

As seen in the figures, the eRASS1 clusters (black stars) are in satisfactory agreement with observational results, which exhibit large scatter in M_⋆, BCG at a fixed halo mass (constrained as σ_BCG ≈ 0.6 or ≈0.25 dex in this work). After removing the redshift dependence (γ_BCG ≈ −0.4) in the left panel of Fig. 9, the BCG stellar mass of the previous studies reveal a mass slope that is consistent with that of the eRASS1 sample. It is noteworthy that the eRASS1 clusters span the widest halo mass range from $\approx {10}^{13}{\mathrm{M}}_{\odot }$ to $\approx {10}^{15}{\mathrm{M}}_{\odot }$, again highlighting the power of eROSITA in finding galaxy clusters and groups. Meanwhile, the agreement in the mass trend between the eRASS1 clusters and the simulation-based results is seen at a similar level (the right panel in Fig. 9). Nevertheless, we observe systematic offset in the absolute scale of the BCG stellar mass. Specifically, the BCG stellar masses from both the hydrodynamic simulations of the IllustrisTNG and Three Hundred projects are systematically higher (by a factor of ≈2 to ≈3) than eRASS1 clusters and other observational results, except for Gonzalez et al. (2013), Kravtsov et al. (2018), and DeMaio et al. (2018). On the other hand, the results from the abundance-matching method (Moster et al. 2013; Behroozi et al. 2013; Cui et al. 2022), which statistically assigns the mass M from N-body simulations to observed halos with the stellar mass estimates, are in better agreement with the eRASS1 clusters and other observational studies. This systematic difference implies possible systematics in the stellar mass produced in the hydrodynamic simulations and/or the photometric measurements in the survey pipelines (see Huang et al. 2018a; Akino et al. 2022, for details in the systematics of the cmodel photometry).

After accounting for the halo-mass dependence (B_BCG ≈ 0.4) in Fig. 10, the eRASS1 measurements similarly show no significant redshift trend, in good agreement with previous studies. This finding strongly suggests that the BCG stellar mass given a fixed halo mass (i.e. the ratio M_⋆, BCG/M) remains relatively stable at the high-mass end, with no clear evidence of evolution out to redshift z ≈ 1. The same picture of the little to no increase in the BCG stellar mass at the pivotal mass ($M_{\mathrm{piv}}=1.4\times {10}^{14}{h}^{-1}{\mathrm{M}}_{\odot }$) is also suggested by the abundance-matching results of Moster et al. (2013) and Behroozi et al. (2013). The increase of the BCG stellar mass at the pivotal halo mass from z = 0.8 to z = 0.1 is estimated to be at a level of ${\left(\frac{1+0.1}{1+0.8}\right)}^{-0.36\pm 0.13}-1\approx (20\pm 8)\%$, informed by the Gaussian prior from Moster et al. (2013). However, Girelli et al. (2020) reported a positive redshift trend, deviating from this pattern. Our finding is consistent with Zhang et al. (2016), where the authors obtained a mildly increasing BCG stellar mass at a fixed halo mass for X-ray-selected clusters with a scaling of ≈(1+z)^{−0.19 ± 0.34}, corresponding to ≈0.13 dex or ≈34% from z ≈ 1 to z ≈ 0. A similar picture is also suggested by Lin et al. (2017, ≈35% increase from z ≈ 1 to z ≈ 0.3).

We note that our results quantify the scaling of the BCG stellar mass statistically as a function of the halo mass and redshift, which is complementary to the method that traces the growth of the BCG stellar mass by paring the progenitor and descendant of a cluster (see Lin et al. 2025). With a sufficiently large sample of clusters covering wide ranges of the halo mass and redshift, our approach with a flexible scaling form is expected to be equivalent to the progenitor-descendant pairing strategy. However, this is not feasible for the eRASS1 sample due to the nature of the X-ray selection, which lacks low-mass groups at high redshift that could serve as the progenitors of high-mass clusters at low redshift.

Using the merging tree in numerical simulations, Lin et al. (2013) found a significant growth of a factor ≈2.3 in the BCG stellar mass from z ≈ 1.5 to z ≈ 0.5 but no growth between z ≈ 0.5 and z ≈ 0. Similarly, Lidman et al. (2012) found a significant growth in the BCG stellar mass (after accounting the halo-mass dependence) of ≈(80±30)% from z ≈ 0.9 to z ≈ 0.1 but no growth between z ≈ 0.4 and z ≈ 0.1. This suggests that BCGs must have assembled the majority of their stellar masses at much earlier epochs (z ≳ 1), before massive clusters formed, and have accreted only small amounts at later times. Moreover, BCGs must grow primarily through an ex-situ process (e.g. merging) at late epochs, as no evidence of recent star formation is found. Consequently, the BCG stellar mass fraction (M_⋆, BCG/M) exhibits a negative mass trend, i.e. B_BCG − 1 ≈ −0.6, while the host cluster continues to grow significantly in the total mass M during the formation. This “rapid-then-slow” growth of the BCG stellar mass is also suggested by the semi-analytical model in De Lucia & Blaizot (2007) and further refined with an updated ICL production in Contini et al. (2014). This two-phase hierarchical formation describes that BCGs acquired their stellar masses in the core via in-situ star formation at very early epochs (z ≳ 2) and gradually assembled additional mass at the large radii via ex-situ mergers (De Lucia & Blaizot 2007; Guo et al. 2011). On the other hand, we find that BCGs at a fixed halo mass ($M_{\mathrm{piv}}=1.4\times {10}^{14}{h}^{-1}{\mathrm{M}}_{\odot }$) assemble little stellar mass since redshift z ≈ 0.8, which is in broad agreement with this picture of the “rapid-then-slow” mass assembly. This picture is also in line with the recent observational result in Lin et al. (2025).

6. Conclusions

Using the data from the HSC survey, we presented the weak-lensing mass calibration of 124 X-ray-selected galaxy clusters selected in the first scan (eRASS1) of the eROSITA All-Sky Survey (eRASS; Merloni et al. 2024) and constrained their BCG stellar-mass-to-halo-mass-and-redshift (M_⋆, BCG–M–z) relation or the stellar-mass-to-halo-mass relation. The cluster sample spans a mass range of $3\times {10}^{13}{h}^{-1}{\mathrm{M}}_{\odot }\lesssim M_{500\mathrm{c}}\lesssim {10}^{15}{h}^{-1}{\mathrm{M}}_{\odot }$ and redshift range of 0.1 < z < 0.8, with the majority of low-mass systems ($M_{500\mathrm{c}}\lesssim 1.5\times {10}^{14}{h}^{-1}{\mathrm{M}}_{\odot }$) at redshift z ≲ 0.35.

The eRASS1 clusters were X-ray-selected (Bulbul et al. 2024) and optically confirmed (Kluge et al. 2024), leading to a clean ICM-based sample with a point-source contamination at a level of ≈5%. We used the observed X-ray count rate C_R as the proxy for the halo mass (M_500c or M), based on which the selection function was parametrized and evaluated at the sky location and redshift of each eRASS1 cluster (Clerc et al. 2024).

For the weak-lensing mass calibration, we measured the tangential reduced shear profile g₊ of each cluster using the latest HSC-Y3 dataset (Li et al. 2022). Following the previous studies (Chiu et al. 2022, 2023), we derived the weak-lensing mass M_WL and inferred the underlying halo mass M using the weak-lensing-mass-to-mass-and-redshift (M_WL–M–z) relation. The M_WL–M–z relation was calibrated by hydrodynamic simulations, accounting for the modelling systematics including the halo triaxiality, the miscentring of X-ray-defined centers, the cluster member contamination to the weak-lensing observable, the bias in the photometric redshift of source galaxies, the calibration uncertainty of the galaxy shape measurement, and the baryonic feedback.

For the BCG stellar mass M_⋆, BCG, we employed the SED fitting technique on the grizY five-band photometry from the third Public Data Release (PDR3) of the HSC survey. The SED template was built from the model of Stellar Population Synthesis (SPS) assuming a passively evolving population, which is supported by our BCG data. We additionally examined the systematics arising from the inclusion of the mid-infrared data (from WISE) in the SED fitting and found no statistically significant change to the final results. We conclude that our BCG stellar mass estimates are homogeneous and robust. We present the HSC cutout images of these BCGs and their SED fitting results in Appendix B.

Among the 124 eRASS1 clusters in this study, 73 clusters have measurements for both the shear profile g₊ and the BCG stellar mass M_⋆, BCG, while 23 have only g₊ measurements and 28 have only M_⋆, BCG measurements. In an approach of forward and population modelling, we simultaneously constrained both the count rate-to-mass-and-redshift (C_R–M–z) relation and the BCG stellar-mass-to-halo-mass (M_⋆, BCG–M–z) relation, accounting for the bias and intrinsic scatter of the weak-lensing mass via the M_WL–M–z relation. We also included the correlated intrinsic scatter among the X-ray count rate C_R, the BCG stellar mass M_⋆, BCG, and the weak-lensing mass M_WL, and found no statistically significant correlation.

Using the weak-lensing shear profiles g₊ alone, we obtain a constraint on the C_R–M–z relation as

$$\begin{array}{cc}& \langle ln(\frac{{\widehat{C}}_{\mathrm{R}}}{\mathrm{counts}/\sec }|M,z)\rangle =ln(0.{093}_{-0.036}^{+0.029})\hfill \\ \hfill & +((1.{50}_{-0.30}^{+0.20})+(0.7\pm 1.5)ln\left(\frac{1+z}{1+z_{\mathrm{piv}}}\right))ln\left(\frac{M}{M_{\mathrm{piv}}}\right)\hfill \\ \hfill & +2ln\left(\frac{E\left(z\right)}{E\left(z_{\mathrm{piv}}\right)}\right)-2ln\left(\frac{D_{\mathrm{L}}\left(z\right)}{D_{\mathrm{L}}\left(z_{\mathrm{piv}}\right)}\right)\hfill \\ \hfill & +(-0.8_{-1.1}^{+1.5})ln\left(\frac{1+z}{1+z_{\mathrm{piv}}}\right),\hfill \end{array}$$

with the intrinsic scatter ${\sigma }_{\mathrm{X}}=0.{56}_{-0.23}^{+0.17}$. The resulting C_R–M–z relation suggests a mass trend significantly steeper than the self-similar prediction of the X-ray luminosity-to-mass relation, but reveals no deviation from the self-similar scaling in redshift. As also found in the companion eROSITA studies (Grandis et al. 2024; Kleinebreil et al. 2025; Ghirardini et al. 2024), we confirm the excessively large intrinsic scatter σ_X compared to the typical value (≈0.3) obtained for the luminosity-to-mass relation in the literature. This implies that significant scatter is introduced in measuring the eROSITA observed count rate converted from the physical flux in X-rays. Our results are in good agreement with those calibrated with the weak-lensing data from the DES (Grandis et al. 2024) and KiDS (Kleinebreil et al. 2025) surveys.

The resulting M_⋆, BCG–M–z relation is constrained as

$$\begin{array}{cc}& \langle ln(\frac{{\widehat{M}}_{\star ,\mathrm{BCG}}}{{\mathrm{M}}_{\odot }}|M,z)\rangle =(11.{547}_{-0.074}^{+0.083})ln\left(10\right)\hfill \\ \hfill & +((0.{184}_{-0.16}^{+0.068})+(0.22\pm 0.29)(\frac{\frac{z}{z_{\mathrm{piv}}}}{\frac{1+z}{1+z_{\mathrm{piv}}}}-1))ln\left(\frac{M}{M_{\mathrm{piv}}}\right)\hfill \\ \hfill & +(1.{49}_{-0.75}^{+0.88})ln\left(\frac{1+z}{1+z_{\mathrm{piv}}}\right),\hfill \end{array}$$(25)

with the intrinsic scatter of ${\sigma }_{\mathrm{BCG}}=0.{554}_{-0.049}^{+0.040}$. Including the modelling of the M_⋆, BCG–M–z relation has no significant impact on that of the C_R–M–z relation. The intrinsic scatter of ${\sigma }_{\mathrm{BCG}}=0.{554}_{-0.049}^{+0.040}$ indicates significant variations in M_⋆, BCG at a fixed halo mass M. Given the large scatter and the shallow mass slope $B_{\mathrm{BCG}}=0.{184}_{-0.16}^{+0.068}$, we conclude that the BCG stellar mass is a highly scattering mass proxy. Given that the X-ray selection favors high-mass clusters at high redshift, we find strong degeneracy between the power-law indices of the mass (B_BCG) and redshift (γ_BCG) trends. To break the degeneracy, we introduced an informative prior on γ_BCG, which results in a steeper mass trend B_BCG = 0.38 ± 0.11.

We compared our M_⋆, BCG–M–z relation with both observational and simulation-based studies in the literature. We find satisfactory agreement in the intrinsic scatter and the scaling in the halo mass and redshift, though some studies exhibit systematic offsets in the absolute scale of the BCG stellar mass. However, these offsets are not significant given the large scatter and systematic uncertainties in M_⋆, BCG among previous studies.

Our M_⋆, BCG–M–z relation suggests that the BCG stellar mass fraction (M_⋆, BCG/M) at a fixed halo mass ($M_{\mathrm{piv}}=1.4\times {10}^{14}{h}^{-1}{\mathrm{M}}_{\odot }$) has remained stable since redshift z ≈ 0.8, with an increase in M_⋆, BCG at a level of ≈(20±8)%. As clusters continue to grow, the BCG stellar mass fraction significantly decreases with increasing total mass. This result aligns well with the picture of the rapid-then-slow BCG formation scenario suggested by previous studies (e.g. De Lucia & Blaizot 2007; Lidman et al. 2012; Lin et al. 2013). A combination of the eRASS1 sample and those sampling the low-mass regime at high redshift, possibly selected in the optical or with the planned next-generation X-ray missions (e.g. Zhang et al. 2020; Reynolds et al. 2023; Cruise et al. 2025), will shed light on the growth of the BCG mass using the statistical approach.

In summary, we presented a study of the galaxy population in the context of the BCG stellar mass based on the cluster sample constructed in the first-year database released by eROSITA-DE, paving a way for future work leveraging the final-year sample with the upcoming Stage IV experiments, such as the Rubin Observatory Legacy Survey of Space and Time (LSST Science Collaboration 2009), Euclid (Euclid Collaboration: Mellier et al. 2025), and the Nancy Grace Roman Telescope (Dore et al. 2019).

Acknowledgments

The authors thank anonymous referee for providing useful comments that led to improvements of this paper. I-Non Chiu acknowledges the support from the National Science and Technology Council in Taiwan (Grant NSTC 111-2112-M-006-037-MY3). This work made use of the computational and storage resources in the Academic Sinica Institute of Astronomy and Astrophysics (ASIAA) and the National Center for High-Performance Computing (NCHC) in Taiwan. I-Non Chiu thanks the hospitality of Ting-Wen Lan and Ji-Jia Tang at National Taiwan University (NTU). This work was supported by JSPS KAKENHI Grant Number JP19KK0076. VG acknowledges the financial contribution from the contracts Prin-MUR 2022 supported by Next Generation EU (M4.C2.1.1, n.20227RNLY3 The concordance cosmological model: stress-tests with galaxy clusters). E. Bulbul, V. Ghirardini, A. Liu, S. Zelmer, and X. Zhang acknowledge financial support from the European Research Council (ERC) Consolidator Grant under the European Union’s Horizon 2020 research and innovation program (grant agreement CoG DarkQuest No 101002585). This work is based on data from eROSITA, the soft X-ray instrument aboard SRG, a joint Russian-German science mission supported by the Russian Space Agency (Roskosmos), in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft- und Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its subcontractors, and is operated by NPOL with support from the Max Planck Institute for Extraterrestrial Physics (MPE). The development and construction of the eROSITA X-ray instrument was led by MPE, with contributions from the Dr. Karl Remeis Observatory Bamberg & ECAP (FAU Erlangen-Nuernberg), the University of Hamburg Observatory, the Leibniz Institute for Astrophysics Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the University of Tübingen, with the support of DLR and the Max Planck Society. The Argelander Institute for Astronomy of the University of Bonn and the Ludwig Maximilians Universität Munich also participated in the science preparation for eROSITA. The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from the Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University. This paper makes use of software developed for Vera C. Rubin Observatory. We thank the Rubin Observatory for making their code available as free software at http://pipelines.lsst.io/. This paper is based on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by the Subaru Telescope and Astronomy Data Center (ADC) at NAOJ. Data analysis was in part carried out with the cooperation of Center for Computational Astrophysics (CfCA), NAOJ. We are honored and grateful for the opportunity of observing the Universe from Maunakea, which has the cultural, historical and natural significance in Hawaii. The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg, and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation. This work is possible because of the efforts in the Vera C. Rubin Observatory (Juric et al. 2017; Ivezic et al. 2019) and PS1 (Chambers et al. 2016; Schlafly et al. 2012; Tonry et al. 2012; Magnier et al. 2013), and in the HSC (Aihara et al. 2018a) developments including the deep imaging of the COSMOS field (Tanaka et al. 2017), the on-site quality-assurance system (Furusawa et al. 2018), the Hyper Suprime-Cam (Miyazaki 2015; Miyazaki et al. 2018; Komiyama et al. 2018), the design of the filters (Kawanomoto et al. 2018), the data pipeline (Bosch et al. 2018), the design of bright-star masks (Coupon et al. 2018), the characterization of the photometry by the code Synpipe (Huang et al. 2018b), the photometric redshift estimation (Tanaka et al. 2018), the shear calibration (Mandelbaum et al. 2018b), and the public data releases (Aihara et al. 2018b, 2019. This work made use of the following software(s): IPython (Pérez & Granger 2007), SciPy (Virtanen et al. 2020), TOPCAT (Taylor 2005, 2006), matplotlib (Hunter 2007), Astropy (Astropy Collaboration 2013), NumPy (Van Der Walt et al. 2011), Pathos (McKerns et al. 2012), pyccl (Chisari et al. 2019), getdist (Lewis 2025), and the conda-forge project (conda-forge community 2015).

References