© The Authors 2024

1 Introduction

The so-called ΛCDM concordance cosmological model is currently the most widely accepted theory on the formation and evolution of our Universe, whereby cold dark matter (CDM) and dark energy parameterized by a cosmological constant (Λ) are the essential ingredients. It is a relatively simple model, where just a handful of cosmological parameters allow for accurate predictions for several observations, for instance: the cosmic microwave background (CMB Smoot et al. 1992; Fixsen et al. 1996), large-scale structure (Peebles 1980), abundance of light elements (Walker et al. 1991; Cooke et al. 2018), and accelerated expansion of the Universe (Riess et al. 1998; Perlmutter et al. 1999). In recent years, state-of-the-art experiments have enabled precise measurement of cosmological parameters. For instance, Planck satellite allowed for detailed constraints on the density field at the CMB epoch (Planck Collaboration VI 2020), while deep wide-area surveys, such as the Dark Energy Survey (DES), Hyper Supreme-Cam Subaru Strategic Program (HSC Y3) Survey, and Kilo Degree Survey (KiDS), have allowed for constraints to be placed on cosmological parameters combination as the parameter S₈ (Dark Energy Survey Collaboration 2022; Asgari et al. 2021; Amon et al. 2022; Secco et al. 2022; van den Busch et al. 2022; Li et al. 2023; Dalal et al. 2023; Dark Energy Survey and Kilo-Degree Survey Collaboration 2023), and the widening of the catalog of known SNe Ia and Cepheid variables allowed accurate calibration of distance ladder method and precise measurement of H₀ (Riess et al. 2022). However, tensions between CMB early-time measurements and large-scale structure late-time experiments have emerged: the so-called H₀-tension (e.g. Di Valentino et al. 2021) and the so-called S₈-tension (e.g. Heymans et al. 2021). These tensions indicate that the currently widely preferred cosmological model, the so-called ΛCDM, might be revised if new generation experiments increase the statistical significance of these cosmic discordances. With the discovery of neutrino oscillations between different flavors, electron, muon, and tau, it has become clear that neutrino eigenstates must have non-zero mass (Fukuda et al. 1998; Ahmad et al. 2002). This leads to several cosmological implications for CMB and the large-scale structure experiments (Bashinsky & Seljak 2004; Lesgourgues et al. 2013; Vagnozzi et al. 2017). Despite their small size and mass, neutrinos influence large-scale structures’ formation by hindering small-scale halos’ formation and leaving imprints on the gravitational collapse (Costanzi et al. 2013).

Clusters of galaxies formed from the highest density fluctuations in the early Universe and represent the largest virialized objects that have undergone a gravitational collapse. Their abundance through the halo mass function and spatial distribution traces the growth of linear density perturbations and, therefore, can be utilized as a probe of cosmology. The cluster abundance, dn/dM, is particularly sensitive to two cosmological parameters: the amplitude of the mass fluctuations on scales of 8h⁻¹ Mpc, namely σ₈¹, and the present-day matter density parameter, Ω_m². It is also dependent on the dark energy equation of state, w, and its evolution because of the growth suppression factor g(z) and the Universe’s expansion history (see Clerc & Finoguenov 2023, for a recent review). Compared to other probes, such as baryon acoustic oscillations (BAO, e.g. Alam et al. 2017), type Ia supernovae (SNe Ia, e.g. Scolnic et al. 2018; Brout et al. 2022), cosmic microwave background (CMB, e.g. Komatsu et al. 2011; Planck Collaboration VI 2020), and galaxy and cosmic shear auto- and cross-correlations (Dark Energy Survey Collaboration 2022; Dark Energy Survey and Kilo-Degree Survey Collaboration 2023), cluster number counts do not display strong degeneracies among the primary late time cosmological parameters discussed above, while retaining a comparable precision. Furthermore, the sum of the left-handed neutrino masses, ∑ m_v, can also be constrained via cluster abundances, particularly in combination with external probes, as the cosmic neutrinos have a non-negligible effect on the growth of matter perturbations during the matter and dark energy-dominated era (Vikhlinin et al. 2009). In fact, the number of low-mass halos in cluster abundance measurements can shed light on the summed masses of three left-handed neutrino species.

Multiwavelength surveys of clusters of galaxies have been widely utilized to place constraints on the cosmological parameters of the late-time Universe, with sample sizes ranging from tens to thousands of clusters. Among these are the surveys in the X-ray band (Vikhlinin et al. 2009; Mantz et al. 2015; Garrel et al. 2022; Chiu et al. 2023), in the optical band (Ider Chitham et al. 2020; Dark Energy Survey Collaboration 2020; Costanzi et al. 2021; Lesci et al. 2022), and mm-wave band through the inverse Compton scattering of CMB photons, called Sunyaev-Zel’dovich (hereafter SZ) effect (Planck Collaboration XX 2014; Planck Collaboration XXIV 2016; de Haan et al. 2016; Bocquet et al. 2019). The challenge for using cluster abundances in cosmological experiments lies in the reconstruction of the cluster mass distribution at a population level, as individual masses are not directly observable through X-ray, optical, and SZ observations. For cluster surveys, many observable mass proxies can be linked to cluster masses, such as richness, for instance, the number of photometrically assigned member galaxies, the temperature, X-ray luminosity, or gas mass of the intracluster medium (ICM), the X-ray luminosity, the SZ decrement, the tangential reduced shear imprinted by gravitational lensing on the shapes of background galaxies. These observables can be used to calibrate the observable-mass distributions and thereby reconstruct the cluster mass function. The observables in the optical, such as richness, suffer from projection effects and contamination, limiting their use in cluster abundance measurements (Dark Energy Survey Collaboration 2020; Ider Chitham et al. 2020; Costanzi et al. 2021). For cluster samples selected in X-rays, the count rate in the soft X-ray band, proportional to the integrated surface brightness, is a convenient observable that correlates tightly with the selection process and scales with halo mass and red-shift. Another advantage of the X-ray detection process is that count rate measurements do not suffer from significant projection effects (Garrel et al. 2022). Additionally, the scatter in count rate at a given mass and redshift is dominated by the X-ray emission mechanism, namely, thermal Bremsstrahlung; thus, it is not more sensitive to the underlying cosmology and environment than the scatter in richness.

A well-calibrated scaling relation between the total mass and the observed signal is critical for capitalizing on the statistical potential of large samples of clusters of galaxies. Historically, in their mass calibration, cluster counts experiments have employed X-ray cluster mass measurements relying upon the hydrostatic equilibrium, which assumes that the gravitational potential provided by the dark matter halo is balanced by thermal pressure (Vikhlinin et al. 2009; Mantz et al. 2010; Planck Collaboration XXIV 2016). However, recent numerical simulations indicate that the contribution of non-thermal pressure due to bulk motions and turbulence can lead to significant biases on the cluster masses measured with these traditional techniques (Lau et al. 2009; Rasia et al. 2012; Scheck et al. 2023). The use of biased X-ray mass measurements is also reflected in the cosmological parameters obtained from cluster counts (Planck Collaboration III 2013). In the last decade, a significant improvement over this approach has been introduced by integrating the gravitational lensing signal in the cluster mass calibration process of the cosmological studies (Mantz et al. 2015; Bocquet et al. 2019; Garrel et al. 2022; Chiu et al. 2023). The weak-lensing (WL) effect directly probes the cluster’s gravitational potential, thus allowing the masses to be accurately measured without any a priori assumption about their hydrostatic equilibrium properties. It is worth noting that this method is not entirely bias-free; the correlated large-scale structure along the line-of-sight, triaxiality, and assumed underlying dark matter distribution, baryonic feedback, and miscentering of the WL signal could introduce additional bias and scatter on the mass scaling relations (Hoekstra 2003; Meneghetti et al. 2010; Becker & Kravtsov 2011; Bahé et al. 2012). However, compared to X-ray measurements, the accuracy of WL measurements is inherently less impacted by hydrodynamical modeling as the majority of the mass is collisionless. Most importantly, since the WL mass measurements are calibrated with numerical simulations, the bias and systematic uncertainties can be reliably and quantitatively estimated and integrated into the final cluster mass calibration effort (Schrabback et al. 2018; Dietrich et al. 2019; Grandis et al. 2021a, 2024; Sommer et al. 2022; Zohren et al. 2022).

Another critical aspect in applying wide-area survey cluster catalogs in statistical analyses is understanding and adequately accounting for the underlying selection function. Surveys in the X-ray band are prone to significant selection effects, such as Eddington and Malmquist biases (Eddington 1913; Malmquist 1922). Additionally, the detected cluster sample often undergoes an optical confirmation process for redshift measurements; optical selection effects may play an essential role in the overall scientific results drawn from the observed cluster population. A good knowledge of the selection function, namely, the interface between an observed sample and the parent cluster population, is a prerequisite for statistical studies of cluster samples, such as the constraints on cosmology. Comparisons between samples detected at different wavelengths have shed light on our understanding of some of the selection biases (Nurgaliev et al. 2017; Rossetti et al. 2017; Willis et al. 2021; Ramos-Ceja et al. 2022). Thus, such comparisons have been successfully applied in the recent applications of the X-ray selected samples of clusters of galaxies in cosmology experiments (e.g., Garrel et al. 2022).

This work presents the constraints on cosmological parameters from the clusters of galaxies detected in the western Galatic hemisphere in the first eROSITA All-Sky Survey. As one of the two X-ray telescopes on board the Spectrum Roentgen Gamma Mission (SRG), eROSITA, was launched on July 13, 2019. The first All-Sky Survey (eRASS1) was successfully executed between December 2019 and June 2020 (Predehl et al. 2021; Sunyaev et al. 2021) and has resulted in the detection of nearly 930 thousand X-ray sources in the western Galactic hemisphere (179.9442 deg < l < 359.9442 deg), where the data and publication rights belong to the eROSITA German Consortium (Predehl et al. 2021; Merloni et al. 2024). In this first All-Sky Survey, we detect a total of 12247 optically confirmed galaxy groups and clusters spanning the redshift range 0.003 < z < 1.32 with a sample purity level of 86% (Bulbul et al. 2024; Kluge et al. 2024). In this work, we use a subsample of this primary eRASS1 cluster catalog, namely, the cosmology subsample constructed with a strict X-ray selection criterion to increase the purity of the sample to about 94%, as described in detail in Bulbul et al. (2024). In the cosmology subsample, we optically confirm 5259 securely detected clusters in the 12 791 deg² common footprint of eRASS1 and the DESI Legacy Survey DR10-South in the photometric redshift range of 0.1 < z_λ < 0.8. This catalog represents the largest ICM-selected cluster sample utilized for cosmological studies.

Besides the statistical power ensured by the extensive eRASS1 cluster catalog, our analysis offers significant improvements over previous works. We employ the observed X-ray count rate, which has the advantage of being coupled with the detection and selection process in the mass calibration phase. The selection effects are modeled using the end-to-end simulations of the eROSITA sky generated based on realistic models of the sources and their distribution in the sky (Comparat et al. 2020; Seppi et al. 2022). We incorporate the selection process into a forward model to represent the entire population of detected clusters in a Bayesian framework. Finally, the residual 6% contamination in our sample is statistically accounted for using a mixture model (see Bulbul et al. 2024; Kluge et al. 2024). For the mass calibration, we utilize a scaling relation between X-ray count rate and cluster mass obtained from the three weak lensing surveys: Dark Energy Survey Year Three (DES Y3), the Hyper Supreme-Cam Subaru Strategic Program Year Three (HSC Y3) Survey, and the Kilo Degree Survey (KiDS) for a total of 2588 clusters in an area of 4968 deg² overlapping with eRASS1 (Chiu et al. 2023; Kleinebreil et al. 2024). The biases due to miscentering, baryonic effects, and correlated large-scale structure along the line of sight are fully accounted for in our weak lensing analysis (see Grandis et al. 2024; Kleinebreil et al. 2024, for detailed analyses and comparisons).

This paper is organized as follows. In Sec. 2, we summarize the datasets in X-ray (eROSITA), optical (DESI Legacy Survey), and weak lensing (DES, KiDS, and HSC) utilized in our analysis. The selection function and Bayesian likelihood modeling method are presented in Sects. 3 and 4. The validation process of the cosmology pipeline and the blinding strategy are described in Sec. 5. In Sects. 6 and 7, we present our results on scaling relations and various models for cosmology. We discuss the systematics in our results and provide an outlook for future eROSITA surveys in Sec. 8. We provide our summary in Sec. 9. Our convention is as follows: the overdensity radius of R₅₀₀ is defined as the radius within which the mean cluster density is 500 times the cosmic critical density at the cluster’s redshift. The mass M₅₀₀ is the cluster mass within a sphere of R₅₀₀. With log, we are referring to the natural logarithm, and when using base 10, we explicitly indicate it with the subscript log₁₀. Throughout this paper, we use the terms consistency, tension, and discrepancy multiple times: we define here that by “consistent” we mean less than 3σ difference, by “in tension” we mean between 3σ and 5σ in difference, and by “discrepant”, we mean above 5σ in difference.

2 Galaxy cluster survey data

2.1 eRASS1 Cluster Catalog

The cluster catalog utilized in this work is based on the soft X-ray band (0.2–2.3 keV) catalog of the first eROSITA All-Sky Survey in the western Galactic hemisphere presented in Merloni et al. (2024). The details of the detection and cleaning processes and compilation of the cluster catalog are described in detail in Bulbul et al. (2024). We restrict ourselves to the cosmology subsample in this work, a subsample of the primary galaxy groups and clusters catalog. To compile the cosmology subsample, Bulbul et al. (2024) applied a conservative extent likelihood cut of ${\mathcal{L}}_{ext}>6$ to the primary sample, resulting in 11 141 cluster candidates detected as extended in the clean eRASS1 cluster catalog.

The confirmation, identification, and redshift measurements of the X-ray-selected cluster candidates are typically performed using the publicly available DESI Legacy Survey Data Release 10 data (LS DR10). LS DR10 is a compilation of three surveys that cover most of the extragalactic sky (Dey et al. 2019; Sevilla-Noarbe et al. 2021). The cosmology subsample is compiled using only the southern area below Dec. ≲ 32.5° (LS DR10-South hereafter) to maximize the homogeneity, as the observations have been taken with the same telescope and set of filters. This ensures that the subset of the eRASS1 clusters has homogeneous systematics on the photometric redshifts $\left({\hat{z}}_{\lambda }\right)$ and richness $(\hat{\lambda })$ measurements. The common footprint between the eROSITA survey in the western Galactic hemisphere and the southern LS DR10-South footprint is 12791 deg², and contains 8129 cluster candidates (Kluge et al. 2024).

We performed the confirmation and identification with the red sequence-based cluster finder algorithm eROMaPPer (Ider Chitham et al. 2020; Kluge et al. 2024). For the catalog used in this work, we combined the g, r, and z filter bands to ensure the homogeneity of the richness, redshifts, and contamination. In the cosmology catalog used in this work, we adopted photometric redshifts for the detected cluster to ensure homogeneity and uniform treatment for all the optically measured properties (see Bulbul et al. 2024, for further details). Of the 7077 candidate clusters in the LS DR10-South footprint, we were able to identify the optical counterparts of 6562 securely detected clusters.

We note that the detection of faint galaxies below the limiting luminosity L > 0.2 L_* is highly uncertain, and the measured richness artificially increases because of Eddington bias. To obtain reliable photometric redshift measurements with well-controlled biases and uncertainties, we applied a further flag of IN_ZVLIM==True (see Kluge et al. 2024, for further details). Additionally, we limited our redshift range to 0.1 < z < 0.8, where the photometric redshifts are the most reliable (see Kluge et al. 2024, for details). As a final cut, we removed all the clusters with measured richness smaller than 3 to decrease the contamination as these would be cases with no optical cluster counterpart (see Kluge et al. 2024, Sect. 3.2.3). This is because they would have a significantly larger chance of being contaminants and, thus, be likely to have the wrong photometric redshifts assigned. The application of richness measurements below 20 has been shown to be problematic in cluster number counts experiments with optical selection only (Dark Energy Survey Collaboration 2020). However in our case, richness is not used as a direct proxy for cosmology; rather, it is mostly used as a clean-up quantity that is very helpful in statistically discriminating between clusters and contaminations in our sample (see Sec. 4.2). Furthermore, richness for an ICM-selected sample (as in our case) is well behaved with small scatter (Saro et al. 2015; Bleem et al. 2020; Grandis et al. 2021b; Giles et al. 2022). The final cosmology subsample comprises 5259 securely confirmed galaxy clusters in the 12791 deg² LS DR10-South area. The projected spatial distribution of the clusters in this sample is shown in Fig. 1, overlaid on the exposure map of the eRASS1 survey. The median redshift in the catalog is 0.29. The redshift and richness histograms of all confirmed clusters are presented in the leftmost panel of Fig. 2. The mixture model, extensively described in Sec. 4.2, implies that the purity of this sample is at the level of 94%, ideal for cosmological studies with cluster counts (see Bulbul et al. 2024; Kluge et al. 2024, for further details).

