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 degSUP2SUP in the western Galactic hemisphere of eROSITAs first AllSky Survey eRASS1. The common footprint region 4968 degSUP2SUP between the eROSITA Survey and Dark Energy Survey DES the KiloDegree Survey KiDS and the Hyper Supreme Camera HSC survey is used for calibration of the scaling between Xray 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 SUBmSUB  0.29SUB0.02SUBSUP0.01SUP and the normalization of the density fluctuations to SUB8SUB  0.88  0.02 and their combination yields SSUB8SUB  SUB8SUBSUBmSUB0.3SUP0.5SUP  0.86  0.01. These results are consistent and achieve at a similar precision with stateoftheart cosmic microwave background CMB measurements. Furthermore the eRASS1 cosmological experiment places a most stringent upper limit on the summed masses of lefthanded light neutrinos to  mSUBvSUB lt 0.43 eV 95 confidence interval from cluster number counts alone. By combining eRASS1 cluster abundance measurements with CMB and groundbased neutrino oscillation experiments we measured the summed neutrino masses to be  mSUBvSUB  0.09SUB0.02SUBSUP0.04SUP eV or  mSUBvSUB  0.12SUB0.02SUBSUP0.03SUP 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  mSUBvSUB and w are left free we find consistent results with the concordance CDM cosmology. Our results from the first AllSky 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 AllSky Surveys as well as for future cluster abundance experiments.
['surveys', 'galaxies: clusters: general', 'galaxies: clusters: intracluster medium', 'cosmology: observations', 'X-rays: galaxies: clusters', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'General Relativity and Quantum Cosmology']
The SRGeROSITA allsky survey Cosmology constraints from cluster abundances in the western Galactic hemisphere
{'Cosmology constraints from cluster abundances in the western Galactic hemisphere': '- V. Ghirardini 1 , E. Bulbul 1 , E. Artis 1 , N. Clerc 2 , C. Garrel 1 , S. Grandis 3 , 5 , M. Kluge 1 , A. Liu 1 , Y. E. Bahar 1 , F. Balzer 1 , I. Chiu 4 , J. Comparat 1 , D. Gruen 5 , F. Kleinebreil 3 , 8 , S. Krippendorf 5 , 6 , A. Merloni 1 , K. Nandra 1 , N. Okabe 7 ,  \nF. Pacaud 8 , P. Predehl 1 , M. E. Ramos-Ceja 1 , T. H. Reiprich 8 , J. S. Sanders 1 , T. Schrabback 3 , 8 , R. Seppi 1 , S. Zelmer 1 , X. Zhang 1 , W. Bornemann 1 , H. Brunner 1 , V. Burwitz 1 , D. Coutinho 1 , K. Dennerl 1 , M. Freyberg 1 , S. Friedrich 1 , 1 1 1 1 9 11 1 1 1 1 ,  \n- R. Gaida , A. Gueguen , F. Haberl , W. Kink , G. Lamer , X. Li , T. Liu , C. Maitra , N. Meidinger , S. Mueller H. Miyatake 13 , 14 , 15 , S. Miyazaki 12 , J. Robrade 10 , A. Schwope 9 , I. Stewart 1\n- 1 Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse 1, 85748 Garching, Germany\n- 2 IRAP, Université de Toulouse, CNRS, UPS, CNES, F-31028 Toulouse, France\n- 3 Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstr. 25 / 8, 6020 Innsbruck, Austria\n- 4 Department of Physics, National Cheng Kung University, 70101 Tainan, Taiwan\n- 5 Universitäts-Sternwarte, Faculty of Physics, LMU Munich, Scheinerstr. 1, 81679 München, Germany\n- 6 Arnold Sommerfeld Center for Theoretical Physics, LMU Munich, Theresienstr. 37, 80333 München, Germany\n- 7 Department of Physical Science, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima,Hiroshima 739-8526, Japan\n- 8 Argelander-Institut für Astronomie (AIfA), Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany\n- 9 Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany\n- 10 Hamburger Sternwarte, Gojenbergsweg 112, 21029 Hamburg, Germany\n- 11 McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA\n- 12 Subaru Telescope, National Astronomical Observatory of Japan, 650 N Aohoku Place Hilo HI 96720 USA\n- 13 Kobayashi-Maskawa Institute for the Origin of Particles and the Universe (KMI), Nagoya University, Nagoya, 464-8602, Japan\n- 14 Institute for Advanced Research, Nagoya University, Nagoya 464-8601, Japan\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  \nJuly 26, 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 2 in the western Galactic hemisphere of eROSITA's first All-Sky Survey (eRASS1). The common footprint region (4968 deg 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 + 0 . 01 -0 . 02 and the normalization of the density fluctuations to σ 8 = 0 . 88 ± 0 . 02 , and their combination yields S 8 = σ 8( Ω m / 0 . 3) 0 . 5 = 0 . 86 ± 0 . 01. These results are consistent and achieve at a similar precision with stateof-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 P 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 P m ν = 0 . 09 + 0 . 04 -0 . 02 eV or P m ν = 0 . 12 + 0 . 03 -0 . 02 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 P m ν and w are left free, we find consistent results with the concordance Λ CDMcosmology. 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.  \nKey words. Surveys - Galaxies: clusters: general - Galaxies: clusters: intracluster medium - X-rays: galaxies: clusters - Cosmology: observations", '1. Introduction': 'The so-called Λ CDM concordance cosmological model is currently the most widely accepted theory on the formation and  \nevolution 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 et al. 2020a), 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 et al. 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 et al. 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 di ff erent 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).  \nClusters 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 8 h -1 Mpc, namely σ 8 1 , and the present-day matter density parameter, Ω 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 et al. 2020a), and galaxy and cosmic shear autoand cross-correlations (Dark Energy Survey Collaboration et al. 2022; Dark Energy Survey and Kilo-Degree Survey Collabora-  \nn et al. 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, P m ν , can also be constrained via cluster abundances, particularly in combination with external probes, as the cosmic neutrinos have a non-negligible e ff ect 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.  \nMultiwavelength 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 et al. 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) e ff ect (Planck Collaboration et al. 2014, 2016b; 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, Xray 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, su ff er from projection e ff ects and contamination, limiting their use in cluster abundance measurements (Dark Energy Survey Collaboration et al. 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 su ff er 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.  \nAwell-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 et al. 2016b). However, recent numerical simulations indicate that the contribution of non-thermal pressure due to bulk motions and turbulence can lead to significant biases on the cluster masses measured with these traditional techniques (Lau et al. 2009; Rasia et al. 2012; Scheck et al. 2023). The use of biased X-ray mass measurements is also reflected in the cosmological parameters obtained from cluster counts (Planck Collabora-  \ntion et al. 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) e ff ect directly probes the cluster\'s gravitational potential, thus allowing the masses to be accurately measured without any a priori assumption about their hydrostatic equilibrium properties. It is worth noting that this method is not entirely bias-free; the correlated large-scale structure along the line-of-sight, triaxiality, and assumed underlying dark matter distribution, baryonic feedback, and miscentering of the WL signal could introduce additional bias and scatter on the mass scaling relations (Hoekstra 2003; Meneghetti et al. 2010; Becker & Kravtsov 2011; Bahé et al. 2012). However, compared to X-ray measurements, the accuracy of WL measurements is inherently less impacted by hydrodynamical modeling as the majority of the mass is collisionless. Most importantly, since the WL mass measurements are calibrated with numerical simulations, the bias and systematic uncertainties can be reliably and quantitatively estimated and integrated into the final cluster mass calibration e ff ort (Schrabback et al. 2018; Dietrich et al. 2019; Grandis et al. 2021a; Sommer et al. 2022; Zohren et al. 2022; Grandis et al. 2024).  \nAnother 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 e ff ects, 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 e ff ects 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 di ff erent 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).  \nThis work presents the constraints on cosmological parameters from the clusters of galaxies detected in the western Galatic hemisphere in the first eROSITA All-Sky Survey. As one of the two X-ray telescopes on board the Spectrum Roentgen Gamma Mission ( SRG ), eROSITA, was launched on July 13, 2019. The first All-Sky Survey (eRASS1) was successfully executed between December 2019 and June 2020 (Predehl et al. 2021; Sunyaev et al. 2021) and has resulted in the detection of nearly 930 thousand X-ray sources in the western Galactic hemisphere (179.9442 deg < l < 359.9442 deg), where the data and publication rights belong to the eROSITA German Consortium (Predehl et al. 2021; Merloni et al. 2024). In this first All-Sky Survey, we detect a total of 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 5,259 securely detected clusters in the 12,791 deg 2 common footprint of eRASS1 and the DESI Legacy Survey DR10-South in the  \nphotometric redshift range of 0 . 1 < z λ < 0 . 8. This catalog represents the largest ICM-selected cluster sample utilized for cosmological studies.  \nBesides the statistical power ensured by the extensive eRASS1 cluster catalog, our analysis o ff ers 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 e ff ects 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 2 overlapping with eRASS1 (Chiu et al. 2023; Kleinebreil et al. 2024). The biases due to miscentering, baryonic e ff ects, 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).  \nThis 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 Sect. 3 and Sect. 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 σ di ff erence, by "in tension" we mean between 3 σ and 5 σ in di ff erence, and by "discrepant," we mean above 5 σ in di ff erence.', '2.1. eRASS1 Cluster Catalog': 'The cluster catalog utilized in this work is based on the soft Xray band (0.2 -2.3 keV) catalog of the first eROSITA All-Sky Survey in the western Galactic hemisphere presented in Merloni et al. (2024). The details of the detection and cleaning processes and compilation of the cluster catalog are described in detail in Bulbul et al. (2024). We restrict ourselves to the cosmology subsample in this work, a subsample of the primary galaxy groups and clusters catalog. To compile the cosmology subsample, Bulbul et al. (2024) applied a conservative extent likelihood cut of L ext > 6 to the primary sample, resulting in 11,141 cluster candidates detected as extended in the clean eRASS1 cluster catalog.  \nFig. 1: Sky distribution of the 5259 confirmed clusters of galaxies in the cosmology subsample, detected in the western Galatic hemisphere of the eRASS1 All-Sky Survey in equatorial coordinates. The common footprint between the Legacy Survey-South and eRASS1 reaches up to 12791 deg 2 . The sample covers a redshift range of 0.1 to 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 2 in total.  \n<!-- image -->  \nThe confirmation, identification, and redshift measurements of the X-ray-selected cluster candidates are typically performed using the publicly available DESI Legacy Survey Data Release 10 data (LS DR10). LS DR10 is a compilation of three surveys that cover most of the extragalactic sky (Dey et al. 2019; Sevilla-Noarbe et al. 2021). The cosmology subsample is compiled using only the southern area below Dec. ≲ 32 . 5 · (LS DR10South 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 (ˆ z λ ) and richness ( ˆ λ ) measurements. The common footprint between the eROSITA survey in the western Galactic hemisphere and the southern LS DR10-South footprint is 12,791 deg 2 , and contains 8,129 cluster candidates (Kluge et al. 2024).  \nWe 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 7,077 candidate clusters in the LS DR10-South footprint, we were able to identify the optical counterparts of 6,562 securely detected clusters.  \nWe 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-  \nntrolled 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 et al. 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 5,259 securely confirmed galaxy clusters in the 12,791 deg 2 LS DR10-South area. The projected spatial distribution of the clusters in this sample is shown in Fig. 1, overlaid on the exposure map of the eRASS1 survey. The median redshift in the catalog is 0.29. The redshift and richness histograms of all confirmed clusters are presented in the leftmost panel of Fig. 2. The mixture model, extensively described in Sect. 4.2, implies that the purity of this sample is at the level of 94%, ideal for cosmological studies with cluster  \ncounts (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 process should be adequately included in the selection model and should show a correlation with cluster mass. Additionally, if the X-ray observable has a negligible variation over the sky and a minimal dependence on cosmology (which are not strictly necessary), the computation cost of the forward hierarchical model significantly decreases. In this regard, we examine extensively several observables, for instance: the extent likelihood, L ext , detection likelihood, L det , and count-rate. In addition, we included commonly used observables in the literature, such as luminosity, gas mass, and the mass proxy YX 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.  \nThe 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 backgroundincluded X-ray images of galaxy clusters to fit cluster and background emission in a total of seven bands between 0.3 keV to 7 keV. We obtain the corrected count rate from the best-fit cluster physical model, including a correction for Galactic hydrogen absorption at photometric redshifts z λ . The measured count rate in the 0.2 -2.3 keV band is used in the further weak lensing mass calibration process and is also shown in the rightmost panel of Fig. 2 and has been released in the cluster catalog (Bulbul et al. 2024).", '2.3. Weak-lensing follow-up observations and analysis': 'In cosmology studies with galaxy cluster abundances, the cluster mass function should be determined with minimal bias. The X-ray observations alone do not provide direct cluster mass measurements, and the usage of hydrostatic mass estimates from X-rays is known to produce significant biases in cosmology results (Planck Collaboration et al. 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 Galac-  \ntic 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 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 M etacalibration pipeline (Hu ff & Mandelbaum 2017; Sheldon & Hu ff 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 et al. 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 e ff ective 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 photoz uncertainties propagated from the calibrations as mentioned above.  \nConsidering the overlap between the DES Y3 and the eRASS1 footprints, we produced tangential shear data for 2201 eRASS1 galaxy clusters, with a total signal-to-noise ratio (S / N) of 65 in the tangential shear profile.', '2.3.2. The Hyper Suprime-Cam (HSC) survey': 'We extract the weak-lensing measurements of eRASS1 clusters using the data from the Hyper Suprime-Cam Subaru Strategic Program (HSC SSP, Aihara et al. 2018), a wide and deep multiband ( g , r , i , z , and Y ) survey in the optical. We make use of the latest three-year weak-lensing data (S19A) (Li et al. 2022) with a total area coverage of ≈ 500 deg 2 , 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 2 , which has been used to deliver robust  \n1  \nFig. 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.  \n<!-- image -->  \nconstraints 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.  \nThe 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 R ∞ z cl + 0 . 2 P ( z ) d z > 0 . 98, where z cl is the cluster redshift and P ( z ) is the photometric redshift distribution of each source. We used the DEmP code (Hsieh & Yee 2014) to estimate the photometric redshift (Tanaka et al. 2018; Nishizawa et al. 2020). The stacked photometric redshift distribution of the sources selected in non-cluster regions was used to estimate the lensing e ffi ciency, which is needed to obtain the weak-lensing mass from the modeling of the shear. With the selected source sample, the tangential shear profile gt ( θ ) as a function of the angular radius θ around the cluster center was calculated by following Sect. 3.4 in Chiu et al. (2022). The best-fit X-ray peak ( RA \\_ XFIT and DEC \\_ XFIT ), the output of the MBProj2D provided in the cosmology catalog of (Bulbul et al. 2024), was used as the cluster center to extract the shear profiles around. The radial binning was defined as ten logarithmic bins in the angular space, corresponding to 0 . 2 Mpc / h and 3 . 5 Mpc / h in the physical coordinate in a fiducial flat Λ CDMcosmology with h ≡ H 0 100 km / s / Mpc = 0 . 7 and Ω m = 0 . 3. We discarded the innermost three bins to mitigate the systematics in the weak-lensing modeling at the cluster core ( ≲ 0 . 5 Mpc / h ). The uncertainty of the shear profile arising from the shape noise and uncorrelated large-scale structures is accounted for in the weak-lensing covariance matrix by following Sect. 3.5 in Chiu et al. (2022).  \nThe 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 di ff erent 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  \nmember contamination to be zero across the radial bins with a 2 σ upper limit of 6% at the innermost bin ( R ≈ 0 . 2 Mpc / h ). It then decays outwards following a projected Navarro-FrenkWhite (NFW) model (Navarro et al. 1996). With a calibration against the COSMOS 30-band photoz sample and a spectroscopic sample, the bias in the photometric redshift of the selected sources can be quantified as ≲ 6% for the cluster redshift z cl ≲ 1 . 2 and accounted for in the weak-lensing modeling. The multiplicative bias, m , in the shape measurements as a function of the source redshift is quantified with a systematic uncertainty, δ m , at 1 . 7% level (Mandelbaum et al. 2018).  \nAs a result, the shear profile gt ( θ ), 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 photoz 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 KiloDegree Survey 5 (Kuijken et al. 2019; Wright et al. 2020; Hildebrandt et al. 2021; Giblin et al. 2021), hereafter referred to as KiDS-1000.  \nWe 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  \nredshift bins. We also used the enhanced redshift calibration presented in van den Busch et al. (2022), which provides the basis for calibrating the shape measurement and photo-z systematic uncertainties.  \nOur 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 e ff ect of blending background galaxies with cluster galaxies on the number density profiles. For this, we simulated KiDSlike galaxy images using GALSIM (Rowe et al. 2015), inject them into the actual r -band detection images of KiDS, and rerun the source detection to extract radial detection probability profiles.  \nWe 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 e ff ects. 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': "Weused 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 3 particles) enable reliable prediction of the large-scale structure and its haloes down to 10 11 M ⊙ . This resolution is su ffi cient 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 e ff ects (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).  \nIn 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).  \nThe model for clusters of galaxies is based on the X-ray observations in the literature (e.g., Reiprich & Böhringer 2002;  \nPierre 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 500 c < 5 × 10 13 M ⊙ ), and it ends up over-predicting the X-ray luminosity of these objects. To overcome this e ff ect, 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.  \nFor 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 di ff erent 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 fivedimensional 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 ) = σ 2 exp( -r 2 / 2), with σ being a free hyperparameter and r as the Euclidian distance between pairs of points, normalized along each dimension by a free hyperparameter (lengthscale). This kernel is the widely-used stationary squared exponential (Rasmussen & Williams 2006) ; lengthscales set the size of "wiggles" in the model while σ controls the average deviation of the latent function from its mean (that is supposed uniform). Hyperparameters are optimized to minimize the marginal likelihood in the nominal implementation of GP classification. The quantity of interest is the posterior of the latent function given the set of training points. For any new set of feature values, the posterior predicts the distribution of the corresponding values of the latent function ; a probabilistic prediction can be obtained by integrating this distribution against the link function.  \nIn 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, CR 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 since 2012) by taking advantage of its numerical scalability properties, that is well-suited to our large training set (of size ≃ 4 . 8 × 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 e ffi ciency while preserving good convergence properties. We found that minibatches of size 2 16 and 12000 iterations provide good results within a tractable time. The optimization of a set of 651 hyperparameters 7  \nFig. 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σ 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.  \n<!-- image -->  \nrelies on the A dadelta 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, di ff ering only by the exact definition of an object that qualifies as "detected" in the simulation. We modified the lower thresholds in selection parameters ( L ext , EXTENT ), choosing (3.6, 16.6), (3.9, 15.4), (4.4, 15.1), (4.8, 13.7), (5.5, 11.6), and (6, 0). By design, the six models predict di ff erent 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 ( L ext > 6, EXTENT > 0). Each model is tested against the remaining third of the simulation sets (each of size ≃ 2 . 5 × 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 modelpredicted 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 AllSky survey to ensure the stability of the results against photon shot noise.  \nWhile the most valuable outcome of the models is the expectation values for the probabilistic predictions P ( I | z , CR , 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 ( CR ), 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  \nFig. 4: Representation of a model P ( I | z , CR , H i ) for fixed values of cluster redshift, z = [0.15, 0.3, 0.5], local galactic absorption n H = 3 . 5 × 10 20 cm -2 , a nominal background level 5 . 2 × 10 -15 erg / s / cm 2 / arcmin 2 , three values of exposure times [50, 150, 500] s, as a function of count-rate.  \n<!-- image -->  \nwith 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 di ff erent 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 measured through eROMaPPer run on the initial cluster catalog as described in Kluge et al. (2024). E ff ectively, these are redshift and richness cuts, 0 . 1 < ˆ z < 0 . 8 and ˆ λ > 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.  \nThis implies that the redshift and richness cuts do not a ff ect 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) a ff ects only the integration limits in likelihood terms where richness or redshift are involved, as discussed in Grandis et al. (2020).  \nWe 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 ( ˆ λ > 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 e ff ect for cluster samples. This test demonstrates that our optical  \nconfirmation 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 ˆ CR indicates the observed n H-corrected count rate, while CR 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; Zubeldia & Challinor 2019; Chiu et al. 2023; Bocquet et al. 2024), we started from a Poissondi ff erential log-likelihood for our set of observables, X-ray n H corrected count rate. ˆ CR , redshift, ˆ z , sky position, ˆ H i , richness, ˆ λ , and weak-lensing tangential shear profile, ˆ gt ;  \nlog L = X j log dNtot d ˆ CRd ˆ zd ˆ λ d ˆ gtd ˆ H i P ( I | ˆ CR , ˆ z , ˆ λ, ˆ H i , ˆ gt ) ! -( GLYPH<20> dNtot d ˆ CRd ˆ zd ˆ λ d ˆ gtd ˆ H i P ( I | ˆ CR , ˆ z , ˆ λ, ˆ H i , ˆ gt ) d ˆ CRd ˆ zd ˆ λ d ˆ gtd ˆ H i GLYPH<21> , (1)  \nwhere 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.  \nWehighlight 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 e ff ect 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:  \nP ( I | ˆ CR , ˆ z , ˆ λ, ˆ H i , ˆ gt ) = P ( I | ˆ CR , ˆ z , ˆ H i ) Θ ( ˆ λ > 3) , (2)  \nwhere Θ ( ˆ λ > 3) is the Heaviside step function, equals to 1 where ˆ λ > 3, and equal to 0 where ˆ λ < 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).  \nWe, therefore, isolate richness and tangential shear profile likelihood terms. In practice, the first term in Equation 1 can be simplified into:  \ndNtot d ˆ CRd ˆ zd ˆ λ d ˆ gtd ˆ H i = dNtot d ˆ CRd ˆ zd ˆ H i P ( ˆ λ, ˆ gt | ˆ CR , ˆ z , ˆ H i ) . (3)  \nGiven these considerations, it is therefore consequent that the measurement of tangential shear profiles does not directly a ff ect 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;  \nlog L = X j log dNtot d ˆ CRd ˆ zd ˆ H i P ( I | ˆ CR , ˆ z , ˆ H i ) Θ ( ˆ λ > 3) ! + X j log P ( ˆ λ, ˆ gt | ˆ CR , ˆ z , ˆ H i , I ) -( GLYPH<20> dNtot d ˆ CRd ˆ zd ˆ λ d ˆ H i P ( I | ˆ CR , ˆ z , ˆ H i ) Θ ( ˆ λ > 3) d ˆ CRd ˆ zd ˆ λ d ˆ H i GLYPH<21> (4)  \nAdditionally, 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-ofsight projection e ff ects 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:  \nP ( ˆ λ, ˆ gt | ˆ CR , ˆ z , ˆ H i , I ) = P ( ˆ λ | ˆ CR , ˆ z , ˆ H i , I ) P (ˆ gt | ˆ CR , ˆ z , ˆ H i , I ) . (5)  \nThus, we can write our final log-likelihood as:  \nlog L = X j log dNtot d ˆ CRd ˆ zd ˆ H i P ( I | ˆ CR , ˆ z , ˆ H i ) ! -( GLYPH<20> dNtot d ˆ CRd ˆ zd ˆ λ d ˆ H i P ( I | ˆ CR , ˆ z , ˆ H i ) Θ ( ˆ λ > 3) d ˆ CRd ˆ zd ˆ λ d ˆ H i GLYPH<21> + X j log P ( ˆ λ | ˆ CR , ˆ z , ˆ H i , I ) + X j log P (ˆ gt | ˆ CR , ˆ z , ˆ H i , I ) . (6)  \nThe 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.  \nIn 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 di ff er from zero. Specifically, redshift from 0 to 0.9, mass from 5 × 10 12 M ⊙ to 5 × 10 15 M ⊙ , richness from 1 to 600, and count rate from 10 -3 cts / s 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 et al. 2016b; 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) 8 . As large cluster samples from wide-area surveys become available, contamination modeling becomes crucial in statistical cosmological studies with cluster abundances.  \nIn eRASS1 cosmology subsample, after the 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 ( √ 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.  \nInformed 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 randomposition 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 cosmologi-  \ncal inference. As a consequence, we need to model the richnessmass relation of the cluster population (see Sect. 4.5).  \nIn 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).  \nThe complete derivation of the mixture model and its e ff ects on each likelihood term in Equation (6) is shown in Appendix A. In short, the probability density functions for every observables ˆ O is written as:  \nP ( ˆ O ) = X i ∈{ C , AGN , RS } f i P ( ˆ O| i) , (7)  \nwhere this equation holds for all our observables ˆ O as ˆ CR , ˆ z , ˆ λ , and ˆ gt . 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 di ff erential cluster number counts term, dNC d ˆ CRd ˆ zd ˆ H i P ( I | ˆ CR , ˆ z , ˆ H i ), which is a significant term in Equation (6), being involved in both terms of the di ff erential 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).  \nWe consider this term as a marginalized likelihood over true count rate, true mass, and true redshift. 9  \ndNC d ˆ CRd ˆ zd ˆ H i P ( I | ˆ CR , ˆ z , ˆ H i ) = P ( ˆ H i ) $ P (ˆ z | z ) P ( I | z , CR , ˆ H i ) × P ( ˆ CR | CR ) dN dV P ( CR | M , z ) dMdCRdz (8)  \ndMdV dz  \nwhere P ( ˆ H i ) is the uniform probability over the sphere that a cluster exists at a specific sky position ˆ H i , P (ˆ z | z ) expresses the uncertainty on our redshift measurements, provided by eROMaPPer (Kluge et al. 2024), P ( I | z , CR , ˆ H i ) is the X-ray selection function (see Sect. 3 and for more details Clerc et al. (2024)), P ( ˆ CR | CR ) is the uncertainty on the X-ray count rate,  \nmeasured by MBProj2D, dN dMdV is the halo mass function 10 , which directly depends on cosmology, dV dz is the fractional comoving volume in an infinitesimal redshift shell (see Equation (28) in Hogg 1999) 11 , P ( CR | M , z ) is the log-normal distribution around the X-ray scaling relation and can be described as:  \nP ( CR | M , z ) = LN ( ⟨ CR | M , z ⟩ , σ X ) . (9)  \nWe 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 di ff erential Poisson term, and can be computed via the equation:  \nN C = & dN C d ˆ CRd ˆ zd ˆ λ d ˆ H i P ( I | ˆ CR , ˆ z , ˆ H i ) Θ ( ˆ λ > 3) d ˆ CRd ˆ zd ˆ λ d ˆ H i (10)  \nwhere several simplifications are applied. First, P (ˆ z | z ), the probability of observing a cluster with measured redshift ˆ z given the true redshift z , is modeled as a mixture of two normal distributions whose parameters are calibrated by eROMaPPer , expressed as:  \nP (ˆ z | z ) = (1 -cz ) N ( bz z , σ z (1 + z )) + cz N ( bz z + c shift , z , σ z (1 + z )) , (11)  \nwhere cz is the fraction of objects for which we measure a redshift shifted by some constant c shift , z , bz is a systematic bias in our redshift estimate, and σ z is the relative error on the measured redshift. We can, therefore, express g ( z , H i ) = R P (ˆ z | z ) w (ˆ z , H i ) d ˆ z as an analytic expression for each sky position, with w (ˆ z , H i ) being the filter function that is redshift and sky position-dependent, as follows:  \nw (ˆ z , H i ) = ( 1 , if (ˆ z > 0 . 1) & (ˆ z < min(0 . 8 , ˆ z max( ˆ H i )) , 0 , otherwise, (12)  \nwhere ˆ z max( ˆ H i ) is the maximum measurable redshift at sky position ˆ H i , given by the ZVLIM map (see Kluge et al. 2024, for further details).  \nThen, P ( ˆ λ | λ ), the probability of observing a cluster with measured richness ˆ λ given the true richness λ , is modeled as in previous cosmological analyses with cluster counts (Dark Energy Survey Collaboration et al. 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 e ff ectively counts galaxy members, and partly demonstrated in Costanzi et al. (2019a). Therefore, we have:  \nP ( ˆ λ | λ ) = N ( λ, √ λ ) . (13)  \nThus, we can express the optical completeness term f ( λ ) = R P ( ˆ λ | λ ) Θ ( ˆ λ > 3) d ˆ λ term analytically. 12 The e ff ects of the optical cut, ˆ λ > 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.  \nFinally, since ˆ CR appears only in the term P ( ˆ CR | CR ), and since probabilities are normalized, the integral over measured count rate is unity and does not a ff ect the calculation of the total number of objects. Therefore, the total number of clusters is expressed as:  \nNC = & P ( I | z , CR , ˆ H i ) g ( z , ˆ H i ) f ( λ ) P ( CR , λ | M , z ) P ( ˆ H i ) × dN dMdV dV dz dMdCRdzd λ d ˆ H i , (14)  \nTo conclude, we now have expressed explicitly all the ingredients required in order to model the cluster number counts likelihood in Equation (8) and to compute the total number of clusters that are detected in Equation (14).', '4.4. Weak-lensing mass calibration likelihood': 'This subsection focuses on the weak lensing likelihood, which is crucial for mass calibration. We make use of Bayes\' theorem to express the weak lensing likelihood as:  \nP (ˆ gt | ˆ CR , ˆ z , ˆ H i , I , C ) = P (ˆ gt , ˆ CR , ˆ z , ˆ H i , I , C ) P ( ˆ CR , ˆ z , ˆ H i , I , C ) = N WL D WL , (15)  \nwhere 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:  \nN WL = P ( ˆ H i ) $ P ( I | CR , ˆ z , ˆ H i ) P ( ˆ CR | CR ) P ( M , ˆ z ) × P (ˆ gt | M WL , ˆ z ) P ( M WL , CR | M , ˆ z ) dMdCRdM WL (16)  \nand  \nD WL = P ( ˆ H i ) " P ( I | CR , ˆ z , ˆ H i ) P ( ˆ CR | CR ) P ( M , ˆ z ) × P ( CR | M , ˆ z ) dMdCR . (17)  \nWe note that the sky position probabilities, P ( ˆ H i ), cancel out in the ratio of N WL / D WL in the weak lensing mass calibration term. P (ˆ gt | M WL , ˆ z ) is the likelihood term that compares the measured weak lensing shear with the one predicted from the weak  \nlensing mass (Grandis et al. 2024). Furthermore, in both the numerator and denominator, we have the error on the measured count rate P ( ˆ CR | CR ), the selection function P ( I | CR , ˆ z , ˆ H i ), the normalized halo mass function, P ( M , ˆ z ). However, since the normalization factor (needed to compute P ( M , z ) from dN dMdV dV dz ) is present both at the numerator and denominator, it cancels itself out, and the halo mass function, dN dMdV dV d ˆ z , can be used directly in the previous equations. All the details regarding the expressions for these terms are described just below Equation (8).  \nFinally, P ( M WL , CR | M , ˆ z ), the scaling relation term between mass and both count rate and weak lensing mass that takes into account the intrinsic correlation between these quantites, is a bivariate log-normal distribution BN ( ¯ µ, ¯ Σ ) , where  \n¯ µ = [ GLYPH<10> log CR | M , z GLYPH<11> , GLYPH<10> log M WL | M , z GLYPH<11> ] , (18)  \n¯ Σ =       σ 2 CR ρ M WL , CR σ CR σ M WL ρ M WL , CR σ CR σ M WL σ 2 M WL       , (19)  \nwhere GLYPH<10> log CR | M , z GLYPH<11> and GLYPH<10> log M WL | M , z GLYPH<11> are the scaling relations between true cluster mass and count rate and weak lensing mass, respectively, see Sect. 6 for further details, and ρ M WL , CR 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:  \nP ( ˆ λ | ˆ CR , ˆ z , ˆ H i , I , C ) = P ( ˆ λ, ˆ CR , ˆ z , ˆ H i , I , C ) P ( ˆ CR , ˆ z , ˆ H i , I , C ) = N opt D opt , (20)  \nwhere (as before) we have factored out the numerator N opt and the denominator D opt. As before, we expand the numerator and denominator as:  \nN opt = P ( ˆ H i ) $ P ( I | CR , ˆ z , ˆ H i ) P ( ˆ CR | CR ) P ( M , ˆ z ) × P ( ˆ λ | λ ) Θ ( ˆ λ > 3) P ( λ, CR | M , ˆ z ) dMdCRd λ, (21)  \nD opt = P ( ˆ H i ) " P ( I | CR , ˆ z , ˆ H i ) P ( ˆ CR | CR ) P ( M , ˆ z ) × P ( CR | M , ˆ z ) dMdCR . (22)  \nHere, in particular, P ( ˆ λ | λ ) 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 | CR , ˆ z , ˆ H i ) is the X-ray selection function, Θ ( ˆ λ > 3) is the optical selection function 13 , P ( M , ˆ z ) is the normalized mass function, P ( ˆ CR | CR ) is the error on the measurement of the count rate,  \nP ( λ, CR | M , ˆ z ) is a bivariate log-normal distribution MN ( ¯ µ, ¯ Σ ) where  \n¯ µ = [ GLYPH<10> log CR | M , ˆ z GLYPH<11> , GLYPH<10> log λ | M , z GLYPH<11> ] , (23)  \n¯ Σ = " σ 2 CR ρλ, CR σ CR σλ ρλ, CR σ CR σλ σ 2 λ + σ 2 P # , (24)  \nwhere GLYPH<10> log CR | M , z GLYPH<11> and GLYPH<10> log λ | M , z GLYPH<11> are the scaling relations between true mass, count rate, and optical richness, details are given in Sect. 6, while ρλ, CR is the intrinsic correlation between richness and count rate and where σλ is richness intrinsic scatter, and σ P is a Poisson term, defined as in Costanzi et al. (2019a) by ⟨ λ | M , z ⟩-1 ⟨ λ | M , z ⟩ 2 . We point out that the denominator for the optical, D opt, and weak lensing, D WL, likelihoods look identical; however, we recall that while the optical term runs over all the objects in the sample, the weak-lensing term runs over only the objects that have tangential shear measurements.', '4.6. Parameters and priors': 'In Table 1, we list all the parameters fitted in our Bayesian statistical model. We include their units, descriptions, and prior functions used. We note that our analysis is not sensitive to some cosmological parameters. For this reason, to avoid the chains reaching unphysical values but still allowing marginalization over these parameters, we limit their range using realistic priors from previous cosmological experiments. Specifically, we limit Ω b in a range allowed by the Big Bang Nucleosynthesis (Beringer et al. 2012; Cooke et al. 2016, 2018; Mossa et al. 2020; Pitrou et al. 2021), ns in a range round Planck CMB results (Planck Collaboration et al. 2020a). For H 0, against which we are not sensitive, we use the Planck CMB Λ CDM results with a tight Gaussian prior set around them (Planck Collaboration et al. 2020a). Given the so-called H 0 tension, we have checked that using a prior from distance ladder method ( H 0 ∼ N (73 . 04 , 1 . 04) km s -1 Mpc -1 , Riess et al. 2022) or a flat prior ( H 0 ∼ U (60 , 80) km s -1 Mpc -1 ) does not a ff ect 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 e ff ects, 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  \nfunction, 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.  \nThe validation of the cosmological framework is performed independently. We developed a new pipeline for generating mock catalogs following the physical formalism (described in Sect. 4 and Sect. 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).  \nWe 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 ⊙ , 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 et al. 2016; Willingale et al. 2013), and exposure time (Predehl et al. 2021) by interpolating the corresponding maps consistent with the survey strategy. To compute the expected observables, such as the count rate and richness, we employed the scaling relations described in 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.  \nThe dark matter surface mass density was calculated by assuming a spherical NFW profile (Navarro et al. 1996) for each cluster. We computed the mean expected shear by integrating this value over the background galaxy distribution from the respective survey. The observed shear was deduced by adding errors similar to those observed in the DES shear catalog(see Sect. 2.3.1 and Grandis et al. 2024). Out of these shear simulations, we selected only those that fall into the optical survey footprint of DES, HSC, and KiDS, while the corresponding shear profiles are assigned to these objects. Finally, by sampling their probability distribution, we included the contaminants, RSs and AGNs (see Sect. 4.2) and used the contamination fractions from the best fit.  \nWe generated ten di ff erent mocks from the same model and used the cosmological pipeline to recover the input parameters, σ 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.  \nTo test the e ff ects of the selection on our results, we applied realistic changes to our selection function. We generated six se-  \nTable 1: Priors used on the free parameters in the fitting, divided into several parameter types for clarity.  \nNotes. With U (min , max) we indicate a uniform distribution between "min" and "max." With N ( µ, σ ) we indicate a normal distribution centered on µ and with standard deviation σ . With TN ( µ, σ, min , max) we indicate a truncated normal distribution centered on µ and with standard deviation σ bounded between "min" and "max."  \nlection functions 14 described in Sect. 3.4. These models di ff er by the thresholds used to construct the simulation catalog. We then generated six mocks, using a di ff erent selection function each time, freezing all the other components in the likelihood. We then ran the cosmology pipeline with di ff erent selection functions to investigate the impact of the uncertainties on the selection function on our results. We find that the contours shift  \nslightly but always within the 68% percentiles of the input of the σ 8 and Ω m parameters.  \nRecent results by Sommer et al. (2023) have suggested that the underestimated miscentering measurements between ICM centers and halo centroids indicated by the numerical simulations may introduce an additional systematic bias on WL mass measurements by 2-6%. However, the exact extent of the bias remains to be investigated and depends on the assumed X-ray miscentering in hydrodynamical simulations. Although the uncertainties in miscentering distribution calibrated on hydrodynam-  \nical simulations are fully accounted for in our WL shear measurements (see Sect. 6.2 and Grandis et al. (2024)), we further investigated the e ff ect of an additional miscentering bias and its influence on our measurements by introducing a bias at the levels as suggested by Sommer et al. (2023) in Equation (28). We find that a decrease in the masses as found by Sommer et al. (2023) 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.  \nThe redshifts in our catalog were measured using the photometric 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 measurements. This potentially could cause a bias on the cosmological parameters mainly as a result of geometrical e ff ects. This is because dN dMdV is a weak function of redshift, especially in the redshift interval of interest, but the volume dV dz and cosmological distances change drastically. We tested this e ff ect 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, σ z ≈ 0 . 1, a ff ects Ω m and σ 8 by about 3% and 2% respectively, and a trend of catastrophic failures on photo-z with occurrence of 10% results in lower values of Ω m by 4% and higher values of σ 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 Equation (11). The intrinsic scatter of the photometric redshift with respect to the spectroscopic redshift is < 0 . 5% , on average, and the fraction of significant catastrophic failures is just 0.13%; thus, it is much smaller than the extreme parameters required to introduce a significant bias. Therefore, the bias introduced by the use of photometric redshifts on Ω m and σ 8 parameters is negligible.  \nGiven the positive results from the robustness tests and verification of the contamination modeling, assessment of the WL mass, and redshift bias on the cosmological constraints, we do not anticipate any significant bug or unaccounted bias in our pipeline and methodology. Therefore, we conclude that the cosmology pipeline and the methods have been validated successfully.', '5.2. Blinding strategy': 'Blinding the results prior to publication minimizes preconceived assumptions and prevents us from introducing bias into our analysis. After the validation phase was successfully completed on our simulations and mock data, we froze the cosmology pipeline and apply the following blinding strategy to ensure that our conclusions would remain as objective and unbiased as possible. We generated six di ff erent 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 shu ffl ed their indices. During this process, the \'true\' selection function model and its index became unknown to the working group and analysis teams. The cosmology analysis for the flat Λ CDM model  \nTable 2: Best-fit parameters of the scaling relations between WL shear and X-ray count-rate and optical richness.Notes. The reported uncertainties are the 68th percentiles. The posteriors of the best-fit parameters are shown in Fig. 5.  \nwas 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 Λ CDManalysis 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 di ff erent 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 Λ CDMmodel 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 countrate, 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 CR and cluster mass, required in Equation (9), can be described with the following equation:  \n* log CR CR , p GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> M , z + = log AX + bX ( z ) log M Mp + ex ( z ) . (25)  \nWeused the CR , p = 0 . 1 cts / s as the pivot value for the count rate, Mp = 2 × 10 14 M ⊙ , as the pivot value for the mass, and zp = 0 . 35 as the pivot value for the redshift. The redshift dependent slope of the scaling relation is expressed by the term bX ( z ), and the redshift evolution of the scaling relation is expressed by ex ( z ). Explicitly, the former is expressed as:  \nbX ( z ) = BX + FX log 1 + z 1 + zp . (26)  \nFig. 5: Posteriors on the CR -M and λ -M scaling relation parameters. A strong correlation between X-ray and optical scaling relation is visible on the corner plots. The best-fit values and their uncertainties are provided in Table 2. The solid and shaded orange contour represents the 68% and 95% contour levels respectively.  \n<!-- image -->  \nHere, BX is the normalization, and FX allows the slope to be redshift dependent. The self-similar model of Kaiser (1986) predicts the FX slope to be zero. In this work, we set the prior for FX parameter centered around zero to avoid any potential bias in our slopes. The redshift evolution term, ex ( z ), is given by:  \nex ( z ) = DX log dL ( z ) dL ( zp ) + EX log E ( z ) E ( zp ) + GX log 1 + z 1 + zp , (27)  \nfor which EX = 2 and DX = -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 GX parameter. The last term, GX , quantifies the deviation from the self-similar model in redshift evolution (e.g., as used in Bulbul et al. 2019).  \nThe intrinsic scatter of our scaling relations is represented by a single value, σ X , independent of mass and redshift, similar to the studies in the literature on scaling relations from X-ray observations, (for example, Mantz et al. 2016; Schellenberger & Reiprich 2017a; Lovisari et al. 2020, 2021; Du ff y 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 Xray 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  \nshows 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 CR -M scaling relations are provided in Table 2, and their posteriors are shown in Fig. 5.  \nThe only other works in the literature that use eROSITA Xray 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 Ax = 0 . 64 + 0 . 04 -0 . 06 . The reported normalization parameters are Ax = 1 . 24 + 0 . 22 -0 . 19 and Ax = 1 . 33 + 0 . 26 -0 . 20 in Chiu et al. (2022, 2023) in ∼ 3 σ tension with our results. In particular, the di ff erence in normalization is due to the observed degeneracy with the scatter, see Fig. B.1 in Appendix B, 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 a ff ecting the X-ray selection function can be safely ignored (Brunner et al. 2022).  \nWe find that the slope of the power-law index of the mass term is Bx = 1 . 38 + 0 . 03 -0 . 04 . Overall, in terms of mass trends, we find a steeper slope than the self-similar mass trend ( Bx , ss = 1), while many scaling relation studies in the literature report even steeper slopes. We find shallower slope than the reported value of Bx = 1 . 58 + 0 . 17 -0 . 14 in Chiu et al. (2022) and that of Bx = 1 . 86 + 0 . 20 -0 . 15 in the count rate mass relations provided in Chiu et al. (2023). The discrepancy is likely due to di ff erent 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 LX -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 BL x = 1 . 45 ± 0 . 07 in Maughan et al. (2007), 1 . 61 ± 0 . 14 in Vikhlinin et al. (2009), 1 . 62 ± 0 . 12 in Pratt et al. (2009), 1 . 34 ± 0 . 18 in Eckmiller et al. (2011), 1 . 66 ± 0 . 22 in Lovisari et al. (2015), and 1 . 35 ± 0 . 07 in Schellenberger & Reiprich (2017b). Our slopes are consistent with these findings at the 1 σ level. Higher slopes are reported in SZ-selected massive cluster samples covering a wider redshift range ( BL x = 1 . 89 ± 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).  \nThe best-fitting scatter σ X , 0 . 98 + 0 . 05 -0 . 04 , in our sample is larger than the reported values from CR -M or LX -Mscaling 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,  \nwe 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:  \n* log M WL Mp GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> M , z + = b ( z ) + bM log M Mp ! . (28)  \nHere, b ( z ) is the redshift-dependent evolution of the weaklensing bias as a function of redshift. It is computed by interpolating simulations snapshots at redshifts: 0 . 24, 0 . 42, 0 . 68, and 0 . 95.  \nThe 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; Naiman et al. 2018; Nelson et al. 2019), 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 (ˆ gt | M WL , ˆ 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 bl for each snapshot, and a global mass trend b M.  \nAs 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 e ff ectively propagating all WL related systematic uncertainties to a distribution on the WL bias values bl and b M. To compress this information and ingest it into the cosmological pipeline as priors, we derive the mean WL bias µ l . The scatter of the Monte Carlo realizations around this mean is well captured by two principal components PC 1 , 2 , l (Grandis et al. 2024). For the mass trend we report mean µ bM , and standard deviation σ bM , i . We similarly fit the scatter in the WL mass, finding that its redshift trend is well described by one principal component.  \nSince 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 bM are survey dependent since the specific properties of each survey a ff ect the reconstructed weak lensing mass slightly di ff erently, 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 µ b , i ( z ), PC 1 , i ( z ), PC 2 , i ( z ), µ bM , i , and σ bM , i with i ∈ { DES, KiDS, HSC } . To emulate the scatter in the WLbias simulations, we introduce three additional random normal parameters, A WL, B WL, and C WL, which are independent of the surveys employed. In practice, b ( z ) and bM can be expressed as:  \nb ( z ) = µ b , i ( z ) + A WL PC 1 , i ( z ) + B WL PC 2 , i ( z ) , bM = µ bM , i + C WL σ bM , i , (29)  \nwhere µ 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 bM , 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 di ff erent mass trends in their WL biases.  \nAs 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:  \nlog σ 2 M WL = s ( z ) + sM log M Mp ! , (30)  \nwhere the simulations have shown that the mass component of the scatter is consistent with zero; we fix sM = 0 in further analysis. For the redshift evolution of the scatter, we follow a similar procedure for the redshift evolution of the bias b ( z ). As before, we take the values at simulation snapshots and interpolate. We use the survey-specific mean value and principal component µ s , i ( z ) and PCi ( z ) and a random normal parameter D WL and express s ( z ) as:  \ns ( z ) = µ s , i ( z ) + D WL PCi ( z ) , (31)  \nwhere the average scatter for each survey and principal component ( µ s , i ( z ) and PCi ( 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.  \nThis method e ff ectively 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 di ff erent 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 e ff ectively 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  \net 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:  \nGLYPH<10> log λ | M , z GLYPH<11> = log A λ + b λ ( z ) log M Mp ! + C λ log 1 + z 1 + zp ! , (32)  \nwhile the slope of the cluster mass term b λ ( z ) can be expanded as:  \nb λ ( z ) = B λ + D λ log 1 + z 1 + zp ! , (33)  \nwhere A λ is the normalization of the scaling relation, B λ is the mass slope, C λ is the redshift evolution, and D λ is the redshift evolution of the mass slope.  \nIn 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-torichness scaling relation and the parameters of the mass-to-count rate scaling relation degenerate with each other.  \nThe best-fit parameters of the richness-M scaling relations are provided in Table 2, and their posteriors are shown in Fig. 5. In the case of the richness-to-mass scaling relation, Equation (32), the best-fit value of the mass slope, B λ = 1 . 01 + 0 . 02 -0 . 03 , is consistent with the majority of reported findings in the literature within 1 σ level, e.g., B λ = 0 . 99 + 0 . 06 -0 . 07 ± 0 . 04 and B λ = 1 . 07 + 0 . 22 -0 . 20 in the CODEX sample (Capasso et al. 2019; Ider Chitham et al. 2020, respectively), B λ = 1 . 02 ± 0 . 08 in SPT sample (Bleem et al. 2020), B λ = 1 . 028 + 0 . 043 -0 . 037 ± 0 . 04 and B λ = 1 . 03 ± 0 . 10 in the DES Y1 and SPT samples Costanzi et al. (2021); Grandis et al. (2021b), B λ = 1 . 05 + 0 . 13 -0 . 12 in Chiu et al. (2023) in the eFEDS sample within the stated uncertainities. Our measurements are 8 σ away from the value reported by the DES collaboration on Y1 cosmology analysis, where they find B λ = 0 . 737 ± 0 . 028 ± 0 . 004 (McClintock et al. 2019; Dark Energy Survey Collaboration et al. 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 et al. 2020).  \nWefindasample lower scatter in the λ -M scaling relation of 0 . 3 ± 0 . 03 compared to the relation between CR -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).  \nThe redshift evolution term has a best-fit value of C λ = -0 . 49 + 0 . 14 -0 . 15 shows a slight redshift trend in our sample. The bestfit value is consistent with the results from the CODEX sample for a sample with spectroscopic member galaxies λ > 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 σ level. In the same CODEX sample, the negative trend in redshift evolution becomes stronger as the authors change the sample with clusters with a larger number of member galaxies, for instance, λ > 20 ( -1 . 00 + 0 . 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 ± 0 . 30 ± 0 . 06 for a sample with richness cut of λ > 40. In contrast, the DES / SPT analyses reported in Costanzi et al. (2021) 0 . 95 ± 0 . 3 is in tension with our results at the > 5 σ level when the SPT number counts are included in the analysis, while consistent with our finding at < 2 σ level, -0 . 02 ± 0 . 34, when only SPT observable-to-mass relation is used (Costanzi et al. 2021). The reported redshift trends or slopes on the richness-mass relations depend on the cuts applied to the optically selected samples because of the rapid change in contamination fraction as a function of richness (Dark Energy Survey Collaboration et al. 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 e ff ect of di ff erent selection effects on the best-fit values of this relation and their comparisons.  \nWe note that the cluster richness is a sum of the color-based membership probability of the galaxies in the field (Ryko ff et al. 2014). In Kluge et al. (2024), we show that the eROMaPPer returns di ff erent richness measurements because of the di ff erent luminosity thresholds adopted for the combination in the different optical filter bands, for instance, grz , griz , grzw 1, due to variations in membership probabilities. Despite these di ff erences, Kluge et al. (2024) demonstrate that there is a simple constant scaling factor between these richnesses, for example, λ grz = S norm λ griz = S ' norm λ grzw 1, without any obvious dependence on redshift. This result is also supported by the reported comparison of the richness measured independently in the literature. For the eRASS1 clusters detected in the common area of the DES Y1 survey (Dark Energy Survey Collaboration et al. 2020), the comparisons between the richness of eRASS1 1 clusters measured in the grz band and of DES Y1 clusters in the grizy band shows a scaling factor that has no dependency on redshift; as such the observed relation is λ grz = 0 . 79 λ DES Y1 Kluge et al. (see 2024, for further details). Therefore, this demonstrates that all the comparisons of richness to mass scaling relations we have performed above are reasonable because richness scaling a ff ects 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 di ff erent 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  \ncombined, the dark sector comprises 95% of the energy density of the Universe (Planck Collaboration et al. 2020a). In the Λ CDMframework, the parameters that describe the Universe are the Hubble constant H 0, the total matter density Ω m (including dark matter and baryons), the spectral index of the primordial power spectrum ns , 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 Universem t 0, the critical density of the Universem ρ 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 16 , H 0, ns , and Ω b.  \nIn 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 w CDM cosmology, we allow the dark energy equation of state to be free and leave summed neutrino masses fixed to 0 eV. In the ν CDMcosmology, we allow summed neutrino masses to be free and leave the dark energy equation of state fixed to -1. We also provide joint constraints on the equation of the state of dark energy and the sum of neutrino masses in the flat ν w CDMcosmological model.', '7.1. Goodness of fit': 'In principle, determining the goodness of fit of the various cosmology models is not a trivial process in our analysis. The cosmology likelihood is not a simple distribution but comprises various components following di ff erent 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 distribution. Adding the likelihood for the mixture model complicates the process even further.  \nTo assess the goodness of fit of the standard Λ CDM cosmological model to the cluster abundance, richness, and count-rate distributions, we performed 1000 weighted 17 realizations from best-fitting parameters extracted from the posterior distribution and calculate the χ 2 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 simulated number of objects in each bin of the observables is relevant in assessing the goodness of fit. In Fig. 6, we provide the distribution of the observables employed in the analysis: count rate, redshift, and richness, and the corresponding best-fit models. Our methodology also allows us to simultaneously distinguish the relative contributions of the clusters and the contamination (AGN and RS) to the overall total model. For illustration purposes, the Pearson residuals are widely used in assessing the goodness of fit in the fits of the data with Poissonian nature statistics, which is also shown in Fig. 6. As evident from the figure, the total model shown in the shaded region in orange produces a good fit to the data with χ 2 = 30 . 3 for the redshift distribution, χ 2 = 24 . 5 for the count rate distribution, and  \nRichness  \nFig. 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.  \n<!-- image -->  \nχ 2 = 71 . 5 for richness distribution, all three with 26 degrees of freedom (d.o.f.).  \nArticle number, page 20 of 41  \nThe 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 χ 2 values of χ 2 DES = 173 . 2 + 43 . 2 -27 . 6 (150 d.o.f.), χ 2 HSC = 140 . 0 + 3 . 2 -6 . 2 (70 d.o.f.), and χ 2 KiDS = 60 . 2 + 9 . 7 -5 . 6 (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 χ 2 distributions are obtained from the higher 1 σ (upper edges of the shaded region in Fig. 7) and the lower 1 σ (lower edge of the shaded region in Fig. 7) of the predicted mean shear profile, respectively. We note that HSC does not yield a good fit to our data. We note from the middle panel in Fig. 7) that this is mostly due to the highest redshift low count rate bins. To assess the robustness of our results and conclusions below, we have checked that running the analysis without the inclusion of HSC does not a ff ect 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).  \nOur strategy in analyzing cluster number counts involves utilizing weak lensing tangential shear information that allows us to simultaneously fit and constrain cosmological parameters and the parameters of the scaling law between cluster mass and observables. We perform extensive checks on the influence of the individual WL survey measurements in our analysis to ensure that these surveys do not introduce a significant bias in the scaling relations or fit for the cluster number counts presented in the companion papers; Grandis et al. (2024) and Kleinebreil et al. (2024). To this end, the scaling relation between count rate to mass is tested on the individual weak lensing datasets by setting the first two terms in the likelihood, in Equation (6), to zero and fitting only the last two terms. We find that the overall model is consistent between the three surveys, and the scaling relation has enough degrees of freedom (see Table 2 for free parameters) to be essentially fixed by the mass calibration. Furthermore, we highlight that using three independent WL surveys simultaneously safeguards us against possible hidden systematics from any individual WL survey.', '7.2. Constraints on the Λ CDM cosmology': "Cluster number counts can provide powerful constraints on the cosmological parameters, especially on Ω m and σ 8, when the selection e ff ects are adequately taken into account, and the bias in the mass calibration process is well-controlled. In our baseline Λ CDMcosmology analysis, we freeze the dark energy equation of state to w = -1, and the remaining best-fit parameters, Ω m and σ 8, and S 8 are:  \nΩ m = 0 . 29 + 0 . 01 -0 . 02 , σ 8 = 0 . 88 ± 0 . 02 , S 8 = 0 . 86 ± 0 . 01 . (34)  \nThe posterior distributions of Ω m and σ 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  \nCR  \nFig. 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).  \n<!-- image -->  \ncluster number count experiments utilizing the WL mass calibration in the literature, e.g., WtG (Mantz et al. 2015), SPTSZ (Bocquet et al. 2019), XXL (Garrel et al. 2022), and eFEDS (Chiu et al. 2023). We find consistent results within 2 -3 σ with the results in the literature, except for the cosmology results from DES cluster abundance measurements Dark Energy Survey Collaboration et al. (2020). In Dark Energy Survey Collaboration et al. (2020), authors report a value of σ 8 = 0 . 85 + 0 . 04 -0 . 06 consistent with our result, whereas their best-fit Ω 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 di ff erence could be due to a biased massrichness 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  \ncan then calibrate our masses with state-of-the-art wide area WL surveys (KiDS, DES, HSC) with a high S / N.  \nWe compare our findings from cluster abundances with the recently published results from other cosmology probes, for instance, CMB measurements from Planck CMB(Planck Collaboration et al. 2020a), the joint cosmic shear constraints from DES and KiDS (Dark Energy Survey and Kilo-Degree Survey Collaboration et al. 2023), the 3 × 2 correlation constraints from DES (Dark Energy Survey Collaboration et al. 2022), lensing measurements from ACT DR6 + BAO (Qu et al. 2024), and lensing measurements from Planck (Planck Collaboration et al. 2020b) in Fig. 9. Both σ 8 and Ω m parameter measurements from the eROSITA cluster abundances agree well within 2 σ confidence interval with CMB (Planck Collaboration et al. 2020b) and lensing measurements (Dark Energy Survey Collaboration et al. 2022; Qu et al. 2024; Planck Collaboration et al. 2020b); however, they are in slight tension with the constraints from the cosmic shear measured from the combination of DES and KiDS data (Dark Energy Survey Collaboration et al. 2022; Dark Energy Survey and Kilo-Degree Survey Collaboration et al. 2023).  \nThe correlation and the correlation coe ffi cients between the free parameters in the fit are shown in Figs. B.2 and B.3 in Appendix B. We find that the terms representing the slopes ( BX and  \nFig. 8. Posterior distributions on the parameters Ω m and σ 8 from the Λ CDM fit on eRASS1 data shown in orange. In blue, we show the Planck CMB 2020 constraints without combination with BAO and SNe Ia Planck Collaboration et al. (2020a). 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 et al. 2020), SPT-SZ (Bocquet et al. 2019), eFEDS (Chiu et al. 2023), XXL (Garrel et al. 2022), and SPT (Bocquet et al. 2024) surveys.  \n<!-- image -->  \nFig. 9. Posterior distributions on Ω m and σ 8 from the Λ CDM fit to eRASS1 data are shown in orange. Similarly, the Planck CMB 2020 best-fit contours and eRASS1 with Planck CMB are shown in blue and red (Planck Collaboration et al. 2020a). 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 et al. 2023), DES 3x2 clustering constraints (Dark Energy Survey Collaboration et al. 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 et al. 2020b), and SDSS BAO RSD (Alam et al. 2021).  \n<!-- image -->  \nB λ ) and the redshift evolution ( GX and C λ ) of the scaling relations between X-ray observable and richness are correlated with Ω m. In contrast, these terms are anti-correlated with σ 8. Interestingly, S 8 shows a strong correlation with the amplitude of the scaling relations ( AX and A λ ), 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 di ff erences between the chains and assume that it is distributed according to a multivariate normal distribution with data  \ninferred mean and covariance. This allows for a calculation of the χ 2 = ¯ µ T C -1 ¯ µ ; and 3) we computed the survival function 18 of the χ 2 distribution.  \nThis method allows for an immediate estimate of the departure from the best fit in terms of its standard deviation. We find that our measurements on Ω m and σ 8 agree with the Planck Collaboration et al. (2020a)'s CMB measurements (TT, TE, EE + low l + lowE combination; we use Planck CMB 2020 hereafter) at the 2 . 3 σ confidence level.  \nTo break the degeneracy between Ω m and σ 8 and have a fair comparison with the Planck CMB2020 constraints, we combine the two probes given the consistency between the two. The combination of eROSITA cluster abundances and Planck CMB2020  \nchains yields the following cosmological parameters:  \nΩ m = 0 . 32 ± 0 . 01 , σ 8 = 0 . 82 ± 0 . 01 , S 8 = 0 . 85 ± 0 . 01 . (35)  \nWeshow the combined results in thick red contours in Figs. 8 and 9. Compared to the previous eROSITA-only results in Equation (34), we see a shift to a higher value of Ω m and to a lower value for σ 8, a significant improvement on our constraints on the Ω m and σ 8 parameter plane. This is because the Planck CMB 2020 constraints on Ω m and σ 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 di ff erent 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 Λ CDManalysis, we set the neutrino mass to 0 eV. Figure B.4 in Appendix B 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).  \nA 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 Ω m and σ 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 α best value that minimizes the standard deviation of σ 8 × ( Ω m / 0 . 3) α along the direction of our degeneracy, is α best = 0 . 42 ± 0 . 02. Our best-fit value is close to 0.5, the canonical value of α used in the definition of S 8 in Planck Collaboration et al. (2020a). The cluster cosmology experiments, e.g., Bocquet et al. (2019) reports α best = 0 . 2 for the SPT-SZ sample, while Chiu et al. (2023) finds α best = 0 . 3 for the eFEDS clusters, a smaller value than what we find in this work. On the other hand, the power spectrum analysis of the shear measurements of the HSC field produces an α best value of 0.45, consistent with our results Hikage et al. (2019). Regarding the S 8 constraints, we obtain a degeneracy between σ 8 and Ω m, which is almost perfectly aligned with the slope in the S 8 definition. This makes the constraints on σ 8 tighter in our experiment than other surveys. Similarly, the recent SPT constraints show the slope α best ∼ 0 . 25 (Bocquet et al. 2024), and indeed their optimal parameter combination σ 8 × ( Ω 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.  \nFigure 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 et al. 2023), clustering (Dark Energy Survey Collaboration et al. 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 et al. 2020a), ACT and Planck lensing constraints (Qu et al. 2024; Madhavacheril et al. 2024), the HSC clustering results (Sugiyama et al. 2023), and the cluster count measurements based on the XXL survey within 2 σ confidence interval (Planck Collaboration et al. 2020a; Garrel et al. 2022). How-  \nver, it is in significant tension ( ≥ 3 σ level) with constraints from other probes, in particular from cosmic shear from DES, HSC, KiDS surveys individually, and the combined KiDS and DES cosmic shear measurements (Amon et al. 2022; Secco et al. 2022; van den Busch et al. 2022; Li et al. 2023), DES 3 × 2 clustering (Sugiyama et al. 2023), and the SDSS BAO and redshift space distortions (RSD) (Alam et al. 2021) measurements. Given the consistency of our results with the CMB measurements, we combine our results with Planck 2020 CMB chains. The results of this combined analysis result in a negligible improvement on S 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.  \nAs 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 e ff ect of the assumption on the halo mass function on our results, we fit the Λ CDM to our data with the Tinker et al. (2010) and Despali et al. (2016) halo mass functions. The results of the comparisons are shown in Fig.B.6. We observe a slight decrease in the parameter σ 8 from 0 . 88 ± 0 . 02 to 0 . 86 ± 0 . 02 when the Tinker et al. (2010) and Despali et al. (2016) halo mass functions are used. The shift is statistically consistent within 1 σ intervals with the Tinker et al. (2008) results. The change in the halo mass function does not significantly impact the parameter Ω 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 ν 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 di ff erent flavors, electron, muon, and tau neutrinos, it has become clear that neutrino eigenstates must have non-zero mass (Fukuda et al. 1998; Ahmad et al. 2002). The theoretical predictions indicate a hierarchical relation between the mass eigenstates of the three species of neutrinos; however, their ordering in mass remains elusive. Based on the normal neutrino mass hierarchy model, the third neutrino mass eigenstate, associated with the tau neutrino, is the heaviest among the three. This is followed by the eigenstate of the muon neutrino, which has an intermediate mass, and that of the electron neutrino, namely, m τ > m µ > me . The ordering is swapped in the inverted mass hierarchy model such that m τ < m µ < me . 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 P m ν > 0 . 059 eV for normal mass hierarchy and P m ν > 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  \nFig. 10. Comparison of the constraints on S 8 = σ 8 × √ Ω m / 0 . 3 with literature results from SPTSZ (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 CMB2020 (Planck Collaboration et al. 2020a), 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 et al. 2020b), from DES 3x2 correlation function (Dark Energy Survey Collaboration et al. 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 et al. 2023), from SDSS BAO and RSD (Alam et al. 2021). All the errors shown are at the 1 σ level. The vertical line and grey shaded areas represent our measurement's location and 1 -2 σ confidence level.  \n<!-- image -->  \ndecay of tritium imply the sum to be P m ν < 1 . 1 eV at 90% confidence level, a narrow range for the allowed summed mass (Aker et al. 2019). The space-based Planck CMBmeasurements provide a similar upper limit of P m ν < 0 . 26 eV at 95% confidence level (Planck Collaboration et al. 2020a).  \nAs 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 U (0 eV , 1 eV). Including massive neutrinos in cosmological analysis requires a fine-tuned interaction with the choice of the halo mass function. As shown in Costanzi et al. (2013), a few changes to the use of the matter power spectrum are fundamental to properly make use of Tinker et al. (2008) halo mass function in the presence of massive neutrinos. As the pyccl 3.0.0 library of Chisari et al. (2019) available on GitHub 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 a ff ecting our posterior distributions. The cosmological constraints with free neutrino mass components are:  \nΩ m = 0 . 28 + 0 . 02 -0 . 02 , σ 8 = 0 . 88 ± 0 . 03 , S 8 = 0 . 85 ± 0 . 02 , X m ν < 0 . 43eV (95% CL) . (36)  \nThese results are visualized in Fig. 11. The upper limit to the sum of neutrino masses is P m ν < 0 . 43 eV. The addition of the  \n(Costanzi et al. 2013) correction to the power spectrum relaxes the upper limits on neutrino masses from 0.22 eV to 0.43 eV, as expected. Our results represent the tightest limits on the sum of neutrino masses from cluster abundance experiments; for instance, the SPT-SZ sample results in an upper limit of < 0 . 74 eV (95% confidence interval Bocquet et al. 2019), while the updated SPT results yield < 0 . 55 eV (95% confidence interval Bocquet et al. 2024). We verify that the values of the Ω m and σ 8 remain statistically consistent when P 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 σ level.  \nGiven the excellent agreement with Planck CMB measurements, the resulting cosmological parameters of two probes can be combined in a statistically meaningful way, enabling a much tighter measurement of the impact of massive neutrinos on both the formation and evolution of the large-scale structure and on the primordial density field. As performed in the Λ CDM analysis, we combine our results with the Planck CMB constraints to break the degeneracy between Ω m and σ 8. We obtain:  \nΩ m = 0 . 31 ± 0 . 01 , σ 8 = 0 . 83 ± 0 . 01 , S 8 = 0 . 84 ± 0 . 01 , X m ν < 0 . 14 eV (95% CL) (37)  \nSimultaneously fitting our measurements with Planck CMB 2020 likelihood chains yield an upper limit of P m ν < 0 . 14 eV, consistent with the results in the literature results with marginally consistent with the inverted mass hierarchy model, which requires P m ν > 0 . 101 eV, where just 18% of the eROSITAPlanck combined samples are above this limit and consistent with the inverted hierarchical model.  \nAs a final step, we combined eRASS1 cluster abundance measurements via an importance sampling with Planck CMB and with the lower limits from ground-based oscillation experiments (Tanabashi et al. 2018). In the case of a normal mass hierarchy scenario, the summed masses of P m ν = 0 . 09 + 0 . 04 -0 . 02 eV; while assuming the inverted mass hierarchy model, we obtain  \nFig. 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 P m ν < 0 . 43eV from cluster number counts alone.  \n<!-- image -->  \nsummed masses of P m ν = 0 . 12 + 0 . 03 -0 . 02 eV. It is interesting to point out (from Fig. C.1 in Appendix C) that the constraints on the mass of the lightest neutrino eigenstates are similar: for both mass hierarchy we obtain m 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 w CDM 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 et al. 2020a), 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 e ff ort invested into constraining the significant contribution of dark energy, its nature remains elusive. One of the main science drivers of eROSITA is to place constraints on its equation of state using the cluster abundances. We define w , the equation of state, as the ratio of the pressure of dark energy to its energy density ( w = p /ρ ) and assume that its density evolves with ρ ∝ a -3(1 + w ) in an isotropic Universe. Previous dark energy equation state measurements point to w ∼ -1 with no evolution with redshift (Astier et al. 2006; Beutler et al. 2011; Alam et al. 2017; Planck Collaboration et al. 2020a). In this section, we consider an extension to the standard Λ CDMbyallowing the equation of state of dark energy to w to vary.  \nTo constrain the value of w , we set the sum of neutrino masses back to zero ( P m ν = 0 eV), and we allow the w to be free with a uniform prior U ( -2 , -0 . 33) in the w CDM model. Naturally, a change in the equation of the state a ff ects 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:  \nΩ m = 0 . 28 ± 0 . 02 , σ 8 = 0 . 88 ± 0 . 02 , S 8 = 0 . 86 ± 0 . 02 , w = -1 . 12 ± 0 . 12 . (38)  \nOur fits appear to prefer w to be -1 . 12 ± 0 . 12, shown in Fig. 12. consistent with the results from Planck CMB (Planck Collaboration et al. 2020a), DES Y3 3x2 pt (Dark Energy Survey Collaboration et al. 2022), Pantheon SNe Ia (Scolnic et al. 2018, 2022; Brout et al. 2022), DES Y5 SNe Ia (DES Collaboration et al. 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 Λ CDMcase.  \nGiven 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  \n0.5  \nFig. 12. Best-fit contours of w and Ω m in the w CDMcosmological 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.  \n<!-- image -->  \nFig. 13: Correlation between w and the X-ray count rate to mass scaling relation normalization and redshift evolution. The correlation for the w CDM case is shown in orange and blue for the ν w CDMcase.  \n<!-- image -->  \nbreak degeneracies and obtain;  \nΩ m = 0 . 30 ± 0 . 01 , σ 8 = 0 . 87 ± 0 . 02 , S 8 = 0 . 87 ± 0 . 01 , w = -0 . 95 + 0 . 05 -0 . 04 . (39)  \nCompared with our results, the parameters Ω m, σ 8, and S 8 are not significantly a ff ected 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  \na significantly increased precision. w remains consistent with the standard Λ CDMvalue of -1 at 1 σ level.  \nDark energy is responsible for a late-time expansion of the Universe and, therefore, has a negligible e ff ect 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 CMB2020 and SDSS BOSS DR16 BAO results, respectively. Given the consistency, we combine eRASS1 chains with both Pantheon SNe, Planck CMB, and SDSS BAO data. Including these two extra constraints helps improve the precision of the parameters significantly. We measure;  \nΩ m = 0 . 31 ± 0 . 01 , σ 8 = 0 . 84 ± 0 . 01 , S 8 = 0 . 85 ± 0 . 01 , w = -1 . 06 ± 0 . 03 . (40)  \nWhen 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 ± 0 . 03 with 3% precision, thus improving the eROSITA only results by a factor of ∼ 3, and consistent with w = -1 at 2 σ level. The shift on the best-fit eROSITA parameters is less than 1 σ . In particular, Ω m is unchanged, σ 8 and S 8 shifts to a slightly lower value.', '7.5. Constraints on the ν w CDM cosmology': 'In connection with the previous sections, to test the e ff ect of the neutrino masses on the dark energy equation of state, we allow both w and P m ν to vary in our fits. The priors are set to uniform distributions for U ( -2 , -0 . 33) on w , and U (0 eV , 1 eV) for P m ν . We find:  \nΩ m = 0 . 27 ± 0 . 02 , σ 8 = 0 . 89 ± 0 . 03 , S 8 = 0 . 84 ± 0 . 02 , w = -1 . 11 ± 0 . 14 , X m ν < 0 . 44 eV (95% CL) .  \n(41)  \nFig. 14. Best-fit values for ν w CDMcosmological 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 Λ CDMmodel.  \n<!-- image -->  \nThis case combines the two previous cosmologies: as for the ν CDM and w CDM case, our results in the Ω m σ 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.  \nDue to the lack of public Planck CMB-only likelihood chains on the simultaneous fits on the w and P m ν in Planck CMB 2020, we combine our results with the readily available CMB + lensing + BAO in Planck Collaboration et al. (2016a) chains. Relying on the 1.6 σ agreement with our results, we find:  \nX (42)  \nΩ m = 0 . 29 ± 0 . 01 , σ 8 = 0 . 85 ± 0 . 02 , S 8 = 0 . 83 ± 0 . 01 , w = -1 . 15 ± 0 . 10 , m ν < 0 . 37 eV (95% CL) .  \nIn this case, the constraints on the w move closer to the value -1, with the summed neutrino masses with an upper limit of P m ν < 0 . 37 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 et al. (2016a) 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 ± 0 . 21 (Bocquet et al. 2019). In the same work, the reported sum of neutrino masses is P m ν = 0 . 50 ± 0 . 24 eV in the same combination of datasets, whereas we find tighter upper limits on P m ν in this work.  \n7.6. Model comparison  \nIn this paper, we fit our eRASS1 dataset with four di ff erent cosmological models, Λ CDM, ν CDM, w CDM, and ν w CDM. We compare the goodness of fit for these di ff erent models to assess whether the data prefer one of the cosmological models in this section. To achieve this, we use the Bayes factor (Je ff reys 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 20 directly as an output. Therefore, it is straightforward to use it to compute the Bayes factor 21 and Je ff reys scale to compare di ff erent cosmology fits (Je ff reys 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.  \nIn 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  \nTable 3: Bayesian evidence for the cosmological models tested.  \nNotes. To interpret this table, we need the Je ff reys scale (Je ff reys 1961): if the di ff erence between the two is greater than 5, then there is significant evidence to prefer the model with the largest Bayesian evidence; if the di ff erence 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.  \nof the standard Λ CDM case, see Equations (34), (36), (38), and (41).', '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 eROSITAICM 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.  \nContamination 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 e ff ort 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 in Appendix B 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 a ff ecting 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 e ff ects in our likelihood and, therefore, marginalize them out to obtain the presented precision cosmological constraints.  \nSelection 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  \nsimulations 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-o ff between modeling complexity and cosmology- and / or modeldependency. 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)).  \nWeak-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-ofsight (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 Nbody 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).  \nForward 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 a ff ects our detection e ffi ciency (see Fig. 16 in Predehl et al. 2021). Additionally, galactic absorption from atomic and molecular hydrogen (HI4PI Collaboration et al. 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 abil-  \nity 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 e ff ectively 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 a ff ected by the challenges in measuring galaxy shapes and photometric redshift. The key di ff erence comes from the theoretical uncertainties that must be accounted for when predicting the WL signals for cosmological inference.  \nWhile cosmic shear can be predicted very accurately at large scales, predictions at intermediate scales of the order of 10 Mpc could be a ff ected 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.  \nThese 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 una ff ected 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 a ff ected by the intermediate scale issues that a ff ect the cosmic shear measurements.  \nWe 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 analytically (Press & Schechter 1974). Detection through X-rays is one of  \nthe most reliable methods for locating clusters and for constructing pure cluster catalogs without any major bias because the Xray 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 e ff ects are minimized. The strict cuts we apply on the X-ray detection chain and the utilization of their follow-up optical observations can eliminate any leftover contamination in the cluster samples, reducing the contamination to < 6%. The remaining contamination is modeled through our mixture method in-situ (Bulbul et al. 2024; Kluge et al. 2024). The selection function is produced from the end-toend 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).  \nWeak 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 a ff ected 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 di ff erence between WL masses calibrated on gravity-only simulations and hydro-dynamical simulations is 2% (Grandis et al. 2021a), while the di ff erence between the matter power spectra between gravity-only and hydrodynamical simulations, but also among di ff erent 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 di ff erent sources of weak lensing signals. In light of these arguments, WL-calibrated number counts of Xray-selected clusters may well be regarded as more robust than cosmic shear analyses.", '8.3. Constraints on the total mass of neutrinos': 'The e ff ect of massive neutrinos on cluster abundance has been explored in this work. Particularly, they a ff ect 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 e ff ect is, therefore, the dumping of structure formation compared to a standard Λ CDM cosmological model with massless neutrinos.  \nIn 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 cluster sample found P m ν < 0 . 74 eV at 95% confidence level, and even CMB results from Planck Collaboration et al. (2020a) yield P m ν < 0 . 26 eV at 95% confidence level. The tightest constraints on the neutrino masses are obtained via combinations of multiple probes. For instance, combining Lymanα forest results and CMB measurements give P m ν < 0 . 15 eV (Rossi et al. 2015), while BAO and CMB together yields P m ν < 0 . 15 eV at 95% confidence level (Vagnozzi et al. 2017). More recent results from the Planck CMBandBOSSBAOmeasurements provide slightly  \ntighter upper limits of P m ν < 0 . 14 eV at 95% confidence level (Tanseri et al. 2022).  \nThe results in our work from eROSITA are competitive by themselves since we obtain P m ν < 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 combination of these two state-of-the-art cosmological probes that we obtain very tight constraints of P m ν < 0 . 14 eV at 95% confidence 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 CMB2020rule out at a 82% confidence level the inverted hierarchy model, which requires P m ν > 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 di ff erent neutrino mass hierarchy models of neutrino eigenstates. Furthermore, measuring the mass of the lightest neutrino eigenstate (combining our results with oscillation 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 e ff ort 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 et al. (2020a), 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).  \nAs 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 CMBanalysis (Planck Collaboration et al. 2020a). Still, the value we find, w = -1 . 12 ± 0 . 12 , is consistent at 1 σ level with the canonical Λ CDM value of -1, and is also consistent with previously reported literature results from DES Y3 3x2 pt ( w = -0 . 98 + 0 . 32 -0 . 20 , Dark Energy Survey Collaboration et al. 2022), SPT clusters ( w = -1 . 55 ± 0 . 41, Bocquet et al. 2019), SDSS eBOSS DR16 BAO ( w = -0 . 69 ± 0 . 15, Alam et al. 2021), Pantheon SNe Ia ( w = -1 . 09 ± 0 . 22, Scolnic et al. 2018), and Pantheon + SNe Ia ( w = -0 . 90 ± 0 . 14, Brout et al. 2022). We note that di ff erently 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 . 04 + 0 . 05 when combining with Pantheon + SNe Ia (Brout et al. 2022), and w = -1 . 06 ± 0 . 03 when Planck CMB(Planck Collaboration et al. 2020a) and SDSS BOSS DR16 BAO (Alam et al. 2017) are included.  \nShown 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  \nthe X-ray count rate to mass scaling relation, GX . 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 ∼ 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 Λ CDM cosmology, the neutrino masses, and the dark energy equation of state. Although we fully account for all known systematics e ff ects 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 semianalytical Tinker et al. (2008) mass function, is constructed for the standard FlatΛ CDM analyses. However, we utilize it even when the total mass of neutrinos P m ν 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.  \nAdditionally, 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 e ff ects 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 ⊙ . 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 ⊙ , 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 e ff ect 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 e ff ect in the WL bias calibration, as discussed in Grandis et al. (2021a) and Grandis et al. (2024).  \nFurther, the e ff ect 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 par-  \nticular, Fumagalli et al. (2021) clearly showed that for large area surveys with a mass limit around few 10 14 M ⊙ , such as SPT and eRASS1, the results from a Poisson likelihood are not a ff ected by the addition of super-sample covariance. In SPT cosmology analyses (Bocquet et al. 2019, 2024), the authors also ignore the e ff ect 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 a ff ected 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 supersample covariance in our cosmology analyses in the near future.  \nCalibrating the point-spread function and mirror e ff ective area of eROSITA is an ongoing e ff ort (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 di ff erence does not a ff ect 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.  \nThe 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.  \nThus 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 σ contour in the Ω m σ 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-SZpolECS-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 et al. 2020). Overall, eROSITA provides a factor of 5 to 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.10.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 sophisti-  \ncated 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 2 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 e ff ects 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 AllSky survey, leads to better-controlled systematic uncertainties than those available in previous literature works. Thus, it allows for reliable and stringent constraints to be placed on the cosmological parameters.  \nWe 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 a ff ecting the mass reconstruction) does not improve our cosmological results significantly because it is the precision (and not the value) that matters.  \nFrom the eRASS1 sample and cosmological analysis presented in this paper, we derived the σ 8 = 0 . 88 ± 0 . 02 with 3% precision and Ω m = 0 . 29 + 0 . 01 -0 . 02 with ∼ 5% precision for the Λ CDM scenario. To remove the degeneracy between σ 8 and Ω m and to have a fair comparison between CMB results, we combined these two parameters and measure the S 8 = 0 . 86 ± 0 . 01 parameter with a 2% precision. Our results, from a late-time probe, are competitive and consistent with Planck CMB2020 measurements (Planck Collaboration et al. 2020a), 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 σ 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).  \nTo test the e ff ect 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 P m ν < 0 . 43 eV (95% upper limit), tighter than the upper limits provided by CMB measurements, for instance, P m ν < 0 . 66 eV by (Komatsu et al. 2011) and P m ν < 0 . 26 eV by Planck Collaboration et al. (2020a) (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 P m ν < 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 P m ν < 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 P m ν = 0 . 09 + 0 . 04 -0 . 02 eV and P m ν = 0 . 12 + 0 . 03 -0 . 02 eV , assuming a normal and inverted mass hierarchy scenario for neutrino eigenstates.  \nWe 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 Ω m, σ 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 P m ν 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.  \nFinally, 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Λ 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 Ω m and σ 8 rely mainly on the mass slope and the redshift evolution of the scaling relation.  \nWith 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 to 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 22 .  \nAcknowledgements. 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 by MPE, with contributions from the Dr. Karl Remeis Observatory Bamberg & ECAP (FAU Erlangen-Nuernberg), the University of Hamburg Observatory, the Leibniz Institute for Astrophysics Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the University of Tübingen, with the support of DLR and the Max Planck Society. The Argelander Institute for Astronomy of  \nthe University of Bonn and the Ludwig Maximilians Universität Munich also participated in the science preparation for eROSITA. The eROSITA data shown here were processed using the eSASS software system developed by the German eROSITA consortium.  \nV. Ghirardini, E. Bulbul, A. Liu, C. Garrel, S. Zelmer, and X. Zhang acknowledge financial support from the European Research Council (ERC) Consolidator Grant under the European Union's Horizon 2020 research and innovation program (grant agreement CoG DarkQuest No 101002585). N. Clerc was financially supported by CNES. T. Schrabback and F. Kleinebreil acknowledge support from the German Federal Ministry for Economic A ff airs and Energy (BMWi) provided through DLR under projects 50OR2002, 50OR2106, and 50OR2302, as well as the support provided by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant 415537506. V.G. also thanks Margherita Ghirardini for her invaluable joyful support.  \nThe Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID #2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop. ID #2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF's NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. Pipeline processing and analyses of the data were supported by NOIRLab and the Lawrence Berkeley National Laboratory (LBNL). The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du'ag (Kitt Peak), a mountain with particular significance to the Tohono O'odham Nation.  \nFunding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology FacilitiesCouncil of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft, and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC / CSIC), the Institut de Física d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the OzDES Membership Consortium, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University. Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 177.A-3016, 177.A-3017, 177.A-3018 and 179.A-2004, and on data products produced by the KiDS consortium. The KiDS production team acknowledges support from: Deutsche Forschungsgemeinschaft, ERC, NOVA and NWO-M grants; Target; the University of Padova, and the University Federico II (Naples). N.O. and S.M. acknowledge JSPS KAKENHI Grant Number JP19KK0076. The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from the Japanese Cabinet O ffi ce, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University. This paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software at http: // dm.lsst.org This paper is based [in part] on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center (ADC) at NAOJ. Data analysis was in part carried out with the cooperation of Center for Computational Astrophysics (CfCA), NAOJ. We are honored and grateful for the opportunity of observing the Universe from Maunakea, which has the cultural, historical and natural significance in Hawaii. I. C. ac-  \npport from the National Science and Technology Council in Taiwan (Grant NSTC 111-2112-M-006-037-MY3).  \nThis work made use of the following Python software packages: SciPy 23 (Virtanen et al. 2020), Matplotlib 24 (Hunter 2007), Astropy 25 (Astropy Collaboration et al. 2022), NumPy 26 (Harris et al. 2020), CAMB (Lewis & Challinor 2011), pyCCL 27 (Chisari et al. 2019), GPy 28 (GPy since 2012), climin 29 (Bayer et al. 2015), ultranest 30 (Buchner 2021)", 'References': "Adami, C., Giles, P., Koulouridis, E., et al. 2018, A&A, 620, A5  \n- Ahmad, Q. R., Allen, R. C., Andersen, T. C., et al. 2002, Phys. Rev. Lett., 89, 011301\n- Aihara, H., Arimoto, N., Armstrong, R., et al. 2018, PASJ, 70, S4\n- Aird, J., Coil, A. L., Georgakakis, A., et al. 2015, MNRAS, 451, 1892\n- Aker, M., Altenmüller, K., Arenz, M., et al. 2019, Phys. Rev. Lett., 123, 221802 Alam, S., Ata, M., Bailey, S., et al. 2017, MNRAS, 470, 2617\n- Alam, S., Aubert, M., Avila, S., et al. 2021, Phys. Rev. D, 103, 083533\n- Amendola, L., Appleby, S., Avgoustidis, A., et al. 2018, Living Reviews in Relativity, 21, 2\n- Amon, A., Gruen, D., Troxel, M. A., et al. 2022, Phys. Rev. D, 105, 023514 Anderson, M. E., Gaspari, M., White, S. D. M., Wang, W., & Dai, X. 2015, MNRAS, 449, 3806\n- Applegate, D. E., von der Linden, A., Kelly, P. L., et al. 2014, MNRAS, 439, 48 Asgari, M., Lin, C.-A., Joachimi, B., et al. 2021, A&A, 645, A104 Astier, P., Guy, J., Regnault, N., et al. 2006, A&A, 447, 31 Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935,\n- 167\n- Athar, M. S., Barwick, S. W., Brunner, T., et al. 2022, Progress in Particle and Nuclear Physics, 124, 103947\n- Bahar, Y. E., Bulbul, E., Clerc, N., et al. 2022, A&A, 661, A7\n- Bahar, Y. E., Bulbul, E., Ghirardini, V., et al. 2024, arXiv e-prints, arXiv:2401.17276\n- Bahé, Y. M., McCarthy, I. G., & King, L. J. 2012, MNRAS, 421, 1073  \nBarnes, D. J., Kay, S. T., Bahé, Y. M., et al. 2017, MNRAS, 471, 1088  \n- Bashinsky, S. & Seljak, U. 2004, Phys. Rev. D, 69, 083002\n- Bayer, J., Osendorfer, C., Diot-Girard, S., Rueckstiess, T., & Urban, S. 2015, TUM, Tech. Rep.\n- Becker, M. R. & Kravtsov, A. V. 2011, ApJ, 740, 25\n- Beringer, J., Arguin, J. F., Barnett, R. M., et al. 2012, Phys. Rev. D, 86, 010001\n- Beutler, F., Blake, C., Colless, M., et al. 2011, MNRAS, 416, 3017\n- Bleem, L. E., Bocquet, S., Stalder, B., et al. 2020, ApJS, 247, 25\n- Bocquet, S., Dietrich, J. P., Schrabback, T., et al. 2019, ApJ, 878, 55\n- Bocquet, S., Grandis, S., Bleem, L. E., et al. 2024, arXiv e-prints, arXiv:2401.02075\n- Bocquet, S., Saro, A., Dolag, K., & Mohr, J. J. 2016, MNRAS, 456, 2361\n- Brout, D., Scolnic, D., Popovic, B., et al. 2022, ApJ, 938, 110\n- Brunner, H., Liu, T., Lamer, G., et al. 2022, A&A, 661, A1\n- Buchner, J. 2021, The Journal of Open Source Software, 6, 3001\n- Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125  \nBulbul, E., Chiu, I. N., Mohr, J. J., et al. 2019, ApJ, 871, 50  \n- Bulbul, E., Liu, A., Kluge, M., et al. 2024, A&A, 685, A106\n- Bulbul, E., Liu, A., Pasini, T., et al. 2022, A&A, 661, A10\n- Burns, J. O., Skillman, S. W., & O'Shea, B. W. 2010, ApJ, 721, 1105\n- Capasso, R., Mohr, J. J., Saro, A., et al. 2019, MNRAS, 486, 1594  \nCapasso, R., Mohr, J. J., Saro, A., et al. 2020, MNRAS, 494, 2736  \n- Chisari, N. E., Alonso, D., Krause, E., et al. 2019, CCL: Core Cosmology Library, Astrophysics Source Code Library, record ascl:1901.003\n- Chiu, I. N., Ghirardini, V., Liu, A., et al. 2022, A&A, 661, A11\n- Chiu, I. N., Klein, M., Mohr, J., & Bocquet, S. 2023, MNRAS, 522, 1601\n- Chuang, C.-H., Yepes, G., Kitaura, F.-S., et al. 2019, MNRAS, 487, 48\n- Clerc, N., Comparat, J., Seppi, R., et al. 2024, A&A, 687, A238  \nClerc, N. & Finoguenov, A. 2023, in Handbook of X-ray and Gamma-ray Astro-  \nphysics. Edited by Cosimo Bambi and Andrea Santangelo, 123  \n- Comparat, J., Eckert, D., Finoguenov, A., et al. 2020, The Open Journal of Astrophysics, 3, 13  \nComparat, J., Luo, W., Merloni, A., et al. 2023, A&A, 673, A122  \n- 24 https: // matplotlib.org /\n- 25 https: // www.astropy.org /\n- 26 https: // numpy.org /\n- 27 https: // github.com / LSSTDESC / CCL\n- 28 https: // github.com / She ffi eldML / GPy\n- 29 https: // github.com / BRML / climin\n- 30 https: // github.com / JohannesBuchner / UltraNest /\n- Comparat, J., Merloni, A., Salvato, M., et al. 2019, MNRAS, 487, 2005  \nComparat, J., Prada, F., Yepes, G., & Klypin, A. 2017, MNRAS, 469, 4157  \n- Cooke, R. J., Pettini, M., Nollett, K. M., & Jorgenson, R. 2016, ApJ, 830, 148\n- Cooke, R. J., Pettini, M., & Steidel, C. C. 2018, ApJ, 855, 102\n- Corless, V. L. & King, L. J. 2007, MNRAS, 380, 149\n- Costanzi, M., Rozo, E., Ryko ff , E. S., et al. 2019a, MNRAS, 482, 490\n- Costanzi, M., Rozo, E., Simet, M., et al. 2019b, MNRAS, 488, 4779  \nCostanzi, M., Saro, A., Bocquet, S., et al. 2021, Phys. Rev. D, 103, 043522  \n- Costanzi, M., Villaescusa-Navarro, F., Viel, M., et al. 2013, J. Cosmology Astropart. Phys., 2013, 012\n- Dalal, R., Li, X., Nicola, A., et al. 2023, Phys. Rev. D, 108, 123519 Dark Energy Survey and Kilo-Degree Survey Collaboration, Abbott, T. M. C., Aguena, M., et al. 2023, The Open Journal of Astrophysics, 6, 36\n- Dark Energy Survey Collaboration, Abbott, T. M. C., Aguena, M., et al. 2020, Phys. Rev. D, 102, 023509\n- Dark Energy Survey Collaboration, Abbott, T. M. C., Aguena, M., et al. 2022, Phys. Rev. D, 105, 023520  \nDauser, T., Falkner, S., Lorenz, M., et al. 2019, A&A, 630, A66  \nde Haan, T., Benson, B. A., Bleem, L. E., et al. 2016, ApJ, 832, 95  \n- DES Collaboration, Abbott, T. M. C., Acevedo, M., et al. 2024, arXiv e-prints, arXiv:2401.02929  \nDespali, G., Giocoli, C., Angulo, R. E., et al. 2016, MNRAS, 456, 2486  \nDey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168  \n- Di Valentino, E., Mena, O., Pan, S., et al. 2021, Classical and Quantum Gravity, 38, 153001\n- Dietrich, J. P., Bocquet, S., Schrabback, T., et al. 2019, MNRAS, 483, 2871 Dolag, K., Bartelmann, M., & Lesch, H. 1999, A&A, 348, 351  \nDolag, K., Vazza, F., Brunetti, G., & Tormen, G. 2005, MNRAS, 364, 753  \n- Du ff y, R. T., Logan, C. H. A., Maughan, B. J., et al. 2022, MNRAS, 512, 2525  \nEckert, D., Ghirardini, V., Ettori, S., et al. 2019, A&A, 621, A40  \n- Eckmiller, H. J., Hudson, D. S., & Reiprich, T. H. 2011, A&A, 535, A105 Eddington, A. S. 1913, MNRAS, 73, 359\n- Esteban, I., Gonzalez-Garcia, M. C., Maltoni, M., Schwetz, T., & Zhou, A. 2020, Journal of High Energy Physics, 2020, 178\n- Fang, W., Hu, W., & Lewis, A. 2008, Phys. Rev. D, 78, 087303\n- Fixsen, D. J., Cheng, E. S., Gales, J. M., et al. 1996, ApJ, 473, 576\n- Freyberg, M., Perinati, E., Pacaud, F., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11444, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 114441O  \nFukuda, Y., Hayakawa, T., Ichihara, E., et al. 1998, Phys. Rev. Lett., 81, 1562 Fumagalli, A., Saro, A., Borgani, S., et al. 2021, A&A, 652, A21 Garrel, C., Pierre, M., Valageas, P., et al. 2022, A&A, 663, A3 Gatti, M., Giannini, G., Bernstein, G. M., et al. 2022, MNRAS, 510, 1223 Giblin, B., Heymans, C., Asgari, M., et al. 2021, A&A, 645, A105 Giles, P. A., Maughan, B. J., Pacaud, F., et al. 2016, A&A, 592, A3 Giles, P. A., Romer, A. K., Wilkinson, R., et al. 2022, MNRAS, 516, 3878 GPy. since 2012, GPy: A Gaussian process framework in python, http:// github.com/SheffieldML/GPy  \n- Grandis, S., Bocquet, S., Mohr, J. J., Klein, M., & Dolag, K. 2021a, MNRAS, 507, 5671\n- Grandis, S., Ghirardini, V., Bocquet, S., et al. 2024, A&A, 687, A178  \nGrandis, S., Klein, M., Mohr, J. J., et al. 2020, MNRAS, 498, 771  \n- Grandis, S., Mohr, J. J., Costanzi, M., et al. 2021b, MNRAS, 504, 1253\n- Grandis, S., Mohr, J. J., Dietrich, J. P., et al. 2019, MNRAS, 488, 2041\n- Gruen, D., Seitz, S., Brimioulle, F., et al. 2014, MNRAS, 442, 1507\n- Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357  \nHensman, J., Fusi, N., & Lawrence, N. D. 2013, ArXiv, abs  \n/  \n1309.6835  \nHensman, J., Matthews, A., & Ghahramani, Z. 2015, in Proceedings of the Eigh-  \nteenth International Conference on Artificial Intelligence and Statistics, ed.  \nG. Guy Lebanon & S. V. N. Vishwanathan, Vol. 38, 351-360  \nHeymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140  \n- HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116 Hikage, C., Oguri, M., Hamana, T., et al. 2019, PASJ, 71, 43\n- Hildebrandt, H., van den Busch, J. L., Wright, A. H., et al. 2021, A&A, 647, A124\n- Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19\n- Hoekstra, H. 2001, A&A, 370, 743\n- Hoekstra, H. 2003, MNRAS, 339, 1155\n- Hoekstra, H., Mahdavi, A., Babul, A., & Bildfell, C. 2012, MNRAS, 427, 1298\n- Hogg, D. W. 1999, arXiv e-prints, astro\n- Hsieh, B. C. & Yee, H. K. C. 2014, ApJ, 792, 102\n- Hu, W. & Kravtsov, A. V. 2003, ApJ, 584, 702\n- Hu ff , E. & Mandelbaum, R. 2017, arXiv e-prints, arXiv:1702.02600\n- Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90  \nHuterer, D. 2023, A&A Rev., 31, 2  \n- Ider Chitham, J., Comparat, J., Finoguenov, A., et al. 2020, MNRAS, 499, 4768 Je ff reys, H. 1961, The Theory of Probability  \nKaiser, N. 1986, MNRAS, 222, 323  \nKleinebreil, F., Grandis, S., Schrabback, T., et al. 2024, arXiv e-prints, arXiv:2402.08456  \nKluge, M., Comparat, J., Liu, A., et al. 2024, arXiv e-prints, arXiv:2402.08453 Klypin, A. & Prada, F. 2019, MNRAS, 489, 1684 Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18 Kowalski, M., Rubin, D., Aldering, G., et al. 2008, ApJ, 686, 749 Kuijken, K., Heymans, C., Dvornik, A., et al. 2019, A&A, 625, A2 Kunz, M. & Sapone, D. 2006, Phys. Rev. D, 74, 123503 Lau, E. T., Kravtsov, A. V., & Nagai, D. 2009, ApJ, 705, 1129 Lesci, G. F., Marulli, F., Moscardini, L., et al. 2022, A&A, 659, A88 Lesgourgues, J., Mangano, G., Miele, G., & Pastor, S. 2013, Neutrino Cosmology  \nLewis, A. 2019, arXiv e-prints, arXiv:1910.13970  \nLewis, A. & Challinor, A. 2011, CAMB: Code for Anisotropies in the Microwave Background, Astrophysics Source Code Library, record ascl:1102.026 Li, X., Miyatake, H., Luo, W., et al. 2022, PASJ, 74, 421 Li, X., Zhang, T., Sugiyama, S., et al. 2023, Phys. Rev. D, 108, 123518 Liu, T., Merloni, A., Comparat, J., et al. 2022, A&A, 661, A27 Lovisari, L., Ettori, S., Gaspari, M., & Giles, P. A. 2021, Universe, 7, 139 Lovisari, L., Reiprich, T. H., & Schellenberger, G. 2015, A&A, 573, A118 Lovisari, L., Schellenberger, G., Sereno, M., et al. 2020, ApJ, 892, 102 Madhavacheril, M. S., Qu, F. J., Sherwin, B. D., et al. 2024, ApJ, 962, 113 Malmquist, K. G. 1922, Meddelanden fran Lunds Astronomiska Observatorium Serie I, 100, 1  \nMandelbaum, R. 2018, ARA&A, 56, 393  \nMandelbaum, R., Lanusse, F., Leauthaud, A., et al. 2018, MNRAS, 481, 3170 Mantz, A., Allen, S. W., Rapetti, D., & Ebeling, H. 2010, MNRAS, 406, 1759 Mantz, A. B., Allen, S. W., Morris, R. G., et al. 2016, MNRAS, 463, 3582 Mantz, A. B., von der Linden, A., Allen, S. W., et al. 2015, MNRAS, 446, 2205 Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113 Maughan, B. J., Jones, C., Jones, L. R., & Van Speybroeck, L. 2007, ApJ, 659, 1125  \nMcClintock, T., Varga, T. N., Gruen, D., et al. 2019, MNRAS, 482, 1352 Melchior, P., Gruen, D., McClintock, T., et al. 2017, MNRAS, 469, 4899 Meneghetti, M., Rasia, E., Merten, J., et al. 2010, A&A, 514, A93 Merloni, A., Lamer, G., Liu, T., et al. 2024, A&A, 682, A34 Miyatake, H., Sugiyama, S., Takada, M., et al. 2023, Phys. Rev. D, 108, 123517 More, S., Sugiyama, S., Miyatake, H., et al. 2023, Phys. Rev. D, 108, 123520 Mossa, V., Stöckel, K., Cavanna, F., et al. 2020, Nature, 587, 210 Murray, S. G., Power, C., & Robotham, A. S. G. 2013, Astronomy and Computing, 3, 23  \nMyles, J., Alarcon, A., Amon, A., et al. 2021, MNRAS, 505, 4249 Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624 Nelson, D., Springel, V., Pillepich, A., et al. 2019, Computational Astrophysics and Cosmology, 6, 2  \nNishizawa, A. J., Hsieh, B.-C., Tanaka, M., & Takata, T. 2020, arXiv e-prints, arXiv:2003.01511  \nNurgaliev, D., McDonald, M., Benson, B. A., et al. 2017, ApJ, 841, 5 Peebles, P. J. E. 1980, The large-scale structure of the universe Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565 Pierre, M., Pacaud, F., Adami, C., et al. 2016, A&A, 592, A1 Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648 Pitrou, C., Coc, A., Uzan, J.-P., & Vangioni, E. 2021, MNRAS, 502, 2474 Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A, 571, A20 Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2013, A&A, 550, A129 Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016a, A&A, 594, A13 Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016b, A&A, 594, A24 Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020a, A&A, 641, A6 Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020b, A&A, 641, A8 Pop, A.-R., Hernquist, L., Nagai, D., et al. 2022, arXiv e-prints,  \n- arXiv:2205.11528  \nPratt, G. W., Croston, J. H., Arnaud, M., & Böhringer, H. 2009, A&A, 498, 361 Predehl, P., Andritschke, R., Arefiev, V., et al. 2021, A&A, 647, A1 Press, W. H. & Schechter, P. 1974, ApJ, 187, 425  \nQian, X. & Vogel, P. 2015, Progress in Particle and Nuclear Physics, 83, 1  \nQu, F. J., Sherwin, B. D., Madhavacheril, M. S., et al. 2024, ApJ, 962, 112  \nRamos-Ceja, M. E., Oguri, M., Miyazaki, S., et al. 2022, A&A, 661, A14  \nRasia, E., Ettori, S., Moscardini, L., et al. 2006, MNRAS, 369, 2013  \nRasia, E., Meneghetti, M., Martino, R., et al. 2012, New Journal of Physics, 14,  \n055018  \nRasmussen, C. E. & Williams, C. K. I. 2006, Gaussian processes for machine learning (MIT Press)  \nReiprich, T. H. & Böhringer, H. 2002, ApJ, 567, 716  \nRiess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009 Riess, A. G., Yuan, W., Macri, L. M., et al. 2022, ApJ, 934, L7 Rossetti, M., Gastaldello, F., Eckert, D., et al. 2017, MNRAS, 468, 1917 Rossi, G., Yèche, C., Palanque-Delabrouille, N., & Lesgourgues, J. 2015, Phys. Rev. D, 92, 063505  \nArticle number, page 34 of 41  \n- Rowe, B. T. P., Jarvis, M., Mandelbaum, R., et al. 2015, Astronomy and Com-  \nputing, 10, 121 Ryko ff , E. S., Rozo, E., Busha, M. T., et al. 2014, ApJ, 785, 104 Salvato, M., Buchner, J., Budavári, T., et al. 2018, MNRAS, 473, 4937 Salvato, M., Wolf, J., Dwelly, T., et al. 2022, A&A, 661, A3 Sanders, J. S., Fabian, A. C., Russell, H. R., & Walker, S. A. 2018, MNRAS, 474, 1065 Saro, A., Bocquet, S., Rozo, E., et al. 2015, MNRAS, 454, 2305 Scheck, D., Sanders, J. S., Bi ffi , V., et al. 2023, A&A, 670, A33 Schellenberger, G. & Reiprich, T. H. 2017a, MNRAS, 469, 3738 Schellenberger, G. & Reiprich, T. H. 2017b, MNRAS, 471, 1370 Schrabback, T., Applegate, D., Dietrich, J. P., et al. 2018, MNRAS, 474, 2635 Scolnic, D., Brout, D., Carr, A., et al. 2022, ApJ, 938, 113 Scolnic, D. M., Jones, D. O., Rest, A., et al. 2018, ApJ, 859, 101 Secco, L. F., Samuro ff , S., Krause, E., et al. 2022, Phys. Rev. D, 105, 023515 Seppi, R., Comparat, J., Bulbul, E., et al. 2022, A&A, 665, A78 Seppi, R., Comparat, J., Nandra, K., et al. 2021, A&A, 652, A155 Sevilla-Noarbe, I., Bechtol, K., Carrasco Kind, M., et al. 2021, ApJS, 254, 24 Sheldon, E. S. & Hu ff , E. M. 2017, ApJ, 841, 24 Smoot, G. F., Bennett, C. L., Kogut, A., et al. 1992, ApJ, 396, L1 Sommer, M. W., Schrabback, T., Applegate, D. E., et al. 2022, MNRAS, 509, 1127 Sommer, M. W., Schrabback, T., Ragagnin, A., & Rockenfeller, R. 2023, arXiv e-prints, arXiv:2306.13187 Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676 Sugiyama, S., Miyatake, H., More, S., et al. 2023, Phys. Rev. D, 108, 123521 Sunyaev, R., Arefiev, V., Babyshkin, V., et al. 2021, A&A, 656, A132 Tanabashi, M., Hagiwara, K., Hikasa, K., et al. 2018, Phys. Rev. D, 98, 030001 Tanaka, M., Coupon, J., Hsieh, B.-C., et al. 2018, PASJ, 70, S9 Tanseri, I., Hagstotz, S., Vagnozzi, S., Giusarma, E., & Freese, K. 2022, Journal of High Energy Astrophysics, 36, 1 Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709 Tinker, J. L., Robertson, B. E., Kravtsov, A. V., et al. 2010, ApJ, 724, 878 Titsias, M. 2009, in Proceedings of Machine Learning Research, Vol. 5, Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics, ed. D. van Dyk & M. Welling (Hilton Clearwater Beach Resort, Clearwater Beach, Florida USA: PMLR), 567-574 Trotta, R. 2008, Contemporary Physics, 49, 71 Vagnozzi, S., Giusarma, E., Mena, O., et al. 2017, Phys. Rev. D, 96, 123503 Valageas, P., Clerc, N., Pacaud, F., & Pierre, M. 2011, A&A, 536, A95 van den Busch, J. L., Wright, A. H., Hildebrandt, H., et al. 2022, A&A, 664, A170 Varga, T. N., DeRose, J., Gruen, D., et al. 2019, MNRAS, 489, 2511 Vazza, F., Brunetti, G., Kritsuk, A., et al. 2009, A&A, 504, 33 Vikhlinin, A., Kravtsov, A. V., Burenin, R. A., et al. 2009, ApJ, 692, 1060 Vikman, A. 2005, Phys. Rev. D, 71, 023515 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261 von der Linden, A., Allen, M. T., Applegate, D. E., et al. 2014, MNRAS, 439, 2 Walker, T. P., Steigman, G., Schramm, D. N., Olive, K. A., & Kang, H.-S. 1991, ApJ, 376, 51 Willingale, R., Starling, R. L. C., Beardmore, A. P., Tanvir, N. R., & O'Brien, P. T. 2013, MNRAS, 431, 394 Willis, J. P., Oguri, M., Ramos-Ceja, M. E., et al. 2021, MNRAS, 503, 5624 Wright, A. H., Hildebrandt, H., van den Busch, J. L., & Heymans, C. 2020, A&A, 637, A100 Zeiler, M. D. 2012, arXiv e-prints, arXiv:1212.5701 Zeng, Z., Yeung, S., & Chu, M.-c. 2019, J. Cosmology Astropart. Phys., 2019, 015  \nZohren, H., Schrabback, T., Bocquet, S., et al. 2022, A&A, 668, A18 Zubeldia, Í. & Challinor, A. 2019, MNRAS, 489, 401", 'Appendix A: Mixture model': 'We start by employing our mixture model on the first term of Equation (6), by explicitly dividing the three populations that comprise our sample; clusters, RSs, and AGNs, as:  \ndNtot d ˆ CRd ˆ zd ˆ H i P ( I | ˆ CR , ˆ z , ˆ H i ) = dNC d ˆ CRd ˆ zd ˆ H i P ( I | ˆ CR , ˆ z , ˆ H i ) + NC (1 -f RS -f AGN) GLYPH<18> f RS P ( ˆ CR , ˆ z , ˆ H i | RS) + f AGN P ( ˆ CR , ˆ z , ˆ H i | AGN) GLYPH<19> (A.1)  \nwhere P ( ˆ CR , ˆ z , ˆ H i | RS) and P ( ˆ CR , ˆ z , ˆ H i | AGN) are the probability density function (PDF) for the contaminants that represent RSs and AGNs, respectively, while the first term, dNC d ˆ CRd ˆ zd ˆ H i P ( I | ˆ CR , ˆ z , ˆ H i ), represents the distribution of our selection observable for the detected clusters.  \nOn the second term of Equation (1), we use the relation between these three populations when computing the total number of detected objects as,  \nNtot = N C 1 -f RS -f AGN (A.2)  \nwhere  \nN C = & dN C d ˆ CRd ˆ zd ˆ λ d ˆ H i P ( I | ˆ CR , ˆ z , ˆ H i ) Θ ( ˆ λ > 3) d ˆ CRd ˆ zd ˆ λ d ˆ H i (A.3)  \nOn the third term in Equation (1), when using the optical information, we write  \nP ( ˆ λ | ˆ CR , ˆ z , ˆ H i , I ) = GLYPH<18> 1 -f RS P ( ˆ CR , ˆ z , ˆ H i | RS , I ) P ( ˆ CR , ˆ z , ˆ H i , I ) -f AGN P ( ˆ CR , ˆ z , ˆ H i | AGN , I ) P ( ˆ CR , ˆ z , ˆ H i , I ) GLYPH<19> × P ( ˆ λ | ˆ CR , ˆ z , ˆ H i , I , C ) + f RS P ( ˆ λ, ˆ z | RS , I ) P ( ˆ CR , ˆ H i | RS , I ) P ( ˆ CR , ˆ z , ˆ H i , I ) + f AGN P ( ˆ λ, ˆ z | AGN , I ) P ( ˆ CR , ˆ H i | AGN , I ) P ( ˆ CR , ˆ z , ˆ H i , I ) (A.4)  \nFinally, on the last term in Equation (1), when using the weak lensing information;  \nP (ˆ gt | ˆ CR , ˆ z , ˆ H i , I ) = GLYPH<18> 1 -f RS P ( ˆ CR , ˆ z , ˆ H i | RS , I ) P ( ˆ CR , ˆ z , ˆ H i , I ) -f AGN P ( ˆ CR , ˆ z , ˆ H i | AGN , I ) P ( ˆ CR , ˆ z , ˆ H i , I ) GLYPH<19> × P (ˆ gt | ˆ CR , ˆ z , ˆ H i , I , C ) + f RS P (ˆ gt | RS , I ) P (ˆ z | RS , I ) P ( ˆ CR , ˆ H i | RS , I ) P ( ˆ CR , ˆ z , ˆ H i , I ) + f AGN P (ˆ gt | ˆ z , AGN , I ) P (ˆ z | AGN , I ) P ( ˆ CR , ˆ H i | AGN , I ) P ( ˆ CR , ˆ z , ˆ H i , I ) . (A.5)  \nTo estimate the overall total contamination in the sample, we use the PDF of count rate and sky position information, that is: P ( ˆ CR , ˆ H i | AGN , I ), and P ( ˆ CR , ˆ H i | RS , I ) of the AGN and RS populations, obtained from our realistic eRASS1 simulations (Seppi et al. 2022; Comparat et al. 2020). We obtain the AGN and RS PDF in redshift and richness (Kluge et al. 2024), i.e., P ( ˆ λ, ˆ z | AGN , I ) from the eROMaPPer run on locations of the eRASS1 point sources (Kluge et al. 2024; Merloni et al. 2024) and P ( ˆ λ, ˆ z | RS , I ) from the eROMaPPer run at random lines of sight.  \nFinally, we use the predicted AGN shear profiles reported in the Comparat et al. (2023). In contrast, we assume that the shear profile is zero at all radii for the RS component. Therefore the shear signal is consistent with a noise fluctuation, as elaborated in Grandis et al. (2024), Sect. 5.1.1.', 'Appendix B.1: Relation between scatter and normalization': 'The reported normalization parameters for the X-ray count rate to mass scaling relation reported in Chiu et al. (2022, 2023) are in tension with our results. In particular, the di ff erence in normalization is due to the observed degeneracy between the normalization and the scatter. Chiu et al. (2022, 2023) assume an informative prior on the count rate mass scatter, σ X = 0 . 30 ± 0 . 08 derived by Grandis et al. (2019) from the luminosity and temperature scaling relations measured by Bulbul et al. (2019). This choice artificially breaks the amplitude scatter degeneracy. We show this in the following Fig. B.1 and the selection of the sample (Bahar et al. 2022). For a more detailed comparison of the scaling relation parameters, we refer to Kleinebreil et al. (2024) and Grandis et al. (2024).  \n3  \nFig. B.1: Degeneracy between scatter and amplitude of the count rate to mass scaling relation. We show our obtained value and the comparison with individual mass calibration from DES, KiDS, and HSC as reported in Kleinebreil et al. (2024); Grandis et al. (2024). The previous measurements from the eFEDS survey presented in Chiu et al. (2022) and Chiu et al. (2023) are shown with crosses.  \n<!-- image -->', 'Appendix B.2: Relation between cosmology and scatter parameters': 'In the following figures, Fig. B.2 and B.3, we show the posterior and correlation between the cosmology parameters in the fit and the X-ray scaling relation parameters.  \nFig. B.2: Degeneracies between cosmological parameters and X-ray count rate scaling relation parameters. It is interesting to point out that the scaling relation parameters a ff ect the degeneracy between Ω m and σ 8.  \n<!-- image -->  \nFig. B.3: Correlation between fitted parameters.  \n<!-- image -->', 'Appendix B.3: Ω m vs σ 8 in the several models fitted to our dataset': 'In this paper, we have fitted for several di ff erent cosmological models. All of these extensions have in common the starting point in the standard vanilla Λ CDM cosmological model. We have fitted the Λ CDM, and then we have allowed several parameters to be free. Particularly dark energy equation of state w and summed neutrino masses P m ν . Finally, for comparison purposes with Λ CDM fit in Planck analysis (Planck Collaboration et al. 2020a), we have also fitted the Λ CDM model with P m ν fixed to 0.06 eV. It is important to show that these extra free parameters do not significantly a ff ect the base robust parameters we are fitting and constraining in the Λ CDM: Ω m vs σ 8. We show this test in Fig. B.4.  \nFig. B.6: Robustness of our main Λ CDMresults (orange) against changes in the halo mass function: Tinker et al. (2008) (blue), and Despali et al. (2016) (green), and 5.5% systematics uncertainty added to the halo mass function of Tinker et al. (2008) following the prescription of (Costanzi et al. 2019b; Dark Energy Survey Collaboration et al. 2020)  \n<!-- image -->  \nm  \nFig. B.4: Posterior on Ω mσ 8 from the Λ CDM, ν CDM, w CDM, ν w CDM, and Λ CDM with P m ν = 0 . 06 eV fits from our eRASS1 cluster dataset. We see the nice overlap in the parameter space that indicates that our fit in the Ω mσ 8 plane is stable against modification to the standard Λ CDM cosmological model.', 'Appendix B.4: Tests of the effect on Ω m vs σ 8 from our modeling choices': 'In this paper, we have applied several cuts to the datasets, for instance, 0 . 1 < z < 0 . 8, EXT\\_LIKE > 6, we have chosen a specific halo mass function in the Tinker et al. (2008), and we have employed the mixture model. The following two figures show the e ff ect of our constraints on Ω m vs σ 8 in Λ CDM cosmological model. In Fig. B.5 in Appendix B, we show the impact of not applying the mixture model, cutting the data in redshift at 0.6 instead of 0.8, and changing the EXT\\_LIKE cut to a slightly more conservative value of 10. We see changes within the 1 σ level for all these cases. In Fig B.6, we show the e ff ect of changing the halo mass function from Tinker et al. (2008) to either Tinker et al. (2010) or Despali et al. (2016). We further tested the e ff ects of systematic uncertainties on the halo mass function by applying the formalism introduced in Costanzi et al. (2019b) and used in the literature (Dark Energy Survey Collaboration et al. 2020; Bocquet et al. 2024) to correct Tinker et al. (2008) model adding two nuisance parameters, q HMF and s HMF, which allow for extra degrees of freedom for both normalization and slope of the halo mass function. We used the following priors on these parameters: q HMF ∼ N (1 , 0 . 035) and s HMF ∼ N (0 , 0 . 02) as in the analysis of  \nBocquet et al. (2024). We can see that Ω m is unchanged, while σ 8 is lowered by around 1 σ . This indicates that change in mass function does not significantly impact our constraints; however, in future eROSITA cosmological analysis, we will need to pay much more care and attention to our modeling choice regarding the halo mass function.  \nFig. B.5: Comparison of our main Λ CDM results (orange) with several tests: switching o ff the mixture model likelihood of Sect. 4.2 (blue), restricting the redshift range z < 0 . 6.  \n<!-- image -->  \n<!-- image -->  \nm  \n.', 'Appendix C: Constraints from neutrino summed masses to the mass of the lightest': 'A constraint on the summed neutrino mass will shed light on the mass of the lightest neutrino eigenstate. Obviously, this also depends on the neutrino mass hierarchy as:  \nX m ν =            mL + q m 2 L + ∆ m 2 21 + q m 2 L + ∆ m 2 31 Normal Hierarchy mL + q m 2 L + ∆ m 2 32 + q m 2 L + ∆ m 2 31 Inverted Hierarchy (C.1)  \nwhere mL is the mass of the lightest neutrino eigenstate. We reproduce Fig. 3 in Qian & Vogel (2015) with updated parameters from NuFIT (Esteban et al. 2020) in Fig. C.1. In the same figure, we also show with horizontal dashed lines the constraints on summed neutrino masses for the combination of eROSITA with Planck CMB and lower allowed limit from neutrino oscillation experiments (Tanabashi et al. 2018).  \nFig. C.1: Current estimates of the relation between the mass of the lightest neutrino eigenstate and the summed neutrino masses (solid lines). The horizontal dashed lines and the shaded area represent the eRASS1 and CMB combined estimate. The vertical dash-dotted line represents the lightest neutrino value corresponding to our summed mass estimate. Finally, the vertical dotted line corresponds to the 1 σ uncertainty around the lightest neutrino mass estimate. We use blue and orange to distinguish between normal and inverted hierarchy. It is remarkable to show that although the di ff erence between the models, the estimation for the lightest neutrino mass is almost identical.  \n<!-- image -->', 'Appendix D: Optical and X-ray completeness': 'Fig. D.1: Optical (blue), X-ray (orange), and total (green) completeness computed as a function of mass from the best-fit values. The three panels show di ff erent redshift, 0.3, 0.5, and 0.15, and in each panel, a solid line corresponds to a exposure of 150 s, dashed to 500 s, and dotted to 50 s. The shaded area around each line represent the 1 σ uncertainty from resampling the best fitting chain 1000 times.  \n<!-- image -->  \nAs we have chosen to pursue a forward modeling approach, the completeness is not pre-determined, but it is modeled on the fly by propagating the selection function in count rate and the selection function in richness to halo mass via the observable mass relations and their correlated scatters. Following the approach proposed in Grandis et al. (2020) for modeling the total selection function of X-ray-selected, optically confirmed clusters, we computed the completeness curves as a function of mass and redshift for optical and X-ray selection functions. For the optical, we have:  \nP ( I | M , z ) = " Θ ( ˆ λ > 3) P ( ˆ λ | λ ) P ( λ | M , z ) d λ d ˆ λ, (D.1)  \nwhile for the X-ray, we have  \nP ( I | M , z , ˆ H i ) = " P ( I | CR , z , ˆ H i ) P ( CR | M , z ) dCR . (D.2)  \nWe computed both of these as a function of mass and displayed them in Figure D.1. In practice, the two incompleteness functions are allowed to correlate by the unknown correlation coefficient ρλ, C R , which we also sample on the fly. As appears from the figure, the optical cut ˆ λ > 3 only plays a minor role for low redshift, high exposure time regions, where it leads to a minimal extra incompleteness for the few 3 × 10 13 M ⊙ halos that up-scatter into our sample. Computing the total number of expected objects with and without completeness indicates that this e ff ect would reduce the total number of expected detections by 1.1%, if not modeled.', 'Appendix E: σ 8': 'σ 8 is defined as the amplitude of linear matter fluctuations in a spherical top hat of radius R = 8 h -1 Mpc at redshift z = 0, i.e.,  \nσ 2 8 = Z ∞ 0 ∆ 2 ( k , z = 0) 3 j 1( kR ) kR ! 2 d log k (E.1)  \nwhere j 1( kR ) is the spherical Bessel function of the first kind of order one, and ∆ 2 ( k , z = 0) is the power of linear matter fluctuation in logarithmic bins of wave number k at redshift z = 0,  \n∆ 2 ( k , z = 0) = k 3 P ( k , z = 0) 2 π 2 = As GLYPH<0> Ω m h 2 GLYPH<1> 2 4 25 k k piv ! ns -1 k 100 ! 4 g 2 ( z = 0) T 2 ( k ) (E.2)  \nwhere g ( z ) is the growth suppression factor (see Equation 10 in Huterer 2023) and T ( k ) is the linear transfer function. It is therefore clear that σ 8 is correlated with the combination of cosmological parameters √ As Ω m h 2 .', 'Appendix F: Systematics on scaling relation modeling': 'The measured intrinsic scatter in count rate is surprisingly large, σ X = 0 . 98 + 0 . 05 -0 . 04 , see in Table 2, compared with the value of ∼ 0.4 reported in WtG (Mantz et al. 2016), or the value of 0.33 reported in eFEDS (Chiu et al. 2022). Given that our cluster sample spans a much larger range in mass and redshift than previous samples, this may indicate mass- and / or redshift-dependence of the scatter. Here, we explore a more complex model of the scatter and its impact on scaling relation and cosmological parameters. We note that this test was performed after unblinding and that we consider the model presented in Sect. 6 as our main analysis.  \nThus, we model the scatter, allowing mass and redshift trends as  \nσ X ( M , z ) = σ X , 0 + σ X , z log 1 + z 1 + zp + σ X , M log M Mp (F.1)  \nfor the X-ray count rate, and  \nσλ ( M , z ) = σλ, 0 + σλ, z log 1 + z 1 + zp + σλ, M log M Mp (F.2)  \nfor the richness. Finally, we further unfreeze sM parameter, see Equation (30), enabling the scatter of the WL bias to depend on mass, as it was already redshift dependent.  \nWe further decided to allow a mass trend in the mass slope of the count rate to mass scaling relation. We thus modify the scaling relations as in Sect. 6, changing Equation (25) as follows:  \n* log CR CR , p GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> M , z + = log AX + bX ( M , z ) log M Mp + ex ( z ) , (F.3)  \nwhere ex ( z ) remains as of Equation (27), while the mass-redshift dependent slope of the scaling relation is expressed by the term bX ( M , z ) as:  \nbX ( z ) = BX + FX log 1 + z 1 + zp + CX log M Mp . (F.4)  \nHere, CX allows the mass-dependent slope.  \nIt is clear that this test introduces six extra free parameters, and therefore, assesses the level of the systematics of the model choices we have made throughout this paper. What we find is quite interesting. In fact, we discover a degeneracy between mixture model parameters, fAGN and fRS , and the scatter evolution and mass dependence parameters. This degeneracy can very well be explained physically. The population that occupies the lowcount-rate, low-richness can be comprised of both contaminants that our detection and confirmation pipeline is not able to distinguish from a galaxy cluster or low-mass groups of galaxies that upscatter in the scaling relation by, for instance, because of the presence of AGNs in its core (Bulbul et al. 2022). Our simulations (Comparat et al. 2020; Seppi et al. 2022) suggest that low-count-rate, low-richness entries in our catalog are due to AGN and RS contaminants. However, only a chance overlap between AGNs and clusters was introduced in eROSITA digital twin, while AGNs inside group and cluster-sized halos were not introduced. However, the situation in the real extended source cosmology catalog might be somewhat di ff erent.  \nAllowing the scatter to change as a function of mass and redshift e ff ectively enables the degeneracy between these two populations mentioned above, and indeed, the chains prefer a nonnegligible fraction of low-mass halos that upscatter in the scaling relations. Even though these changes do not a ff ect the value of S 8 = 0 . 866 ± 0 . 012, and the Ω m and σ 8 change by about 1 σ but still remain fully consistent with Planck CMB (Planck Collaboration et al. 2020a), we remain very cautious regarding the physical interpretation of the corresponding best fitting parameters in this model. In fact, despite the several improvements over the consistency with literature results on the value of the scatter at the pivot mass, now 0 . 66 ± 0 . 04, we have several indications that this model physically fails in determining some parameters:  \n- -The slope of the X-ray count rate scaling relation becomes significantly higher, 2 . 0 ± 0 . 04, much higher than previously reported in the literature for X-ray luminosity to mass scaling relation, see Sect. 6.1.\n- -The mass trend of the scatter is very strong, -0 . 37 ± 0 . 02, with suspiciously high precision, quite higher than what is expected from numerical simulations ( ∼ 0.1 Pop et al. 2022), and with significant implication for high-mass clusters ( > 10 15 M ⊙ ) having below 10% scatter in their luminosity to mass scaling relation, much lower than what was previously reported in dedicated studies of high-mass clusters, e.g. 40% in WtG (Mantz et al. 2016).\n- -Mixture model parameters fAGN and fRS are forced to 0 value. This would imply that our cleaning procedure has reduced the fraction of contaminants to nearly zero. However, we know from our simulations (Comparat et al. 2020; Seppi et al. 2022) and optical confirmation (Kluge et al. 2024) that a remaining contamination fraction of the order of 5% should be left in the cosmology sample.  \nIt is thus clear that despite the intriguing results of this run, further dedicated studies are needed to physically understand what is causing the observed parameters, currently not consistent with our best understanding of the eeROSITA sample. Further multi-wavelength studies of the eROSITA cluster population are clearly necessary; for instance, if these populations are indeed AGN, they should have characteristic infrared emission (Salvato et al. 2018, 2022). In our upcoming scaling relation paper (Pacaud et al. in preparation), we will investigate several aspects of the population modeling of galaxy clusters through scaling relation between mass and observables, with particular attention to the aspects discussed in this appendix. In our future eROSITA cosmological analysis, we will use these results to update our population modeling.'}