# The SRG/eROSITA All-Sky Survey: Cosmology Constraints from Cluster Abundances in the Western Galactic Hemisphere

V Ghirardini, E Bulbul, E Artis, N Clerc, C Garrel, S Grandis, M Kluge, A Liu, Y.E Bahar, F Balzer, et al.

► To cite this version:

V Ghirardini, E Bulbul, E Artis, N Clerc, C Garrel, et al. The SRG/eROSITA All-Sky Survey: Cosmology Constraints from Cluster Abundances in the Western Galactic Hemisphere. Astronomy & Astrophysics - A&A, 2024, 689, pp.A298. 10.1051/0004-6361/202348852. hal-04480641

HAL Id: hal-04480641

https://hal.science/hal-04480641v1

Submitted on 26 Jun 2025

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# The SRG/eROSITA all-sky survey

## Cosmology constraints from cluster abundances in the western Galactic hemisphere

V. Ghirardini<fnref n="1" />,*, E. Bulbul<fnref n="1" />, E. Artis<fnref n="1" />, N. Clerc<fnref n="2" />, C. Garrel<fnref n="1" />, S. Grandis<fnref n="3" />,<fnref n="5" />, M. Kluge<fnref n="1" />, A. Liu<fnref n="1" />, Y. E. Bahar<fnref n="1" />, F. Balzer<fnref n="1" />, I. Chiu<fnref n="4" />, J. Comparat<fnref n="1" />, D. Gruen<fnref n="5" />, F. Kleinebreil<fnref n="3" />,<fnref n="8" />, S. Krippendorf<fnref n="5" />,<fnref n="6" />, A. Merloni<fnref n="1" />, K. Nandra<fnref n="1" />, N. Okabe<fnref n="7" />, F. Pacaud<fnref n="8" />, P. Predehl<fnref n="1" />, M. E. Ramos-Ceja<fnref n="1" />, T. H. Reiprich<fnref n="8" />, J. S. Sanders<fnref n="1" />, T. Schrabback<fnref n="3" />,<fnref n="8" />, R. Seppi<fnref n="1" />, S. Zelmer<fnref n="1" />, X. Zhang<fnref n="1" />, W. Bornemann<fnref n="1" />, H. Brunner<fnref n="1" />, V. Burwitz<fnref n="1" />, D. Coutinho<fnref n="1" />, K. Dennerl<fnref n="1" />, M. Freyberg<fnref n="1" />, S. Friedrich<fnref n="1" />, R. Gaida<fnref n="1" />, A. Gueguen<fnref n="1" />, F. Haberl<fnref n="1" />, W. Kink<fnref n="1" />, G. Lamer<fnref n="9" />, X. Li<fnref n="11" />, T. Liu<fnref n="1" />, C. Maitra<fnref n="1" />, N. Meidinger<fnref n="1" />, S. Mueller<fnref n="1" />, H. Miyatake<fnref n="13" />,<fnref n="14" />,<fnref n="15" />, S. Miyazaki<fnref n="12" />, J. Robrade<fnref n="10" />, A. Schwope<fnref n="9" />, and I. Stewart<fnref n="1" />

Received 5 December 2023 / Accepted 17 June 2024

## ABSTRACT

The evolution of the cluster mass function traces the growth of linear density perturbations, providing valuable insights into the growth of structures, the nature of dark matter, and the cosmological parameters governing the Universe. The primary science goal of eROSITA, on board the Spectrum Roentgen Gamma (SRG) mission, is to constrain cosmology through the evolution of the cluster mass function. In this paper, we present a set of cosmological constraints obtained from 5259 clusters of galaxies detected over an area of 12791 deg<fnref n="2" /> in the western Galactic hemisphere of eROSITA’s first All-Sky Survey (eRASS1). The common footprint region (4968 deg<fnref n="2" />) between the eROSITA Survey and Dark Energy Survey (DES), the Kilo-Degree Survey (KiDS), and the Hyper Supreme Camera (HSC) survey is used for calibration of the scaling between X-ray count rate of the clusters and their total mass through measurements of their weak gravitational lensing signal. The eRASS1 cluster abundances constrain the ΛCDM parameters, namely, the energy density of the total matter to Ω_m = 0.29<sup>+0.01</sup><sub>-0.02</sub> and the normalization of the density fluctuations to σ_8 = 0.88 ± 0.02, and their combination yields S_8 = σ_8(Ω_m/0.3)<sup>0.5</sup> = 0.86 ± 0.01. These results are consistent and achieve at a similar precision with state-of-the-art cosmic microwave background (CMB) measurements. Furthermore, the eRASS1 cosmological experiment places a most stringent upper limit on the summed masses of left-handed light neutrinos to Σ m_ν < 0.43 eV (95% confidence interval) from cluster number counts alone. By combining eRASS1 cluster abundance measurements with CMB- and ground-based neutrino oscillation experiments, we measured the summed neutrino masses to be Σ m_ν = 0.09<sup>+0.04</sup><sub>-0.02</sub> eV or Σ m_ν = 0.12<sup>+0.03</sup><sub>-0.02</sub> eV, assuming a normal or inverted mass hierarchy scenario for neutrino eigenstates. The eRASS1 cluster abundances significantly improve the constraints on the dark energy equation of state parameter to w = -1.12 ± 0.12. When Σ m_ν and w are left free, we find consistent results with the concordance ΛCDM cosmology. Our results from the first All-Sky Survey improve the cosmological constraints by over a factor of 5 to 9 over the previous cluster surveys, establishing cluster abundance measurements for precision cosmology and setting the stage for deeper eROSITA All-Sky Surveys, as well as for future cluster abundance experiments.

**Key words.** surveys – galaxies: clusters: general – galaxies: clusters: intracluster medium – cosmology: observations – X-rays: galaxies: clusters

<footnote n="1">¹ Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse 1, 85748 Garching, Germany</footnote>

<footnote n="2">² IRAP, Université de Toulouse, CNRS, UPS, CNES, 31028 Toulouse, France</footnote>

<footnote n="3">³ Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstr. 25/8, 6020 Innsbruck, Austria</footnote>

<footnote n="4">⁴ Department of Physics, National Cheng Kung University, 70101 Tainan, Taiwan</footnote>

<footnote n="5">⁵ Universitäts-Sternwarte, Faculty of Physics, LMU Munich, Scheinerstr. 1, 81679 München, Germany</footnote>

<footnote n="6">⁶ Arnold Sommerfeld Center for Theoretical Physics, LMU Munich, Theresienstr. 37, 80333 München, Germany</footnote>

<footnote n="7">⁷ Department of Physical Science, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan</footnote>

<footnote n="8">⁸ Argelander-Institut für Astronomie (AIfA), Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany</footnote>

<footnote n="9">⁹ Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany</footnote>

<footnote n="10">¹⁰ Hamburger Sternwarte, Gojenbergsweg 112, 21029 Hamburg, Germany</footnote>

<footnote n="11">¹¹ McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA</footnote>

<footnote n="12">¹² Subaru Telescope, National Astronomical Observatory of Japan, 650 N Aohoku Place Hilo HI 96720 USA</footnote>

<footnote n="13">¹³ Kobayashi-Maskawa Institute for the Origin of Particles and the Universe (KMI), Nagoya University, Nagoya 464-8602, Japan</footnote>

<footnote n="14">¹⁴ Institute for Advanced Research, Nagoya University, Nagoya 464-8601, Japan</footnote>

<footnote n="15">¹⁵ Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study (UTIAS), The University of Tokyo, Chiba 277-8583, Japan</footnote>

<footnote>* Corresponding author; vittorio@mpe.mpg.de</footnote>
# 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_8$ (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_0$ (Riess et al. 2022). However, tensions between CMB early-time measurements and large-scale structure late-time experiments have emerged: the so-called $H_0$-tension (e.g. Di Valentino et al. 2021) and the so-called $S_8$-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^{-1}$ Mpc, namely $\sigma_8^1$, and the present-day matter density parameter, $\Omega_m^2$. 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, $\sum m_\nu$, 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 redshift. 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

<footnote n="1">¹ See Appendix E for detailed definition of $\sigma_8$.</footnote>

<footnote n="2">² $\Omega_m = \rho_m/\rho_c$, where $\rho_m$ is the matter density and $\rho_c$ is the critical density of the Universe.</footnote>
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 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 Galactic 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 12 247 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*<sub>l</sub> < 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 Sect. 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 Sect. 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 Sect. 8. We provide our summary in Sect. 9. Our convention is as follows: the overdensity radius of $R_{500}$ 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_{500}$ is the cluster mass within a sphere of $R_{500}$. With log, we are referring to the natural logarithm, and when using base 10, we explicitly indicate it with the subscript $log_{10}$. Throughout this paper, we use the terms consistency, tension, and discrepancy multiple times: we define here that by “consistent” we mean less than $3\sigma$ difference, by “in tension” we mean between $3\sigma$ and $5\sigma$ in difference, and by “discrepant”, we mean above $5\sigma$ 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
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Sky distribution of the 5259 confirmed clusters of galaxies in the cosmology subsample, detected in the western Galactic 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 deg². 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 deg² in total.</figcaption></figure>

groups and clusters catalog. To compile the cosmology subsample, Bulbul et al. (2024) applied a conservative extent likelihood cut of $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. $\le 32.5^\circ$ (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 ($\hat{z}_\lambda$) and richness ($\lambda$) measurements. The common footprint between the eROSITA survey in the western Galactic hemisphere and the southern LS DR10-South footprint is 12 791 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 Sect. 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 12 791 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 left-most panel of Fig. 2. The mixture model, extensively described in Sect. 4.2, implies that the purity of this sample is at the level
<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>

of 94%, ideal for cosmological studies with cluster counts (see Bulbul et al. 2024; Kluge et al. 2024, for further details).

## 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 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_\lambda$. 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$^2$, are shown in Fig. 1. In particular, eRASS1 has a 4060 deg$^2$ common area with DES, 679 deg$^2$ with HSC, and 1108 deg$^2$ with KiDS$^3$. 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

<footnote n="3">³ The reported areas include the entire footprint, without taking into account any masking or regions where shape catalog was not available.</footnote>
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^4 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 multi-band (*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_{cl}+0.2}^{\infty} P(z) dz > 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(\theta)$ as a function of the

angular radius $\theta$ 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 $\Lambda$CDM cosmology with $h = \frac{H_0}{100 \text{ km/s/Mpc}} = 0.7$ and $\Omega_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\sigma$ upper limit of 6% at the innermost bin ($R \approx 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} \le 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, $\delta m$, at 1.7% level (Mandelbaum et al. 2018).

As a result, the shear profile $g_t(\theta)$, 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<fnref n="5" /> (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

<footnote n="4">⁴ DES Y3 catalogs are publicly available at https://des.ncsa.illinois.edu/releases/y3a2/Y3key-catalogs</footnote>

<footnote n="5">⁵ KiDS catalogs are publicly available at https://kids.strw.leidenuniv.nl/DR4/KiDS-1000_shearcatalogue.php</footnote>
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^{-1}$ Gpc) and resolution (4096³ particles) enable reliable prediction of the large-scale structure and its haloes down to $10^{11} M_{\odot}$. 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 \times 10^{13} M_{\odot}$), 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) = \sigma^2 \exp(-r^2/2)$, with $\sigma$ 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 $\sigma$ 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^6$, and the sky position, $H_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 $\approx 4.8 \times 10^5$). 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^{16}$ and 12 000 iterations provide good results within a tractable time. The optimization of a set of 651 hyperparameters$^7$ 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 Sect. 5), six models were computed, differing only by the exact definition of

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Reliability of the selection function model, as estimated from a test sample comprising 2 566 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\sigma$ 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.</figcaption></figure>