Fig. 1 Sky distribution of the 5259 confirmed clusters of galaxies in the cosmology subsample, detected in the western Galatic hemisphere of the eRASS1 All-Sky Survey in equatorial coordinates. The common footprint between the Legacy Survey-South and eRASS1 reaches up to 12 791 deg2. The sample covers a redshift range of 0.1–0.8. The redshifts of the sources are color-coded, and the sizes of the marks are proportional to the observed eROSITA X-ray count rate. The overlaid image in the background shows the exposure map of the eRASS1 Survey in the 0.2–2.3 keV energy band. The regions plotted in solid lines show the footprints of optical surveys employed in the week lensing mass calibration DES in blue, KiDS in red, and HSC Y3 in green, covering 4968 deg2 in total.

Fig. 2 Photometric redshift, optical richness, and X-ray count rate histogram of the 5259 confirmed eRASS1 clusters in the cosmology subsample (top). The median redshift of the eRASS1 cosmology catalog is 0.29, the median richness is 35, and the median count rate is 0.33 cts/s.

2.2 Selection of the X-ray observable

The choice of an adequate X-ray observable employed in the internal cluster mass calibration process is vital in cluster cosmology studies. Ideally, the X-ray observable employed in the selection process should be adequately included in the selection model and should show a correlation with cluster mass. Additionally, if the X-ray observable has a negligible variation over the sky and a minimal dependence on cosmology (which are not strictly necessary), the computation cost of the forward hierarchical model significantly decreases. In this regard, we examine extensively several observables, for instance: the extent likelihood, ${\mathcal{L}}_{ext}$, detection likelihood, ${\mathcal{L}}_{det}$, and count-rate. In addition, we included commonly used observables in the literature, such as luminosity, gas mass, and the mass proxy Y_X in the observer’s frame, along with their correlation with mass using the eROSITA’s digital twin (Seppi et al. 2022). Among these observables, we find that the X-ray count rate satisfies all the conditions mentioned above, especially with the clear advantage that it has minimal dependence on the underlying cosmological model. Furthermore, the sky variation can be further decreased, and, similarly to how we have divided the counts by exposure map to correct for exposure time variation over the sky, we correct the count rate for line of sight absorption due to hydrogen in the so-called n_H-corrected count rate.

The raw count rate that is the output of the detection chain presented in Merloni et al. (2024) is only a rough estimate of the true clusters’ count rate and shows a large scatter in correlations with cluster mass. We improve these pipeline measurements through a rigorous X-ray analysis correcting for the Galactic absorption and a more accurate ICM and background modeling to reduce the scatter with the total mass. The details of the X-ray processing are provided in Bulbul et al. (2024); we briefly summarize the process here. X-ray image fitting is performed using the MultiBand Projector in 2D (MBProj2D) tool (Sanders et al. 2018), which forward-models background-included X-ray images of galaxy clusters to fit cluster and background emission in a total of seven bands between 0.3 and 7 keV. We obtain the corrected count rate from the best-fit cluster physical model, including a correction for Galactic hydrogen absorption at photometric redshifts z_λ. The measured count rate in the 0.2–2.3 keV band is used in the further weak lensing mass calibration process and is also shown in the rightmost panel of Fig. 2 and has been released in the cluster catalog (Bulbul et al. 2024).

2.3 Weak-lensing follow-up observations and analysis

In cosmology studies with galaxy cluster abundances, the cluster mass function should be determined with minimal bias. The X-ray observations alone do not provide direct cluster mass measurements, and the usage of hydrostatic mass estimates from X-rays is known to produce significant biases in cosmology results (Planck Collaboration III 2013). Therefore, we employed weak-lensing follow-up observations of the eROSITA-selected clusters to obtain reliable mass measurements with minimal bias in our mass calibration and cosmology analysis. This work uses dedicated analysis of the weak lensing signal around the clusters in three major optical surveys, HSC, DES, and KiDS, which have overlapping footprints with eRASS1 in the western Galactic hemisphere. The common footprints of these surveys, covering a total area of 4968 deg², are shown in Fig. 1. In particular, eRASS1 has a 4060 deg² common area with DES, 679 deg² with HSC, and 1108 deg² with KiDS³. The detailed analysis methods and the mass calibration with these surveys are covered in our companion papers by Grandis et al. (2024) and Kleinebreil et al. (2024), and are similar to Chiu et al. (2022) for HSC. A thorough comparison of the results from these independent surveys is presented in Kleinebreil et al. (2024). In the following sections, we summarize the shear extraction analysis and mass calibration strategy used in these three surveys.

2.3.1 Dark Energy Survey

Here, summarize here the measurements and calibration strategy presented in Grandis et al. (2024), aimed at extracting the tangential shear profile around eRASS1 clusters in the DES footprint. In this analysis, we utilize observations from the first three years of observations (DES Y3), covering the full DES survey footprint. The DES Y3⁴ shape catalog (Gatti et al. 2022) is built from the r, i, z-bands using the METACALIBRATION pipeline (Huff & Mandelbaum 2017; Sheldon & Huff 2017). Our analysis uses the same selection of lensing source galaxies in tomographic bins as the DES Y3 3x2 pt analysis (Dark Energy Survey Collaboration 2022). This selection is defined and calibrated in Myles et al. (2021) and Gatti et al. (2022). We perform a background selection for each eRASS1 cluster and group by weighting the respective tomographic bins and excluding tomographic source redshift bins whose median redshift is smaller than the cluster redshift. In the calibration steps, we provide (1) a fitting formula for the effective shape noise as a function of cluster redshift, (2) a measure of the contamination of our background selection by cluster sources, validated with three independent measurement techniques, (3) a statistical test that the cross shear is consistent with zero, (4) tangential shear profiles around individual clusters binned in angular annuli, with their respective source redshift distributions and shape noise induced uncertainty, and (5) a percent-level accurate calibration of the shape and photo-z uncertainties propagated from the calibrations as mentioned above.

Considering the overlap between the DES Y3 and the eRASS1 footprints, we produced tangential shear data for 2201 eRASS1 galaxy clusters, with a total signal-to-noise ratio (S/N) of 65 in the tangential shear profile.

2.3.2 The Hyper Suprime-Cam (HSC) survey

We extract the weak-lensing measurements of eRASS1 clusters using the data from the Hyper Suprime-Cam Subaru Strategic Program (HSC SSP, Aihara et al. 2018), a wide and deep multiband (g, r, i, z, and Y) survey in the optical. We make use of the latest three-year weak-lensing data (S19A) (Li et al. 2022) with a total area coverage of ≈500 deg², by which 96 eRASS1 clusters from the western Galactic half of the eROSITA sky are covered. As detailed in Li et al. (2022), the shape measurements in the HSC weak-lensing data are carried out in the i-band imaging with a limiting magnitude of i ≤ 24.5 mag, and have been extensively calibrated through image simulations (Mandelbaum et al. 2018). This results in a shape catalog with a density of ≈ 20 galaxies/arcmin², which has been used to deliver robust constraints on cosmology using cosmic shears (Dalal et al. 2023; Li et al. 2023; Sugiyama et al. 2023; More et al. 2023; Miyatake et al. 2023). For the weak-lensing analysis of eRASS1 clusters, we follow the same procedure as presented in Chiu et al. (2022), which we refer to for more details.

The weak-lensing measurements were extracted based on individual clusters. For each eRASS1 cluster in the HSC S19A footprint, lensing sources were selected on the basis that they were securely located at the background of the cluster. This was done by requiring ${\int }_{z_{\mathrm{cl}}+0.2}^{\mathrm{\infty }}P(z)\mathrm{d}z>0.98$, where z_cl is the cluster redshift and P(z) is the photometric redshift distribution of each source. We used the DEmP code (Hsieh & Yee 2014) to estimate the photometric redshift (Tanaka et al. 2018; Nishizawa et al. 2020). The stacked photometric redshift distribution of the sources selected in non-cluster regions was used to estimate the lensing efficiency, which is needed to obtain the weak-lensing mass from the modeling of the shear. With the selected source sample, the tangential shear profile g_t (θ) as a function of the angular radius θ around the cluster center was calculated by following Sect. 3.4 in Chiu et al. (2022). The best-fit X-ray peak (RA_XFIT and DEC_XFIT), the output of the MBProj2D provided in the cosmology catalog of Bulbul et al. (2024), was used as the cluster center to extract the shear profiles around. The radial binning was defined as ten logarithmic bins in the angular space, corresponding to 0.2 Mpc/h and 3.5 Mpc/h in the physical coordinate in a fiducial flat ΛCDM cosmology with $h\equiv \frac{H_0}{100\mathrm{km}/\mathrm{s}/\mathrm{Mpc}}=0.7$ and Ω_m = 0.3. We discarded the innermost three bins to mitigate the systematics in the weak-lensing modeling at the cluster core (≲0.5 Mpc/h). The uncertainty of the shear profile arising from the shape noise and uncorrelated large-scale structures is accounted for in the weak-lensing covariance matrix by following Sect. 3.5 in Chiu et al. (2022).

The contamination of cluster members to the shear signal is quantified by a P (z)-decomposition method (Gruen et al. 2014; Melchior et al. 2017; Varga et al. 2019). We find no sign of cluster member contamination in the stacked shear profiles in different redshift and richness bins. Conservatively, and as in the previous analysis of the eROSITA Final Equatorial Depth Survey (eFEDS) clusters (Chiu et al. 2022), we set the cluster member contamination to be zero across the radial bins with a 2σ upper limit of 6% at the innermost bin (R ≈ 0.2 Mpc/h). It then decays outwards following a projected Navarro-Frenk-White (NFW) model (Navarro et al. 1996). With a calibration against the COSMOS 30-band photo-z sample and a spectroscopic sample, the bias in the photometric redshift of the selected sources can be quantified as ≲6% for the cluster redshift z_cl ≲ 1.2 and accounted for in the weak-lensing modeling. The multiplicative bias, m, in the shape measurements as a function of the source redshift is quantified with a systematic uncertainty, δm, at 1.7% level (Mandelbaum et al. 2018).

As a result, the shear profile g_t (θ), the lensing covariance matrix that serves the measurement uncertainty and the photometric redshift distribution of the selected source sample are obtained as the HSC weak-lensing data products for 96 eRASS1 clusters, with a total S/N of 40. The systematic uncertainties, including the cluster member contamination, the photo-z bias, and the multiplicative bias associated with these weak-lensing measurements, are quantified and taken into account in deriving the weak-lensing mass bias b_WL(M, z) (see Kleinebreil et al. 2024).

2.3.3 Kilo Degree Survey

We used the gold sample of weak lensing and photometric redshift measurements from the fourth data release of the Kilo-Degree Survey⁵ (Kuijken et al. 2019; Wright et al. 2020; Hildebrandt et al. 2021; Giblin et al. 2021), hereafter referred to as KiDS-1000.

We extracted individual reduced tangential shear profiles for a total of 236 eRASS1 galaxy clusters in both the KiDS-North field (101 clusters) and the KiDS-South field (136 clusters), as both have overlap with the eRASS1 footprint with a total S/N of 19. A detailed presentation of the method is given in Kleinebreil et al. (2024). Our analysis was run separately in the tomographic redshift bins of KiDS-1000, which provide spectroscopically calibrated redshift distributions for five photometric galaxy redshift bins. We also used the enhanced redshift calibration presented in van den Busch et al. (2022), which provides the basis for calibrating the shape measurement and photo-z systematic uncertainties.

Our analysis includes a model for the cluster member contamination that is based on radial background galaxy number density profiles (e.g. Applegate et al. 2014). The model accounts for the effect of blending background galaxies with cluster galaxies on the number density profiles. For this, we simulated KiDS-like galaxy images using GALSIM (Rowe et al. 2015), inject them into the actual r-band detection images of KiDS, and re-run the source detection to extract radial detection probability profiles.

We point out the importance of KiDS because it provides valuable overlap with HSC and DES, thus allowing us to perform several consistency checks and thorough comparisons between three surveys employed in the mass calibration in this work (Kleinebreil et al. 2024). In particular, there are 45 clusters in the overlap region between HSC and KiDS and 125 clusters in the overlap region between KiDS and DES. We also remark that there is no overlap between DES and HSC in the sky area relevant to the cluster cosmology catalog; therefore, these consistency checks are vital to show the robustness of our analysis.

3 The X-ray selection function

Determining the mass distribution and the consequent cosmological analyses with cluster surveys heavily depend on the accurate knowledge of the survey’s incompleteness, contamination, and selection effects. To characterize the survey selection function, we use the high fidelity end-to-end simulations of eRASS1 from the eROSITA’s digital twin (Seppi et al. 2022). In this section, we provide brief descriptions of the numerical simulations and source modeling employed in creating the selection function since the detailed modeling and testing of it are already published in Comparat et al. (2019, 2020), Liu et al. (2022), Seppi et al. (2022), and Clerc et al. (2024).

3.1 Simulations

We used all the snapshots of UNIT1i dark matter-only simulation with redshift z < 6.1 to generate a full sky light cone (Chuang et al. 2019). Its box size (1 h⁻¹ Gpc) and resolution (4096³ particles) enable reliable prediction of the large-scale structure and its haloes down to 10¹¹ M_⊙. This resolution is sufficient to model the X-ray active galactic nucleus (AGN) population at the depth of eRASS1 (Comparat et al. 2019). The box size is large enough to contain a fair census of massive clusters and to suppress finite volume effects (Comparat et al. 2017; Klypin & Prada 2019). The light cone is populated using models of the X-ray emission from clusters and AGNs (Comparat et al. 2019, 2020; Liu et al. 2022; Seppi et al. 2022).

In the digital twin, the AGNs are simulated with an abundance matching relation between X-ray luminosity and stellar mass (Comparat et al. 2019). By construction, this reproduces the AGN luminosity function as modeled by Aird et al. (2015). We find that the AGN model performs well as it reproduces the general statistics of the AGN population: the source number with respect to the detection flux distribution (log N-log S), luminosity function, and two-point correlation function (Liu et al. 2022; Comparat et al. 2023).

The model for clusters of galaxies is based on the X-ray observations in the literature (e.g., Reiprich & Böhringer 2002; Pierre et al. 2016; Adami et al. 2018; Sanders et al. 2018; Eckert et al. 2019). A Gaussian process randomly draws temperature, mass, redshift, and emissivity profiles from a combined covariance matrix. The resulting model from Comparat et al. (2020) reproduces the scaling relations between mass, temperature, and luminosity (e.g., Giles et al. 2016) at the cluster scale. However, since the model relies on observed clusters with a high S/N it is not trained at the low-mass group’s scale (M_500c < 5 × 10¹³ M_⊙), and it ends up over-predicting the X-ray luminosity of these objects. To overcome this effect, the fluxes are scaled to match the expected scaling relation between luminosity and stellar mass (Anderson et al. 2015) in the low-mass regime (see Seppi et al. 2022, for details). The profiles are assigned to dark matter halos by a near-neighbor process, where the central emissivity is related to the dynamical state of the dark matter halo (see Seppi et al. 2021, for a detailed discussion on dark matter halo property and dynamical state). In total, we simulated 1 116 758 clusters and groups.

3.2 Event generation and catalog inference

To mimic the eROSITA survey, the next step is to generate X-ray events with the SIXTE simulator (Dauser et al. 2019). SIXTE uses the ancillary response file, the redistribution matrix file, PSF, background, and the eRASS1 attitude file of the spacecraft in the most up-to-date calibration database in a self-consistent manner with the observations. This method allows the simulations to accurately account for the instrumental response and exposure variations across the sky due to the scanning strategy of SRG. The output is a collection of event lists with a realistic angular position, energy, and arrival time. These event lists are processed through the eSASS software (Brunner et al. 2022) as described in Merloni et al. (2024). We matched the resulting source catalog to the input catalog from simulations by tracing the origin of each photon. The sources are grouped into five classes, described as follows: (1) PNT is a primary counterpart of a simulated point source. The detected source is uniquely assigned to an AGN or star. (2) EXT is a primary counterpart of a simulated extended source. The detected source is uniquely assigned to a cluster. (3) PNT2 is a secondary counterpart of a simulated point source. The detected source is assigned to an AGN or star, but the primary detection containing a large number of events is already classified as PNT or EXT. This is a split source corresponding to an AGN or star. (4) EXT2 is a secondary counterpart of a simulated extended source. The detected source is assigned to a cluster, but the primary detection containing a large number of events is already classified as PNT or EXT. This is a split source corresponding to a cluster. (5) BKG is a false detection due to random fluctuations in the background. The entry in the source catalog is not associated with any simulated source.

For AGNs and stars, we created a single realization of the events. The sources and the realization scheme are described in detail in Seppi et al. (2022). We create 100 realizations for the clusters initiated with a different random seed. The remaining 99 other realizations are used to determine accurately the selection function (see Clerc et al. 2024, for further details). Compared to the original simulation (see Sect. 3.2 in Seppi et al. 2022), the generation of cluster events was successful in the whole western Galactic half of the All-Sky survey.

3.3 Assigning a detection probability

We build a detection probability model based on the simulated eROSITA survey as described in the previous section. We select Gaussian process (GP) classifiers, which satisfy the requirements of short computational response, reduced prior assumptions, and can issue meaningful uncertainties. Given a set of five selected labels (features) from the simulated catalogs, we train a real-valued latent function over the corresponding five-dimensional parameter space, under the assumption of a GP prior (e.g., Rasmussen & Williams 2006). We map its values to the interval [0, 1] through a “probit” bijection (link function). A Bernoulli likelihood describes the distribution of binary detection flags associated with the simulated clusters, given the values of the mapped functional. Our GP prior assumes a pairwise covariance function of the form k(r) = σ² exp(−r²/2), with σ being a free hyperparameter and r as the Euclidian distance between pairs of points, normalized along each dimension by a free hyperparameter (lengthscale). This kernel is the widely-used stationary squared exponential (Rasmussen & Williams 2006); lengthscales set the size of “wiggles” in the model while σ controls the average deviation of the latent function from its mean (that is supposed uniform). Hyperparameters are optimized to minimize the marginal likelihood in the nominal implementation of GP classification. The quantity of interest is the posterior of the latent function given the set of training points. For any new set of feature values, the posterior predicts the distribution of the corresponding values of the latent function; a probabilistic prediction can be obtained by integrating this distribution against the link function.

In practice, we select two-thirds of the simulated clusters in the eRASS1 sky to train a classifier with features being the true redshift, z, the Galactic column density (n_H)-corrected count rate, C_R⁶, and the sky position, ℋ_i (this indicates a combination of three sky-dependent quantities, uniquely determined by its sky coordinates: the local background surface brightness, the local n_H, and the exposure time) A sky mask is applied to reproduce the sample selection, and we only considered objects simulated at 0.05 < z < 0.85, slightly larger than the cosmology subsample to take into account leaking in and out of detected clusters due to redshift uncertainty. After normalizing the feature values into an adequate parameter space, we resort to using stochastic variational Gaussian processes (SVGP Hensman et al. 2015) implementation in the GPy library (GPy 2012) by taking advantage of its numerical scalability properties, that is well-suited to our large training set (of size ≃4.8 × 10⁵). The sparse GP classification scheme (Titsias 2009) embedded in the algorithm relies on a small set of inducing points (30 in our case). These inducing points “summarize” the whole training dataset and are optimized with other hyperparameters. SVGP implements stochastic optimization on several mini-batches and variational inference (Hensman et al. 2013) to approximate the marginal likelihood, further increasing efficiency while preserving good convergence properties. We found that minibatches of size 2¹⁶ and 12 000 iterations provide good results within a tractable time. The optimization of a set of 651 hyperparameters⁷ relies on the ADADELTA method (Zeiler 2012) as implemented in the climin library.

3.4 Selection models for validation and blinding

For validation and blinding purposes (see details in Sec. 5), six models were computed, differing only by the exact definition of an object that qualifies as “detected” in the simulation. We modified the lower thresholds in selection parameters (${\mathcal{L}}_{ext}$, EXTENT), choosing (3.6, 16.6), (3.9, 15.4), (4.4, 15.1), (4.8, 13.7), (5.5, 11.6), and (6, 0). By design, the six models predict different populations of detected clusters, yet we have tuned the thresholds to ensure the total number of objects is similar. Among these models, only one corresponds to the true set of selection cuts, which is described in Sec. 2 (${\mathcal{L}}_{ext}>6$, EXTENT > 0). Each model is tested against the remaining third of the simulation sets (each of size ≃2.5 × 10⁵). Figure 3 shows the outcome of testing the true selection model against a concatenation of 10 test samples. The figure highlights the satisfactory agreement between the model-predicted detection probability and the actual detection rate per bin to reflect the internal consistency of our selection model. In addition to this test, the total number of predicted detections is checked to be consistent with the number of actual detections in the test samples. We repeat this series of tests using a collection of 10 independent Poissonian realizations of the eRASS1 All-Sky survey to ensure the stability of the results against photon shot noise.

While the most valuable outcome of the models is the expectation values for the probabilistic predictions P(I|z, C_R, ℋ_i), the probability of detecting a cluster (represented by I, standing for inclusion or detection, as in Mantz et al. 2010) at a given sky position (ℋ_i) a given redshift (z), and a specific intrinsic count rate (C_R), the very nature of GP enables deriving confidence levels around the expectation, conditioned to the training set and the GP prior; in general, the width of confidence levels correlates with the sparsity of the training sample, that depends on location in the feature parameter space. This assigned probability is fully integrated into the likelihood for modeling the selection function. In Figure 4, we show a summary plot showing example outputs of the selection function at different exposure times and redshifts as a function of count rate.

Fig. 3 Reliability of the selection function model, as estimated from a test sample comprising 2566 057 simulated galaxy clusters, none are used in the training phase. The x-axis represents the detection probability predicted by the selection function model. Top panel: y-axis representing the fraction of detected objects among all simulated objects in each bin. The shaded area corresponds to the 1σ binomial uncertainty. The dotted line represents the 1:1 relation. Bottom panel: histogram of clusters in the test sample per bin of predicted detection probability.

3.5 Optical selection function

As described in Sec. 2.1, the cosmology cluster catalog is compiled after optical measurements of richness and redshift measured through eROMaPPer run on the initial cluster catalog as described in Kluge et al. (2024). Effectively, these are redshift and richness cuts, $0.1<\hat{z}<0.8$ and $\hat{\lambda }>3$, which are applied to the primary X-ray catalog presented in Bulbul et al. (2024). This means our selection is a two-step process, where X-ray selection is applied first and optical selection later.

This implies that the redshift and richness cuts do not affect the X-ray selection function as these cuts are applied after the first initial X-ray selection. The modeling of this optical selection through the optical selection function is given by the Heaviside step function, which (in practice) affects only the integration limits in likelihood terms where richness or redshift are involved, as discussed in Grandis et al. (2020).

We stress that the role of the optical selection function is to model the contamination in the cluster sample. The cluster detection and completeness are purely governed by the X-ray detection pipeline and modeled by the X-ray selection function through the eROSITA’s digital twin (Seppi et al. 2022). The optical identification process does not remove a significant number of clusters in the sample owing to the low richness limit $(\hat{\lambda }>3)$ we employ in this process. To demonstrate this, eROMaPPer was run on SZ-based cluster samples, which are known to be complete. In the companion paper (Kluge et al. 2024) Sect. 5, we show that more than 97% (optical completeness) of SZ-detected clusters are identified by eROMaPPer. The remaining few percent are ambiguous cases where clusters overlap along the line of sight (and therefore also the X-ray signals overlap), a known effect for cluster samples. This test demonstrates that our optical confirmation process does not yield to optical incompleteness and the method presented in this section successfully models both purity and completeness of the sample. The small impact of optical selection is further illustrated in Appendix D (particularly in Fig. D.1).

Fig. 4 Representation of a model P(I|z, CR, ℋi) for fixed values of cluster redshift, z = [0.15, 0.3, 0.5], local galactic absorption nH = 3.5 × 1020 cm−2, a nominal background level 5.2 × 10−15erg/s/cm2/arcmin2, three values of exposure times [50, 150, 500] s, as a function of count-rate.

4 Likelihood

This section derives the global likelihood function used to fit the data, sample the contamination, and the scaling relations to produce the cosmological inference. In the following, a variable with a “hat” indicates an observed quantity, while without it indicates an intrinsic quantity, for instance ${\hat{C}}_R$ indicates the observed n_H-corrected count rate, while C_R indicates the intrinsic n_H-corrected count rate.

4.1 Global likelihood function

The likelihood is the probability of the data set modeled as a function of the parameters of the statistical model. Similarly to the previous studies in the literature (e.g., Mantz et al. 2010, 2015; Bocquet et al. 2019, 2024; Zubeldia & Challinor 2019; Chiu et al. 2023), we started from a Poisson-differential log-likelihood for our set of observables, X-ray n_H corrected count rate. ${\hat{C}}_R$ redshift, $\hat{z}$, sky position, ${\hat{\mathcal{H}}}_i$, richness, $\hat{\lambda }$, and weak-lensing tangential shear profile, ${\hat{g}}_t$; $$\begin{array}{c}\log \mathcal{L}=\sum _j\log \left(\frac{\mathrm{d}{N}_{tot}}{\mathrm{d}{\hat{C}}_R\mathrm{d}\hat{z}\mathrm{d}\hat{\lambda }\mathrm{d}{\hat{g}}_t\mathrm{d}{\hat{\mathcal{H}}}_i}P(I\mid {\hat{C}}_R,\hat{z},\hat{\lambda },{\hat{\mathcal{H}}}_i,{\hat{g}}_t)\right)\\ -\int \cdots \int [\frac{\mathrm{d}{N}_{tot}}{\mathrm{d}{\hat{C}}_R\mathrm{d}\hat{z}\mathrm{d}\hat{\lambda }\mathrm{d}{\hat{g}}_t\mathrm{d}{\hat{\mathcal{H}}}_i}P(I\mid {\hat{C}}_R,\hat{z},\hat{\lambda },{\hat{\mathcal{H}}}_i,{\hat{g}}_t)\\ \mathrm{d}{\hat{C}}_R\mathrm{d}\hat{z}\mathrm{d}\hat{\lambda }\mathrm{d}{\hat{g}}_t\mathrm{d}{\hat{\mathcal{H}}}_i],\end{array}$$(1)

where the first term runs over all detected clusters, while the second one is a computation representing the total number of objects detected in the eRASS1 survey.

We highlight that (as tends to happen in the process of detecting and including clusters in our dataset) our selection function does not depend on the measured tangential shear profile. Further, it is crucial that we do not discard the measured shear profiles whose signal is consistent with noise from the fitting process to avoid having an effect on the selection function. Furthermore, since the likelihood mimics the detection and measuring process, the selection function can be separated into the X-ray selection followed by the optical selection, as follows: $$P(I\mid {\hat{C}}_R,\hat{z},\hat{\lambda },{\hat{\mathcal{H}}}_i,{\hat{g}}_t)=P(I\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i)\mathrm{\Theta }(\hat{\lambda }>3),$$(2)

where $\mathrm{\Theta }(\hat{\lambda }>3)$ is the Heaviside step function, equals to 1 where $\hat{\lambda }>3$, and equal to 0 where $\hat{\lambda }<3$. This cut is necessary to reproduce the richness cut applied by the optical confirmation process with eROMaPPer code runs. For more information, we refer to Sect. 2.1 and Kluge et al. (2024).

We, therefore, isolate richness and tangential shear profile likelihood terms. In practice, the first term in Equation (1) can be simplified into: $$\frac{\mathrm{d}{N}_{tot}}{\mathrm{d}{\hat{C}}_R\mathrm{d}\hat{z}\mathrm{d}\hat{\lambda }\mathrm{d}{\hat{g}}_t\mathrm{d}{\hat{\mathcal{H}}}_i}=\frac{\mathrm{d}{N}_{tot}}{\mathrm{d}{\hat{C}}_R\mathrm{d}\hat{z}\mathrm{d}{\hat{\mathcal{H}}}_i}P(\hat{\lambda },{\hat{g}}_t\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i).$$(3)

Given these considerations, it is therefore consequent that the measurement of tangential shear profiles does not directly affect the total number of detected objects. This means that in the likelihood term that computes the total number of detected objects, the second of Equation (1), the dependence on tangential shear disappears, and Equation (1) therefore becomes; $$\begin{array}{rl}& \log \mathcal{L}=\sum _j\log (\frac{\mathrm{d}{N}_{tot}}{\mathrm{d}{\hat{C}}_R\mathrm{d}\hat{z}\mathrm{d}{\hat{\mathcal{H}}}_i}P(I\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i)\mathrm{\Theta }(\hat{\lambda }>3))\\ & +\sum _j\log P(\hat{\lambda },{\hat{g}}_t\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i,I)\\ & -\int \cdots \int [\frac{\mathrm{d}{N}_{tot}}{\mathrm{d}{\hat{C}}_R\mathrm{d}\hat{z}\mathrm{d}\hat{\lambda }\mathrm{d}{\hat{\mathcal{H}}}_i}P(I\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i)\mathrm{\Theta }(\hat{\lambda }>3)\\ & \mathrm{d}{\hat{C}}_R\mathrm{d}\hat{z}\mathrm{d}\hat{\lambda }\mathrm{d}{\hat{\mathcal{H}}}_i]\end{array}$$(4)

Additionally, although optical richness and weak lensing shear profiles are correlated signals, we can safely assume that the correlation between the two quantities is a subdominant factor. The richness used in this work involves a level of X-ray information which makes the richness better behaved than if it was measured from a purely optical survey without the X-ray or SZe information (Saro et al. 2015; Bleem et al. 2020; Grandis et al. 2021b; Giles et al. 2022). In fact, the extra prior regarding the location of the center of the halo and the presence of an X-ray emitting halo allow for more accurate richness measurement. Additionally, a significant reduction in the number of halos with low mass and spurious sources not associated with X-ray emitting halos that pollute the mass-observable space because of line-of-sight projection effects is elevated in the ICM-selected samples. As a consequence, the intrinsic correlation between weak lensing mass and richness is expected to be small as it is absorbed in the other two correlations we consider in this work, specifically the correlation between WL mass and count rate and the correlation between richness and count rate. Therefore, the likelihood of observing the measured tangential shear profile and the measured richness are uncorrelated and can be simplified as: $$P(\hat{\lambda },{\hat{g}}_t\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i,I)=P(\hat{\lambda }\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i,I)P({\hat{g}}_t\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i,I).$$(5)

Thus, we can write our final log-likelihood as: $$\begin{array}{rl}& \log \mathcal{L}=\sum _j\log \left(\frac{\mathrm{d}{N}_{tot}}{\mathrm{d}{\hat{C}}_{R}d\hat{z}\mathrm{d}{\hat{\mathcal{H}}}_i}P(I\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i)\right)\\ & -\int \cdots \int [\frac{\mathrm{d}{N}_{tot}}{\mathrm{d}{\hat{C}}_R\mathrm{d}\hat{z}\mathrm{d}\hat{\lambda }\mathrm{d}{\hat{\mathcal{H}}}_i}P(I\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i)\mathrm{\Theta }(\hat{\lambda }>3)\\ & \mathrm{d}{\hat{C}}_R\mathrm{d}\hat{z}\mathrm{d}\hat{\lambda }\mathrm{d}{\hat{\mathcal{H}}}_i]\\ & +\sum _j\log P(\hat{\lambda }\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i,I)\\ & +\sum _j\log P({\hat{g}}_t\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i,I).\end{array}$$(6)

The first term in Equation (6) represents the number count likelihood coming from X-ray information, and the second term is the expected total number of detected objects, which depends on the selection function. In contrast, the latter two terms describe the information we obtain from our follow-up observables: optical richness and tangential shear.