an object that qualifies as “detected” in the simulation. We modified the lower thresholds in selection parameters ($\mathcal{L}_{\text{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 Sect. 2 ($\mathcal{L}_{\text{ext}} > 6$, EXTENT > 0). Each model is tested against the remaining third of the simulation sets (each of size $\approx 2.5 \times 10^5$). 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, H_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 ($H_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.

## 3.5. Optical selection function

As described in Sect. 2.1, the cosmology cluster catalog is compiled after optical measurements of richness and redshift

<footnote>$^6$ We note that by correcting the observed counts using exposure map to get count rate and correcting for the absorption due to hydrogen we effectively make the X-ray observable sky independent.</footnote>

<footnote>$^7$ With $\sigma$, the lengthscales along the five dimensions, 30 × 5 coordinates for the inducing points, 30 numbers for the mean values at the inducing points, and 30 × 31/2 values for the Cholesky decomposition of the covariance matrix between those values.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Representation of a model $P(I|z, C_R, \mathcal{H}_i)$ for fixed values of cluster redshift, $z = [0.15, 0.3, 0.5]$, local galactic absorption $n_H = 3.5 \times 10^{20} \text{cm}^{-2}$, a nominal background level $5.2 \times 10^{-15} \text{erg/s/cm}^2/\text{arcmin}^2$, three values of exposure times [50, 150, 500] s, as a function of count-rate.</figcaption></figure>

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 $\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 ($\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).

## 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, $\lambda$, and weak-lensing tangential shear profile, $\hat{g}_i$;

$$ \log \mathcal{L} = \sum_j \log \left( \frac{dN_{tot}}{d\hat{C}_R d\hat{z} d\hat{\lambda} d\hat{g}_i d\hat{\mathcal{H}}_i} P(I|\hat{C}_R, \hat{z}, \hat{\lambda}, \hat{\mathcal{H}}_i, \hat{g}_i) \right) - \int \cdots \int \left[ \frac{dN_{tot}}{d\hat{C}_R d\hat{z} d\hat{\lambda} d\hat{g}_i d\hat{\mathcal{H}}_i} P(I|\hat{C}_R, \hat{z}, \hat{\lambda}, \hat{\mathcal{H}}_i, \hat{g}_i) \right. \\ \left. \qquad \qquad \qquad \qquad \qquad d\hat{C}_R d\hat{z} d\hat{\lambda} d\hat{g}_i d\hat{\mathcal{H}}_i \right], \quad (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|\hat{C}_R, \hat{z}, \hat{\lambda}, \hat{\mathcal{H}}_i, \hat{g}_i) = P(I|\hat{C}_R, \hat{z}, \hat{\mathcal{H}}_i) \Theta(\hat{\lambda} > 3), \quad (2) $$

where $\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{dN_{tot}}{d\hat{C}_R d\hat{z} d\hat{\lambda} d\hat{g}_i d\hat{\mathcal{H}}_i} = \frac{dN_{tot}}{d\hat{C}_R d\hat{z} d\hat{\mathcal{H}}_i} P(\hat{\lambda}, \hat{g}_i | \hat{C}_R, \hat{z}, \hat{\mathcal{H}}_i). \quad (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;

$$ \log \mathcal{L} = \sum_j \log \left( \frac{dN_{tot}}{d\hat{C}_R d\hat{z} d\hat{\lambda} d\hat{g}_i d\hat{\mathcal{H}}_i} P(I|\hat{C}_R, \hat{z}, \hat{\mathcal{H}}_i)\Theta(\hat{\lambda} > 3) \right) \\ + \sum_j \log P(\hat{\lambda}, \hat{g}_i | \hat{C}_R, \hat{z}, \hat{\mathcal{H}}_i, I) \\ - \int \cdots \int \left[ \frac{dN_{tot}}{d\hat{C}_R d\hat{z} d\hat{\lambda} d\hat{g}_i d\hat{\mathcal{H}}_i} P(I|\hat{C}_R, \hat{z}, \hat{\mathcal{H}}_i)\Theta(\hat{\lambda} > 3) \right. \\ \left. \qquad \qquad \qquad d\hat{C}_R d\hat{z} d\hat{\lambda} d\hat{g}_i d\hat{\mathcal{H}}_i \right] \quad (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(\lambda, \hat{g}_i | \hat{C}_R, \hat{z}, \hat{H}_i, I) = P(\lambda | \hat{C}_R, \hat{z}, \hat{H}_i, I) P(\hat{g}_i | \hat{C}_R, \hat{z}, \hat{H}_i, I). \quad (5)$$

Thus, we can write our final log-likelihood as:

$$\begin{aligned} \log \mathcal{L} = & \sum_j \log \left( \frac{dN_{tot}}{d\hat{C}_R d\hat{z} d\hat{H}_i} P(I|\hat{C}_R, \hat{z}, \hat{H}_i) \right) \\ & - \int \cdots \int \left[ \frac{dN_{tot}}{d\hat{C}_R d\hat{z} d\hat{\lambda} d\hat{H}_i} P(I|\hat{C}_R, \hat{z}, \hat{H}_i) \Theta(\lambda > 3) \right. \\ & \left. \qquad d\hat{C}_R d\hat{z} d\hat{\lambda} d\hat{H}_i \right] \\ & + \sum_j \log P(\lambda | \hat{C}_R, \hat{z}, \hat{H}_i, I) \\ & + \sum_j \log P(\hat{g}_i | \hat{C}_R, \hat{z}, \hat{H}_i, I). \end{aligned} \quad (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 \times 10^{12} M_{\odot}$ to $5 \times 10^{15} M_{\odot}$, richness from 1 to 600, and count rate from $10^{-3}$ 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)<fnref n="8" />. 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 sub-populations 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 Sect. 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 Sect. 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.

<footnote n="8">⁸ In Bocquet et al. (2019), the authors compute the corresponding likelihood only for a fraction of systems that are not optically confirmed. This is different than our implementation presented here.</footnote>
In short, the probability density functions for every observables $\hat{O}$ is written as:

$$P(\hat{O}) = \sum_{i \in \{C, \text{AGN}, \text{RS}\}} f_i P(\hat{O}|i), \quad (7)$$

where this equation holds for all our observables $\hat{O}$ as $\hat{C}_R$, $\hat{z}$, $\hat{\lambda}$, and $\hat{g}_i$. $f_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{dN_c}{d\hat{C}_R d\hat{z} d\hat{H}_i} P(I|\hat{C}_R, \hat{z}, \hat{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<fnref n="9" />.

$$\frac{dN_c}{d\hat{C}_R d\hat{z} d\hat{H}_i} P(I|\hat{C}_R, \hat{z}, \hat{H}_i) = P(\hat{H}_i) \iiint P(\hat{z}|z) P(I|z, C_R, \hat{H}_i) \times P(\hat{C}_R|C_R) \frac{dN}{dM dV} \frac{dV}{dz} P(C_R|M, z) dM dC_R dz \quad (8)$$

where $P(\hat{H}_i)$ is the uniform probability over the sphere that a cluster exists at a specific sky position $\hat{H}_i$, $P(\hat{z}|z)$ expresses the uncertainty on our redshift measurements, provided by eROMaPPer (Kluge et al. 2024), $P(I|z, C_R, \hat{H}_i)$ is the X-ray selection function (see Sect. 3 and for more details Clerc et al. 2024), $P(\hat{C}_R|C_R)$ is the uncertainty on the X-ray count rate, measured by MBProj2D, $\frac{dN}{dM dV}$ is the halo mass function<fnref n="10" />, which directly depends on cosmology, $\frac{dV}{dz}$ is the fractional comoving volume in an infinitesimal redshift shell (see Equation (28) in Hogg 1999)<fnref n="11" />, $P(C_R|M, z)$ is the log-normal distribution around the X-ray scaling relation and can be described as:

$$P(C_R|M, z) = \mathcal{L}N(\langle C_R|M, z \rangle, \sigma_X). \quad (9)$$

We refer to Sect. 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_c = \iiint\int \frac{dN_c}{d\hat{C}_R d\hat{z} d\hat{\lambda} d\hat{H}_i} P(I|\hat{C}_R, \hat{z}, \hat{H}_i) \Theta(\hat{\lambda} > 3) d\hat{C}_R d\hat{z} d\hat{\lambda} d\hat{H}_i \quad (10)$$

where several simplifications are applied. First, $P(\hat{z}|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:

$$P(\hat{z}|z) = (1 - c_z) N(b_z z, \sigma_z(1 + z)) + c_z N(b_z z + c_{\text{shift},z}, \sigma_z(1 + z)), \quad (11)$$

where $c_z$ is the fraction of objects for which we measure a redshift shifted by some constant $c_{\text{shift},z}$, $b_z$ is a systematic bias in our redshift estimate, and $\sigma_z$ is the relative error on the measured redshift. We can, therefore, express $g(z, \hat{H}_i) = \int P(\hat{z}|z)w(\hat{z}, \hat{H}_i)d\hat{z}$ as an analytic expression for each sky position, with $w(\hat{z}, \hat{H}_i)$ being the filter function that is redshift and sky position-dependent, as follows:

$$w(\hat{z}, \mathcal{H}_i) = \begin{cases} 1, & \text{if } (\hat{z} > 0.1) \& (\hat{z} < \min(0.8, \hat{z}_{\text{max}}(\hat{H}_i))), \\ 0, & \text{otherwise,} \end{cases} \quad (12)$$

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

Then, $P(\hat{\lambda}|\lambda)$, the probability of observing a cluster with measured richness $\hat{\lambda}$ given the true richness $\lambda$, 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}|\lambda) = N(\lambda, \sqrt{\lambda}). \quad (13)$$

Thus, we can express the optical completeness term $f(\lambda) = \int P(\hat{\lambda}|\lambda)\Theta(\hat{\lambda} > 3)d\hat{\lambda}$ term analytically<fnref n="12" />. 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|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:

$$N_c = \iiint\int P(I|z, C_R, \hat{H}_i)g(z, \hat{H}_i)f(\lambda)P(C_R, \lambda|M, z)P(\hat{H}_i) \times \frac{dN}{dM dV} \frac{dV}{dz} dM dC_R dz d\lambda d\hat{H}_i. \quad (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).

<footnote n="12">¹² In short, for redshift and richness, for which the distributions are normal distributions integrated in some bounded range, the integral is expressed as an error function.</footnote>

<footnote n="9">⁹ Technically also over the true cluster position; however, we assume that cluster positions are measured precisely, e.g., $P(\hat{H}_i|\hat{H}_i) = \delta(\hat{H}_i - \hat{H}_i)$, where $\delta$ represents the Dirac's delta function. Indeed, the relevant observing conditions do not vary within the astrometric uncertainty on the cluster position.</footnote>

<footnote n="10">¹⁰ $\frac{dN}{dM dV} = f(\sigma) \frac{\rho}{M} \frac{d \log(\sigma^{-1})}{d \log M}$ where $f(\sigma)$ is a function of $\sigma(R, z)$, which is defined as $\sigma^2(R, z) = \frac{1}{2\pi^2} \int k^2 P(k, z) |W(kR)|^2 dk$</footnote>

<footnote n="11">¹¹ We point out that $\frac{dN}{dM dV}$ is the halo mass function as a function of comoving volume, while $\frac{dN}{dM dV dz}$ is the halo mass function as a function of redshift. Cosmology calculators, like pyccl, give the user the former, and the latter is left to calculate for the user.</footnote>
### 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}_i | \hat{C}_R, \hat{z}, \hat{H}_i, I, C) = \frac{P(\hat{g}_i, \hat{C}_R, \hat{z}, \hat{H}_i, I, C)}{P(\hat{C}_R, \hat{z}, \hat{H}_i, I, C)} = \frac{N_{WL}}{D_{WL}}, \quad (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 Sect. 4.3, we marginalize the expressions $N_{WL}$ and $D_{WL}$, similarly to Equation (8), as follows:

$$N_{WL} = P(\hat{H}_i) \iint \int P(I | C_R, \hat{z}, \hat{H}_i) P(\hat{C}_R | C_R) P(M, \hat{z}) \times P(\hat{g}_i | M_{WL}, \hat{z}) P(M_{WL}, C_R | M, \hat{z}) dM dC_R dM_{WL} \quad (16)$$

and

$$D_{WL} = P(\hat{H}_i) \iint P(I | C_R, \hat{z}, \hat{H}_i) P(\hat{C}_R | C_R) P(M, \hat{z}) \times P(C_R | M, \hat{z}) dM dC_R. \quad (17)$$

We note that the sky position probabilities, $P(\hat{H}_i)$, cancel out in the ratio of $N_{WL} / D_{WL}$ in the weak lensing mass calibration term. $P(\hat{g}_i | M_{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 | C_R)$, the selection function $P(I | C_R, \hat{z}, \hat{H}_i)$, the normalized halo mass function, $P(M, \hat{z})$. However, since the normalization factor (needed to compute $P(M, z)$ from $\frac{dN}{dM d\hat{z}}$) is present both at the numerator and denominator, it cancels itself out, and the halo mass function, $\frac{dN}{dM 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_{WL}, C_R | 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 quantities, is a bi-variate log-normal distribution $BN(\bar{\mu}, \bar{\Sigma})$, where

$$\bar{\mu} = [\langle \log C_R | M, z \rangle, \langle \log M_{WL} | M, z \rangle], \quad (18)$$

$$\bar{\Sigma} = \begin{bmatrix} \sigma_{C_R}^2 & \rho_{M_{WL}, C_R} \sigma_{C_R} \sigma_{M_{WL}} \\ \rho_{M_{WL}, C_R} \sigma_{C_R} \sigma_{M_{WL}} & \sigma_{M_{WL}}^2 \end{bmatrix}, \quad (19)$$

where $\langle \log C_R | M, z \rangle$ and $\langle \log M_{WL} | M, z \rangle$ are the scaling relations between true cluster mass and count rate and weak lensing mass, respectively, see Sect. 6 for further details, and $\rho_{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} | \hat{C}_R, \hat{z}, \hat{H}_i, I, C) = \frac{P(\hat{\lambda}, \hat{C}_R, \hat{z}, \hat{H}_i, I, C)}{P(\hat{C}_R, \hat{z}, \hat{H}_i, I, C)} = \frac{N_{opt}}{D_{opt}}, \quad (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:

$$N_{opt} = P(\hat{H}_i) \iiint P(I | C_R, \hat{z}, \hat{H}_i) P(\hat{C}_R | C_R) P(M, \hat{z}) \times P(\hat{\lambda} | \lambda) \Theta(\hat{\lambda} > 3) P(\lambda, C_R | M, \hat{z}) dM dC_R d\lambda, \quad (21)$$

$$D_{opt} = P(\hat{H}_i) \iint P(I | C_R, \hat{z}, \hat{H}_i) P(\hat{C}_R | C_R) P(M, \hat{z}) \times P(C_R | M, \hat{z}) dM dC_R. \quad (22)$$

Here, in particular, $P(\hat{\lambda} | \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 | C_R, \hat{z}, \hat{H}_i)$ is the X-ray selection function, $\Theta(\hat{\lambda} > 3)$ is the optical selection function$^{13}$, $P(M, \hat{z})$ is the normalized mass function, $P(\hat{C}_R | C_R)$ is the error on the measurement of the count rate, $P(\lambda, C_R | M, \hat{z})$ is a bivariate log-normal distribution $MN(\bar{\mu}, \bar{\Sigma})$ where

$$\bar{\mu} = [\langle \log C_R | M, \hat{z} \rangle, \langle \log \lambda | M, z \rangle], \quad (23)$$

$$\bar{\Sigma} = \begin{bmatrix} \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{bmatrix}, \quad (24)$$

where $\langle \log C_R | M, z \rangle$ and $\langle \log \lambda | M, z \rangle$ are the scaling relations between true mass, count rate, and optical richness, details are given in Sect. 6, while $\rho_{\lambda, C_R}$ is the intrinsic correlation between richness and count rate and where $\sigma_{\lambda}$ is richness intrinsic scatter, and $\sigma_{P}$ is a Poisson term, defined as in Costanzi et al. (2019a) by $\langle \lambda | M, z \rangle^{-1}$. 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 $\Omega_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 around Planck CMB results (Planck Collaboration VI 2020). For $H_0$, against which we are not sensitive, we use the Planck CMB $\Lambda$CDM results with a tight Gaussian prior set around them (Planck Collaboration VI 2020). Given the so-called $H_0$ tension, we have checked that using a prior from distance ladder method

<footnote>$^{13}$ Even though the optical selection function has no effect since this term is computed only for clusters that are detected, we still leave it inside the equation to recall that the optical uncertainty term $P(\hat{\lambda} | \lambda)$ is not normalized from $-\infty$ to $+\infty$, but from 3 to $\infty$, where 3 is our implicit optical cut.</footnote>
<table><caption>**Table 1.** Priors used on the free parameters in the fitting, divided into several parameter types for clarity.</caption><thead><tr><th>Parameter</th><th>Units</th><th>Description</th><th>Prior</th></tr></thead><tbody><tr><td colspan="4"><strong>Cosmology</strong></td></tr><tr><td>Ω<sub>m</sub></td><td>–</td><td>Mean matter density at present time</td><td>U(0.05, 0.95)</td></tr><tr><td>log<sub>10</sub> A<sub>s</sub></td><td>–</td><td>Amplitude of the primordial power spectrum</td><td>U(-10, -8)</td></tr><tr><td>H<sub>0</sub></td><td>km/s/Mpc</td><td>Hubble expansion rate at present time</td><td>N(67.77, 0.6)</td></tr><tr><td>Ω<sub>b</sub></td><td>–</td><td>Mean baryon density at present time</td><td>U(0.046, 0.052)</td></tr><tr><td>n<sub>s</sub></td><td>–</td><td>Spectra index of the primordial power spectrum</td><td>U(0.92, 1.0)</td></tr><tr><td>w<sub>0</sub></td><td>–</td><td>Dark energy equation of state. Fixed to -1 in ΛCDM</td><td>U(-2.5, -0.33)</td></tr><tr><td>Σ m<sub>v</sub></td><td>eV</td><td>Summed neutrino masses. Fixed to 0 eV in ΛCDM</td><td>U(0, 1)</td></tr><tr><td colspan="4"><strong>X-ray scaling relation</strong></td></tr><tr><td>A<sub>X</sub></td><td>–</td><td>Normalization of the M - C<sub>R</sub> scaling relation</td><td>U(0.01, 3)</td></tr><tr><td>B<sub>X</sub></td><td>–</td><td>Mass slope of the M - C<sub>R</sub> scaling relation</td><td>U(0.1, 5)</td></tr><tr><td>D<sub>X</sub></td><td>–</td><td>Luminosity distance evolution of the M - C<sub>R</sub> scaling relation</td><td>Fixed to -2</td></tr><tr><td>E<sub>X</sub></td><td>–</td><td>Scale factor evolution of the M - C<sub>R</sub> scaling relation</td><td>Fixed to 2</td></tr><tr><td>F<sub>X</sub></td><td>–</td><td>Redshift evolution of the mass slope of the M - C<sub>R</sub> scaling relation</td><td>U(-5, 5)</td></tr><tr><td>G<sub>X</sub></td><td>–</td><td>Redshift evolution of the normalization of the M - C<sub>R</sub> scaling relation</td><td>U(-5, 5)</td></tr><tr><td>σ<sub>X</sub></td><td>–</td><td>Intrinsic scatter of the M - C<sub>R</sub> scaling relation</td><td>U(0.05, 2)</td></tr><tr><td colspan="4"><strong>Weak lensing mass calibration</strong></td></tr><tr><td>A<sub>WL</sub></td><td>–</td><td>Scatter in the weak lensing bias from the first principal component</td><td>N(0, 1)</td></tr><tr><td>B<sub>WL</sub></td><td>–</td><td>Scatter in the weak lensing bias from the second principal component</td><td>N(0, 1)</td></tr><tr><td>C<sub>WL</sub></td><td>–</td><td>Standardize mass slope of the weak lensing bias</td><td>N(0, 1)</td></tr><tr><td>D<sub>WL</sub></td><td>–</td><td>Redshift dependent intrinsic scatter in the weak lensing bias</td><td>N(0, 1)</td></tr><tr><td>ρ<sub>MWL,C<sub>R</sub></sub></td><td>–</td><td>Intrinsic correlation between weak lensing mass and count rate</td><td>U(-0.9, 0.9)</td></tr><tr><td colspan="4"><strong>Richness mass calibration</strong></td></tr><tr><td>log A<sub>λ</sub></td><td>–</td><td>Normalization of the M - λ scaling relation</td><td>U(1.4, 6)</td></tr><tr><td>B<sub>λ</sub></td><td>–</td><td>Mass slope of the M - λ scaling relation</td><td>U(0, 2)</td></tr><tr><td>C<sub>λ</sub></td><td>–</td><td>Redshift evolution of the normalization of the M - λ scaling relation</td><td>U(-2, 2)</td></tr><tr><td>D<sub>λ</sub></td><td>–</td><td>Redshift evolution of the mass slope of the M - λ scaling relation</td><td>U(-2, 2)</td></tr><tr><td>σ<sub>λ</sub></td><td>–</td><td>Intrinsic scatter of the M - λ scaling relation</td><td>U(0.05, 2)</td></tr><tr><td>ρ<sub>λ,C<sub>R</sub></sub></td><td>–</td><td>Intrinsic correlation between richness and count rate</td><td>U(-0.9, 0.9)</td></tr><tr><td colspan="4"><strong>Contamination modeling</strong></td></tr><tr><td>f<sub>AGN</sub></td><td>–</td><td>Fraction of AGN contaminants in the extended source sample</td><td>U(0, 0.1)</td></tr><tr><td>f<sub>RS</sub></td><td>–</td><td>Fraction of RS contaminants in the extended source sample</td><td>U(0, 0.15)</td></tr><tr><td colspan="4"><strong>Redshift uncertainty</strong></td></tr><tr><td>σ<sub>z</sub></td><td>–</td><td>Relative error on the measured redshift</td><td>TN(0.0050, 0.0011, 0, 1)</td></tr><tr><td>b<sub>z</sub></td><td>–</td><td>Systematic bias in our redshift estimate</td><td>N(1.005, 0.037)</td></tr><tr><td>c<sub>z</sub></td><td>–</td><td>Fraction of objects for which we measure a shifted redshift</td><td>TN(0.0013, 0.0010, 0, 1)</td></tr><tr><td>c<sub>shift,z</sub></td><td>–</td><td>Amount of redshift shift for c<sub>z</sub> fraction of objects</td><td>U(-0.1, 0.1)</td></tr></tbody></table>

**Notes.** With $U(\text{min}, \text{max})$ we indicate a uniform distribution between “min” and “max”. With $N(\mu, \sigma)$ we indicate a normal distribution centered on $\mu$ and with standard deviation $\sigma$. With $TN(\mu, \sigma, \text{min}, \text{max})$ we indicate a truncated normal distribution centered on $\mu$ and with standard deviation $\sigma$ bounded between “min” and “max”.

($H_0 \sim N(73.04, 1.04)$ km s⁻¹ Mpc⁻¹, Riess et al. 2022) or a flat prior ($H_0 \sim 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_{500} > 10^{12} M_{\odot}$, 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 rela-
tions described in Sect. 6 with fiducial values from our fits.
Observable quantities were then drawn according to the expected
statistical distribution. As described in Sect. 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 respec-
tive 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 probabil-
ity distribution, we included the contaminants, RSs and AGNs
(see Sect. 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 param-
eters, σ8 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<fnref n="14" /> described in Sect. 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 σ8 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 simula-
tions 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 Sect. 6.2 and Grandis et al. 2024), we fur-
ther 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 Equa-
tion (28). We find that a decrease in the masses as found by
Sommer et al. (2024) indeed decreases the values of Ωm and σ8
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 pho-
tometric information from the LS DR10-South grz bands; see
Sect. 2.1 (Kluge et al. 2024). This means that their uncertainties
are significantly larger than spectroscopic redshift measure-
ments. This potentially could cause a bias on the cosmological
parameters mainly as a result of geometrical effects. This is
because $\frac{dN}{dM\,dV}$ is a weak function of redshift, especially in the
redshift interval of interest, but the volume $\frac{dV}{dz}$ and cosmologi-
cal 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 σ8 measurements are unrealistic. For example, a systematic
error of 10% on photo-zs, namely, $\sigma_z \approx 0.1$, affects Ωm and σ8
by about 3 and 2% respectively, and a trend of catastrophic fail-
ures on photo-z with occurrence of 10% results in lower values
of Ωm by 4% and higher values of σ8 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 Equa-
tion (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 σ8 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

<footnote n="14">¹⁴ We note that at this point, the “true” selection function is still
unknown to us as part of our blinding strategy.</footnote>
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 Sect. 3.4, only one corresponding to the true selection model, namely, with the $L_{ext}$ thresholds that match the construction of the cosmology subsample. These selection functions were tested on the mocks described in Sect. 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 $\Lambda$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 $\Lambda$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 Sect. 5.1). We then unblinded the results and completed the rest of the draft of the paper, extending our analysis beyond the $\Lambda$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:

$$ \left\langle \log \frac{C_R}{C_{R,p}} \middle| M, z \right\rangle = \log A_X + b_X(z) \log \frac{M}{M_p} + e_X(z). \quad (25) $$

We used the $C_{R,p} = 0.1$ cts / s as the pivot value for the count rate, $M_p = 2 \times 10^{14}$ M$_\odot$, 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

<table><caption>**Table 2.** Best-fit parameters of the scaling relations between WL shear and X-ray count-rate and optical richness.</caption><thead><tr><th>Scaling relation</th><th>Parameter</th><th>Best-fit value</th></tr></thead><tbody><tr><td rowspan="5"><i>C</i><sub>R</sub> − <i>M</i></td><td><i>A</i><sub>X</sub></td><td>0.64<sup>+0.04</sup><sub>−0.06</sub></td></tr><tr><td><i>B</i><sub>X</sub></td><td>1.38<sup>+0.03</sup><sub>−0.04</sub></td></tr><tr><td><i>F</i><sub>X</sub></td><td>−0.33 ± 0.12</td></tr><tr><td><i>G</i><sub>X</sub></td><td>0.29<sup>+0.12</sup><sub>−0.13</sub></td></tr><tr><td>σ<sub>X</sub></td><td>0.98<sup>+0.05</sup><sub>−0.04</sub></td></tr><tr><td rowspan="5"><i>λ</i> − <i>M</i></td><td><i>A</i><sub>λ</sub></td><td>3.08 ± 0.05</td></tr><tr><td><i>B</i><sub>λ</sub></td><td>1.01<sup>+0.02</sup><sub>−0.03</sub></td></tr><tr><td><i>C</i><sub>λ</sub></td><td>−0.49<sup>+0.14</sup><sub>−0.15</sub></td></tr><tr><td><i>D</i><sub>λ</sub></td><td>0.30 ± 0.11</td></tr><tr><td>σ<sub>λ</sub></td><td>0.32<sup>+0.04</sup><sub>−0.03</sub></td></tr></tbody></table>

**Notes.** The reported uncertainties are the 68th percentiles. The posteriors of the best-fit parameters are shown in Fig. 5.

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}. \quad (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(z_p)} + E_X \log \frac{E(z)}{E(z_p)} + G_X \log \frac{1+z}{1+z_p}, \quad (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, $\sigma_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.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Posteriors on the $C_R - M$ and $\lambda - 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.</figcaption></figure>

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.04}_{-0.06}$. The reported normalization parameters are $A_x = 1.24^{+0.22}_{-0.19}$ and $A_x = 1.33^{+0.26}_{-0.20}$ in Chiu et al. (2022, 2023) in $\sim 3\sigma$ 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$^2$, 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.03}_{-0.04}$. 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.17}_{-0.14}$ in Chiu et al. (2022) and that of $B_x = 1.86^{+0.20}_{-0.15}$ 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 \pm 0.14$ in Vikhlinin et al. (2009), $1.62 \pm 0.12$ in Pratt et al. (2009), $1.34 \pm 0.18$ in Eckmiller et al. (2011), $1.66 \pm 0.22$ in Lovisari et al. (2015), and $1.35 \pm 0.07$ in Schellenberger & Reiprich (2017b). Our slopes are consistent with these findings at the $1\sigma$ level. Higher slopes are reported in SZ-selected massive cluster samples covering a wider redshift range ($B_{L_x} = 1.89 \pm 0.18$, Bulbul et al. 2019). When the authors excised the core, they found shallower slopes $1.60^{+0.15}_{-0.16}$, 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.05}_{-0.04}$, in our sample is larger than the reported values from $C_R - \bar{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 Sect. 7).

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

$$ \left\langle \log \frac{M_{\mathrm{WL}}}{M_p} \middle| M, z \right\rangle = b(z) + b_M \log \left( \frac{M}{M_p} \right). \quad (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}_i|M_{\mathrm{WL}}, \hat{z})$ described in Sect. 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 $\mu_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}$. 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)$, $PC_{1,i}(z)$, $PC_{2,i}(z)$, $\mu_{bM,i}$, and $\sigma_{bM,i}$ with $i \in \{\mathrm{DES}, \mathrm{KiDS}, \mathrm{HSC}\}$. To emulate the scatter in the WL bias simulations, we introduce three additional random normal parameters, $A_{\mathrm{WL}}$, $B_{\mathrm{WL}}$, and $C_{\mathrm{WL}}$, which are independent of the surveys employed. In practice, $b(z)$ and $b_M$ can be expressed as:

$$ \begin{aligned} b(z) &= \mu_{b,i}(z) + A_{\mathrm{WL}} PC_{1,i}(z) + B_{\mathrm{WL}} PC_{2,i}(z), \\ b_M &= \mu_{bM,i} + C_{\mathrm{WL}} \sigma_{bM,i}, \end{aligned} \quad (29) $$

where $\mu_{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{WL}}}^2 = s(z) + s_M \log \left( \frac{M}{M_p} \right), \quad (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 $\mu_{s,i}(z)$ and $PC_i(z)$ and a random normal parameter $D_{\mathrm{WL}}$ and express $s(z)$ as:

$$ s(z) = \mu_{s,i}(z) + D_{\mathrm{WL}} PC_i(z), \quad (31) $$

where the average scatter for each survey and principal component ($\mu_{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:

$$ \langle \log \lambda | M, z \rangle = \log A_{\lambda} + b_{\lambda}(z) \log \left( \frac{M}{M_{p}} \right) + C_{\lambda} \log \left( \frac{1+z}{1+z_{p}} \right), \quad (32) $$

while the slope of the cluster mass term $b_\lambda(z)$ can be expanded as:

$$ b_{\lambda}(z) = B_{\lambda} + D_{\lambda} \log \left( \frac{1+z}{1+z_{p}} \right), \quad (33) $$

where $A_\lambda$ is the normalization of the scaling relation, $B_\lambda$ is the mass slope, $C_\lambda$ is the redshift evolution, and $D_\lambda$ 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_\lambda = 1.01^{+0.02}_{-0.03}$, is consistent with the majority of reported findings in the literature within $1\sigma$ level, e.g., $B_\lambda = 0.99^{+0.06}_{-0.07} \pm 0.04$ and $B_\lambda = 1.07^{+0.22}_{-0.20}$ in the CODEX sample (Capasso et al. 2019;

Ider Chitham et al. 2020, respectively), $B_\lambda = 1.02 \pm 0.08$ in SPT sample (Bleem et al. 2020), $B_\lambda = 1.028^{+0.043}_{-0.037} \pm 0.04$ and $B_\lambda = 1.03 \pm 0.10$ in the DES Y1 and SPT samples Costanzi et al. (2021); Grandis et al. (2021b), $B_\lambda = 1.05^{+0.13}_{-0.12}$ in Chiu et al. (2023) in the eFEDS sample within the stated uncertainties. Our measurements are $8\sigma$ away from the value reported by the DES collaboration on Y1 cosmology analysis, where they find $B_\lambda = 0.737 \pm 0.028 \pm 0.004$ (McClintock et al. 2019; Dark Energy Survey Collaboration 2020)$^{15}$. 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 $\lambda - M$ scaling relation of $0.3 \pm 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_\lambda = -0.49^{+0.14}_{-0.15}$ 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 $\lambda > 3$ ($-0.26^{+0.23}_{-0.24}$, Capasso et al. 2019); however, their results are consistent with no-redshift evolution, while ours indicate a redshift evolution at $3\sigma$ 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, $\lambda > 20$ ($-1.00^{+0.49}_{-0.56}$; 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 \pm 0.30 \pm 0.06$ for a sample with richness cut of $\lambda > 40$. In contrast, the DES/SPT analyses reported in Costanzi et al. (2021) $0.95 \pm 0.3$ is in tension with our results at the $>5\sigma$ level when the SPT number counts are included in the analysis, while consistent with our finding at $<2\sigma$ level, $-0.02 \pm 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, $\lambda_{grz} = S_{\text{norm}}\lambda_{griz} = S'_{\text{norm}}\lambda_{grzw1}$, without any obvious dependence

<footnote>$^{15}$ In the paper they reported the inverse of this number: $\frac{1}{B_\lambda} = 1.356 \pm 0.051 \pm 0.008$.</footnote>
on redshift. This result is also supported by the reported com-
parison 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 mea-
sured 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 $\lambda_{grz} = 0.79\lambda_{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 σ_8, and the reionization
optical depth τ. The other cosmological parameters can be
derived from this baseline set of parameters, for instance, the
age of the Universe m₀, the critical density of the Universe m_c,
and a combination of Ω_m with σ_8 along the principal direction
of the degeneracy between the two, S_8. In our baseline model,
the flat ΛCDM, we fit for Ω_m, σ_8 <fnref n="16" />, 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 con-
straints 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 cos-
mology models is not a trivial process in our analysis. The
cosmology likelihood is not a simple distribution but com-
prises various components following different distributions. For
instance, the component for the cluster number counts in the
global likelihood (Sect. 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 distribu-
tion. 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<fnref n="17" /> 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 Sect. 5.1 for further details. Overall, the total simu-
lated number of objects in each bin of the observables is relevant
in assessing the goodness of fit. In Fig. 6, we provide the distri-
bution 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 pur-
poses, 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 dis-
tribution, χ² = 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 χ²_DES = 173.2⁺⁴³.²⁻²⁷.²
(150 d.o.f.), χ²_HSC = 140.0⁺³.²⁻⁶.² (70 d.o.f.), and χ²_KiDS = 60.2⁺⁹.⁷⁻⁵.⁶
(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 pro-
file, 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 scal-
ing 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

<footnote n="16">¹⁶ To be more precise, in the code, we do not fit σ₈ but Aₛ, the normal-
ization of the power spectrum, and σ₈ is a derived from the constraints
on Aₛ, see Appendix E.</footnote>

<footnote n="17">¹⁷ We use the ultranest (Buchner 2021) provided weights to ensure that the points sampled from the posterior distribution provide a non-negligible weight to the sampled posterior.</footnote>
<figure><img src="image_7.png" /><figcaption>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.</figcaption></figure>

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.

## 7.2. Constraints on the ΛCDM cosmology

Cluster number counts can provide powerful constraints on the cosmological parameters, especially on $\Omega_m$ and $\sigma_8$, 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, $\Omega_m$ and $\sigma_8$, and $S_8$ are:

$$ \begin{aligned} \Omega_m &= 0.29^{+0.01}_{-0.02}, \\ \sigma_8 &= 0.88 \pm 0.02, \\ S_8 &= 0.86 \pm 0.01. \end{aligned} \qquad (34) $$

The posterior distributions of $\Omega_m$ and $\sigma_8$ 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\sigma$ 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.04}_{-0.06}$ consistent with our result, whereas their best-fit $\Omega_m = 0.179^{+0.031}_{-0.038}$ 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 $\sigma_8$ and $\Omega_m$ parameter measurements from the eROSITA cluster abundances agree well within $2\sigma$ 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
<figure><img src="image_8.png" /><figcaption>**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).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Posterior distributions on the parameters $\Omega_m$ and $\sigma_8$ from the $\Lambda$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.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 9. Posterior distributions on $\Omega_m$ and $\sigma_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).</figcaption></figure>

(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_\lambda$) and the redshift evolution ($G_X$ and $C_\lambda$) of the scaling relations between X-ray observable and richness are correlated with $\Omega_m$. In contrast, these terms are anti-correlated with $\sigma_8$. Interestingly, $S_8$ shows a strong correlation with the amplitude of the scaling relations ($A_X$ and $A_\lambda$), but we do not observe any correlation between $S_8$ 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 = \bar{\mu}^T C^{-1} \bar{\mu}$; and 3) we computed the survival function<fnref n="18" /> of the $\chi^2$ 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 $\Omega_m$ and $\sigma_8$ 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\sigma$ confidence level.

To break the degeneracy between $\Omega_m$ and $\sigma_8$ 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{aligned} \Omega_m &= 0.32 \pm 0.01, \\ \sigma_8 &= 0.82 \pm 0.01, \\ S_8 &= 0.85 \pm 0.01. \end{aligned} \qquad (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 $\Omega_m$ and to a lower value for $\sigma_8$, a significant improvement on our constraints on the $\Omega_m$ and $\sigma_8$ parameter plane. This is because the Planck CMB 2020 constraints on $\Omega_m$ and $\sigma_8$ are tight and degenerate along a direction orthogonal to the $S_8$ 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_8$, which combines the constraints on parameters $\Omega_m$ and $\sigma_8$. Our data constrains the parameter $S_8$ with high precision due to the direction and orientation of the degeneracy between these two parameters. The best-fit $\alpha_{best}$ value that minimizes the standard deviation of $\sigma_8 \times (\Omega_m/0.3)^\alpha$ along the direction of our degeneracy, is $\alpha_{best} = 0.42 \pm 0.02$. Our best-fit value is close to 0.5, the canonical value of $\alpha$ used in the definition of $S_8$ in Planck Collaboration VI (2020). The cluster cosmology experiments, e.g., Bocquet et al. (2019) reports $\alpha_{best} = 0.2$ for the SPT-SZ sample, while Chiu et al. (2023) finds $\alpha_{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 $\alpha_{best}$ value of 0.45, consistent with our results Hikage et al. (2019). Regarding the $S_8$ constraints, we obtain a degeneracy between $\sigma_8$ and $\Omega_m$, which is almost perfectly aligned with the slope in the $S_8$ definition. This makes the

<footnote n="18">¹⁸ The survival function complements one of the cumulative distribution functions.</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Comparison of the constraints on $S_8 = \sigma_8 \times \sqrt{\Omega_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\sigma$ level. The vertical line and grey shaded areas represent our measurement's location and $1-2\sigma$ confidence level.</figcaption></figure>

constraints on $\sigma_8$ tighter in our experiment than other surveys. Similarly, the recent SPT constraints show the slope $\alpha_{best} \sim 0.25$ (Bocquet et al. 2024), and indeed their optimal parameter combination $\sigma_8 \times (\Omega_m/0.3)^{0.25}$ measure a relative uncertainty very close to the uncertainty we obtain on $S_8$ in the eRASS1 sample. The degeneracy reported in other experiments does not have a direction in line with the $S_8$ definition.

Figure 10 shows our constraint on the $S_8$ 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\sigma$ confidence interval (Planck Collaboration VI 2020; Garrel et al. 2022). However, it is in significant tension ($\ge 3\sigma$ 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_8$ 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_8$ 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 $\Lambda$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 $\sigma_8$ from $0.88 \pm 0.02$ to $0.86 \pm 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\sigma$ intervals with the Tinker et al. (2008) results. The change in the halo mass function does not significantly impact the parameter $\Omega_m = 0.30^{+0.01}_{-0.02}$ 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).

### 7.3. Constraints on the $\nu$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_\tau > m_\mu > m_e$. The ordering is swapped in the inverted mass hierarchy model such that $m_\tau < m_\mu < 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 $\sum m_\nu > 0.059$ eV for normal mass hierarchy and $\sum m_\nu > 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 $\sum m_\nu < 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 $\sum m_\nu < 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 \text{ eV}, 1 \text{ eV})$. Including massive neutrinos in cosmological analysis requires a fine-tuned interaction
<figure><img src="image_12.png" /><figcaption>Fig. 11. Best-fit values for νCDM 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 $\sum m_ν < 0.43$eV from cluster number counts alone.</figcaption></figure>

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<fnref n="19" /> 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{aligned}
\Omega_m &= 0.28^{+0.02}_{-0.02}, \\
\sigma_8 &= 0.88 \pm 0.03, \\
S_8 &= 0.85 \pm 0.02, \\
\sum m_ν < 0.43\text{eV} \quad (95\% \text{ CL}).
\end{aligned} 
\qquad (36) 
$$

These results are visualized in Fig. 11. The upper limit to the sum of neutrino masses is $\sum m_ν < 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 $\Omega_m$ and $\sigma_8$ remain statistically consistent when $\sum m_ν$ is allowed to be non-zero, see Fig. B.4. We find that eRASS1 and Planck CMB νCDM parameters are consistent at the 2.0$\sigma$ 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 $\Lambda$CDM analysis, we combine our results with the Planck CMB constraints to break the degeneracy between $\Omega_m$ and $\sigma_8$. We obtain:

$$ 
\begin{aligned}
\Omega_m &= 0.31 \pm 0.01, \\
\sigma_8 &= 0.83 \pm 0.01, \\
S_8 &= 0.84 \pm 0.01, \\
\sum m_ν < 0.14 \text{ eV} \quad (95\% \text{ CL}).
\end{aligned} 
\qquad (37) 
$$

Simultaneously fitting our measurements with Planck CMB 2020 likelihood chains yield an upper limit of $\sum m_ν < 0.14$ eV, consistent with the results in the literature results with marginally consistent with the inverted mass hierarchy model, which requires $\sum 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

<footnote n="19">¹⁹ https://github.com/LSSTDESC/CCL/tree/master</footnote>
<figure><img src="image_13.png" /><figcaption>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 3x2 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.</figcaption></figure>

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_\nu = 0.09^{+0.04}_{-0.02}$ eV; while assuming the inverted mass hierarchy model, we obtain summed masses of $\sum m_\nu = 0.12^{+0.03}_{-0.02}$ 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.016}_{-0.010}$ eV (68% confidence intervals) or with an upper limit of 0.044 eV.

## 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/\rho$) and assume that its density evolves with $\rho \propto a^{-3(1+w)}$ in an isotropic Universe. Previous dark energy equation state measurements point to $w \sim -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 ($\sum m_\nu = 0$ eV), and we allow the w to be free with a uniform prior $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{aligned}
\Omega_m &= 0.28 \pm 0.02, \\
\sigma_8 &= 0.88 \pm 0.02, \\
S_8 &= 0.86 \pm 0.02, \\
w &= -1.12 \pm 0.12.
\end{aligned} 
\qquad (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 3x2 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{aligned}
\Omega_m &= 0.30 \pm 0.01, \\
\sigma_8 &= 0.87 \pm 0.02, \\
S_8 &= 0.87 \pm 0.01, \\
w &= -0.95^{+0.05}_{-0.04}.
\end{aligned} 
\qquad (39) 
$$

Compared with our results, the parameters Ω_m, σ_8, and S_8 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.05-0.04, 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
<figure><img src="image_14.png" /><figcaption>Fig. 13. Correlation between w and the X-ray count rate to mass scaling relation normalization and redshift evolution. The correlation for the wCDM case is shown in orange and blue for the vwCDM case.</figcaption></figure>

parameters significantly. We measure;

$$ 
\begin{aligned}
\Omega_m &= 0.31 \pm 0.01, \\
\sigma_8 &= 0.84 \pm 0.01, \\
S_8 &= 0.85 \pm 0.01, \\
w &= -1.06 \pm 0.03.
\end{aligned} 
\qquad (40) 
$$

When combined with the geometrical probes, such as BAO, SNe Ia, and CMB, our data constrain the best-fit value of the dark energy equation of state to be $w = -1.06 \pm 0.03$ with 3% precision, thus improving the eROSITA only results by a factor of ~3, and consistent with $w = -1$ at $2\sigma$ level. The shift on the best-fit eROSITA parameters is less than $1\sigma$. In particular, $\Omega_m$ is unchanged, $\sigma_8$ and $S_8$ shifts to a slightly lower value.

## 7.5. Constraints on the vwCDM cosmology

In connection with the previous sections, to test the effect of the neutrino masses on the dark energy equation of state, we allow both w and $\Sigma m_v$ to vary in our fits. The priors are set to uniform distributions for $U(-2, -0.33)$ on w, and $U(0 \text{ eV}, 1 \text{ eV})$ for $\Sigma m_v$. We find:

$$ 
\begin{aligned}
\Omega_m &= 0.27 \pm 0.02, \\
\sigma_8 &= 0.89 \pm 0.03, \\
S_8 &= 0.84 \pm 0.02, \\
w &= -1.11 \pm 0.14, \\
\Sigma m_v < 0.44 \text{ eV} \quad (95\% \text{ CL}).
\end{aligned} 
\qquad (41) 
$$

This case combines the two previous cosmologies: as for the vwCDM and wCDM case, our results in the $\Omega_m - \sigma_8$ parameter space are unchanged. The degeneracy that w has with normalization and redshift evolution of the scaling relation remains the

same as that seen in Fig. 13. The posterior distribution of these cosmological parameters for the vwCDM model are shown in Fig. 14.

Due to the lack of public Planck CMB-only likelihood chains on the simultaneous fits on the w and $\Sigma m_v$ in Planck CMB 2020, we combine our results with the readily available CMB+lensing+BAO in Planck Collaboration XIII (2016) chains. Relying on the $1.6\sigma$ agreement with our results, we find:

$$ 
\begin{aligned}
\Omega_m &= 0.29 \pm 0.01, \\
\sigma_8 &= 0.85 \pm 0.02, \\
S_8 &= 0.83 \pm 0.01, \\
w &= -1.15 \pm 0.10, \\
\Sigma m_v < 0.37 \text{ eV} \quad (95\% \text{ CL}).
\end{aligned} 
\qquad (42) 
$$

In this case, the constraints on the w move closer to the value -1, with the summed neutrino masses with an upper limit of $\Sigma m_v < 0.37 \text{ eV}$. Compared with the other combinations with CMB results, we notice that the improvements are modest, which is expected since the cosmological constraints in Planck Collaboration XIII (2016) are weaker than what is published in the latest Planck CMB 2020 work. Our best-fit result on w is closer to the best-fit value in the SPT-SZ+Planck combination, where the authors report $w = -1.12 \pm 0.21$ (Bocquet et al. 2019). In the same work, the reported sum of neutrino masses is $\Sigma m_v = 0.50 \pm 0.24 \text{ eV}$ in the same combination of datasets, whereas we find tighter upper limits on $\Sigma m_v$ in this work.

## 7.6. Model comparison

In this paper, we fit our eRASS1 dataset with four different cosmological models, ΛCDM, vCDM, wCDM, and vwCDM. We compare the goodness of fit for these different models to assess whether the data prefer one of the cosmological models in this section. To achieve this, we use the Bayes factor (Jeffreys 1961; Trotta 2008; Buchner et al. 2014), which, in our case, becomes available with the fitting sampler we use in our analysis. The ultranest method gives the Bayesian evidence<fnref n="20" /> directly as an output. Therefore, it is straightforward to use it to compute the Bayes factor<fnref n="21" /> and Jeffreys scale to compare different cosmology fits (Jeffreys 1961). More complex models tend to be flexible since they can accommodate a broader range of possible data with a larger number of free parameters; however, that flexibility is automatically penalized by Bayesian inference in this framework and ends up giving us an objective estimator. Specifically, since a more complex model predicts that the observed data will be drawn from a more extensive range of possibilities, a model with extra free parameters will have a larger volume over which to integrate to compute the Bayesian evidence and, therefore, it may be disfavoured.

In Table 3, we show the Bayesian evidence obtained in the various models we test with the eRASS1 dataset. We find no evidence for a clear preference among the models we test since all the measured bets-fit parameters have values similar to that of the standard ΛCDM case, see Equations (34), (36), (38), and (41).

<footnote n="20">²⁰ The Bayesian evidence is defined as $Z = \int L(D|\theta)\pi(\theta)d\theta$, where $L(D|\theta)$ is the likelihood and $\pi(\theta)$ is the prior.</footnote>

<footnote n="21">²¹ The Bayes factor is the ratio of the Bayesian evidence of two different models, $B = Z_{M1}/Z_{M2}$.</footnote>
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Best-fit values for *vu*CDM cosmological model from eROSITA experiment (orange), *Planck* CMB 2015 data (blue), and combination of eROSITA and *Planck* CMB 2015 (red). The best-fit parameters are consistent with the canonical ΛCDM model.</figcaption></figure>

<table><caption>**Table 3.** Bayesian evidence for the cosmological models tested.</caption><thead><tr><th>Cosmology</th><th>Bayesian evidence</th></tr></thead><tbody><tr><td>ΛCDM</td><td>-11419.6 ± 0.7</td></tr><tr><td>ΛCDM with Σ m<sub>v</sub> = 0.06 eV</td><td>-11418.3 ± 0.5</td></tr><tr><td>vCDM</td><td>-11418.8 ± 0.8</td></tr><tr><td>wCDM</td><td>-11419.5 ± 0.5</td></tr><tr><td>vwCDM</td><td>-11421.2 ± 0.6</td></tr></tbody></table>

**Notes.** To interpret this table, we need the Jeffreys scale (Jeffreys 1961): if the difference between the two is greater than 5, then there is significant evidence to prefer the model with the largest Bayesian evidence; if the difference between the two is between 2.5 and 5, then there is weak preference for the model with the largest Bayesian evidence; otherwise, the models are statistically indistinguishable.

# 8. Discussion

## 8.1. Methodology

This paper presents the cosmological parameters obtained from the first All-Sky Survey of eROSITA by studying the abundance of clusters of galaxies detected in the survey. The tight constraints on cosmology we present in this work result from statistical power with the employment of the largest pure eROSITA-ICM selected cluster sample and the well-controlled systematics related to mass calibration and selection modeling in our analysis. Below, we discuss several significant improvements and critical aspects that allow us to constrain cosmology with cluster number counts with well-understood and characterized systematics.

**Contamination modeling.** the purity level of the cosmology subsample can reach up to 94% level after the X-ray extent likelihood $L_{ext}$ cut on the X-ray selection on the cluster sample and optical confirmation process performed as described in Bulbul et al. (2024); Kluge et al. (2024). The remainder of the contamination is modeled then through the novel mixture method (see Sect. 4.2), a forward model for eliminating the contamination in statistical samples, considering the X-ray and optical properties of the contaminants and clusters modeled through the simulations in our analysis. This effort made it possible to maintain the maximally acceptable $L_{ext}$ cut, namely, >6, without reducing the number of clusters in our sample. Figure B.5 shows the impact on our analysis of the contamination modeling. It is worth noting that the contamination also impacts the mass calibration process. When considering a population of detected clusters at a fixed count rate, the weak lensing signal from the source is partly due to contaminants along the line of sight. The contamination reduces the expected mass since they have, on average, lower tangential shear than actual halos. This depends on the count rate affecting the overall induced mass bias. This means that non-negligible contamination will also bias the slope and scatter of the scaling relations. In the present and future eROSITA cosmology analyses, the mixture model is necessary since it allows us to account for and correct these effects in our likelihood and, therefore, marginalize them out to obtain the presented precision cosmological constraints.

**Selection function.** the completeness of the eROSITA survey is estimated from the state-of-the-art realistic twin simulations of the sky (Comparat et al. 2020; Seppi et al. 2022). The simulations include the most relevant components to an
X-ray survey, namely, the detector properties, instrumental noise, and both resolved and unresolved astrophysical sources. Our selection function model constructed from these realistic simulations depends on a limited number of parameters due to a trade-off between modeling complexity and cosmology- and/or model-dependency. The reformulation of the selection function modeling into a classification problem enabled us to leverage the efficiency of Gaussian process classifiers; their robust mathematical framework and quick computational response make them well-suited for our analysis. The full construction of the selection function is presented in an accompanying paper by Clerc et al. (2024). In this analysis, we model the detection probability of clusters as a function of count rate, redshift, and sky position through the selection function. Although robust, our selection function is calibrated with a single cosmological setup using realistic eRASS1 simulations of the sky. This makes our overall analysis marginally dependent on the underlying cosmological setup. However, the cluster count rate profiles at fixed global count rate value ingested in the selection function are representative enough (Comparat et al. 2020; Seppi et al. 2022), and the choice of the X-ray observable makes this marginal dependence quite negligible. Additionally, the tests performed on the selection function indicate that the assumptions of cosmology do not significantly alter our probability distribution and our results on cosmology (see Sect. 5.1 and Clerc et al. 2024).

**Weak-lensing mass calibration and bias.** in the past decade, a cutting-edge approach has been developed for cosmology with cluster abundances utilizing weak lensing measurements in the mass calibration process to minimize the bias on cosmological parameters due to hydrostatic equilibrium assumption. The most crucial drawback of hydrostatic masses is the large impact of unaccounted nonthermal pressure support due to baryonic physics in the form of shocks, bulk motion, turbulence, and magnetic field (Dolag et al. 1999, 2005; Rasia et al. 2006; Vazza et al. 2009; Burns et al. 2010; Lau et al. 2009). On the other hand, the weak lensing signal, although not bias-free, depends on the total enclosed mass within an aperture; therefore, it is a more reliable measurement. For this reason, mass calibration using weak lensing shear data is the preferred choice in cosmology once the small biases are calibrated via numerical simulations (Hoekstra et al. 2012; von der Linden et al. 2014; McClintock et al. 2019). Triaxiality (Corless & King 2007; Meneghetti et al. 2010), substructures along the line-of-sight (Hoekstra 2001; Becker & Kravtsov 2011), and miscentering can bias weak lensing shear estimation in the three independent surveys (DES, KiDS, and HSC) we employed in this analysis. However, we obtain a strong and realistic proxy for the weak lensing bias by calibrating the shear measurements with the N-body or hydrodynamic simulations. We account for the total bias in our cosmology framework, as described in our supporting papers for weak lensing follow-up by Kleinebreil et al. (2024) and Grandis et al. (2024).

**Forward modeling approach.** in our analysis, we build a Bayesian forward modeling approach to account for the position-dependent detection chain of eROSITA. eRASS1 survey strategy is designed such that the exposure time varies with the location across the sky with significantly higher exposures close to the north and south ecliptic poles, which affects our detection efficiency (see Fig. 16 in Predehl et al. 2021). Additionally, galactic absorption from atomic and molecular hydrogen (HI4PI Collaboration 2016; Willingale et al. 2013) and eROSITA's detector and sky background (Freyberg et al. 2020) also varies

significantly with sky position. This means our ability to detect clusters and accurately measure their count rate is not uniform but position-dependent. In our work, we integrate the variation in sky position in the implementation in all components of our analysis, for instance, in the selection function while modeling the uncertainty of the observables, predicting our contamination level via mixture modeling, and the forward-modeled likelihood. Finally, Bayesian hierarchical forward modeling is the approach of choice when trying to understand a complex system based on robust theoretical considerations, which has strong predictive power in understanding the dependence of the data on input parameters. Cosmology studies with cluster abundances are well suited to employ modern advanced statistical approaches, such as the Bayesian hierarchical forward modeling approach, especially with state-of-the-art computing systems that allow for fast evaluation of complex models. One possible drawback of such modeling, compared to machine learning, for instance is the somewhat fixed relation between mass and observable through the scaling relation and its scatter around it. We have tested the impact of changing a single scatter value model with a mass and redshift-dependent scatter in Appendix F, showing that our conclusions are effectively unchanged.

## 8.2. $S_8$ Measurements

The discrepancy in the inferred value of the amplitude of matter fluctuation between this analysis and the results from cosmic shear, see Fig. 10, merits special attention. We use WL data to reconstruct our cluster sample's mass distribution accurately. Comparing this with cosmological predictions of halo mass function (Tinker et al. 2008) provides constraints on the cosmological parameters. Cosmic shear experiments use the same WL data but study its correlation function of large and intermediate scales, which is also measured with percent accuracy from the same optical data. Both approaches are thus similarly affected by the challenges in measuring galaxy shapes and photometric redshift. The key difference comes from the theoretical uncertainties that must be accounted for when predicting the WL signals for cosmological inference.

While cosmic shear can be predicted very accurately at large scales, predictions at intermediate scales of the order of 10 Mpc could be affected by several astrophysical systematics, as discussed in detail in Mandelbaum (2018) and Huterer (2023). Among these is the modeling of the intrinsic alignment of source galaxy shapes with the lensing potential, the non-linear matter power spectrum calculation, and the impact that baryon feedback and its associated hydro-dynamical modeling uncertainties have on the cosmic shear correlation functions.

These concerns do not impact the halo mass function derived from the linear fluctuation matter power spectrum, as collapsed halos retain information on the linear fluctuation at the epoch when detached from background expansion. Weak lensing measurements of galaxy clusters are also unaffected by these: the halos dominate the line of sight potential and are thus the major source of the weak lensing signal. We deliberately limit our analysis to scales smaller than $3.2/(1+z)$ Mpc/h to ensure that we measure only the clusters' potential and are thus not affected by the intermediate scale issues that affect the cosmic shear measurements.

We stress that the systematics related to the theoretical foundation of cluster number counts, cluster detection, selection effects, and mass calibration are well understood and accounted for in our analysis. Cluster number counts rely on limited
assumptions on mass function, which can be computed analyt-
ically (Press & Schechter 1974). Detection through X-rays is one
of the most reliable methods for locating clusters and for con-
structing pure cluster catalogs without any major bias because
the X-ray signal from the hot ICM is proportional to the square of
the electron number density and, therefore, the density contrast
of a galaxy cluster is enhanced, and projection effects are mini-
mized. The strict cuts we apply on the X-ray detection chain and
the utilization of their follow-up optical observations can elimi-
nate any leftover contamination in the cluster samples, reducing
the contamination to <6%. The remaining contamination is mod-
eled through our mixture method in-situ (Bulbul et al. 2024;
Kluge et al. 2024). The selection function is produced from the
end-to-end simulations of the eROSITA footprint and is well
understood and tested against the mocks (Comparat et al. 2019;
Seppi et al. 2022; Clerc et al. 2024).

Weak lensing measurements required for mass calibration of galaxy clusters share with cosmic shear analysis some observational systematics, like photo-z calibration and shape measurement (Kleinebreil et al. 2024; Grandis et al. 2024). However, cluster mass calibration through weak lensing observations is significantly less affected by the physical systematics that impact cosmic shear (intrinsic alignment of the source galaxies, the non-linear calculation of the power spectrum). For instance, the difference between WL masses calibrated on gravity-only simulations and hydro-dynamical simulations is 2% (Grandis et al. 2021a), while the difference between the matter power spectra between gravity-only and hydrodynamical simulations, but also among different hydro-dynamical simulations, is the order of 10%. We have further tested in the validation process that the effect on the mock observations from additional biases is small (<1σ on the cosmological constraints) and does not significantly impact our results from cluster counts (see Sect. 5.1). In short, weak lensing results are remarkably robust for small scales around massive halos but start to be more challenging to model accurately for intermediate scales averaged over large sky areas because of the different sources of weak lensing signals. In light of these arguments, WL-calibrated number counts of X-ray-selected clusters may well be regarded as more robust than cosmic shear analyses.

### 8.3. Constraints on the total mass of neutrinos

The effect of massive neutrinos on cluster abundance has been explored in this work. Particularly, they affect the primordial power spectrum and the assembly of large-scale structures (Costanzi et al. 2013; Zeng et al. 2019). Massive neutrinos will not cluster on scales smaller than their free-streaming length due to their high thermal velocities. With non-zero masses, neutrinos are part of the total matter field and contribute to its smoothing. The global effect is, therefore, the dumping of structure formation compared to a standard ΛCDM cosmological model with massless neutrinos.

In the literature, measurements of summed neutrino masses
from individual cosmological experiments are generally not very
tight. For example, Bocquet et al. (2019) using the SPT clus-
ter sample found $\sum m_\nu < 0.74$ eV at 95% confidence level, and
even CMB results from Planck Collaboration VI (2020) yield
$\sum m_\nu < 0.26$ eV at 95% confidence level. The tightest constraints
on the neutrino masses are obtained via combinations of multi-
ple probes. For instance, combining Lyman-$\alpha$ forest results and
CMB measurements give $\sum m_\nu < 0.15$ eV (Rossi et al. 2015),
while BAO and CMB together yields $\sum m_\nu < 0.15$ eV at 95%
confidence level (Vagnozzi et al. 2017). More recent results from

the Planck CMB and BOSS BAO measurements provide slightly
tighter upper limits of $\sum m_\nu < 0.14$ eV at 95% confidence level
(Tanseri et al. 2022).

The results in our work from eROSITA are competitive by
themselves since we obtain $\sum m_\nu < 0.43$ eV at 95% confidence
level, close to a level that is competitive with the constraints from
Planck CMB 2020. However, as in the literature, it is the com-
bination of these two state-of-the-art cosmological probes that
we obtain very tight constraints of $\sum m_\nu < 0.14$ eV at 95% confi-
dence level. We stress that although the neutrino hierarchy mass
problem remains unsolved, our findings bring us remarkably
close to a breakthrough. In fact, our result combined with Planck
CMB 2020 rule out at a 82% confidence level the inverted hier-
archy model, which requires $\sum m_\nu > 0.101$ eV. As we move
forward to the future cosmology with the exploitation of deeper
eROSITA surveys, we will significantly improve the constraints
on neutrino summed masses, possibly with higher confidence in
discriminating against different neutrino mass hierarchy models
of neutrino eigenstates. Furthermore, measuring the mass of the
lightest neutrino eigenstate (combining our results with oscilla-
tion experiments) with an order of magnitude better precision
will guide the upcoming ground- and space-based experiments
to finally pin down the neutrino masses.

## 8.4. Constraints on the dark energy equation of state

Accurate measurements of the dark energy equation of state $w$ as the ratio between pressure and density are extremely important in our effort to understand the nature of this elusive Universe component responsible for the late-time accelerated expansion of the Universe (Riess et al. 1998; Perlmutter et al. 1999). As discussed in Planck Collaboration VI (2020), we also adopt the parameterized post-Friedmann (PPF) model (Fang et al. 2008) to allow the dark energy equation of state to cross the ΛCDM value of -1. When this crossing happens (i.e., $w < -1$), dark energy is said to be considered “phantom”, and many mechanisms can explain the crossing (see, for example, Vikman 2005; Kunz & Sapone 2006; Amendola et al. 2018, for extensive theoretical background on the topic).

As shown in Sect. 7.4 eROSITA cluster abundances slightly prefer a phantom equation of state, i.e., $w < -1$, as in the case of the Planck CMB analysis (Planck Collaboration VI 2020). Still, the value we find, $w = -1.12 \pm 0.12$, is consistent at $1\sigma$ level with the canonical ΛCDM value of $-1$, and is also consistent with previously reported literature results from DES Y3 $3 \times 2$ pt ($w = -0.98_{-0.20}^{+0.32}$, Dark Energy Survey Collaboration 2022), SPT clusters ($w = -1.55 \pm 0.41$, Bocquet et al. 2019), SDSS eBOSS DR16 BAO ($w = -0.69 \pm 0.15$, Alam et al. 2021), Pantheon SNe Ia ($w = -1.09 \pm 0.22$, Scolnic et al. 2018), and Pantheon+ SNe Ia ($w = -0.90 \pm 0.14$, Brout et al. 2022). We note that differently from these other cosmological experiments, our analysis strategy using eROSITA selected clusters enables the constrain of all the main late-time parameters with high precision at the 10% percent level without the need to combine with additional cosmological experiments to break the degeneracies between the measured parameters. The addition of external cosmological experiments enables much tighter constraints on the dark energy equation of state, as we find $w = -0.95_{-0.05}^{+0.04}$ when combining with Pantheon+ SNe Ia (Brout et al. 2022), and $w = -1.06 \pm 0.03$ when Planck CMB (Planck Collaboration VI 2020) and SDSS BOSS DR16 BAO (Alam et al. 2017) are included.

Shown clearly in Fig. 13, $w$ is slightly correlated with the
amplitude of the X-ray count rate to mass scaling relation, and
it is significantly anti-correlated with the redshift evolution of
the X-ray count rate to mass scaling relation, $G_X$. This highlights the need for a new cluster sample spanning a larger redshift span, with corresponding WL mass measurements, to accurately calibrate scaling relations out to high redshifts $z \sim 1$. Deeper eROSITA surveys will expand our redshift coverage and improve count rate measurement with higher S/N data. The cosmology analyses performed with future eROSITA-based cluster catalogs will allow us to constrain the dark energy equation of state with higher precision.

## 8.5. Future improvements

The critical aspects summarized above allow us to place stringent constraints on the standard $\Lambda$CDM cosmology, the neutrino masses, and the dark energy equation of state. Although we fully account for all known systematics effects that could lead to significant uncertainties compared to the statistical precision of eRASS1 sample in our analysis, we will improve our methodology by considering the following additional systematics to realize precision cosmology with the future deeper eROSITA surveys. The halo mass function employed in this work, the semi-analytical Tinker et al. (2008) mass function, is constructed for the standard Flat-$\Lambda$CDM analyses. However, we utilize it even when the total mass of neutrinos $\Sigma m_\nu$ or dark energy equation of state $w$ parameters are varied. Our future work will use the halo mass function emulators with consistent physics to overcome this caveat.

Additionally, the Tinker et al. (2008) halo mass function is calibrated on the dark matter-only simulations and, therefore, does not include any baryonic physics. It has been shown that baryons can impact the halo mass function, especially when the low-mass and low-redshift haloes, such as galaxy groups, are included in cosmology analyses (Bocquet et al. 2016). In our analysis, the baryonic effects are negligible due to our sample's limited number of low-mass haloes. We estimate that roughly 12% of the sources in the cosmology catalog have a mass smaller than $10^{14} M_\odot$. Following the method established in the cleaning of the group catalog in Bahar et al. (2024) to identify the group size haloes, we find that our sample has a 30% spurious fraction below masses $M< 10^{14} M_\odot$, which will not be included in deeper eROSITA data. Despite including these spurious sources, we do not apply further cuts to avoid modifications in the selection function and allow our mixture model to model the contaminants. This effect is completely absorbed in the weak lensing mass bias calibrated in the hydrodynamical simulations. However, the baryonic contribution will become a critical bias as we include low-mass haloes in cosmology analysis as the depth of the eROSITA survey increases. Indeed, at low redshifts, the limiting factor for our WL mass calibration is the residual uncertainty from baryonic physics. In the future, the inclusion of low mass halos in the samples will be handled either with a modification to the halo mass function or by further absorbing this effect in the WL bias calibration, as discussed in Grandis et al. (2021a) and Grandis et al. (2024).

Further, the effect of the super-sample covariance (Hu & Kravtsov 2003) on the recovered contour levels of cosmological parameters are negligible in our analysis with eRASS1 sample. Following the methods of Valageas et al. (2011) and Clerc et al. (2024), we estimate the ratio between the sample variance and the Poisson shot noise and found that the sample variance is just a fifth of the Poisson shot noise, whereas an SPT-like sample has a ratio of a tenth. The eRASS1 super-sample covariance is small and can be neglected in this analysis. This analysis is also consistent with the reported values in the literature. In particular,

Fumagalli et al. (2021) clearly showed that for large area surveys with a mass limit around few $10^{14} M_\odot$, such as SPT and eRASS1, the results from a Poisson likelihood are not affected by the addition of super-sample covariance. In SPT cosmology analyses (Bocquet et al. 2019, 2024), the authors also ignore the effect of sample variance. Fumagalli et al. (2021) also showed that the deeper eROSITA surveys with a pseudo-flux-limited selection function detecting 200k objects are significantly affected by super-sample covariance. Therefore, as the eROSITA survey depth increases and the mass limit of the cluster sample decreases, super-sample covariance will play an essential role in the total error budget in precision cosmology achieved by future eROSITA surveys. We will, therefore, fully account for super-sample covariance in our cosmology analyses in the near future.

Calibrating the point-spread function and mirror effective area of eROSITA is an ongoing effort (see Merloni et al. 2024). In a companion paper by Bulbul et al. (2024), we observe a flux discrepancy in comparisons of eROSITA and Chandra data for extended sources. We stress that this observed difference does not affect our cosmology results. Increasing or decreasing the count rate measurements with a constant factor for all clusters will only shift the normalization of the scaling relations between the count rate and cluster mass and does not alter the slope; therefore, it does not influence our results in a significant way.

The covariant intrinsic correlation between three observables, count rate, weak lensing mass, and richness, and the uncertainties on the selection function (Clerc et al. 2024) will be considered in future eROSITA cosmological analysis. Significant improvements in both memory usage and likelihood computation time are needed to achieve both of these goals.

Thus far, eROSITA completed four All-Sky surveys and 40% of the fifth survey before it was placed in a safe mode (see Merloni et al. 2024, for further details). To project the results from this work to future eROSITA cosmology work with deeper surveys and eRASS1's impact compared to other X-ray cluster surveys, utilizing weak lensing observations in their mass calibration, we compute the figure of merit (FOM), as presented in Garrel et al. (2022). We calculate the area inside the $1\sigma$ contour in the $\Omega_m - \sigma_8$ plane using Green's theorem. We find that eRASS1 provides an improvement on the cosmological parameters by a factor of 9.4, 8.1, 6.6, 2.7, 5.4, and 5.6 over the latest results from the XXL survey, eFEDS, SPT-SZ, SPT-SZ-poleECS-pol500, WtG, DES surveys, respectively (Garrel et al. 2022; Chiu et al. 2023; Bocquet et al. 2019, 2024; Mantz et al. 2015; Dark Energy Survey Collaboration 2020). Overall, eROSITA provides a factor of 5–10 over the published cosmology results from cluster abundances in the literature. This exercise demonstrates the revolutionary power of eROSITA in cluster cosmology.

## 9. Summary

In this work, we present the constraints on the cosmological parameters from the clusters of galaxies detected in the first eROSITA's All-Sky survey in the western Galactic hemisphere through galaxy cluster mass function. We employ the X-ray count rate of the 5259 confirmed clusters of galaxies, the largest pure ICM-selected sample to date, in the redshift range of 0.1–0.8. The photometric redshifts are measured self-consistently using the Legacy Survey South DR10 data covering a 12 791 deg$^2$ area. The mass indicator X-ray observable, count-rate, is measured from the eRASS1 survey by
reprocessing the X-ray data with advanced analysis, including
rigorous background subtraction, the correction for the Galactic
absorption, and the sophisticated ICM modeling. The mass
calibration between the count rate and weak lensing shear signal
from the three large area surveys (DES, HSC, and KiDS) with
overlapping footprints of 4968 deg² is utilized in this work.
In this common area, we use the shear information of a total
of 2348 clusters of galaxies in mass calibration likelihood to
fine-tune the mass to X-ray observable relation. The sample's
selection effects and left-over contamination are fully accounted
for in the probability distribution obtained from the eROSITA's
digital twin in a Bayesian inference forward modeling procedure.
Our strategy, combined with higher statistical precision ensured
by the eROSITA All-Sky survey, leads to better-controlled sys-
tematic uncertainties than those available in previous literature
works. Thus, it allows for reliable and stringent constraints to be
placed on the cosmological parameters.

We built scaling relations between the X-ray count rate and the weak lensing mass for the clusters in the common footprint with DES, KiDS, and HSC. We find that the slope and normalization of the relation are consistent with the studies in the literature utilizing the count rate and luminosity at the 1σ level. Finally, we highlight that the scaling relation scatter is orthogonal to all constrained cosmological parameters, indicating that reducing the mass-to-count rate scatter (while affecting the mass reconstruction) does not improve our cosmological results significantly because it is the precision (and not the value) that matters.

From the eRASS1 sample and cosmological analysis presented in this paper, we derived the $\sigma_8 = 0.88 \pm 0.02$ with 3% precision and $\Omega_m = 0.29_{-0.02}^{+0.01}$ with ~5% precision for the $\Lambda$CDM scenario. To remove the degeneracy between $\sigma_8$ and $\Omega_m$ and to have a fair comparison between CMB results, we combined these two parameters and measure the $S_8 = 0.86 \pm 0.01$ parameter with a 2% precision. Our results, from a late-time probe, are competitive and consistent with Planck CMB 2020 measurements (Planck Collaboration VI 2020), an early-time probe, not confirming the previously reported $S_8$ tension between early and late-time probes (Heymans et al. 2021; Huterer 2023). However, we observe a >3$\sigma$ discrepancy on the $S_8$ parameter measured in some of the cosmic shear, clustering, and BAO/RSD experiments in the literature (Amon et al. 2022; Secco et al. 2022; van den Busch et al. 2022; Li et al. 2023; Sugiyama et al. 2023; Alam et al. 2021).

To test the effect of left-handed neutrinos, we fit an additional set of models to eRASS1 data where the total neutrino mass was allowed to vary. The upper limit to the sum of neutrino masses is $\sum m_\nu < 0.43$ eV (95% upper limit), tighter than the upper limits provided by CMB measurements, for instance, $\sum m_\nu < 0.66$ eV by (Komatsu et al. 2011) and $\sum m_\nu < 0.26$ eV by Planck Collaboration VI (2020) (all 95% confidence level), and much tighter than previous cluster count experiments, e.g., SPT-SZ sample (Bocquet et al. 2019), place an upper limit of $\sum m_\nu < 0.74$ eV (95% upper limit) on the sum of the neutrino mass. We finally combine eRASS1 cluster abundance measurements with Planck CMB 2020 chains and find that the constraints on the total mass of neutrinos improve to $\sum m_\nu < 0.11$ eV. This result represents the tightest bound on the total mass of the left-handed neutrinos to date, and it excludes the inverted hierarchy mass model with a 93% level. When we combine the eRASS1 cluster abundance measurements with Planck CMB and the neutrino oscillation experiments (Tanabashi et al. 2018), we measured the summed masses to be $\sum m_\nu = 0.09_{-0.02}^{+0.04}$ eV

and $\sum m_\nu = 0.12_{-0.02}^{+0.03}$ eV, assuming a normal and inverted mass hierarchy scenario for neutrino eigenstates.

We also explored several other cosmological models with eRASS1 cluster abundances and their combination with Planck CMB or Pantheon+ SNe Ia measurements. Allowing the dark energy equation of state parameter *w* to be free, we find a best-fit measurement of *w* = -1.12 ± 0.12 with 10% precision, consistent at 1σ with the canonical value of -1 of the ΛCDM model. This measurement is a significant improvement over the constraints placed by the previous cluster abundance experiments with the SPT-SZ survey (w = -1.55 ± 0.41, Bocquet et al. 2019). This measurement highlights the remarkable constraining power of eROSITA of all the late-time parameters, particularly $\Omega_m$, $\sigma_8$, and *w*, without needing an external cosmology experiment to break the degeneracies. When the likelihood is combined with Planck CMB and other geometrical probes, such as BAO, SNe Ia, to break parameter degeneracies, the best-fit value does not change significantly; however, the uncertainties improve significantly (w = -1.06 ± 0.03) to a 3% level. Leaving both *w* and $\sum m_\nu$ free in our fits does not change the cosmological parameters significantly; the upper limit on the summed neutrino masses becomes 0.23 eV, while the equation of state parameter remains consistent with the canonical value (w = -1.09 ± 0.15). Including the information from the Planck CMB 2016 data does not significantly improve our constraints on this case.

Finally, we compare all the tested cosmology models against each other in a Bayesian formalism to assess if a preferred model by the eRASS1 data exists. However, we have found no clear preference for a model more complex than the canonical-$\Lambda$CDM. Furthermore, we studied the correlation between cosmological and scaling relation parameters, finding that $S_8$ depends only on the normalization of the scaling relation, while $\Omega_m$ and $\sigma_8$ rely mainly on the mass slope and the redshift evolution of the scaling relation.

With this work, we have established the reliability of the cluster abundance experiments in precision cosmology. The cosmological constraints obtained from the first eROSITA All-Sky Survey performed within the first six months of operations provide a factor of 5–9 improvements over the previous results from similar cluster surveys in the literature. Future eRASS1 cosmology analyses, enabled by deeper eROSITA surveys, will allow for higher precision measurements to be drawn for the fundamental cosmological parameters over eRASS1 as the number of clusters increases, and the redshift span of the sample widens. The robustness of the cluster cosmology experiment demonstrated here corroborates and lends support to many future large experiments, such as *Euclid* and *Vera C Rubin Observatory* studies. To enable a fair comparison of our cosmological results from cluster abundances with other probes, we provide the chains for all our fits online in *GetDist* (Lewis 2019) readable format<fnref n="22" />.

*Acknowledgements.* The authors thank the referee for helpful and constructive comments on the draft. The authors thank Prof. Catherine Heymans for her help and availability for the implementation of the blinding strategy and her useful comments on the draft and Prof. Anja von der Linden for her comments on the population modeling. This work is based on data from eROSITA, the soft X-ray instrument aboard SRG, a joint Russian-German science mission supported by the Russian Space Agency (Roskosmos), in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft und Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its subcontractors and is operated by NPOL with support from the Max Planck Institute for Extraterrestrial Physics (MPE). The development and construction of the eROSITA X-ray instrument was led

<footnote n="22">²² The link has become available after acceptance of the paper.</footnote>