In the following subsections, we expand on the terms in the likelihood: first, we show how our mixture model allows us to distinguish the clusters from the contamination, for instance, AGNs and background fluctuations, in the sample. We then examine the X-ray likelihood term, the joint X-ray/optical cluster total number of predicted clusters, the WL mass calibration, and the richness mass calibration. The integrals are computed to obtain the likelihood value for a set of parameters in a range where the integrands differ from zero. Specifically, redshift from 0 to 0.9, mass from 5 × 10¹² M_⊙ to 5 × 10¹⁵ M_⊙, richness from 1 to 600, and count rate from 10⁻³ to 100 cts/s.

4.2 Mixture method for the contamination modeling

Statistical studies previously performed with limited samples of clusters often had a negligible fraction of contaminants. They did not require extensive modeling as the fraction of expected contaminants was sub-dominant compared to the Poisson noise (Planck Collaboration XXIV 2016; Garrel et al. 2022). The effects of contamination become relevant only when the contamination fraction exceeds the Poisson noise of the data. It has recently become a routine to include the contamination modeling in cluster abundance studies (see the implementation in the SPT-SZ sample Bocquet et al. 2019)⁸. As large cluster samples from wide-area surveys become available, contamination modeling becomes crucial in statistical cosmological studies with cluster abundances.

In eRASS1 cosmology subsample, after the ${\mathcal{L}}_{ext}>6$ cut is applied, the expected contamination level is estimated to be 20% from our realistic simulations (Seppi et al. 2022). The optical confirmation process is expected to reduce the contamination to ~6% after removing the unconfirmed sources in the sample (Kluge et al. 2024; Bulbul et al. 2024). The Poissonian noise of our data $(\sqrt{N}/N)$ is at the 1.4% level; therefore, in our analysis, contamination modeling is crucial to derive unbiased results from statistical fits. To this end, we develop a mixture model, which is a probabilistic model that represents multiple subpopulations of sources within a given sample without requiring exact knowledge of the nature of each object.

Informed by the prior work on the X-ray selection on simulation (Seppi et al. 2022), we model our cluster sample as the mixture of three populations: clusters, AGNs, and random sources (RSs) from background fluctuations. Motivated by the fact that for clusters, both richness and count rate correlate with halo mass at a given redshift, we model the cluster population with multivariate observable mass scaling relations (see Sec. 4.3). For the other two populations, the X-ray emission is modeled as uncorrelated with the richness-redshift distribution, with the count rate distribution as a function of sky position coming from our X-ray simulations and the richness–redshift distribution coming from optical follow-up runs of the point source catalog and random position for the AGN and RS populations, respectively. This serves two related purposes. Firstly, we can fit for the contamination of our sample on the fly. Secondly, we ensure that for the cluster population, the optical signal and the derived redshift are linked to the same host halo and its X-ray signature, preventing chance superpositions. In the context of X-ray-selected cluster samples, erroneously associating a high redshift to a lower count rate X-ray signature would mislead us into interpreting that object as a high mass, high redshift cluster, biasing the cosmological inference. As a consequence, we need to model the richness-mass relation of the cluster population (see Sec. 4.5).

In practice, we can a-posteriori assign a probability to each source belonging to one or the other sub-population based on its measured properties (Kluge et al. 2024). In detail, our realistic simulations suggest that the three classes of X-ray sources comprise the extended source catalog: clusters, AGNs, and RSs from background fluctuations. Clusters are the population of interest and form the majority of the cleaned extended source sample (~90–95%). The global population modeling can be represented as a sum of cluster distribution, which is cosmology dependent, and the distribution of RSs and AGNs, which are assumed to have known probability density functions but unknown overall numbers, thus allowing for a simultaneous fit for cosmology as well as the overall number of contaminants in the likelihood Equation (6).

The complete derivation of the mixture model and its effects on each likelihood term in Equation (6) is shown in Appendix A. In short, the probability density functions for every observables $\hat{\mathcal{O}}$ is written as: $$P(\hat{\mathcal{O}})={\displaystyle \sum _{\mathrm{i}\in \{\mathrm{C},\mathrm{A}\mathrm{G}\mathrm{N},\mathrm{R}\mathrm{S}\}}}{f}_{\mathrm{i}}P(\hat{\mathcal{O}}\mid \mathrm{i}),$$(7)

where this equation holds for all our observables $\hat{\mathcal{O}}$ as ${\hat{C}}_R,\hat{z}$, $\hat{\lambda }$, and ${\hat{g}}_t.f_{\mathrm{i}}$ are the fraction of each population in the eRASS1 sample and sum up to unity. Therefore, the mixture model adds two extra parameters to the fit, which respectively describe the amount of AGNs and RSs in our cluster cosmology sample, while the constraint of unity summed fraction fixes the number of clusters.

4.3 Cluster number counts

This subsection focuses on the differential cluster number counts term, $\frac{\mathrm{d}{N}_C}{\mathrm{d}{\hat{C}}_R\mathrm{d}\hat{z}d{\hat{\mathcal{H}}}_i}P(I\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i)$, which is a significant term in Equation (6), being involved in both terms of the differential Poissonian likelihood required to compute the number counts likelihood. Note that the corresponding terms for AGNs and RSs from the mixture model are directly recoverable from the mixture model; see Equation (A.1).

We consider this term as a marginalized likelihood over true count rate, true mass, and true redshift⁹.

$$\begin{array}{rl}& \frac{\mathrm{d}{N}_C}{\mathrm{d}{\hat{C}}_R\mathrm{d}\hat{z}\mathrm{d}{\hat{\mathcal{H}}}_i}P(I\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i)=P\left({\hat{\mathcal{H}}}_i\right)\iiint P(\hat{z}\mid z)P(I\mid z,C_R,{\hat{\mathcal{H}}}_i)\\ & \times P\left({\hat{C}}_R\mid C_R\right)\frac{\mathrm{d}N}{\mathrm{d}M\mathrm{d}V}\frac{\mathrm{d}V}{\mathrm{d}z}P(C_R\mid M,z)\mathrm{d}M\mathrm{d}{C}_R\mathrm{d}z\end{array}$$(8)

where $P\left({\hat{\mathcal{H}}}_i\right)$ is the uniform probability over the sphere that a cluster exists at a specific sky position ${\hat{\mathcal{H}}}_i,P(\hat{z}\mid z)$ expresses the uncertainty on our redshift measurements, provided by eROMaPPer (Kluge et al. 2024), $P(I\mid z,C_R,{\hat{\mathcal{H}}}_i)$ is the X-ray selection function (see Sec. 3 and for more details Clerc et al. 2024), $P({\hat{C}}_R\mid C_R)$ is the uncertainty on the X-ray count rate, measured by MBProj2D, $\frac{\mathrm{d}N}{\mathrm{d}M\mathrm{d}V}$ is the halo mass function¹⁰, which directly depends on cosmology, $\frac{\mathrm{d}V}{\mathrm{d}z}$ is the fractional comoving volume in an infinitesimal redshift shell (see Equation (28) in Hogg 1999)¹¹, P(C_R|M, z) is the log-normal distribution around the X-ray scaling relation and can be described as: $$P(C_R\mid M,z)=\mathcal{L}\mathcal{N}(⟨C_R\mid M,z⟩,{\sigma }_X).$$(9)

We refer to Sec. 6 for details on the parameterized scaling relation. Furthermore, the total number of objects detected in eRASS1, depends on the calculation of the number of clusters and on the fraction of contaminants, as shown in Equation (A.2). The total number of detected clusters is simply the integration over the observable space of the differential Poisson term, and can be computed via the equation: $$N_{\mathrm{C}}=⨌\frac{\mathrm{d}{N}_{\mathrm{C}}}{\mathrm{d}{\hat{C}}_R\mathrm{d}\hat{z}\mathrm{d}\hat{\lambda }\mathrm{d}{\hat{\mathcal{H}}}_i}P(I\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i)\mathrm{\Theta }(\hat{\lambda }>3)\mathrm{d}{\hat{C}}_R\mathrm{d}\hat{z}\mathrm{d}\hat{\lambda }\mathrm{d}{\hat{\mathcal{H}}}_i$$(10)

where several simplifications are applied. First, $P(\hat{z}\mid z)$, the probability of observing a cluster with measured redshift $\hat{z}$ given the true redshift z, is modeled as a mixture of two normal distributions whose parameters are calibrated by eROMaPPer, expressed as: $$\begin{array}{rl}P(\hat{z}\mid z)=& (1-c_z)\mathcal{N}(b_{z}z,{\sigma }_z(1+z))\\ & +c_z\mathcal{N}(b_{z}z+c_{\text{shift },z},{\sigma }_z(1+z)),\end{array}$$(11)

where c_z is the fraction of objects for which we measure a redshift shifted by some constant c_shift,z, b_z is a systematic bias in our redshift estimate, and σ_z is the relative error on the measured redshift. We can, therefore, express $g(z,{\mathcal{H}}_i)=$ $\int P(\hat{z}\mid z)w(\hat{z},{\mathcal{H}}_i)\mathrm{d}\hat{z}$ as an analytic expression for each sky position, with $w(\hat{z},{\mathcal{H}}_i)$ being the filter function that is redshift and sky position-dependent, as follows: $$w(\hat{z},{\mathcal{H}}_i)=\{\begin{array}{ll}1,& \text{ if }(\hat{z}>0.1)\mathrm{\&}(\hat{z}<min(0.8,{\hat{z}}_{max}\left({\hat{\mathcal{H}}}_i\right)),\\ 0,& \text{ otherwise },\end{array}$$(12)

where ${\hat{z}}_{\text{max}}\left({\hat{\mathcal{H}}}_i\right)$ is the maximum measurable redshift at sky position ${\hat{\mathcal{H}}}_i$, given by the ZVLIM map (see Kluge et al. 2024, for further details).

Then, $P(\hat{\lambda }\mid \lambda )$, the probability of observing a cluster with measured richness $\hat{\lambda }$ given the true richness λ, is modeled as in previous cosmological analyses with cluster counts (Dark Energy Survey Collaboration 2020; Chiu et al. 2023; Bocquet et al. 2024) as a Poisson distribution in its Gaussian limit. A Poisson distribution is a reasonable theoretical assumption as richness effectively counts galaxy members, and partly demonstrated in Costanzi et al. (2019a). Therefore, we have: $$P(\hat{\lambda }\mid \lambda )=\mathcal{N}(\lambda ,\sqrt{\lambda })\text{. }$$(13)

Thus, we can express the optical completeness term $f(\lambda )=$ $\int P(\hat{\lambda }\mid \lambda )\mathrm{\Theta }(\hat{\lambda }>3)\mathrm{d}\hat{\lambda }$ term analytically¹². The effects of the optical cut, $\hat{\lambda }>3$, on the overall sample completeness is very small, decreasing the overall number of clusters computed from Equation (10), and thus the completeness by 1.1%. Furthermore, photometric redshifts from optical structures with 3 or fewer member galaxies are arguably quite unreliable.

Finally, since ${\hat{C}}_R$ appears only in the term $P({\hat{C}}_R\mid C_R)$, and since probabilities are normalized, the integral over measured count rate is unity and does not affect the calculation of the total number of objects. Therefore, the total number of clusters is expressed as: $$\begin{array}{rl}{N}_C=& ⨌P(I\mid z,C_R,{\hat{\mathcal{H}}}_i)g(z,{\hat{\mathcal{H}}}_i)f(\lambda )P(C_R,\lambda \mid M,z)P\left({\hat{\mathcal{H}}}_i\right)\\ & \times \frac{\mathrm{d}N}{\mathrm{d}M\mathrm{d}V}\frac{\mathrm{d}V}{\mathrm{d}z}\mathrm{d}M\mathrm{d}{C}_R\mathrm{d}z\mathrm{d}\lambda \mathrm{d}{\hat{\mathcal{H}}}_i.\end{array}$$(14)

To conclude, we now have expressed explicitly all the ingredients required in order to model the cluster number counts likelihood in Equation (8) and to compute the total number of clusters that are detected in Equation (14).

4.4 Weak-lensing mass calibration likelihood

This subsection focuses on the weak lensing likelihood, which is crucial for mass calibration. We make use of Bayes’ theorem to express the weak lensing likelihood as: $$P({\hat{g}}_t\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i,I,C)=\frac{P({\hat{g}}_t,{\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i,I,C)}{P({\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i,I,C)}=\frac{N_{\mathrm{W}\mathrm{L}}}{D_{\mathrm{W}\mathrm{L}}},$$(15)

where we have factored out the numerator N_WL and the denominator D_WL; otherwise, the expression will not fit a double column article. Similarly to what we have done in Sec. 4.3, we marginalize the expressions N_WL and D_WL, similarly to Equation (8), as follows: $$\begin{array}{rl}{N}_{\mathrm{W}\mathrm{L}}=& P\left({\hat{\mathcal{H}}}_i\right)\iiint P(I\mid C_R,\hat{z},{\hat{\mathcal{H}}}_i)P\left({\hat{C}}_R\mid C_R\right)P(M,\hat{z})\\ & \times P({\hat{g}}_t\mid M_{\mathrm{W}\mathrm{L}},\hat{z})P(M_{\mathrm{W}\mathrm{L}},C_R\mid M,\hat{z})\mathrm{d}M\mathrm{d}{C}_R\mathrm{d}{M}_{\mathrm{W}\mathrm{L}}\end{array}$$(16)

and $$\begin{array}{rl}{D}_{\mathrm{W}\mathrm{L}}=& P\left({\hat{\mathcal{H}}}_i\right)\iint P(I\mid C_R,\hat{z},{\hat{\mathcal{H}}}_i)P\left({\hat{C}}_R\mid C_R\right)P(M,\hat{z})\\ & \times P(C_R\mid M,\hat{z})\mathrm{d}M\mathrm{d}{C}_R.\end{array}$$(17)

We note that the sky position probabilities, $P\left({\hat{\mathcal{H}}}_i\right)$, cancel out in the ratio of N_WL / D_WL in the weak lensing mass calibration term. $P({\hat{g}}_t\mid M_{\mathrm{WL}},\hat{z})$ is the likelihood term that compares the measured weak lensing shear with the one predicted from the weak lensing mass (Grandis et al. 2024). Furthermore, in both the numerator and denominator, we have the error on the measured count rate $P({\hat{C}}_R\mid C_R)$, the selection function $P(I\mid C_R,\hat{z},{\hat{\mathcal{H}}}_i)$, the normalized halo mass function, $P(M,\hat{z})$. However, since the normalization factor (needed to compute $P(M,z)$ from $\frac{\mathrm{d}N}{\mathrm{d}M\mathrm{d}V}\frac{\mathrm{d}V}{\mathrm{d}z}$) is present both at the numerator and denominator, it cancels itself out, and the halo mass function, $\frac{\mathrm{d}N}{\mathrm{d}M\mathrm{d}V}\frac{\mathrm{d}V}{\mathrm{d}\hat{z}}$, can be used directly in the previous equations. All the details regarding the expressions for these terms are described just below Equation (8).

Finally, $P(M_{\mathrm{WL}},C_R\mid M,\hat{z})$, the scaling relation term between mass and both count rate and weak lensing mass that takes into account the intrinsic correlation between these quantites, is a bi-variate log-normal distribution $\mathcal{B}\mathcal{N}(\overline{\mu },\overline{\mathrm{\Sigma }})$, where $$\overline{\mu }=[⟨\log C_R\mid M,z⟩,⟨\log M_{\mathrm{W}\mathrm{L}}\mid M,z⟩],$$(18) $$\overline{\mathrm{\Sigma }}=\left[\begin{array}{cc}{\sigma }_{C_R}^2& {\rho }_{M_{\mathrm{W}\mathrm{L}},C_R}{\sigma }_{C_R}{\sigma }_{M_{\mathrm{W}\mathrm{L}}}\\ {\rho }_{M_{\mathrm{W}\mathrm{L}},C_R}{\sigma }_{C_R}{\sigma }_{M_{\mathrm{W}\mathrm{L}}}& {\sigma }_{M_{\mathrm{W}\mathrm{L}}}^2\end{array}\right],$$(19)

where ⟨log C_R|M, z⟩ and ⟨log M_WL|M, z⟩ are the scaling relations between true cluster mass and count rate and weak lensing mass, respectively, see Sec. 6 for further details, and ρM_WL,C_R is the intrinsic correlation between weak lensing mass and count rate.

4.5 Optical likelihood

Similarly to the previous subsection, we apply Bayes’s theorem to express the optical likelihood, meaning the probability of measuring the observed richness for a cluster that has an observed count rate and redshift at a specific sky position, as: $$P(\hat{\lambda }\mid {\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i,I,C)=\frac{P(\hat{\lambda },{\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i,I,C)}{P({\hat{C}}_R,\hat{z},{\hat{\mathcal{H}}}_i,I,C)}=\frac{N_{\mathrm{o}\mathrm{p}\mathrm{t}}}{D_{\mathrm{o}\mathrm{p}\mathrm{t}}},$$(20)

where (as before) we have factored out the numerator N_opt and the denominator D_opt. As before, we expand the numerator and denominator as: $$\begin{array}{rl}{N}_{\mathrm{o}\mathrm{p}\mathrm{t}}=& P\left({\hat{\mathcal{H}}}_i\right)\iiint P(I\mid C_R,\hat{z},{\hat{\mathcal{H}}}_i)P\left({\hat{C}}_R\mid C_R\right)P(M,\hat{z})\\ & \times P(\hat{\lambda }\mid \lambda )\mathrm{\Theta }(\hat{\lambda }>3)P(\lambda ,C_R\mid M,\hat{z})\mathrm{d}M\mathrm{d}{C}_R\mathrm{d}\lambda ,\end{array}$$(21) $$\begin{array}{rl}{D}_{\mathrm{o}\mathrm{p}\mathrm{t}}=& P\left({\hat{\mathcal{H}}}_i\right)\iint P(I\mid C_R,\hat{z},{\hat{\mathcal{H}}}_i)P\left({\hat{C}}_R\mid C_R\right)P(M,\hat{z})\\ & \times P(C_R\mid M,\hat{z})\mathrm{d}M\mathrm{d}{C}_R.\end{array}$$(22)

Here, in particular, $P(\hat{\lambda }\mid \lambda )$ is the likelihood term that evaluates the likelihood of observing a richness given the true richness, given by the measurement uncertainty on the richness given by eROMapper. Similarly to the weak lensing likelihood, we have that $P(I\mid C_R,\hat{z},{\hat{\mathcal{H}}}_i)$ is the X-ray selection function, $\mathrm{\Theta }(\hat{\lambda }>3)$ is the optical selection function¹³, $,P(M,\hat{z})$ is the normalized mass function, $P({\hat{C}}_R\mid C_R)$ is the error on the measurement of the count rate, $P(\lambda ,C_R\mid M,\hat{z})$ is a bivariate log-normal distribution $\mathcal{M}\mathcal{N}(\overline{\mu },\overline{\mathrm{\Sigma }})$ where $$\overline{\mu }=[⟨\log C_R\mid M,\hat{z}⟩,⟨\log \lambda \mid M,z⟩],$$(23) $$\overline{\mathrm{\Sigma }}=\left[\begin{array}{cc}{\sigma }_{C_R}^2& {\rho }_{\lambda ,C_R}{\sigma }_{C_R}{\sigma }_{\lambda }\\ {\rho }_{\lambda ,C_R}{\sigma }_{C_R}{\sigma }_{\lambda }& {\sigma }_{\lambda }^2+{\sigma }_P^2\end{array}\right],$$(24)

where ⟨log C_R|M, z⟩ and ⟨log λ|M, z⟩ are the scaling relations between true mass, count rate, and optical richness, details are given in Sec. 6, while ρ_λ,C_R is the intrinsic correlation between richness and count rate and where σ_λ is richness intrinsic scatter, and is a Poisson term, defined as in Costanzi et al. (2019a) by $\frac{⟨\lambda \mid M,z⟩-1}{⟨\lambda \mid M,z{⟩}^2}$. We point out that the denominator for the optical, D_opt, and weak lensing, D_WL, likelihoods look identical; however, we recall that while the optical term runs over all the objects in the sample, the weak-lensing term runs over only the objects that have tangential shear measurements.

4.6 Parameters and priors

In Table 1, we list all the parameters fitted in our Bayesian statistical model. We include their units, descriptions, and prior functions used. We note that our analysis is not sensitive to some cosmological parameters. For this reason, to avoid the chains reaching unphysical values but still allowing marginalization over these parameters, we limit their range using realistic priors from previous cosmological experiments. Specifically, we limit Ω_b in a range allowed by the Big Bang Nucleosynthesis (Beringer et al. 2012; Cooke et al. 2016, 2018; Mossa et al. 2020; Pitrou et al. 2021), n_s in a range round Planck CMB results (Planck Collaboration VI 2020). For H₀, against which we are not sensitive, we use the Planck CMB ΛCDM results with a tight Gaussian prior set around them (Planck Collaboration VI 2020). Given the so-called H₀ tension, we have checked that using a prior from distance ladder method (H₀ ~ $\mathcal{N}$(73.04, 1.04) km s⁻¹ Mpc⁻¹, Riess et al. 2022) or a flat prior (H₀ ~ $\mathcal{U}$(60, 80) km s⁻¹ Mpc⁻¹) does not affect our results and conclusions.

5 Validation of the cosmology pipeline and the blinding strategy

To ensure the results from the run of the cosmology pipeline on the eROSITA cosmology subsample are robust against any additional unaccounted systematics, we extensively test our pipeline against the realistic simulations of the eROSITA’s digital twin (Seppi et al. 2022) and the mock observations that reflect the characteristics of our modeling. We perform our analysis blindly to avoid any unintentional unconscious bias toward the already published results in the literature on cosmology. Specifically, we investigate the influence of uncertainties on photometric redshifts, selection function, WL miscentering, other WL-related systematic effects, and contamination modeling on our results. This section describes in detail our testing of the cosmology pipeline, assessment of potential biases, and our blinding strategy.

5.1 Mocks generation and validation of the pipeline

As cosmology studies with cluster counts rely on specific modeling assumptions, such as modeling of the mass bias, selection function, and redshift uncertainty, we tested the robustness of the cosmological framework against the parameters used in this analysis. To verify the cosmology framework we made the cosmology pipeline and theoretical modeling available to the collaboration members for an independent review of each component.

The validation of the cosmological framework is performed independently. We developed a new pipeline for generating mock catalogs following the physical formalism (described in Sects. 4 and 6), without relying on the eROSITA cosmology pipeline. We ensure that the observables generated by these mocks follow the distribution of the eROSITA X-ray count rate and optical information from LS DR10-South, and weak lensing follow-up data sets from the DES, KiDS, and HSC surveys. Then the cosmology pipeline is tested against these mocks to make sure the pipeline reproduces the input parameters in the assumed theoretical framework. This approach allows us to independently test the statistical robustness of all the contributions to the total likelihood, shown in Equation (6).

We generated mock cluster catalogs with multi-wavelength information using the dark matter halo mass function package hmf (Murray et al. 2013). In the mocks, we utilized halos with M₅₀₀ > 10¹²M_⊙, and distributed them homogeneously across the celestial sphere to reflect the eROSITA clusters sample. At the position of each halo, we assigned a local background (Freyberg et al. 2020), absorption from hydrogen (HI4PI Collaboration 2016; Willingale et al. 2013), and exposure time (Predehl et al. 2021) by interpolating the corresponding maps consistent with the survey strategy. To compute the expected observables, such as the count rate and richness, we employed the scaling relations described in Sec. 6 with fiducial values from our fits. Observable quantities were then drawn according to the expected statistical distribution. As described in Sec. 3.3, we computed the detection probability for each halo and select edclusters accordingly based on eROSITA’s selection function. Finally, we tested the components of the global cosmology likelihood to ensure that systematics trends in the data are recovered well and estimated accurately by the pipeline.

The dark matter surface mass density was calculated by assuming a spherical NFW profile (Navarro et al. 1996) for each cluster. We computed the mean expected shear by integrating this value over the background galaxy distribution from the respective survey. The observed shear was deduced by adding errors similar to those observed in the DES shear catalog(see Sect. 2.3.1 and Grandis et al. 2024). Out of these shear simulations, we selected only those that fall into the optical survey footprint of DES, HSC, and KiDS, while the corresponding shear profiles are assigned to these objects. Finally, by sampling their probability distribution, we included the contaminants, RSs and AGNs (see Sec. 4.2) and used the contamination fractions from the best fit.

We generated ten different mocks from the same model and used the cosmological pipeline to recover the input parameters, σ₈ and Ω_m. Most of the output parameter values fall well within our 68% confidence level of the input values. We note that the rest of the parameters also fall well within our 68% confidence level of the input values; in particular, we also recovered the input number of contamination level accurately in these tests. This validation step ensures that the components of the global likelihood, mass calibration, cluster counts, and contamination modeling work as expected free of additional systematics. As a result of this test, we demonstrate that the pipeline is free of bugs and passes the statistical robustness test successfully.

To test the effects of the selection on our results, we applied realistic changes to our selection function. We generated six selection functions¹⁴ described in Sec. 3.4. These models differ by the thresholds used to construct the simulation catalog. We then generated six mocks, using a different selection function each time, freezing all the other components in the likelihood. We then ran the cosmology pipeline with different selection functions to investigate the impact of the uncertainties on the selection function on our results. We find that the contours shift slightly but always within the 68% percentiles of the input of the σ₈ and Ω_m parameters.

Recent results by Sommer et al. (2024) have suggested that the underestimated miscentering measurements between ICM centers and halo centroids indicated by the numerical simulations may introduce an additional systematic bias on WL mass measurements by 2–6%. However, the exact extent of the bias remains to be investigated and depends on the assumed X-ray miscentering in hydrodynamical simulations. Although the uncertainties in miscentering distribution calibrated on hydro-dynamical simulations are fully accounted for in our WL shear measurements (see Sec. 6.2 and Grandis et al. 2024), we further investigated the effect of an additional miscentering bias and its influence on our measurements by introducing a bias at the levels as suggested by Sommer et al. (2024) in Equation (28). We find that a decrease in the masses as found by Sommer et al. (2024) indeed decreases the values of Ω_m and σ₈ as anticipated; however, the best-fit values remain statistically consistent within the 68% contour level of the input values for both parameters.

The redshifts in our catalog were measured using the photometric information from the LS DR10-South grz bands; see Sec. 2.1 (Kluge et al. 2024). This means that their uncertainties are significantly larger than spectroscopic redshift measurements. This potentially could cause a bias on the cosmological parameters mainly as a result of geometrical effects. This is because $\frac{\mathrm{d}N}{\mathrm{d}M\mathrm{d}V}$ is a weak function of redshift, especially in the redshift interval of interest, but the volume $\frac{\mathrm{d}V}{\mathrm{d}z}$ and cosmological distances change drastically. We tested this effect for various parameters in Equation (11) when generating mocks. We find that the parameters that would result in significant bias in Ω_m and σ₈ measurements are unrealistic. For example, a systematic error of 10% on photo-zs, namely, σ_z ≈ 0.1, affects Ω_m and σ₈ by about 3 and 2% respectively, and a trend of catastrophic failures on photo-z with occurrence of 10% results in lower values of Ω_m by 4% and higher values of σ₈ by 2%. The comparisons of the photometric measurements against spectroscopic redshifts in the redshift range 0.1 to 0.8, given in Kluge et al. (2024), allows us to assess the actual values for the parameters in Equation (11). The intrinsic scatter of the photometric redshift with respect to the spectroscopic redshift is <0.5%, on average, and the fraction of significant catastrophic failures is just 0.13%; thus, it is much smaller than the extreme parameters required to introduce a significant bias. Therefore, the bias introduced by the use of photometric redshifts on Ω_m and σ₈ parameters is negligible.

Given the positive results from the robustness tests and verification of the contamination modeling, assessment of the WL mass, and redshift bias on the cosmological constraints, we do not anticipate any significant bug or unaccounted bias in our pipeline and methodology. Therefore, we conclude that the cosmology pipeline and the methods have been validated successfully.

5.2 Blinding strategy

Blinding the results prior to publication minimizes preconceived assumptions and prevents us from introducing bias into our analysis. After the validation phase was successfully completed on our simulations and mock data, we froze the cosmology pipeline and apply the following blinding strategy to ensure that our conclusions would remain as objective and unbiased as possible. We generated six different selection function models, as presented in Sec. 3.4, only one corresponding to the true selection model, namely, with the ${\mathcal{L}}_{ext}$ thresholds that match the construction of the cosmology subsample. These selection functions were tested on the mocks described in Sec. 5.1 to ensure that they do not display any unphysical behavior that will lead to unintentional blinding. A total of six selection functions were then sent to our external blinding expert, Professor Catherine Heymans. The eROSITA-DE cluster working group received four of these six selection models after the blinding lead shuffled their indices. During this process, the “true” selection function model and its index became unknown to the working group and analysis teams. The cosmology analysis for the flat ΛCDM model was carried out using the four independent selection functions simultaneously without any prior knowledge of the “true” selection function. This strategy resulted in four independent cosmology analyses and contours. A draft of the paper was prepared on the four distinct cosmology results from the ΛCDM analysis and circulated for the working group’s review. During this phase, we satisfied the tests described in the validation phase; the cosmology pipeline was extensively tested against any bugs or biases and completed the internal checks on the different components of the analysis (see Sec. 5.1). We then unblinded the results and completed the rest of the draft of the paper, extending our analysis beyond the ΛCDM model without introducing any modifications to the analysis method or results after unblinding.

6 Scaling relations

This section describes our formalism for the scaling relations between mass and our observables accounting for the redshift evolution. This step is vital for constructing the mass function from the clusters detected in the eRASS1 survey with X-ray count-rate, optical richness, redshift, and shear profile information. We also compare our findings and the reported best-fit values in the literature from similar studies.

6.1 X-ray Count Rate – Mass Scaling Relation

The general form of the scaling relation between the X-ray observable C_R and cluster mass, required in Equation (9), can be described with the following equation: $$⟨\log \frac{C_R}{C_{R,p}}\mid M,z⟩=\log A_X+b_X(z)\log \frac{M}{M_p}+e_x(z).$$(25)

We used the $C_{R,p}=0.1\mathrm{cts}/\mathrm{s}$ as the pivot value for the count rate, M_p = 2 × 10¹⁴ M_⊙, as the pivot value for the mass, and z_p = 0.35 as the pivot value for the redshift. The redshift dependent slope of the scaling relation is expressed by the term b_X(z), and the redshift evolution of the scaling relation is expressed by e_x(z). Explicitly, the former is expressed as: $$b_X(z)=B_X+F_X\log \frac{1+z}{1+z_p}.$$(26)

Here, B_X is the normalization, and F_X allows the slope to be redshift dependent. The self-similar model of Kaiser (1986) predicts the F_X slope to be zero. In this work, we set the prior for F_X parameter centered around zero to avoid any potential bias in our slopes. The redshift evolution term, e_x(z), is given by: $$e_x(z)=D_X\log \frac{d_L(z)}{d_L\left(z_p\right)}+E_X\log \frac{E(z)}{E\left(z_p\right)}+G_X\log \frac{1+z}{1+z_p},$$(27)

for which E_X = 2 and D_X = −2 are the values adopted by the self-similar model (Kaiser 1986) for the case of luminosity; our n_H-corrected count rate has in principle the same dependence except for the K-correction, which is captured by the G_X parameter. The last term, G_X, quantifies the deviation from the self-similar model in redshift evolution (e.g., as used in Bulbul et al. 2019). The intrinsic scatter of our scaling relations is represented by a single value, σ_X, independent of mass and redshift, similar to the studies in the literature on scaling relations from X-ray observations, (for example, Mantz et al. 2016; Schellenberger & Reiprich 2017a; Lovisari et al. 2020, 2021; Duffy et al. 2022; Giles et al. 2022), where the limited number of clusters and mass range does not allow to constrain mass or redshift trend for the scatter. This assumption is also adopted in the cosmological analyses utilizing galaxy cluster number counts (Mantz et al. 2015; Schellenberger & Reiprich 2017b; Bocquet et al. 2019; Costanzi et al. 2021; Bocquet et al. 2024). Additionally, the X-ray luminosity’s scatter as a function of mass has been studied for the mass range similar to the sample used in this work in IllustrisTNG numerical simulations of Pop et al. (2022); see bottom panel of their Figure 6. The authors show that the scatter shows a subdominant mass trend. As we extend the mass and redshift range of the eROSITA samples in the deeper surveys, the proper modeling of the intrinsic scatter in cosmology analyses becomes important. The best-fit values of the variables in the C_R − M scaling relations are provided in Table 2, and their posteriors are shown in Fig. 5.

The only other works in the literature that use eROSITA X-ray count-rate in mass scaling relations are the ones by Chiu et al. (2022, 2023) on the eFEDS sample. Our best-fit normalization parameter is $A_x={0.64}_{-0.06}^{+0.04}$. The reported normalization parameters are $A_x={1.24}_{-0.19}^{+0.22}$ and $A_x={1.33}_{-0.20}^{+0.26}$ in Chiu et al. (2022, 2023) in ~3σ tension with our results. In particular, the difference in normalization is due to the observed degeneracy with the scatter, see Fig. B.1, the selection of the sample (Bahar et al. 2022), and the selection function parameters. In fact, in Chiu et al. (2022, 2023), the authors have significantly simplified the selection function, particularly in its sky variation. eFEDS covers an area of 140 deg², albeit much smaller in area compared to eRASS1, the eFEDS area is not small enough such that fundamental sky properties affecting the X-ray selection function can be safely ignored (Brunner et al. 2022).

We find that the slope of the power-law index of the mass term is $B_x={1.38}_{-0.04}^{+0.03}$. Overall, in terms of mass trends, we find a steeper slope than the self-similar mass trend (B_x,ss = 1), while many scaling relation studies in the literature report even steeper slopes. We find shallower slope than the reported value of $B_x={1.58}_{-0.14}^{+0.17}$ in Chiu et al. (2022) and that of $B_x={1.86}_{-0.15}^{+0.20}$ in the count rate mass relations provided in Chiu et al. (2023). The discrepancy is likely due to different treatments of the selection function between these two sets of results. In this work, we accurately model our selection with the one-to-one simulations of the eROSITA X-ray sky, while in Chiu et al. (2023), an empirical model for the selection is used. In principle, as the X-ray count rate is closely related to the X-ray luminosity, the mass slope of the L_X − M scaling relations reported in the literature can be compared to our slope measurements. The slopes reported in the literature range between 1.3 and 1.9. The studies resulting in shallow slopes employ X-ray hydrostatic masses $B_{L_x}=$ $1.45\pm 0.07$ in Maughan et al. (2007), 1.61 ± 0.14 in Vikhlinin et al. (2009), 1.62 ± 0.12 in Pratt et al. (2009), 1.34 ± 0.18 in Eckmiller et al. (2011), 1.66 ± 0.22 in Lovisari et al. (2015), and 1.35 ± 0.07 in Schellenberger & Reiprich (2017b). Our slopes are consistent with these findings at the 1σ level. Higher slopes are reported in SZ-selected massive cluster samples covering a wider redshift range (B_Lx = 1.89 ± 0.18, Bulbul et al. 2019). When the authors excised the core, they found shallower slopes ${1.60}_{-0.16}^{+0.15}$, consistent with our finding within uncertainties. Another sample with a higher reported mass slope is the CODEX with the best-fit value of ${2.01}_{-0.09}^{+0.09}$ (Capasso et al. 2020).

The best-fitting scatter ${\sigma }_X,{0.98}_{-0.04}^{+0.05}$, in our sample is larger than the reported values from C_R − M or L_X-M scaling in the literature, e.g., 0.17 and 0.25 in the X-ray selected REXCESS (Pratt et al. 2009), HIFLUGCS samples (Lovisari et al. 2015), and 0.27 in the SZ selected samples (Bulbul et al. 2019). The log-normal scatter of the soft-band luminosity is constrained to be 0.12 in Chiu et al. (2022) for eFEDS clusters. We also note that the scatter in the C-Eagle cosmological hydrodynamical simulations of clusters is 0.30, which is also smaller than our measurements (Barnes et al. 2017). This is partially because our mass range and redshift are more extensive than those employed in previous studies. We do not apply any cuts on the cluster morphology and do not excise core emissions. This strategy might lead to a larger spread than previously observed in scaling relations. However, we note that this large scatter has a negligible impact on the accuracy of the cosmological parameters (see Sec. 7).

Fig. 5 Posteriors on the CR − M and λ − M scaling relation parameters. A strong correlation between X-ray and optical scaling relation is visible on the corner plots. The best-fit values and their uncertainties are provided in Table 2. The solid and shaded orange contour represents the 68 and 95% contour levels respectively.

6.2 Weak-lensing inferred mass-mass scaling relation

Measuring the tangential shear of galaxies lying behind galaxy clusters is a powerful tool for estimating the cluster mass. As with all mass proxies, the weak-lensing inferred mass is also a biased mass proxy; however, the factors contributing to the bias are under control and can be accurately calibrated using simulations with knowledge of the survey-specific background source distribution. The general form of the scaling relation that relates the weak-lensing inferred mass with the true mass (see Grandis et al. 2021a) required in Equation (18), is described by the following equation: $$⟨\log \frac{M_{\mathrm{W}\mathrm{L}}}{M_p}\mid M,z⟩=b(z)+b_M\log \left(\frac{M}{M_p}\right).$$(28)

Here, b(z) is the redshift-dependent evolution of the weak-lensing bias as a function of redshift. It is computed by interpolating simulations snapshots at redshifts: 0.24, 0.42, 0.68, and 0.95.

The calibration approach is first presented in Grandis et al. (2021a), and consists of Monte Carlo simulations, where each single iteration consists of three steps: (1) starting from surface mass density maps of massive halos from the cosmological hydro-dynamical TNG300 simulations (Pillepich et al. 2018; Marinacci et al. 2018; Springel et al. 2018; Nelson et al. 2018, 2019; Naiman et al. 2018), and the properties of the cluster sample and the WL surveys, we predict a catalog with realistic synthetic shear profiles; (2) we use the same expression $P({\hat{g}}_t\mid M_{\mathrm{WL}},\hat{z})$ described in Sec. 4.4, which depends on the dataset used, in our case DES, KiDS, and HSC, evaluated on the synthetic shear profiles, to assign WL masses in the simulation, and (3) we finally fit the relation given above between WL mass and halo mass, recovering a WL bias values b_l for each snapshot, and a global mass trend b_M.

As introduced in Grandis et al. (2021a), the true novelty of this approach is that the simulation above can be repeated order 1000 times while varying the inputs within their uncertainties, thus effectively propagating all WL related systematic uncertainties to a distribution on the WL bias values b_l and b_M. To compress this information and ingest it into the cosmological pipeline as priors, we derive the mean WL bias μ_l. The scatter of the Monte Carlo realizations around this mean is well captured by two principal components PC_1,2,l (Grandis et al. 2024). For the mass trend we report mean ${\mu }_{b_{M,}}$ and standard deviation ${\sigma }_{b_{M,i}}$. We similarly fit the scatter in the Wl mass, finding that its redshift trend is well described by one principal component.

Since we use the weak lensing shear in our mass calibration from the three independent surveys, DES, KiDS, HSC, we calibrate the count rate to mass scaling relations jointly. We note that the parameters b(z) and b_M are survey dependent since the specific properties of each survey affect the reconstructed weak lensing mass slightly differently, see Grandis et al. (2024) and Kleinebreil et al. (2024). However, since only the mean relation and uncertainty around it are survey-dependent, we express b(z) using the set of parameters ${\mu }_{b,i}(z),P{C}_{1,i}(z),P{C}_{2,i}(z),{\mu }_{b_{M,i}}$, and with i ∈ {DES, KiDS, HSC}. To emulate the scatter in the WL bias simulations, we introduce three additional random normal parameters, A_WL, B_WL, and C_WL, which are independent of the surveys employed. In practice, b(z) and b_M can be expressed as: $$\begin{array}{rl}b(z)& ={\mu }_{b,i}(z)+A_{\mathrm{W}\mathrm{L}}P{C}_{1,i}(z)+B_{\mathrm{W}\mathrm{L}}P{C}_{2,i}(z)\\ b_M& ={\mu }_{b_{M,i}}+C_{\mathrm{W}\mathrm{L}}{\sigma }_{b_{M,i}},\end{array}$$(29)

where μ_b, i(z), PC_1,i(z), and PC_2,i(z), the survey-dependent average relation. The two principal components, PC_1,i(z), and PC_2,i(z), are the values computed in the simulations used for interpolation at the specific cluster redshift. For b_M, no interpolation is required, but written in this form allows us to have only one free parameter for the three surveys, as the three surveys have slightly different mass trends in their WL biases.

As mentioned above, in the same simulation setup, we can also derive the redshift and mass trend of the scatter in WL mass, which is given by: $$\log {\sigma }_{M_{\mathrm{W}\mathrm{L}}}^2=s(z)+s_M\log \left(\frac{M}{M_p}\right),$$(30)

where the simulations have shown that the mass component of the scatter is consistent with zero; we fix s_M = 0 in further analysis. For the redshift evolution of the scatter, we follow a similar procedure for the redshift evolution of the bias b(z). As before, we take the values at simulation snapshots and interpolate. We use the survey-specific mean value and principal component μ_s,i(z) and PC_i(z) and a random normal parameter D_WL and express s(z) as: $$s(z)={\mu }_{s,i}(z)+D_{\mathrm{W}\mathrm{L}}P{C}_i(z),$$(31)

where the average scatter for each survey and principal component (μ_s,i(z) and PC_i(z)) are computed in the simulation snapshots and used in interpolation. This formalism allows distinguishing the average survey-dependent components from the simulation-dependent scatter around them, dominated by simulation uncertainty. It leaves us with four free parameters related to WL to fully describe the systematic and statistical uncertainty of the bias.

This method effectively assumes that the WL mass measurement systematics are correlated between the three surveys. This is true for the low redshift cases, where the hydrodynamical modeling uncertainties represent almost the entirety of the systematic uncertainty on the estimated weak lensing masses. At higher redshifts, each survey is dominated by a different systematic: DES by the photometric source redshift uncertainty (see Grandis et al. 2024, Sect. 4.5), HSC by the shape measurement uncertainty (see Kleinebreil et al. 2024), and KiDS by the cluster member contamination model (see Kleinebreil et al. 2024, Appendix A). All these systematics might be partially correlated. Assuming that the WL mass measurement systematics are correlated effectively increases the systematics error budget. Thus, our analysis strategy provides a conservative upper bound on the systematics.

6.3 Richness-mass scaling relation

Although the X-ray count rate is our main observable, richness measurements of the detected sources are required to distinguish the contamination (from an AGN or RS) from the cluster number counts. Similarly to previous studies in the literature, optical confirmation and richness determination are useful for assessing the contamination level of the sample and improving its purity to produce state-of-the-art cosmological analysis (Bocquet et al. 2024). In Kluge et al. (2024), we develop a concept of a redshift-dependent richness cut by combining optical and X-ray properties of the ICM to statistically model both clusters and contamination simultaneously. This method comes at the cost of introducing richness as another observable but with the benefit of self-consistent modeling of clusters and contaminants in the sample. In the likelihood term in Equation (6), the relation between richness and count rate is explicitly expressed in the third term. We modify the optical richness to mass scaling relation in Saro et al. (2015); Bleem et al. (2020); Grandis et al. (2020, 2021b), which we require as Equation (23), as follows: $$⟨\log \lambda \mid M,z⟩=\log A_{\lambda }+b_{\lambda }(z)\log \left(\frac{M}{M_p}\right)+C_{\lambda }\log \left(\frac{1+z}{1+z_p}\right),$$(32)

while the slope of the cluster mass term b_λ(z) can be expanded as: $$b_{\lambda }(z)=B_{\lambda }+D_{\lambda }\log \left(\frac{1+z}{1+z_p}\right),$$(33)

where A_λ is the normalization of the scaling relation, B_λ is the mass slope, C_λ is the redshift evolution, and D_λ is the redshift evolution of the mass slope.

In this formulation, marginalizing over mass allows us to account for the correlation between count rate and richness and their uncertainties. It is crucial to calibrate the scaling relation between richness and count rate well to achieve precise and accurate cosmological results. However, we expect from our mathematical framework that the best-fit parameters of the mass-to-richness scaling relation and the parameters of the mass-to-count rate scaling relation degenerate with each other.

The best-fit parameters of the richness-M scaling relations are provided in Table 2, and their posteriors are shown in Fig. 5. In the case of the richness-to-mass scaling relation, Equation (32), the best-fit value of the mass slope, B_λ ${1.01}_{-0.03}^{+0.02}$ is consistent with the majority of reported findings in the literature within 1σ level, e.g., $B_{\lambda }={0.99}_{-0.07}^{+0.06}\pm 0.04$ and $B_{\lambda }={1.07}_{-0.20}^{+0.22}$ in the CODEX sample (Capasso et al. 2019; Ider Chitham et al. 2020, respectively), B_λ = 1.02 ± 0.08 in SPT sample (Bleem et al. 2020), $B_{\lambda }={1.028}_{-0.037}^{+0.043}\pm 0.04$ and B_λ = 1.03 ± 0.10 in the DES Y1 and SPT samples Costanzi et al. (2021); Grandis et al. (2021b), $B_{\lambda }={1.05}_{-0.12}^{+0.13}$ in Chiu et al. (2023) in the eFEDS sample within the stated uncertainities. Our measurements are 8σ away from the value reported by the DES collaboration on Y1 cosmology analysis, where they find B_λ = 0.737 ± 0.028 ± 0.004 (McClintock et al. 2019; Dark Energy Survey Collaboration 2020)¹⁵. The tension in the richness slope is identified as the driving factor for the discrepancy between cosmological results obtained from DES Y1 WL calibrated cluster number counts and cosmic shear and galaxy auto-and cross-correlations of the same data (Dark Energy Survey Collaboration 2020).

We find a sample lower scatter in the λ − M scaling relation of 0.3 ± 0.03 compared to the relation between C_R − M. This should somehow not mislead since, as mentioned, all the parameters of the mass-to-richness and the mass-to-count scaling relations are degenerate, see Fig. 5. Furthermore, it is reported that the scatter of the mass-to-richness scaling relation is better behaved in an X-ray-selected sample than in an optically-selected sample (Saro et al. 2015; Bleem et al. 2020; Grandis et al. 2021b; Giles et al. 2022).

The redshift evolution term has a best-fit value of C_λ = $-{0.49}_{-0.15}^{+0.14}$ shows a slight redshift trend in our sample. The best-fit value is consistent with the results from the CODEX sample for a sample with spectroscopic member galaxies λ > 3 ($-{0.26}_{-0.24}^{+0.23}$; Capasso et al. 2019); however, their results are consistent with no-redshift evolution, while ours indicate a redshift evolution at 3σ level. In the same CODEX sample, the negative trend in redshift evolution becomes stronger as the authors change the sample with clusters with a larger number of member galaxies, for instance, λ > 20 ($-{1.00}_{-0.56}^{+0.49}$; Capasso et al. 2019). In the DES sample(McClintock et al. 2019), the evolution trend is consistent with our measurements within the uncertainties. In their sample, no significant redshift evolution was observed −0 30 ± 0 30 ± 0 06 for a sample with richness cut of λ > 40. In contrast, the DES/SPT analyses reported in Costanzi et al. (2021) 0 95 ± 0 3 is in tension with our results at the >5σ level when the SPT number counts are included in the analysis, while consistent with our finding at <2σ level, −0.02 ± 0.34, when only SPT observable-to-mass relation is used (Costanzi et al. 2021). The reported redshift trends or slopes on the richness-mass relations depend on the cuts applied to the optically selected samples because of the rapid change in contamination fraction as a function of richness (Dark Energy Survey Collaboration 2020). Our sample relies mostly on X-ray selection, which, after optical cleaning, remains very little contaminated; further study is needed to understand better the effect of different selection effects on the best-fit values of this relation and their comparisons.

We note that the cluster richness is a sum of the color-based membership probability of the galaxies in the field (Rykoff et al. 2014). In Kluge et al. (2024), we show that the eROMaPPer returns different richness measurements because of the different luminosity thresholds adopted for the combination in the different optical filter bands, for instance, grz, griz, grzw1, due to variations in membership probabilities. Despite these differences, Kluge et al. (2024) demonstrate that there is a simple constant scaling factor between these richnesses, for example, λ_grz = S_normλ_griz = S′_normλ_grzw1, without any obvious dependence on redshift. This result is also supported by the reported comparison of the richness measured independently in the literature. For the eRASS1 clusters detected in the common area of the DES Y1 survey (Dark Energy Survey Collaboration 2020), the comparisons between the richness of eRASS1 1 clusters measured in the grz band and of DES Y1 clusters in the grizy band shows a scaling factor that has no dependency on redshift; as such the observed relation is λ_grz = 0.79λ_{DES Y1} Kluge et al. (see 2024, for further details). Therefore, this demonstrates that all the comparisons of richness to mass scaling relations we have performed above are reasonable because richness scaling affects only the normalization of the scaling relations, not the slope or redshift dependence.

7 Cosmological constraints

This section discusses the cosmology parameters constrained by eROSITA cluster abundance measurements compared to and combined with data from other cluster surveys and cosmological probes. In addition, our pipeline allows us to test different cosmological models beyond the ΛCDM concordance model. This current state-of-the-art cosmological model describes the formation and evolution of our Universe. In its name, ΛCDM, indicates the two primary ingredients of this cosmological model; Λ implies the dark energy introduced after the discovery of the accelerated expansion of the Universe (Riess et al. 1998; Perlmutter et al. 1999), while CDM stands for cold dark matter. When combined, the dark sector comprises 95% of the energy density of the Universe (Planck Collaboration VI 2020). In the ΛCDM framework, the parameters that describe the Universe are the Hubble constant H₀, the total matter density Ω_m (including dark matter and baryons), the spectral index of the primordial power spectrum n_s, the baryon density Ω_b, the amplitude of the matter fluctuations σ₈, and the reionization optical depth τ. The other cosmological parameters can be derived from this baseline set of parameters, for instance, the age of the Universem t₀, the critical density of the Universem ρ_c, and a combination of Ω_m with σ₈ along the principal direction of the degeneracy between the two, S₈. In our baseline model, the flat ΛCDM, we fit for Ω_m, σ₈¹⁶, H₀, n_s, and Ω_b.

In the fits providing the constraints beyond the ΛCDM model, where the dark energy equation of state or the summed neutrino masses are allowed to vary, namely, are not fixed to −1 and 0 as in the ΛCDM model. In the wCDM cosmology, we allow the dark energy equation of state to be free and leave summed neutrino masses fixed to 0 eV. In the vCDM cosmology, we allow summed neutrino masses to be free and leave the dark energy equation of state fixed to −1. We also provide joint constraints on the equation of the state of dark energy and the sum of neutrino masses in the flat vwCDM cosmological model.

7.1 Goodness of fit

In principle, determining the goodness of fit of the various cosmology models is not a trivial process in our analysis. The cosmology likelihood is not a simple distribution but comprises various components following different distributions. For instance, the component for the cluster number counts in the global likelihood (Sec. 4) should follow a Poisson distribution, while the likelihoods of the follow-up observables, such as weak lensing shear and optical richness, follow a log-normal distribution. Adding the likelihood for the mixture model complicates the process even further.

To assess the goodness of fit of the standard ΛCDM cosmo-logical model to the cluster abundance, richness, and count-rate distributions, we performed 1000 weighted¹⁷ realizations from best-fitting parameters extracted from the posterior distribution and calculate the χ² value. This simulation step is identical to the generation of the mock observations described in the validation phase; see Sec. 5.1 for further details. Overall, the total simulated number of objects in each bin of the observables is relevant in assessing the goodness of fit. In Fig. 6, we provide the distribution of the observables employed in the analysis: count rate, redshift, and richness, and the corresponding best-fit models. Our methodology also allows us to simultaneously distinguish the relative contributions of the clusters and the contamination (AGN and RS) to the overall total model. For illustration purposes, the Pearson residuals are widely used in assessing the goodness of fit in the fits of the data with Poissonian nature statistics, which is also shown in Fig. 6. As evident from the figure, the total model shown in the shaded region in orange produces a good fit to the data with χ² = 30.3 for the redshift distribution, χ² = 24.5 for the count rate distribution, and χ² = 71.5 for richness distribution, all three with 26 degrees of freedom (d.o.f.).

The best-fit shear model derived from the posterior of the scaling relation between count rate and cluster mass is shown on the observed stacked shear data from the DES (top panel), HSC (middle panel), and the KiDS (bottom panel) of Fig. 7. The stacking procedure is described in detail in Sect. 5.2 of Grandis et al. (2024) for DES, and Sect. 4.3 of Kleinebreil et al. (2024) for KiDS and HSC. We find χ² values of ${\chi }_{\text{DES}}^2={173.2}_{-27.6}^{+43.2}$ (150 d.o.f.), ${\chi }_{\mathrm{HSC}}^2={140.0}_{-6.2}^{+3.2}$ ( 70 d.o.f.), and ${\chi }_{\mathrm{KiDS}}^2={60.2}_{-5.6}^{+9.7}$ (63 d.o.f.) for fits to the shear data from DES, HSC, and KiDS data, respectively. Note that the upper and lower limits in the χ² distributions are obtained from the higher 1σ (upper edges of the shaded region in Fig. 7) and the lower 1σ (lower edge of the shaded region in Fig. 7) of the predicted mean shear profile, respectively. We note that HSC does not yield a good fit to our data. We note from the middle panel in Fig. 7) that this is mostly due to the highest redshift low count rate bins. To assess the robustness of our results and conclusions below, we have checked that running the analysis without the inclusion of HSC does not affect our results. We note that we find good agreement between shear measurements and WL-only mass calibration of eRASS1 clusters in the lensing analyses from three surveys we utilize (Kleinebreil et al. 2024, Sects. 4.1, and 4.2). In short, we show that HSC, DES, and KiDS WL are consistent both at the shear level and at the WL-calibrated scaling relation level (Kleinebreil et al. 2024).

Our strategy in analyzing cluster number counts involves utilizing weak lensing tangential shear information that allows us to simultaneously fit and constrain cosmological parameters and the parameters of the scaling law between cluster mass and observables. We perform extensive checks on the influence of the individual WL survey measurements in our analysis to ensure that these surveys do not introduce a significant bias in the scaling relations or fit for the cluster number counts presented in the companion papers; Grandis et al. (2024) and Kleinebreil et al. (2024). To this end, the scaling relation between count rate to mass is tested on the individual weak lensing datasets by setting the first two terms in the likelihood, in Equation (6), to zero and fitting only the last two terms. We find that the overall model is consistent between the three surveys, and the scaling relation has enough degrees of freedom (see Table 2 for free parameters) to be essentially fixed by the mass calibration. Furthermore, we highlight that using three independent WL surveys simultaneously safeguards us against possible hidden systematics from any individual WL survey.

Fig. 6 Distribution of clusters (black) as a function of redshift (top panel), count rate (middle panel), and richness (bottom panel) compared with the best fitting model prediction of the same quantities (blue) and the contribution to the total model from clusters (orange), AGN contamination (green), and RS contamination (purple). The black solid lines represent the data distribution, the shaded grey area the uncertainty due to measurement errors, and the black dotted lines indicate the Poisson uncertainty on the number of objects in each bin. In all figures, the bottom panel shows the Poisson Pearson residuals and the 2σ contours of the best-fit model.

7.2 Constraints on the ΛCDM cosmology

Cluster number counts can provide powerful constraints on the cosmological parameters, especially on Ω_m and σ₈, when the selection effects are adequately taken into account, and the bias in the mass calibration process is well-controlled. In our baseline ΛCDM cosmology analysis, we freeze the dark energy equation of state to w = −1, and the remaining best-fit parameters, Ω_m and σ₈, and S₈ are: $$\begin{array}{rl}{\mathrm{\Omega }}_{\mathrm{m}}& ={0.29}_{-0.02}^{+0.01},\\ {\sigma }_8& =0.88\pm 0.02,\\ S_8& =0.86\pm \mathrm{0.01.}\end{array}$$(34)

The posterior distributions of Ω_m and σ₈ from eROSITA only analysis are shown in contours in orange in Fig. 8. In the same figure, we also provide a comparison of our results with the other cluster number count experiments utilizing the WL mass calibration in the literature, e.g., WtG (Mantz et al. 2015), SPT-SZ (Bocquet et al. 2019), XXL (Garrel et al. 2022), and eFEDS (Chiu et al. 2023). We find consistent results within 2–3σ with the results in the literature, except for the cosmology results from DES cluster abundance measurements Dark Energy Survey Collaboration (2020). In Dark Energy Survey Collaboration (2020), authors report a value of ${\sigma }_8={0.85}_{-0.06}^{+0.04}$ consistent with our result, whereas their best-fit ${\mathrm{\Omega }}_{\mathrm{m}}={0.179}_{-0.038}^{+0.031}$ measurement is significantly lower than our value. A similar pattern is observed in comparisons of DES with the results from other experiments based on cluster counts. Costanzi et al. (2021) argues that the observed difference could be due to a biased mass-richness scaling relation. Compared to the previous cluster number count experiments mentioned above eRASS1 cosmology has the advantage of having a significantly larger number of clusters, an overall better understanding of the selection function, and, compared with optical surveys, an observable that we expect to be more closely linked to the actual virialized halos instead of line-of-sight projections. Furthermore, we have a mass reconstruction accuracy of up to a few percent, robustly established using the method proposed in Grandis et al. (2021a). We can then calibrate our masses with state-of-the-art wide area WL surveys (KiDS, DES, HSC) with a high S/N.

We compare our findings from cluster abundances with the recently published results from other cosmology probes, for instance, CMB measurements from Planck CMB (Planck Collaboration VI 2020), the joint cosmic shear constraints from DES and KiDS (Dark Energy Survey and Kilo-Degree Survey Collaboration 2023), the 3 × 2 correlation constraints from DES (Dark Energy Survey Collaboration 2022), lensing measurements from ACT DR6 + BAO (Qu et al. 2024), and lensing measurements from Planck (Planck Collaboration VIII 2020) in Fig. 9. Both σ₈ and Ω_m parameter measurements from the eROSITA cluster abundances agree well within 2σ confidence interval with CMB (Planck Collaboration VIII 2020) and lensing measurements (Dark Energy Survey Collaboration 2022; Qu et al. 2024; Planck Collaboration VIII 2020); however, they are in slight tension with the constraints from the cosmic shear measured from the combination of DES and KiDS data (Dark Energy Survey Collaboration 2022; Dark Energy Survey and Kilo-Degree Survey Collaboration 2023).

The correlation and the correlation coefficients between the free parameters in the fit are shown in Figs. B.2 and B.3. We find that the terms representing the slopes (B_X and B_λ) and the redshift evolution (G_X and C_λ) of the scaling relations between X-ray observable and richness are correlated with Ω_m. In contrast, these terms are anti-correlated with σ₈. Interestingly, S₈ shows a strong correlation with the amplitude of the scaling relations (A_X and A_λ), but we do not observe any correlation between S₈ and the mass slope and the redshift evolution. Notably, none of the best-fit cosmological parameters are degenerate with the scatter terms in scaling relations. This indicates that the observables employed in this work do not need to be tightly correlated with cluster mass, as only the mass distribution at the population level plays a determining role in cosmology. In our forward modeling approach, it is crucial to statistically reproduce the observables at a population level; therefore, the X-ray count rate, despite being a high-scatter mass proxy, is a powerful observable because of its role in the selection process. To calculate the consistency of our result with CMB measurements, we (1) resampled our MCMC chains and the Planck CMB 2020 chains to obtain the same number of points; (2) took the differences between the chains and assume that it is distributed according to a multivariate normal distribution with data inferred mean and covariance. This allows for a calculation of the ${\chi }^2={\overline{\mu }}^T{C}^{-1}\overline{\mu }$; and 3) we computed the survival function¹⁸ of the χ² distribution.

This method allows for an immediate estimate of the departure from the best fit in terms of its standard deviation. We find that our measurements on Ω_m and σ₈ agree with the Planck Collaboration VI (2020)’s CMB measurements (TT, TE, EE+lowl+lowE combination; we use Planck CMB 2020 hereafter) at the 2.3σ confidence level.

To break the degeneracy between Ω_m and σ₈ and have a fair comparison with the Planck CMB 2020 constraints, we combine the two probes given the consistency between the two. The combination of eROSITA cluster abundances and Planck CMB 2020 chains yields the following cosmological parameters: $$\begin{array}{rl}{\mathrm{\Omega }}_{\mathrm{m}}& =0.32\pm 0.01,\\ {\sigma }_8& =0.82\pm 0.01,\\ S_8& =0.85\pm \mathrm{0.01.}\end{array}$$(35)

We show the combined results in thick red contours in Figs. 8 and 9. Compared to the previous eROSITA-only results in Equation (34), we see a shift to a higher value of Ω_m and to a lower value for σ₈, a significant improvement on our constraints on the Ω_m and σ₈ parameter plane. This is because the Planck CMB 2020 constraints on Ω_m and σ₈ are tight and degenerate along a direction orthogonal to the S₈ direction. As a result, the combined best-fit parameters are at the parameter space where these two constraints intersect. It is worth noting that Planck CMB 2020 ΛCDM constraints are obtained with a different set of assumptions compared to ours, namely, the sum of neutrino masses are fixed in their fits to 0.06 eV, the lower bound on the neutrino masses from oscillation experiments (Tanabashi et al. 2018). In our ΛCDM analysis, we set the neutrino mass to 0 eV. Figure B.4 shows the negligible impact of fixing the summed neutrino on our constraints to a lower limit of 0.06 eV obtained from the ground-based experiments (Tanabashi et al. 2018).

A more direct way to compare our constraints with the results from the CMB analysis is to measure the S₈, which combines the constraints on parameters Ω_m and σ₈. Our data constrains the parameter S₈ with high precision due to the direction and orientation of the degeneracy between these two parameters. The best-fit α_best value that minimizes the standard deviation of σ₈ × (Ω_m/0.3)^α along the direction of our degeneracy, is α_best = 0.42 ± 0.02. Our best-fit value is close to 0.5, the canonical value of α used in the definition of S₈ in Planck Collaboration VI (2020). The cluster cosmology experiments, e.g., Bocquet et al. (2019) reports α_best = 0.2 for the SPT-SZ sample, while Chiu et al. (2023) finds α_best = 0.3 for the eFEDS clusters, a smaller value than what we find in this work. On the other hand, the power spectrum analysis of the shear measurements of the HSC field produces an α_best value of 0.45, consistent with our results Hikage et al. (2019). Regarding the S₈ constraints, we obtain a degeneracy between σ₈ and Ω_m, which is almost perfectly aligned with the slope in the S₈ definition. This makes the constraints on σ₈ tighter in our experiment than other surveys. Similarly, the recent SPT constraints show the slope α_best ~ 0.25 (Bocquet et al. 2024), and indeed their optimal parameter combination σ₈ × (Ω_m/0.3)^0.25 measure a relative uncertainty very close to the uncertainty we obtain on S₈ in the eRASS1 sample. The degeneracy reported in other experiments does not have a direction in line with the S₈ definition.

Figure 10 shows our constraint on the S₈ parameter and its comparisons with the selected results from experiments utilizing cluster abundances and weak lensing mass calibration in the literature (Bocquet et al. 2019; Garrel et al. 2022), cosmic shear (Amon et al. 2022; Secco et al. 2022; Li et al. 2023; van den Busch et al. 2022; Dark Energy Survey and Kilo-Degree Survey Collaboration 2023), clustering (Dark Energy Survey Collaboration 2022; Sugiyama et al. 2023; Qu et al. 2024). Our measurement is consistent with results from CMB measurements with Planck and WMAP (Hinshaw et al. 2013; Planck Collaboration VI 2020), ACT and Planck lensing constraints (Qu et al. 2024; Madhavacheril et al. 2024), the HSC clustering results (Sugiyama et al. 2023), and the cluster count measurements based on the XXL survey within 2σ confidence interval (Planck Collaboration VI 2020; Garrel et al. 2022). However, it is in significant tension (≥3σ level) with constraints from other probes, in particular from cosmic shear from DES, HSC, KiDS surveys individually, and the combined KiDS and DES cosmic shear measurements (Amon et al. 2022; Secco et al. 2022; van den Busch et al. 2022; Li et al. 2023), DES 3 × 2 clustering (Sugiyama et al. 2023), and the SDSS BAO and redshift space distortions (RSD) (Alam et al. 2021) measurements. Given the consistency of our results with the CMB measurements, we combine our results with Planck 2020 CMB chains. The results of this combined analysis result in a negligible improvement on S ₈ constraints, which can be seen in the second data point from the top on Fig. 10. This behavior can be easily explained that the constraints on the S₈ from eRASS1 are already tight; therefore, the combination of the parameters does not significantly improve precision.

As the error budget on our measurements is smaller than the systematic uncertainties due to the assumption of the halo mass function, we employ the Tinker et al. (2008) halo mass function because of its simple semi-analytical form. To test the effect of the assumption on the halo mass function on our results, we fit the ΛCDM to our data with the Tinker et al. (2010) and Despali et al. (2016) halo mass functions. The results of the comparisons are shown in Fig. B.6. We observe a slight decrease in the parameter σ₈ from 0.88 ± 0.02 to 0.86 ± 0.02 when the Tinker et al. (2010) and Despali et al. (2016) halo mass functions are used. The shift is statistically consistent within 1σ intervals with the Tinker et al. (2008) results. The change in the halo mass function does not significantly impact the parameter ${\mathrm{\Omega }}_{\mathrm{m}}={0.30}_{-0.02}^{+0.01}$ respectively. This shows that the choice of the halo mass function does not yet limit our cosmological constraints in this work. As the eROSITA survey gets deeper and the constraints on the cosmological parameters get tighter, the uncertainties related to the halo mass function will become a dominant source of systematics and should therefore be treated with caution (see tests in Appendix B.4 and Figure B.6).

Fig. 7 Goodness of fit for the three WL surveys, DES (top panel), HSC (middle panel), and KiDS (bottom panel). Each panel shows a stacked observed profile in the count rate and redshift bin (stars), compared with the best-fitting weak lensing shear profiles according to mass calibration, i.e., the weak lensing masses based on the count rates of the clusters and calibrated scaling relations. The color bar indicates the count rate binning applied to DES and HSC, while tomographic bins for KiDS. For a more extensive description of the weak lensing goodness of fit, see Grandis et al. (2024), and Kleinebreil et al. (2024).

Fig. 8 Posterior distributions on the parameters Ωm and σ8 from the ΛCDM fit on eRASS1 data shown in orange. In blue, we show the Planck CMB 2020 constraints without combination with BAO and SNe Ia Planck Collaboration VI (2020). In red, we show the combination of eRASS1 with Planck CMB. As a comparison, we also present previous results from similar cluster surveys that employ weak lensing shear data in their mass calibration, e.g., WtG (Mantz et al. 2015), DES (Dark Energy Survey Collaboration 2020), SPT-SZ (Bocquet et al. 2019), eFEDS (Chiu et al. 2023), XXL (Garrel et al. 2022), and SPT (Bocquet et al. 2024) surveys.

Fig. 9 Posterior distributions on Ωm and σ8 from the ΛCDM fit to eRASS1 data are shown in orange. Similarly, the Planck CMB 2020 best-fit contours and eRASS1 with Planck CMB are shown in blue and red (Planck Collaboration VI 2020). This figure provides a comparison of eRASS1 with other cosmology probes, e.g., the cosmic shear constraints from the joint analysis of DES and KiDS data (Dark Energy Survey and Kilo-Degree Survey Collaboration 2023), DES 3x2 clustering constraints (Dark Energy Survey Collaboration 2022), KiDS 3x2 constraints (Heymans et al. 2021), results from the combined ACT lensing and BAO (Qu et al. 2024), Planck lensing analyses (Planck Collaboration VIII 2020), and SDSS BAO RSD (Alam et al. 2021).

Fig. 10 Comparison of the constraints on $S_8={\sigma }_8\times \sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$ with literature results from SPT-SZ (Bocquet et al. 2019), from SPT-SZ, SPTpol ECS, and SPTpol 500d cluster counts (Bocquet et al. 2024), from XXL cluster number counts and clustering (Garrel et al. 2022), from Planck CMB 2020 (Planck Collaboration VI 2020), from WMAP CMB (Hinshaw et al. 2013), from ACT lensing (Qu et al. 2024), from ACT lensing (Madhavacheril et al. 2024), from Planck lensing (Planck Collaboration VIII 2020), from DES 3x2 correlation function (Dark Energy Survey Collaboration 2022), from HSC 3x2 correlation function (Sugiyama et al. 2023), from DES cosmic shear (Amon et al. 2022; Secco et al. 2022), from HSC cosmic shear (Li et al. 2023), from KiDS cosmic shear (van den Busch et al. 2022), from joint DES and KiDS cosmic shear (Dark Energy Survey and Kilo-Degree Survey Collaboration 2023), from SDSS BAO and RSD (Alam et al. 2021). All the errors shown are at the 1σ level. The vertical line and grey shaded areas represent our measurement’s location and 1–2σ confidence level.

7.3 Constraints on the νCDM cosmology

Despite their small size and mass, the left-handed neutrinos influence the formation of large-scale structure by hindering the formation of small-scale haloes and leaving their imprints on the gravitational collapse process. The relative number of low-mass haloes in cluster abundance measurements can shed light on the summed masses of three left-handed neutrino species. With the discovery of neutrino oscillations between different flavors, electron, muon, and tau neutrinos, it has become clear that neutrino eigenstates must have non-zero mass (Fukuda et al. 1998; Ahmad et al. 2002). The theoretical predictions indicate a hierarchical relation between the mass eigenstates of the three species of neutrinos; however, their ordering in mass remains elusive. Based on the normal neutrino mass hierarchy model, the third neutrino mass eigenstate, associated with the tau neutrino, is the heaviest among the three. This is followed by the eigenstate of the muon neutrino, which has an intermediate mass, and that of the electron neutrino, namely, m_τ > m_μ > m_e. The ordering is swapped in the inverted mass hierarchy model such that m_τ < m_μ <m_e. Recent constraints of their summed mass show that the lower limit provided by the oscillation experiments depends on the assumed underlying hierarchy model (Qian & Vogel 2015; Esteban et al. 2020), with a lower limit of ∑ m_v > 0.059 eV for normal mass hierarchy and ∑ m_v > 0.101 eV for the inverted hierarchy models (Tanabashi et al. 2018; Athar et al. 2022). Constraining the summed mass to <0.1 eV implies that the inverted model is excluded. The ground-based experiments through beta decay of tritium imply the sum to be ∑ m_v < 1.1 eV at 90% confidence level, a narrow range for the allowed summed mass (Aker et al. 2019). The space-based Planck CMB measurements provide a similar upper limit of ∑ m_v < 0.26 eV at 95% confidence level (Planck Collaboration VI 2020).

As we probe the largest collapsed objects in the Universe, galaxy clusters, the cluster number counts can be used to constrain the summed masses of the neutrinos. To do so, we allow the sum of neutrino masses to be free in the cosmology pipeline, with a uniform prior of $\mathcal{U}$(0 eV, 1 eV). Including massive neutrinos in cosmological analysis requires a fine-tuned interaction with the choice of the halo mass function. As shown in Costanzi et al. (2013), a few changes to the use of the matter power spectrum are fundamental to properly make use of Tinker et al. (2008) halo mass function in the presence of massive neutrinos. As the pyccl 3.0.0 library of Chisari et al. (2019) available on GitHub¹⁹ does not implement the prescriptions of Costanzi et al. (2013), therefore accounting for corrections to the power spectrum, matter density, and, finally, to the halo mass function, we modify the local version of pyccl by editing several hardcoded parts of code to allow the user to employ Costanzi et al. (2013) prescription. We stress that our lower limit is set to 0 eV instead of a value of 0.059 eV (Tanabashi et al. 2018) adopted by the Planck CMB 2020 analysis. As the lower limit on the neutrino mass depends on the assumptions of the neutrino mass model, we do not make any prior assumptions about its value and adopt uniform priors on its value to avoid this parameter from affecting our posterior distributions. The cosmological constraints with free neutrino mass components are: $$\begin{array}{rl}{\mathrm{\Omega }}_{\mathrm{m}}& ={0.28}_{-0.02}^{+0.02},\\ {\sigma }_8& =0.88\pm 0.03,\\ S_8& =0.85\pm 0.02,\\ \sum m_{\nu }& <0.43\mathrm{e}\mathrm{V}(95\mathrm{\%}\mathrm{C}\mathrm{L}).\end{array}$$(36)

These results are visualized in Fig. 11. The upper limit to the sum of neutrino masses is ∑ m_v < 0.43 eV. The addition of the Costanzi et al. (2013) correction to the power spectrum relaxes the upper limits on neutrino masses from 0.22 eV to 0.43 eV, as expected. Our results represent the tightest limits on the sum of neutrino masses from cluster abundance experiments; for instance, the SPT-SZ sample results in an upper limit of <0.74 eV (95% confidence interval Bocquet et al. 2019), while the updated SPT results yield < 0.55 eV (95% confidence interval Bocquet et al. 2024). We verify that the values of the Ω_m and cr₈ remain statistically consistent when ∑ m_v is allowed to be non-zero, see Fig. B.4. We find that eRASS1 and Planck CMB vCDM parameters are consistent at the 2.0σ level.

Given the excellent agreement with Planck CMB measurements, the resulting cosmological parameters of two probes can be combined in a statistically meaningful way, enabling a much tighter measurement of the impact of massive neutrinos on both the formation and evolution of the large-scale structure and on the primordial density field. As performed in the ΛCDM analysis, we combine our results with the Planck CMB constraints to break the degeneracy between Ω_m and σ₈. We obtain: $$\begin{array}{rl}{\mathrm{\Omega }}_{\mathrm{m}}& =0.31\pm 0.01,\\ {\sigma }_8& =0.83\pm 0.01,\\ S_8& =0.84\pm 0.01,\\ \sum m_v& <0.14\mathrm{e}\mathrm{V}(95\mathrm{\%}\mathrm{C}\mathrm{L}).\end{array}$$(37)

Simultaneously fitting our measurements with Planck CMB 2020 likelihood chains yield an upper limit of ∑ m_v < 0.14 eV, consistent with the results in the literature results with marginally consistent with the inverted mass hierarchy model, which requires ∑ m_ν > 0.101 eV, where just 18% of the eROSITA- Planck combined samples are above this limit and consistent with the inverted hierarchical model.

As a final step, we combined eRASS1 cluster abundance measurements via an importance sampling with Planck CMB and with the lower limits from ground-based oscillation experiments (Tanabashi et al. 2018). In the case of a normal mass hierarchy scenario, the summed masses of $\sum m_v={0.09}_{-0.02}^{+0.04}\mathrm{eV}$; while assuming the inverted mass hierarchy model, we obtain summed masses of $\sum m_v={0.12}_{-0.02}^{+0.03}\mathrm{eV}$. It is interesting to point out (from Fig. C.1) that the constraints on the mass of the lightest neutrino eigenstates are similar: for both mass hierarchy we obtain $m_{\text{light}}={0.017}_{-0.010}^{+0.016}\mathrm{eV}$ (68% confidence intervals) or with an upper limit of 0.044 eV.

Fig. 11 Best-fit values for vCDM cosmological model from eRASS1 cosmological experiment (orange), Planck CMB 2020 data (blue), and combination of eROSITA and Planck CMB (red) are shown. The SPT-SZ cluster abundances and SDSS BOSS DR16 BAO + CMB lensing results are shown in green and brown (Bocquet et al. 2019; Alam et al. 2017). SPT recent constraints from Bocquet et al. (2024) are added in purple. eRASS1 cosmological experiment places the very tight constraints on the sum of the neutrino masses to be ∑ mv < 0.43eV from cluster number counts alone.

Fig. 12 Best-fit contours of w and Ωm in the wCDM cosmological model from eROSITA experiment (orange), Planck CMB data (blue), SDSS BOSS DR16 BAO constraints (green), Pantheon SNe Ia constraints (brown), DES 3×2 constraints (grey) are shown. We also show the resulting combination of eROSITA and Pantheon SNe Ia (red), the combination of Planck CMB 2020 with Pantheon SNe Ia and BAO (purple), and finally, the combination of eROSITA with Planck CMB 2020, Pantheon SNe Ia and SDSS BOSS DR16 BAO (black). eRASS1 experiment can place tight constraints on the dark energy equation of state parameter without any need to combine with other cosmological experiments.

7.4 Constraints on the wCDM cosmology

The late-time acceleration of the expansion of the Universe is one of the major discoveries of modern physics for the standard ΛCDM cosmology model, where the cosmological constant Λ drives the late-time acceleration (Riess et al. 1998; Perlmutter et al. 1999). Several independent experiments from Supernovae Ia (Riess et al. 1998; Astier et al. 2006; Kowalski et al. 2008), CMB (Komatsu et al. 2011; Planck Collaboration VI 2020), and BAO (Beutler et al. 2011; Alam et al. 2017, 2021) produced constraints on fractional energy density of the dark energy to the total energy of the Universe to be ~67–70%. Despite the intensive experimental effort invested into constraining the significant contribution of dark energy, its nature remains elusive. One of the main science drivers of eROSITA is to place constraints on its equation of state using the cluster abundances. We define w, the equation of state, as the ratio of the pressure of dark energy to its energy density (w = p/ρ) and assume that its density evolves with ρ ∝ a^−3(1+w) in an isotropic Universe. Previous dark energy equation state measurements point to w ~ −1 with no evolution with redshift (Astier et al. 2006; Beutler et al. 2011; Alam et al. 2017; Planck Collaboration VI 2020). In this section, we consider an extension to the standard ΛCDM by allowing the equation of state of dark energy to w to vary.

To constrain the value of w, we set the sum of neutrino masses back to zero (∑ m_v = 0 eV), and we allow the w to be free with a uniform prior $\mathcal{U}$(−2, −0.33) in the wCDM model. Naturally, a change in the equation of the state affects both the background and the hierarchical growth of structure as it changes the expansion history of the Universe and the cosmological distances converted from the redshift measurements. We obtain the best-fit parameters as: $$\begin{array}{rl}{\mathrm{\Omega }}_{\mathrm{m}}& =0.28\pm 0.02,\\ {\sigma }_8& =0.88\pm 0.02,\\ S_8& =0.86\pm 0.02,\\ w& =-1.12\pm \mathrm{0.12.}\end{array}$$(38)

Our fits appear to prefer w to be −1.12 ± 0.12, shown in Fig. 12. consistent with the results from Planck CMB (Planck Collaboration VI 2020), DES Y3 3×2 pt (Dark Energy Survey Collaboration 2022), Pantheon SNe Ia (Scolnic et al. 2018, 2022; Brout et al. 2022), DES Y5 SNe Ia (DES Collaboration 2024), and SDSS BOSS DR16 BAO (Alam et al. 2017), at 1 − 2σ level. The best-fit parameters of the scaling relations, the mixture model, and the correlation between parameters of this model are consistent with the baseline ΛCDM case.

Given the consistency between eROSITA and Pantheon SNe Ia results on w at 0.5σ level and the orthogonal direction of the degeneracy between Ω_m and w, we combine two datasets to break degeneracies and obtain; $$\begin{array}{rl}{\mathrm{\Omega }}_{\mathrm{m}}& =0.30\pm 0.01,\\ {\sigma }_8& =0.87\pm 0.02,\\ S_8& =0.87\pm 0.01,\\ w& =-{0.95}_{-0.04}^{+0.05}.\end{array}$$(39)

Compared with our results, the parameters Ω_m, σ₈, and S₈ are not significantly affected by the inclusion of Pantheon SNe Ia data. However, our best fitting w shifts to a slightly higher value of $-{0.95}_{-0.04}^{+0.05}$, consistent with the eROSITA measurements with a significantly increased precision. w remains consistent with the standard ΛCDM value of −1 at 1σ level.

Dark energy is responsible for a late-time expansion of the Universe and, therefore, has a negligible effect on the CMB signal. The constraining power of Planck CMB data on the dark energy equation of state is quite limited. We find that eRASS1 constraints on w are consistent w at the 1.4σ and 1.8σ level with Planck CMB 2020 and SDSS BOSS DR16 BAO results, respectively. Given the consistency, we combine eRASS1 chains with both Pantheon SNe, Planck CMB, and SDSS BAO data. Including these two extra constraints helps improve the precision of the parameters significantly. We measure; $$\begin{array}{rl}{\mathrm{\Omega }}_{\mathrm{m}}& =0.31\pm 0.01,\\ {\sigma }_8& =0.84\pm 0.01,\\ S_8& =0.85\pm 0.01,\\ w& =-1.06\pm \mathrm{0.03.}\end{array}$$(40)