# The SRG-eROSITA All-Sky Survey: Constraints on $f(R)$ Gravity from Cluster Abundance

E Artis, V Ghirardini, E Bulbul, S Grandis, C Garrel, N Clerc, R Seppi, J

Comparat, M Cataneo, Y.E Bahar, et al.

► To cite this version:

E Artis, V Ghirardini, E Bulbul, S Grandis, C Garrel, et al. The SRG-eROSITA All-Sky Survey: Constraints on $f(R)$ Gravity from Cluster Abundance. Astronomy & Astrophysics - A&A, 2024, 691, pp.A301. 10.1051/0004-6361/202449587. hal-04483514

HAL Id: hal-04483514

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

Submitted on 26 Jun 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# The SRG/eROSITA All-Sky Survey

## Constraints on $f(R)$ gravity from cluster abundances

E. Artis<fnref n="1" />,*, D. V. Ghirardini<fnref n="1" />, D. E. Bulbul<fnref n="1" />, D. S. Grandis<fnref n="5" />,<fnref n="2" />, D. C. Garrel<fnref n="1" />, D. N. Clerc<fnref n="3" />, D. R. Seppi<fnref n="1" />,<fnref n="7" />, J. Comparat<fnref n="1" />, M. Cataneo<fnref n="4" />, Y. E. Bahar<fnref n="1" />, F. Balzer<fnref n="1" />, I. Chiu<fnref n="8" />, D. Gruen<fnref n="5" />, F. Kleinebreil<fnref n="2" />,<fnref n="4" />, M. Kluge<fnref n="1" />, S. Krippendorf<fnref n="5" />,<fnref n="6" />, X. Li<fnref n="10" />, A. Liu<fnref n="1" />, A. Merloni<fnref n="1" />, H. Miyatake<fnref n="12" />,<fnref n="13" />,<fnref n="14" />, S. Miyazaki<fnref n="11" />, K. Nandra<fnref n="1" />, N. Okabe<fnref n="9" />, F. Pacaud<fnref n="4" />, P. Predehl<fnref n="1" />, M. E. Ramos-Ceja<fnref n="1" />, T. H. Reiprich<fnref n="4" />, J. S. Sanders<fnref n="1" />, T. Schrabback<fnref n="2" />,<fnref n="4" />, S. Zelmer<fnref n="1" />, and X. Zhang<fnref n="1" />

Received 13 February 2024 / Accepted 30 September 2024

## ABSTRACT

The evolution of the cluster mass function traces the growth of the linear density perturbations and can be utilized to constrain the parameters of cosmological and alternative gravity models. In this context, we present new constraints on potential deviations from general relativity by investigating the Hu-Sawicki parametrization of the $f(R)$ gravity with the first Spectrum Roentgen Gamma (SRG)/eROSITA All-Sky Survey (eRASS1) cluster catalog in the western Galactic hemisphere in combination with the overlapping Dark Energy Survey Year-3, KiloDegree Survey, and Hyper Suprime-Cam data for weak lensing mass calibration. For the first time, we present constraints obtained from cluster abundances only. When we consider massless neutrinos, we find a strict upper limit of $\log|f_{R0}| < -4.31$ at a 95% confidence level. Massive neutrinos suppress structure growth at small scales, and thus have the opposite effect of $f(R)$ gravity. We consequently investigate the joint fit of the mass of the neutrinos with the modified gravity parameter. We obtain $\log|f_{R0}| < -4.08$ jointly with $\sum m_\nu < 0.49$ eV at a 95% confidence level, which is tighter than the limits in the literature utilizing cluster counts only. At $\log|f_{R0}| = -6$, the number of clusters is not significantly changed by the theory. Consequently, we do not find any statistical deviation from general relativity in the study of eRASS1 cluster abundance. Deeper surveys with eROSITA, increasing the number of detected clusters, will further improve constraints on $\log|f_{R0}|$ and investigate alternative gravity theories.

**Key words.** gravitation – galaxies: clusters: general – cosmological parameters – large-scale structure of Universe

## 1. Introduction

Our understanding of the Universe has been altered significantly with the discovery of its accelerated expansion, first demonstrated by supernovae observations (Riess et al. 1998; Perlmutter et al. 1999). The paradigm of the standard cosmological model widely represents this unknown dark energy with a constant added in Einstein's gravitational equations, $\Lambda$, the so-called cosmological constant. Other ingredients necessary to describe our cosmological model include a nonrelativistic (cold) dark matter (DM) component that represents most of the mass fraction of the large-scale structures (LSSs). These two components constitute the major fraction of the energy content of the present-day Universe and form the backbone of the concordance $\Lambda$CDM model.

So far, when considered together with inflation, this parametrization is relatively successful at describing the main observed properties of the Universe (Planck Collaboration VI 2020).

Galaxy clusters represent the most massive virialized structures in the Universe. The number density of clusters of galaxies per unit mass and redshift – the cluster mass function – and their spatial distribution – the two-point correlation function and higher-order statistics – have great potential to constrain the nature of dark energy and test the robustness of the underlying cosmological models (see Clerc & Finoguenov 2023, for a recent review). Cluster counts are currently the most efficient at placing constraints on the matter density parameter ($\Omega_m$), the root mean square of the density fluctuations in spheres of 8 Mpc $h^{-1}$ ($\sigma_8$), and the sum of the masses of the neutrinos ($\sum m_\nu$). Additionally, they also break the degeneracies between

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

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

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

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

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

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

<footnote n="7">⁷ Department of Astronomy, University of Geneva, Ch. d'Ecogia 16, 1290 Versoix, Switzerland</footnote>

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

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

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

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

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

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

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

<footnote>* Corresponding author; eartis@mpe.mpg.de</footnote>
the cosmological parameters of other probes, such as baryon acoustic oscillations (BAO; Eisenstein et al. 2005; Alam et al. 2017), redshift space distortions (RSD; Percival & White 2009), the weak-lensing (WL) cosmic shear (Troxel et al. 2018; Asgari et al. 2021; Amon et al. 2022; Li et al. 2023), type Ia supernovae and other distance ladders (Riess et al. 2019; Freedman et al. 2019), and the cosmic microwave background (CMB; Bennett et al. 2013; Planck Collaboration VI 2020). Previous and ongoing cluster surveys across a wide wavelength range demonstrate the potential of cluster number counts as a powerful cosmological probe when the biases related to mass calibration are reduced through the utilization of weak gravitational lensing observations (Mantz et al. 2015; Bocquet et al. 2019; Zubeldia & Challinor 2019; Ider Chitham et al. 2020; Garrel et al. 2022; Lesci et al. 2022; Sunayama et al. 2023; Fumagalli et al. 2024; Bocquet et al. 2024a,b, among others). Modern cluster surveys, such as the eROSITA All-Sky cluster survey (see Bulbul et al. 2024; Ghirardini et al. 2024), are revolutionizing the field of cosmology, promising to reach the precision of the CMB experiments on cosmological parameters, and have the potential to measure departures from Einstein's theory of general relativity (GR) and test theories of dark energy (Wolf & Ferreira 2023). In addition to cluster abundance, eROSITA can probe cosmology through the spatial distribution of galaxy clusters (Seppi et al. 2024) and the X-ray power spectrum.

It also provides an avenue for testing alternative gravity theories, which could, in return, explain the late-time accelerated expansion of the Universe. Indeed, a way of explaining the Universe's accelerated expansion is to modify the theory of gravitation on large scales, while retaining GR in the high-density Solar System and on galactic scales. Potential departures from GR can be modeled by a set of theories in which the Ricci scalar (R), the quantity describing the infinitesimal volumes in curved space-time, is replaced in the Einstein-Hilbert action by a generic function, f(R). One of the models introducing a physically motivated function is the one proposed by Hu & Sawicki (2007) (HS-f(R) hereafter). This model has direct phenomenological implications. For example, Tsujikawa et al. (2009) and Gannouji et al. (2009) showed that it would significantly increase the structure formation rate, and thus its predictions can be verified. In the context of cluster abundance studies, it increases the expected number density of massive DM halos. Many recent attempts have been made to provide tools for computing summary statistics measuring departures from GR, with the goal of constraining this model observationally. First, numerous codes adapt the standard Boltzmann solvers; for example, CAMB (Lewis & Challinor 2011) and CLASS (Lesgourgues 2011) for modified gravity (MG hereafter) studies with alternative software such as EFTCAMB (Hu et al. 2014), hi_class (Zumalacárregui et al. 2017), MGCAMB (Zucca et al. 2019; Wang et al. 2023), or FRCAMB Xu (2015). These extensions can compute the power spectrum in different MG scenarios.

However, the hierarchical structure formation scenario assumes that galaxy clusters originate from the gravitational collapse of matter overdensities. The cluster abundance can thus be obtained by computing the characteristics of this collapse in the corresponding theory (Kopp et al. 2013; Lombriser et al. 2013), while considering the standard power spectrum. Starting from this fact, several studies have attempted to constrain the f(R) parameters using clusters in the literature. In the first of the studies, Schmidt et al. (2009) used a sample of 49 low-redshift clusters detected by ROSAT and re-observed by Chandra (Vikhlinin et al. 2009) to constrain the HS-f(R) model. This approach was combined with other cosmological probes by

Lombriser et al. (2012), showing that the best constraints are obtained when cluster abundance is considered. Cataneo et al. (2015) report a robust upper limit for the value of the background scalar field, $log|f_{R0}| < -4.79$, by utilizing a sample of 224 ROSAT-detected clusters combined with CMB data. Additionally, Lombriser et al. (2013) and later Cataneo et al. (2016) developed a halo mass function (HMF) formalism to predict accurate cluster abundances, which was followed by the formalisms of von Braun-Bates et al. (2017) and Gupta et al. (2022). These studies explore the implication of f(R) gravity in the halo formation process. In parallel, Hagstotz et al. (2019) developed an alternative method whereby the HMF modeling accounts for the massive neutrinos and MG effects. The modifications of gravity are then encapsulated in the collapse model and the HMF parametrization. In this work, we have used this approach to predict the cluster counts. We emphasize again that all the cosmological quantities related to the matter power spectrum were computed with the standard version of CAMB, without any need for an application of the modified Boltzmann solver.

With 5259 clusters of galaxies selected in its cosmology sample in a wide redshift range, $0.1 < z < 0.8$, the western Galactic half of the first eROSITA All-Sky Survey (eRASS1) has the statistical power to reach far beyond constraints available from any cluster survey in the literature. Consequently, this work utilizes the largest intracluster medium (ICM) selected cluster sample to date, in combination with WL follow-up measurements for cluster masses and an accurate selection function to constrain deviations from GR by placing robust limits on the HS parametrization of f(R) gravity.

The paper is organized as follows. In Sect. 2, we describe the datasets used. Section 3 introduces the framework used to describe departures from GR in the HS-f(R) gravity. Section 4 provides a detailed description of the statistical methods. Finally, Sect. 5 presents the constraints on the HS-f(R) gravity and a comparison with the ΛCDM concordance model. Throughout this paper, we use the notation $log ≡ log_{10}$ and $ln ≡ log_e$.

## 2. Survey data

The data used in this work are identical to the survey data presented in Bulbul et al. (2024); Kluge et al. (2024); Ghirardini et al. (2024); Grandis et al. (2024). We summarize the multi-wavelength data utilized in this work, including X-ray, optical, and WL observations, in the following section.

### 2.1. eROSITA All-Sky Survey

The soft X-ray telescope eROSITA (Predehl et al. 2021) on board the Spectrum Roentgen Gamma Mission completed its All-Sky Survey program on June 11, 2020, 184 days after the start of the survey. In this work, we have used the data of the western Galactic half of the eROSITA All-Sky survey (359.9442 deg $> l > 179.9442$ deg) belonging to the German eROSITA consortium. The raw data were processed through the eSASS software (Brunner et al. 2022). Bulbul et al. (2024); Kluge et al. (2024) compiled two catalogs of clusters detected in this region: the primary galaxy group and cluster catalog and the cosmology sample based on the catalog of the X-ray sources in the soft band (0.2–2.3 keV) provided in Merloni et al. (2024). We utilized the eRASS1 cosmology sample in this work, which was compiled adopting a selection cut of the extent likelihood, $L_{ext} > 6$ (see Bulbul et al. 2024, for further details). The optical identification of the clusters selected in the LS DR10-South common
footprint produces a final cosmology sample of 5259 clusters of galaxies reaching purity levels of 95%, an ideal sample for the MG studies. Unlike the primary cluster catalog, the redshifts in the cosmology subsample are purely photometric measurements, enabling a consistent assessment of the systematic effects in our analysis. The total survey area covers a 12 791 deg² region in the western Galactic hemisphere. The count rates were extracted with the 2D-image fitting tool MBProj2D described in Bulbul et al. (2024) and were employed as a mass proxy in the WL mass calibration likelihood (Ghirardini et al. 2024).

## 2.2. Weak lensing survey data

To perform our mass calibration in a minimally biased way, we utilized the deep and wide-area optical surveys for lensing measurements with the overlapping footprint with eROSITA in the western Galactic hemisphere. The three wide-area surveys used in this work for mass calibration include the Dark Energy Survey (DES), Kilo Degree Survey (KiDS), and the Hyper Suprime-Cam Survey (HSC). We briefly describe the WL survey data here.

The three-year WL data (S19A) covered by the HSC Sub-aru Strategic Program (Aihara et al. 2018; Li et al. 2022) WL measurements were made around the location of eRASS1 clusters using in the g, r, i, z, and Y bands in the optical. The total area coverage of HSC is ≈500 deg². The shear profile, $g_t(\theta)$, the lensing covariance matrix that serves the measurement uncertainty, and the photometric redshift distribution of the selected source sample were obtained as the HSC WL data products for 96 eRASS1 clusters, with a total signal-to-noise of 40.

We utilized data from the first three years of observations of the Dark Energy Survey (DES Y3). The DES Y3 shape catalog (Gatti et al. 2021) was built from the r, i, z bands using the METACALIBRATION pipeline (Huff & Mandelbaum 2017; Sheldon & Huff 2017). Considering the overlap between the DES Y3 footprint and the eRASS1 footprint, we produced tangential shear data for 2201 eRASS1 galaxy clusters, with a total signal-to-noise of 65 in the tangential shear profile. The analysis and shear profile extraction details are presented in Grandis et al. (2024).

We used the gold sample of WL and photometric redshift measurements from the fourth data release of the Kilo-Degree Survey (Kuijken et al. 2019; Wright et al. 2020; Hildebrandt et al. 2021; Giblin et al. 2021), hereafter referred to as KiDS-1000. We extracted individual reduced tangential shear profiles for a total of 236 eRASS1 galaxy clusters in both the KiDS-North field (101 clusters) and the KiDS-South field (136 clusters), as both overlap with the eRASS1 footprint with a total signal-to-noise of 19.

The WL data in the mass calibration process is detailed and published in Ghirardini et al. (2024); Grandis et al. (2024). We did not reprocess the shear profiles and used the ones used in Ghirardini et al. (2024). This is justified in Sect. 3.3.

## 3. Modified gravity framework

The evolution of the Universe's LSSs is sensitive to modification to the standard frameworks of GR and $\Lambda$CDM. To investigate potential discrepancies, we need to probe different scales and epochs. Figure 1 shows the redshift and scale sensitivity of the most common cosmological probes. Galaxy cluster number counts efficiently place constraints on the low-redshift or large-scale regime not explored by other probes. This section

<figure><img src="image_2.png" /><figcaption>Fig. 1. Estimation of the scale and redshift sensitivity of the different cosmological probes, adapted from Preston et al. (2023). The background colors, from blue to green, represent the transition from the linear to the nonlinear regime. We represent cluster abundance with the red rectangle. The cosmology sample of eRASS1 spans the redshift range 0.1 < z < 0.8, as is represented by the solid red rectangle. Future data releases and experiments will increase the upper limit of this interval to the hatched red rectangle. Cluster abundance experiments thus probe a specific regime at LSSs and the nearby Universe, offering a complementary window on structure formation. The constraints from cluster clustering (Seppi et al. 2024) are omitted here. Additionally, the X-ray power spectrum might probe smaller scales.</figcaption></figure>

describes the framework explored in this paper; namely, the HS-f(R) model.

### 3.1. HS-f(R) gravity in the context of cluster counts

We chose to test the MG model that consists of a modification of the Einstein-Hilbert action in the following form:

$$S = \int dx^4 \sqrt{-g} \left( \frac{R + f(R)}{16\pi G} + \mathcal{L}_m \right), \qquad (1)$$

where $\mathcal{L}_m$ represents the Lagrangian of the matter field, R is the Ricci scalar, and $f(R)$ represents the extension to GR, for which $f(R) = -2\Lambda$. In the context of $f(R)$ gravity, many functional forms are proposed (see De Felice & Tsujikawa 2010, for a review). Throughout this paper, we focus on the Hu & Sawicki (2007) model, which formalizes $f(R)$ as

$$f(R) = -2\Lambda \frac{R^n}{R^n + m^{2n}}, \qquad (2)$$

where $m^2$ represents the curvature scale, and $\Lambda$ and $n$ (the scaling index) are constants. In the high-curvature regime ($R \gg m^2$), Eq. (2) becomes

$$f(R) = -2\Lambda \frac{R^n}{R^n + m^{2n}} = -2\Lambda \frac{1}{1 + \left(\frac{m^2}{R}\right)^n} \sim -2\Lambda \left(1 - \left(\frac{m^2}{R}\right)^n\right).$$

To express this function in terms of $f_R(R) = df/dR$, we rewrote it, and obtained

$$f(R) = -2\Lambda - \frac{f_{R0}}{n} \frac{\overline{R_0}^{n+1}}{R^n},$$
where we have defined the value of the field, $f_{R0} = df/dR(R_0) = -2n\Lambda m^{2n}/\bar{R}_0^{n+1}$, at $z=0$, $\bar{R}_0$ is the Ricci scalar at $z=0$, and overbars denote background quantities. State-of-the-art numerical simulations are required to build a precise HMF suitable for cluster abundance cosmology. In these simulations, it is common to use $n=1$ (e.g. Oyaizu 2008; Zhao et al. 2011; Giocoli et al. 2018). For consistency with the HMF calibration method, we thus froze the scaling index to unity. We finally obtained

$$f(R) \approx -2\Lambda - f_{R0} \frac{\bar{R}_0^2}{R}. \qquad (3)$$

In this model, we note that when $|f_{R0}| \to 0$, we recover the usual GR action with a cosmological constant in Eq. (1). Thus, we interpret $\Lambda$ as the cosmological constant. The only parameter quantifying the deviations from GR is $f_{R0}$. When $|f_{R0}| \ll 1$, the background expansion is not distinguishable from GR. In this paper, our main goal is to measure $f_{R0}$ or to set an upper limit on its value. In the rest of this section, we briefly describe the behavior of this model.

The variational principle applied to the action in Eq. (1) leads to the modified Einstein's equation:

$$G_{\mu\nu} - f_R R_{\mu\nu} - \left(\frac{f}{2} - \Box f_R\right) g_{\mu\nu} - \nabla_{\mu} \nabla_{\nu} f_R = 8\pi G T_{\mu\nu}. \qquad (4)$$

As usual, Eq. (4) provides the equation of motion of the scalar field, $f_R$:

$$\nabla^2 \delta f_R = \frac{a^2}{3} (\delta R(f_R) - 8\pi G \delta \rho_m), \qquad (5)$$

where $\delta$ denotes small perturbations of the related quantities. This also provides a Poisson-like equation given by

$$\nabla^2 \psi = \frac{16\pi G}{3} a^2 \rho_m - \frac{a^2}{6} \delta R(f_R). \qquad (6)$$

For a better understanding of the behavior in this model, we provide the thin-shell condition example first introduced by Khoury & Weltman (2004) and adapted in Li et al. (2012) and Lombriser et al. (2013). Since GR is well constrained in high-density environments like the Solar System, we want to recover it at small scales: this is the so-called “screening mechanism.” As an illustration, we shall consider a top-hat overdensity of radius $R_{in}$. Then, the gravitational pull exerted on a particle of mass, $m$, outside of the overdensity can be expressed as

$$F = G \frac{M_{in}m}{r^2} \left( 1 + \frac{1}{3} \frac{\Delta R}{R_{in}} \right), \qquad (7)$$

where $M_{in}$ is the mass inside the top-hat overdensity, $r$ is the distance to the center, and $\Delta R/R_{in}$ is the screening factor defined as

$$\frac{\Delta R}{R_{in}} = \min \left\{ \frac{3|f_R^{\text{in}} - f_R^{\text{out}}|}{2\phi_N}, 1 \right\}, \qquad (8)$$

where $\phi_N$ is the Newtonian gravitational potential, $\phi_N = GM_{in}/R_{in}$, and $f_R^{\text{in}}$ and $f_R^{\text{out}}$ are the values of the scalar field inside and outside of the top-hat overdensity. Considering that the background expansion is considered in $\Lambda$CDM, we can approximate the background Ricci scalar to the one expected in this model. Inside the top-hat overdensity, we consider that the screening

mechanism is active; that is, we recover GR. We can then express the scalar field as

$$f_R^{\text{in}} = f_{R0} \left( \frac{1 + 4\Omega_{\Lambda 0}/\Omega_{m 0}}{\rho_m^{\text{in}}/\rho_{m 0} + 4\Omega_{\Lambda 0}/\Omega_{m 0}} \right), \qquad (9)$$

where $\rho_m^{\text{in}}$ is the density inside the considered region. In Eq. (9), we note that the scalar field values become smaller in high-density regions ($\rho_{m0} \ll \rho_m^{\text{in}}$). Thus, in Eq. (8), the screening factor becomes

$$\frac{\Delta R}{R_{\text{in}}} \sim \frac{3|f_R^{\text{out}}|}{2\phi_N}.$$

The condition for the screening to be activated is $\Delta R/R_{in} \ll 1$. Thus, we must have $|f_R^{\text{out}}| \ll \phi_N$. Since $f_R$ is a monotonic function of time, a relevant condition would be to have

$$|f_{R0}| < \phi_N.$$

For massive galaxy clusters, the Newtonian potential (in natural units) is on the order of $10^{-5}$. Thus, for these gravitationally bound systems, the deviation to GR is screened (i.e., suppressed) when $|f_{R0}| < 10^{-5}$ (see Sect. 3.3). Overall, the HS-f(R) gravity model is defined with only one additional parameter: $f_{R0}$. This formalism allows the recovery of GR in the necessary regimes, identified from observational constraints. In the next section, we describe the formalism predicting the DM halo abundance.

### 3.2. The f(R) halo mass function

The density of clusters per unit mass and redshift was modeled as follows:

$$\frac{dn}{d \ln M} = \frac{\rho_{m,0}}{M} \frac{d \ln \sigma^{-1}}{d \ln M} f(\sigma), \qquad (10)$$

where $\rho_{m,0}$ is the matter density at present, $f(\sigma)$ is the multiplicity function, and $\sigma(R, z)$ is defined by

$$\sigma^2(R, z) = \frac{1}{2\pi^2} \int_0^\infty k^2 P(k, z) \left( \frac{3j_1(kR)}{kR} \right)^2 dk, \qquad (11)$$

where $j_1(x) = (\sin(x) - x \cos(x))/x^2$ is the spherical Bessel function of the first kind of order one, and $P(k)$ is the matter power spectrum. We followed the approach developed in Hagstotz et al. (2019) to describe the specific behavior of the HMF in the HS-f(R) framework when compared to the standard fitting functions and emulators (Tinker et al. 2008; Despali et al. 2016; Bocquet et al. 2020; Euclid Collaboration 2023) obtained from simulations assuming standard GR. The main advantage of this formalism is that the entire framework is derived from the collapse of the structures in $f(R)$ gravity, with standard computations for the power spectrum. Deviations from the GR HMF were modeled as follows: the multiplicity function can be written as in Eq. (10). The goal is to correct the multiplicity function by a factor encapsulating the deviations from GR; therefore,

$$f^{(R)}(\sigma) = f^{\text{GR}}(\sigma) \cdot \eta(f_{R0}, \sigma), \qquad (12)$$

where the correcting factor is given by

$$\eta(f_{R0}, \sigma) = \frac{f_k^{(R)}(f_{R0}, \sigma)}{f_k^{\text{GR}}(\sigma)}. \qquad (13)$$
The multiplicity functions, including the subscript *k* in the correcting factor, were defined as follows. In the formalism of the multiplicity function, the density field, $\delta$, realizes a random walk (i.e., a Markov process) when smoothed by a simple top-hat filter in *k* space, with the radius, *R*, defined in Eq. (11). This means that in the trajectories, $\delta(R)$, the next positions depend only on the current position. The multiplicity functions in Eq. (13) are obtained that way and are noted with the subscript *k*. They are not realistic enough to describe the halo abundance. The $f(R)$ properties are encoded on the collapse of structures and not on the statistics of the density field, and their ratio is well suited to quantify a departure from GR. However, more realistic filters should be considered to derive a proper fitting function, which causes departures from the uncorrelated random walk. Any multiplicity function, including these non-Markovian corrections found in the literature, can be applied for $f^{GR}$. In this work, we used $f^{GR}(\sigma)$ from Tinker et al. (2008). We emphasize that the resulting fitting function for HMF we have used in this work was initially calibrated for masses defined as $M_{200m}$. We converted the HMF to $M_{500c}$ as eRASS1 mass measurements are calibrated for $M_{500,c}$. Additionally, doing so allows us to stay consistent with the standard $\Lambda$CDM cosmology analysis presented in Ghirardini et al. (2024).

In the case of GR, the multiplicity function is given by

$$f_k^{GR}(\sigma) = \sqrt{\frac{2a}{\pi}} \frac{\delta_c}{\sigma} \exp\left(-a\frac{(\delta_c + \beta\sigma^2)^2}{2\sigma^2}\right), \quad (14)$$

where $a$ and $\beta$ were fixed to the values found in Hagstotz et al. (2019) and provided in Table 1, and $\delta_c$ is the critical density for the collapse in GR given in appendix C of Nakamura & Suto (1997),

$$\delta_c^{GR}(z) = \frac{3(12\pi)^{2/3}}{20} \left(1 - 0.012299 \log\left(1 + \frac{\Omega_m^{-1} - 1}{(1+z)^3}\right)\right). \quad (15)$$

The multiplicity function in the case of HS-$f(R)$ gravity for non-Markovian correction follows

$$f_k^{f(R)}(\sigma) = \sqrt{\frac{2a}{\pi}} \frac{1}{\sigma} e^{-a\bar{B}^2/(2\sigma^2)} \left( \delta_c^{f(R)} - \frac{3M}{2} \frac{\partial \delta_c^{f(R)}}{\partial M} \frac{\partial \ln \sigma}{\partial \ln R} \right), \quad (16)$$

where $a$ is a constant. $\bar{B}$ is the barrier, the threshold beyond which collapse structures are formed. This threshold follows

$$\bar{B} = \delta_c^{f(R)}(f_{R0}, M, z) + \beta\sigma^2,$$

with $\beta$ being another constant. The key of this formalism is thus to compute the spherical collapse threshold in $f(R)$ gravity, noted as $\delta_c^{f(R)}$ hereafter. We used the numerical solution provided by Kopp et al. (2013) and given in Eq. (17). For a review, we point the reader to the later reference and the previous ones in this section. The numerical solution for the spherical collapse threshold follows

$$\delta_c^{f(R)}(f_{R0}, M, z) = \delta_c^{GR}(z) \times \Delta(f_{R0}, M, z), \quad (17)$$

where $\Delta(f_{R0}, M, z)$ is given by

$$\Delta(f_{R0}, M, z) = 1 + b_2 (1+z)^{-a_3} \left(m_b - \sqrt{m_b^2 + 1}\right) + b_3 \left(\tanh(m_b) - 1\right). \quad (18)$$

<table><caption>**Table 1.** Parameters of the HMF used in this analysis.</caption><thead><tr><th>Parameter</th><th></th><th>Value</th></tr></thead><tbody><tr><td>a</td><td>Equations (14) and (16)</td><td>0.73</td></tr><tr><td>β</td><td>Equations (14) and (16)</td><td>0.11</td></tr><tr><td>μ<sub>1</sub></td><td>Equation (19)</td><td>1.38</td></tr><tr><td>μ<sub>2</sub></td><td>Equation (19)</td><td>21.32</td></tr><tr><td>α<sub>4</sub></td><td>Equation (19)</td><td>0.067</td></tr><tr><td>β<sub>3</sub></td><td>Equation (19)</td><td>0.0056</td></tr></tbody></table>

Notes. The parameters reported here are found in Hagstotz et al. (2019).

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Relative difference (in percentages) of the HMF at $z = 0.5$ when we increase $f_{R0}$. We note that $\Delta n = dn/d \ln M(f_{R0}) - dn/d \ln M$, where $dn/d \ln M(f_{R0})$ is the HMF described with Eq. (20) and $dn/d \ln M$ is obtained with the GR fitting function from Tinker et al. (2008). Increased values of $f_{R0}$ create more massive halos.</figcaption></figure>

The different parameters of the collapse model are defined as follows, as they were originally fit:

$$\left\{ \begin{array}{l} b_2 = 0.0166 \\ a_3(f_{R0}) = 1 + \exp(-2.08 (\log|f_{R0}| + 5.56)^2) \\ m_b(f_{R0}, M, z) = (1+z)^{a_3} (\log M - m_1(1+z)^{-a_4}) \\ m_1(f_{R0}) = \mu_1 \log|f_{R0}| + \mu_2 \\ a_4(f_{R0}) = \alpha_4 (\tanh(0.69 (\log|f_{R0}| + 6.65)) + 1) \\ b_3(f_{R0}) = \beta_3 (2.41 - \log|f_{R0}|) \end{array} \right. \quad (19)$$

Three values are then required to describe the collapse threshold in the HS-$f(R)$ parametrization. Those are represented by the parameters $\mu_1, \mu_2, \alpha_4$, and $\beta_3$. In this work, we have used the best-fit values given by Hagstotz et al. (2019), which are summarized in Table 1.

Figure 2 shows the expected distribution of DM halos when we vary $\log f_{R0}$. As was expected, we recover the GR prediction when $|f_{R0}|$ takes small values. At $\log f_{R0} < -5$, the number of clusters is changed at most by 10%, for objects of a mass of $10^{14} M_{\odot}$. For objects with $\log f_{R0} < -6$, the expected distribution clusters cannot be distinguished from the GR prediction. We emphasize that for $\log f_{R0} = -8$, we observe an excess of
very massive halos. However, this is related to modeling inaccu-
racies. Given the small counts of clusters in this mass range, this
has no impact on our cosmological inference pipeline.

The HMF derived here was computed for M<sub>200m</sub>, while the eRASS1 cosmology pipeline considers M<sub>500c</sub>. Thus, we followed (Hu & Kravtsov 2003) to convert the mass: we first computed the HMF from Tinker et al. (2008), dn/dln M<sub>500c</sub>, for a broad range of halo masses, M<sub>500c</sub>, and converted these masses assuming the mass concentration relation from Duffy et al. (2008) and the Navarro, Frenk, and White (NFW) density profile (Navarro et al. 1996). The resultant f(R) HMF thus follows

$$
\frac{dn^{f(R)}}{d \ln M_{500c}} = \frac{\rho_{m,0}}{M_{500c}} \left| \frac{d \ln \sigma}{d \ln M_{500c}} \right| f_{\text{Tinker2008}}(\sigma(M_{500c})) \\
\times \frac{f_{\text{k}}^{f(R)}(\sigma(M_{200m}(M_{500c})))}{f_{\text{k}}^{\text{GR}}(\sigma(M_{200m}(M_{500c})))}. \qquad (20)
$$

This equation represents the HMF used in this work, although
we note that the uncertainties related to the HMF can impact our
constraints (Salvati et al. 2020; Artis et al. 2021).

Indeed, due to different numerical simulations, halo finders,
various fitting functions used in previous work, and baryonic
effects, we expect the departure between the standard GR HMFs
in the literature to be as high as 10% (Knebe et al. 2013;
Euclid Collaboration 2023). Additionally, Cataneo et al. (2015)
derived constraints using a ratio of HMFs computed in GR and
f(R) gravity. Their approach was conservative, as an accurate,
robust, and computationally efficient modeling of the chameleon
screening was lacking at the time. Therefore, it was impossible
to take advantage of the full constraining power of the data. The
main difference is that the authors used the multiplicity func-
tion proposed by Sheth & Tormen (2002). In the ratio of their
multiplicity function, they also computed the variance, σ(R), in
the respective considered theory. In other words, the term related
to the f(R) multiplicity function uses the power spectrum com-
puted in f(R) gravity. Finally, the critical density for the collapse
in f(R) gravity is computed if the theory is always unscreened.
These modeling choices lead to differences of up to 20% with
the predictions considered in this work. As was described ear-
lier, the correction proposed by Hagstotz et al. (2019) is more
recent and takes into account more physical effects, while the
model of Cataneo et al. (2016) is computationally expensive in
the context of a Bayesian analysis. We thus chose to focus on the
latest HMF predictions. However, these systematic uncertainties
are accounted for following Costanzi et al. (2019); that is, we
computed the HMF as

$$
\frac{d\tilde{n}^{f(R)}}{d \ln M_{500c}} = \frac{dn^{f(R)}}{d \ln M_{500c}} (s \ln(M/10^{14}M_{\odot}) + q), \quad (21)
$$

where *q* and *s* represent the errors on the amplitude and slope
of the HMF, respectively. Adopting a conservative approach,
we used the priors *q* ~ *N*(1, 0.5) and *s* ~ *N*(0, 0.5). Annex A
investigates the choice of these priors and their impact on the
HMF.

3.3. Weak lensing mass calibration in modified gravity

In this work, we aim to constrain modifications of the theory of
gravity with the cluster HMF calibrated with the scaling relation
between the X-ray count rate and the WL shear data. Therefore,
the predictions for the WL signal of galaxy clusters in the MG
theories should be needed for self-consistent analysis, but we

show that GR approximation can be valid in this scenario for
the eRASS1 sample. As is discussed in the following, we used
the prediction from GR for the WL signal of massive clusters
by simply utilizing the calibration performed by Grandis et al.
(2024). For readability, we provide the basis of the analysis and
critical aspects (see Bartelmann & Schneider 2001, for detailed
description and review). In the framework of GR, the gravita-
tional potential of a lens cluster deflects the light from distant
background sources. The surface brightness, S_s, at position θ of
a source with negligible dimensions compared to the lens size is
as follows:

$$
S_{\text{obs}}(\theta) = S_s(A\theta), \qquad (22)
$$

where **A** is the lensing distortion matrix. It is written as

$$
\mathbf{A} = \begin{pmatrix} 1 - \kappa - \gamma_1 & \gamma_2 \\ \gamma_2 & 1 - \kappa + \gamma_1 \end{pmatrix}. \qquad (23)
$$

The convergence, $\kappa$, relates to isotropic transformations of the image, while the shear components ($\gamma_1, \gamma_2$) are linked to the image distortion. These two quantities describe the main transformations of the image, and they are computed as follows: the total distribution of mass in clusters is assumed to follow the NFW profile introduced in Navarro et al. (1996), which is well-fit simulated DM halos, and defined as

$$
\rho(r) = \frac{\delta(c_{\Delta})\rho_{c}(z)}{\frac{r}{r_{\Delta}/c_{\Delta}}\left(1 + \frac{r}{r_{\Delta}/c_{\Delta}}\right)^{2}}, \quad (24)
$$

where $r_Δ$ is the radius defined such that the mean density of the
enclosed volume is $Δ$ times the mean density of the Universe.
$\rho_c = 3H^2/8\pi G$ is the critical density and $c_Δ$ is the concentration
parameter. $\delta$ is a function that results from the definition of the
radius, $r_Δ$, and reads:

$$
\delta(c_{\Delta}) = \frac{\Delta}{3} \frac{c_{\Delta}^{3}}{\ln(1 + c_{\Delta}) - c_{\Delta}(1 + c_{\Delta})}. \qquad (25)
$$

This mass profile creates a gravitational potential noted as
$\phi$. We note the lens cluster, $l$, and background sources (typically
galaxies). In the thin lens approximation (the typical size of the
lens is negligible compared to the distance), we can write the
lensing potential of a cluster as

$$
\psi(\theta) = \frac{2}{c^2} \frac{D_{ls}}{D_1 D_s} \int \phi(D_1 \theta, z) \, dz, \qquad (26)
$$

where $\theta = (\theta_1, \theta_2)$ is the position on the sky. $D_{ls}$, $D_1$, and $D_s$ are, respectively, the distance between the lens cluster and the source galaxy, between the observer and the lens, and between the observer and the galaxy. The convergence and shear component can be expressed as a function of the derivatives of the lensing potential through

$$
\left\{
\begin{aligned}
\kappa(\theta) &= \frac{1}{2} \left( \frac{\partial^2 \psi}{\partial \theta_1^2} + \frac{\partial^2 \psi}{\partial \theta_1^2} \right) \\
\gamma_1(\theta) &= \frac{1}{2} \left( \frac{\partial^2 \psi}{\partial \theta_1^2} - \frac{\partial^2 \psi}{\partial \theta_1^2} \right) . \\
\gamma_2(\theta) &= \frac{\partial^2 \psi}{\partial \theta_1 \partial \theta_2}
\end{aligned}
\right.
\qquad (27)
$$
These expressions, combined with the definition of the lensing potential (26), lead to the following: if we define the surface mass density as

$$ \Sigma(R) = 2 \int_{0}^{\infty} \rho(R, z) dz, \qquad (28) $$

where the matter density follows the NFW profile (Eq. (24)), the convergence profile is expressed as

$$ \kappa(R) = \frac{\Sigma(R)}{\Sigma_c}, \qquad (29) $$

and the radial dependence of the tangential shear follows

$$ \gamma_t(R) = \frac{\Sigma(<R) - \Sigma(R)}{\Sigma_c}, \qquad (30) $$

where $\Sigma(<R)$ is the mean surface mass density in the radius, $R$, and $\Sigma_c = c^2/4\pi G \cdot D_s/D_1 D_{ls}$. These quantities were derived in Wright & Brainerd (2000). The WL mass calibration was performed through the “reduced tangential shear,”

$$ g_t(R) = \frac{\gamma(R)}{1 - \kappa(R)}. \qquad (31) $$

The reduced tangential shear profiles are related to the underlying mass through $P(\hat{g}_t|M_{WL}, \hat{z})$ (see Sect. 4.2). Since it is the consequence of the gravitational potential, it carries the necessary information on the underlying mass. To summarize, we assumed a standard lens equation resulting from GR to constrain our masses.

In this work, we consider one specific way of testing potential deviations: we probe the Hu-Sawicki parameterization of the $f(R)$ gravity.

The theory of gravity impacts cluster WL in two ways:

1. The paths of light through the cluster's gravitational potential might be altered by the modified theory, leading to the same lensing potential (Bekenstein & Sanders 1994) causing a different distortion in the shape of background galaxies.

2. The effective Newtonian potential of galaxy clusters in a modified theory might deviate from the one predicted in GR. The first point is of no concern in the modification we explore. In the case of $f(R)$ gravity, the cosmic expansion remains the same (see Sect. 3.1), but the geodesics are potentially affected. However, Schmidt (2010) showed that they remain virtually unchanged, up to a negligible factor on the order of $\le 10^{-4}$. We can thus use the generic lensing equations in HS-$f(R)$ gravity.

Halo matter profiles in $f(R)$ gravity have been extensively studied by Mitchell et al. (2018); Ruan et al. (2023), focusing on the mass concentration relation and the 3D density profiles. Neither of these two observables exactly matches the 2D projected matter density maps needed to anchor the WL mass calibration in a cluster cosmological context (Grandis et al. 2021a). Nonetheless, the general findings of Mitchell et al. (2018) about the maximal halo mass at which halos become fully screened can still be applied. That work finds that halos with a mass larger than $10^{p_2} h^{-1} M_{\odot}$ are screened, and thus behave like they would in GR. The parameter is given by

$$ p_2 = 1.50 (\log |f_{R0}| - \log(1 + z)) + 21.6, \qquad (32) $$

where we rounded before the numerical errors reported in Mitchell et al. (2018). We evaluated this expression at $z = 0.1$, the lower redshift limit of our cluster sample. Given that we

are working with an X-ray-selected sample, this is also the redshift at which we reach the lowest limiting mass. For the values $\log_{10} |f_{R0}| = 5.25, 5.5, 6$, we find $p_2 = 13.69, 13.31, 12.56$. The first value of $p_2$ corresponds to a mass of $10^{13.85} M_{\odot}$. Considering that 95.4% of the objects in the cosmology sample have an estimated mass greater than this value, we can argue that the cluster one halo regime that we analyze is also unaltered in the $f(R)$ gravity parameter space we explore.

## 4. Statistical inference

This section provides a brief description of the cluster cosmology inference pipeline. A detailed review can be found in the ΛCDM eRASS1 cosmology results (Ghirardini et al. 2024). The main element is that the cluster abundance statistic follows a Poisson statistic (i.e., the cluster abundance as a function of their observable is a Poisson process), coupled with WL shear measurements, calibrating the scaling relation between mass and X-ray observables. A mixture model for eliminating the contamination due to active galactic nuclei (AGNs) and background fluctuations misclassified as clusters in the eRASS1 cosmology subsample was employed to perform a precision cosmology experiment.

We express $\lambda$ as the intensity of the Poisson process (the number density of objects per unit of observables) and $x$ as the vector of the observables. This intensity depends on the model, and we note as $\Theta$ the set of cosmological parameters. By definition, the expected number of objects whose observable properties belong to the subsample, $\Omega$, of the observable parameter space follows

$$ N_{\{x \in \Omega\}}(\Theta) = \int_{\Omega} \lambda(x|\Theta) dx. \qquad (33) $$

In our model, the global observable vector is defined as

$$ x = \{ \tilde{C}_R, \tilde{z}, \tilde{\lambda}, \mathcal{H}, g_+ \}, $$

where $\tilde{C}_R$ are the observed count rate, $\tilde{z}$ are the observed photometric redshifts, $\tilde{\lambda}$ is the observed richness, $\mathcal{H}$ is the sky positions, which need to be considered due to the uneven exposure of the survey, and $g_+$ is the shear profile obtained from the different WL surveys described in Sect. 2.2. The intensity of the process thus describes the number density of the detected clusters in the observables parameter space. The observed value is the Poisson realization of the expected value, $N_{\{x \in \Omega\}}(\Theta)$. In the framework we adopt in this paper, we need to describe the intensity of the different species we are considering in order to predict the abundance of the objects. Then, the general shape of the Poisson log-likelihood will follow

$$ \ln \mathcal{L}(\Theta) = \sum_{i} \ln(\lambda(x_i|\Theta)) - \int_{x} \lambda(x|\Theta) dx, \qquad (34) $$

where the index, $i$, runs over the observed objects. Using the Poisson likelihood, we neglect the cosmic variance (Hu & Kravtsov 2003) and the super sample variance (Lacasa & Grain 2019) as these effects are subdominant given the area covered by eRASS1 and the limited number of clusters included in this analysis (see Ghirardini et al. 2024, for discussion). Our catalog contains three classes of objects: galaxy clusters (C), which are of interest, and contaminants, such as AGNs, and background fluctuations misclassified as clusters (NC). Our model simultaneously accounts for the cluster counts and the contaminant
fractions through the Poisson mixture model. The total density is the sum of the three-component model:

$$ \lambda_{\text{tot}}(x|\Theta) = \lambda_{\text{C}}(x|\Theta) + \lambda_{\text{AGN}}(x) + \lambda_{\text{NC}}(x). \quad (35) $$

In terms of the total number of objects in each class, we obtain

$$ N_{\text{tot}}(\Theta) = N_{\text{C}}(\Theta) + f_{\text{AGN}}N_{\text{tot}}(\Theta) + f_{\text{NC}}N_{\text{tot}}(\Theta), \quad (36) $$

where $f_{\text{AGN}}$ and $f_{\text{NC}}$ are the respective fractions of contaminants. Consequently, we can express the number of AGNs, false detections, and the total number of objects in the catalog as a function of the cosmology-dependent number of the cluster, $N_{\text{C}}$:

$$ \left\{ \begin{array}{l} N_{\text{tot}}(\Theta) = (1/(1 - f_{\text{AGN}} - f_{\text{NC}}))N_{\text{C}}(\Theta) \\ N_{\text{AGN}}(\Theta) = (f_{\text{AGN}}/(1 - f_{\text{AGN}} - f_{\text{C}}))N_{\text{C}}(\Theta) . \\ N_{\text{NC}}(\Theta) = (f_{\text{NC}}/(1 - f_{\text{AGN}} - f_{\text{NC}}))N_{\text{C}}(\Theta) \end{array} \right. \quad (37) $$

Starting from the formalism above, we can write the total number of objects as

$$ N_{\text{tot}}(\Theta) = \frac{1}{1 - f_{\text{AGN}} - f_{\text{NC}}} \int_{x} \lambda_{\text{C}}(x|\Theta) \, dx. \quad (38) $$

The number density of contaminants follows $\lambda_{\text{AGN}}(x|\Theta) = N_{\text{AGN}}(\Theta)P_{\text{AGN}}(x)$, and $\lambda_{\text{NC}}(x|\Theta) = N_{\text{NC}}(\Theta)P_{\text{NC}}(x)$. In this equation, $P$ describes the probability distribution function of the respective object, depending on the observables. These terms are described in Sect. 4.3. Finally, the likelihood (Eq. (34)) becomes

$$ \ln \mathcal{L}(\Theta) = \sum_i \ln \left( \lambda_C(x_i|\Theta) + N_{\text{AGN}}(\Theta)P_{\text{AGN}}(x_i) + N_{\text{NC}}(\Theta)P_{\text{NC}}(x_i) \right) - \frac{1}{1 - f_{\text{AGN}} - f_{\text{NC}}} \int_x \lambda_C(x|\Theta) \, dx. \quad (39) $$

The method described in the appendix of Bocquet et al. (2015) separates the mass calibration likelihood from the number count likelihood. For the clusters that belong to the footprints of the WL surveys, the number density becomes

$$ \lambda_C(x_{\backslash \{g_+\}}, g_+|\Theta) = \lambda_C(x_{\backslash \{g_+\}}|\Theta)P(g_+|x_{\backslash \{g_+\}}, \Theta). \quad (40) $$

This formalism is applied in Eq. (39) and allows for the calibration of X-ray and optical scaling relations (see Sect. 4.2 and Ghirardini et al. (2024) for details of the mass calibration).

In the following sections, we briefly describe the components of Eq. (39). The cluster abundance term, $\lambda_C$, is described in Sect. 4.1, the WL likelihood is given in Sect. 4.2, and the mixture model accounting for contamination is provided in Sect. 4.3.

### 4.1. Cluster abundance likelihood

The term $\lambda_C$ in Eq. (39) represents the expected number density of clusters. It is obtained through

$$ \lambda_C = \frac{dN_C}{d\hat{C}_R d\hat{z} d\hat{H}} = \int_{z} \int_{M} \int_{C_R} \mathcal{P}(I|C_R, z, \hat{H}) \mathcal{P}(\hat{C}_R|C_R) \mathcal{P}(C_R|M, z) \mathcal{P}(\hat{z}|z) \frac{dn}{dM} \frac{dV}{dz} dM dz dC_R. \quad (41) $$

$dn/dM$ is the HMF, $dV/dz$ is the comoving volume, and $\mathcal{P}(I|C_R, z, \hat{H})$ is the probability of detecting a cluster; that is,

the selection function, the fraction of detected objects at a given observed sky position, redshift, and count rates, shown in Clerc et al. (2024). We note that the fraction of detected objects is estimated for a ΛCDM cosmology. We assume that this quantity is not significantly modified by the models considered in this work. The other terms are explained below. Consistently with (Ghirardini et al. 2024), the count rates are assumed to be linked to the underlying mass of the objects through (Ghirardini et al. 2024)

$$ \left\langle \ln \frac{\overline{C_R}}{C_{R,p}} | M, z \right\rangle = \ln A_X + b_X(M, z) \cdot \ln \frac{M}{M_p} + e_X(z), \quad (42) $$

where $C_{R,p} = 0.1 \text{ cts/s}$ and $M_p = 2 \times 10^{14} M_\odot$ are the pivot count rates and mass. All the masses are expressed in solar masses, $M_\odot$. The other terms are

$$ b_X(M, z) = \left( B_X + C_X \cdot \ln \frac{M}{M_p} + F_X \cdot \ln \frac{1+z}{1+z_p} \right), \quad (43) $$

where $z_p = 0.35$ is the pivot redshift, and

$$ e_X(z) = D_X \cdot \ln \frac{d_L(z)}{d_L(z_p)} + E_X \cdot \ln \frac{E(z)}{E(z_p)} + G_X \cdot \ln \frac{1+z}{1+z_p}, \quad (44) $$

where $d_L$ is the luminosity distance and $E(z) = H(z)/H_0$ is the evolution function. The parameters {$A_X, B_X, C_X, D_X, E_X, F_X$} were fit jointly with the cosmological parameters to provide consistency and avoid astrophysical biases. We assume that the true count rates are the following,

$$ C_R \sim \text{Log-}\mathcal{N}(\overline{C_R}(M, z), \sigma_X), $$

where Log-N characterizes the log-normal distribution, $\overline{C_R}(M, z)$ is the output of the scaling relation in Eq. (42), and $\sigma_X$ is the dispersion on the selection function, fit together with the cosmological parameters. We also include the measurement errors on count rates supplied by the MultiBand Projector in 2D (MBProj2D) tool (Sanders et al. 2018), and the observed count rates are

$$ \hat{C}_R \sim \text{Log-}\mathcal{N}(C_R, \hat{\sigma}), $$

where $\hat{\sigma}$ is fixed at the value given by the X-ray processing pipeline, MBProj2D (see Bulbul et al. 2024, for details). The observed redshifts were obtained through eROMaPPer, the optical identification tool (see Kluge et al. 2024), as in the following:

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

where $c_z, b_z, \sigma_z$, and $c_{\text{shift},z}$ are the fit parameters of the distribution. The richness, a proxy for the number of galaxies belonging to a given cluster, was computed starting from Grandis et al. (2021b) as

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

where all the masses are expressed in solar mass ($M_\odot$), and the redshift dependence of the mass slope follows

$$ b_{\lambda}(z) = B_{\lambda} + D_{\lambda} \ln \left( \frac{1+z}{1+z_{p}} \right). \quad (47) $$

The optical richness is not one of the main observables used in the mass calibration. It is only used to eliminate the contamination in the sample through the mixture method. We kept the richness limit low ($\lambda > 3$) to ensure that no additional optical selection effects were introduced in the analysis.
<table><caption>**Table 2.** Priors on the cosmological parameters used in this analysis when they are involved in the corresponding model.</caption><thead><tr><th>Parameter</th><th>Units</th><th>Description</th><th>Prior</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>–</td><td>Mean matter density at present time</td><td>U(0.2, 0.95)</td></tr><tr><td>log A<sub>s</sub></td><td>–</td><td>Amplitude of the primordial power spectrum</td><td>U(-10, -8)</td></tr><tr><td>H<sub>0</sub></td><td>km/s/Mpc</td><td>Hubble expansion rate at present time</td><td>N(67.77, 0.6)</td></tr><tr><td>Ω<sub>b</sub></td><td>–</td><td>Mean baryon density at present time</td><td>U(0.046, 0.052)</td></tr><tr><td>n<sub>s</sub></td><td>–</td><td>Spectra index of the primordial power spectrum</td><td>U(0.92, 1.0)</td></tr><tr><td>Σ m<sub>v</sub></td><td>eV</td><td>Summed neutrino masses</td><td>U(0, 1)</td></tr><tr><td>log |f<sub>RO</sub>|</td><td>–</td><td>Logarithm of the derivative of f(R) with respect to the Ricci scalar taken at present-day</td><td>U(-8, -3)</td></tr></tbody></table>

**Notes.** For the full set of priors used for the other parameters, including scaling relation and nuisance parameters, see Ghirardini et al. (2024). With $U(\text{min}, \text{max})$, we indicate a uniform distribution between “min” and “max”. With $N(\mu, \sigma)$, we indicate a normal distribution centered on $\mu$ and with the standard deviation, $\sigma$.

## 4.2. Weak-lensing likelihood

An essential step between the cluster number counts and X-ray observables obtained from the eROSITA survey and the HMF is calibrating the cluster mass scaling relations. We utilized the WL shear measurements in calibrating cluster scaling relations as it is currently the most reliable method with minimal bias available (Euclid Collaboration 2024). We used the surveys presented in Sect. 2 – DES, KiDS, and HSC – to achieve this goal. Following Grandis et al. (2024), we assumed a relation between WL inferred masses and true mass in the following form:

$$ \left\langle \ln \frac{M_{\mathrm{WL}}}{M_{\mathrm{P}}} \middle| M, z \right\rangle = b(z) + b_{\mathrm{M}} \ln \left( \frac{M}{M_{\mathrm{P}}} \right). \qquad (48) $$

The DES, KiDS, and HSC surveys were used simultaneously to calibrate the count rate to mass (Eq. (42)) and richness to mass (Eq. (46)) relations. The scatter around this scaling is

$$ \log \sigma_{M_{\mathrm{WL}}}^2 = s(z) + s_M \log \left( \frac{M}{M_{\mathrm{P}}} \right), $$

where the mass component is consistent with zero, and we fixed $s_M = 0$ in the rest of the analysis, as is prescribed by Grandis et al. (2024).

## 4.3. Contamination fraction

Although the eRASS1 cosmology subsample is a high-fidelity catalog with a purity level reaching 95%, the remaining 5% of contaminants should be eliminated from our analysis. We applied a mixture model fitting for the fraction of AGN and random sources. We refer the reader to Ghirardini et al. (2024) for a detailed method review.

# 5. Results

This section provides our results for the HS-f(R) gravity based on cluster counts from the eRASS1 cosmology catalog with WL mass calibration. We first show that our mass scale is consistent with the base ΛCDM analysis from Ghirardini et al. (2024). The priors on the cosmological parameters are shown in Table 2.

## 5.1. Comparison of the mass scale with GR-ΛCDM

One of the key ingredients of cluster abundance cosmology is the relation between the observables used to detect the objects and the underlying cluster mass. In this work, the main observables are the X-ray count rates of the clusters observed in

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Comparison of the posteriors on the scaling relation obtained in the f(R) framework (the mass of the neutrinos is fixed) and the standard ΛCDM analysis presented in Ghirardini et al. (2024) (gray). The different models are fully consistent with one another. We overplot the contours obtained from WL only from Grandis et al. (2024) (green).</figcaption></figure>

the first eROSITA All-Sky Survey. The modeling of the scaling relations between observables and the cluster mass is described in Sect. 4. The parameters of this count rate to mass relation are fit jointly with the cosmological parameters. This means that when considering extensions to GR that do not affect cluster masses due to screening mechanisms, we should find consistent scaling relations, as the fit cosmology is similar to the one found in ΛCDM.

The redshift dependence of the mass slope, $C_X$, was fixed to zero, the luminosity distance dependence, $D_X$, to -2, and normalized Hubble parameter dependence, $E_X$, to 2. Figure 3 shows the fit parameters, $A_X$, $B_X$, $F_X$, and $G_X$, when compared to the standard ΛCDM analysis presented in Ghirardini et al. (2024) and the WL mass calibration part of the likelihood presented in Grandis et al. (2024). The latter represents the information from the WL mass calibration of DES clusters.

Overall, there is good agreement between our different models whenever we consider deviation from GR. The mass calibration parameters are indeed not affecting our results. This means that the mass calibration is robust and that changes in the cosmological model primarily affect the underlying distribution of clusters represented by the HMF.
<table><caption>**Table 3.** Constraints on the HS-f(R) gravity model in different scenarios.</caption><thead><tr><th></th><th>log |f<sub>R0</sub>|</th><th>Ω<sub>m</sub></th><th>σ<sub>8</sub></th><th>Σ m<sub>ν</sub></th><th>S<sub>8</sub></th></tr></thead><tbody><tr><td>HS-f(R)</td><td>&lt;-4.31</td><td>0.28 ± 0.02</td><td>0.89 ± 0.04</td><td>—</td><td>0.85 ± 0.04</td></tr><tr><td>HS-f(R) + Σ m<sub>ν</sub></td><td>&lt;-4.08</td><td>0.31 ± 0.03</td><td>0.89 ± 0.05</td><td>&lt;0.49 eV</td><td>0.84 ± 0.04</td></tr></tbody></table>

**Notes.** The errors provided for the constrained parameters are the 68% confidence intervals. The upper limits are given at the 95% confidence level.

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Comparison of the posteriors obtained on the $\Omega - \sigma_8$ plane. We represent the $\Lambda$CDM analysis introduced in Ghirardini et al. (2024) in gray, while the posteriors of HS-f(R) case are shown in red.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Constraints obtained on log $f_{R0}$ with eRASS1 clusters. No external cosmological datasets are combined with our data, which is shown in this figure. Constraints with a fixed neutrino mass are shown in red, while constraints for which the neutrino mass is free are shown in magenta.</figcaption></figure>

## 5.2. Constraints on f(R) gravity

We applied two different models to constrain the HS-f(R) gravity parameterization using the Hagstotz et al. (2019) formalism. In the first one, we assumed massless neutrinos, while in the second one, we simultaneously constrained the sum of the neutrino masses and the parameter log |f<sub>R0</sub>| (Eq. (3)). We summarize our results in Table 3. In both cases, eRASS1 cluster number counts provide upper limits on log |f<sub>R0</sub>|. In the first case, we find the following constraints, summarized in Figure 4:

$$ 
\begin{aligned}
\Omega_m &= 0.28 \pm 0.02 \\
\sigma_8 &= 0.89 \pm 0.04, \\
S_8 &= 0.85 \pm 0.04
\end{aligned} 
\qquad (49) 
$$

at a 68% confidence level. These constraints are in good agreement with the results of the standard $\Lambda$CDM analysis presented in Ghirardini et al. (2024). For the first time, we have obtained constraints on the f(R) parameters with cluster abundance only, obtaining

$$ \log |f_{R0}| < -4.31 \text{ at 95\% confidence } (\sum m_{\nu} = 0 \text{ eV}). \qquad (50) $$

Figure 5 presents the measurements in the $\Omega_m - \log |f_{R0}|$ plane. Our constraints show a consistency with GR, as has been suggested by previous surveys using cluster abundance to probe HS-f(R) gravity (see Sect. 5.3). Our results suggest that large unscreened departures from GR can be ruled out by cluster counts alone in the near future, providing an independent test of gravity on scales ~10 Mpc. Additionally, Hagstotz et al. (2019) show that massive neutrinos can counteract the effect of MG by reducing the abundance of massive structures. For instance,

massive neutrinos increase the collapse barrier presented in Eq. (17). For $\sum m_{\nu} = 0.3$ eV, the collapse barrier is not distinguishable from the GR barrier if $\log |f_{R0}| = -5$. It is thus necessary to investigate the case in which the sum of the mass of the neutrinos is fit jointly with $\log |f_{R0}|$. We used the prescription provided in Costanzi et al. (2013) and derived the linear quantities from the cold DM power spectrum instead of the total matter power spectrum. When we freed the masses of the neutrinos, we obtained

$$ \log |f_{R0}| < -4.08 \text{ at 95\% confidence } (\sum m_{\nu} \text{ free}). \qquad (51) $$

We additionally find the following constraints on the cosmological parameters, including the neutrino masses:

$$ 
\begin{aligned}
\Omega_m &= 0.31 \pm 0.03 \\
\sigma_8 &= 0.87 \pm 0.05 \\
S_8 &= 0.84 \pm 0.04 \\
\sum m_{\nu} &< 0.49 \text{ eV at 95\% confidence}
\end{aligned} 
\qquad (52) 
$$

These constraints are again fully consistent with those obtained in the base $\Lambda$CDM (Ghirardini et al. 2024) (see Table 3). We investigated the consistency between our $S_8$ measurements and the one given in the standard $\Lambda$CDM analysis in Ghirardini et al. (2024), who finds

$$ S_8 = 0.86 \pm 0.01. \qquad (53) $$

Figure 6 shows the different models, compared to the results obtained by Planck Collaboration VI (2020). In all cases, our $S_8$
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Comparison of the $S_8$ parameter obtained using the different models. We show the $1\sigma$ uncertainties. The case introducing massless neutrinos is in red, while the one considering massive neutrinos is in pink. Results from the CMB are shown in green. The gray area represents the 68% and 95% confidence regions of the ΛCDM analysis from Ghirardini et al. (2024).</figcaption></figure>

value is in good agreement with the CMB. Additionally, the three models are fully compatible. However, the upper limits on the neutrino mass are higher due to the increased free parameters in the fits. It is worth noting that we solely used cluster abundances with no priors on the cosmological parameters when constraining $f_{R0}$. The ability to find upper limits with cluster counts alone can primarily be explained by the increased statistics provided by the eRASS1 cluster sample. This allows us to provide a self-consistent physical framework. Moreover, the constraints benefit from the shear measurement from high-quality DES, KiDS, and HSC mass calibration employed in this work (Grandis et al. 2024). It is therefore clear that in both cases, whether or not the neutrino masses are allowed to be free, the constraints from the eRASS1 cluster count measurements remain consistent with GR with a cosmological constant, for which $f_{R0} = 0$, as is shown in Eq. (3).

### 5.3. Comparison with other constraints

Hereafter, all the constraints are shown at the 95% confidence level. In this work, we provide constraints on $f(R)$ gravity with cluster abundance only for the first time. This allows us to test gravity with a self-consistent physical framework at a given physical scale and redshift range. Indeed, recent forecasts (Mitchell et al. 2021; Vogt et al. 2024) predict that future cluster abundance studies alone will be able to distinguish GR from alternative scenarios. Similarly, an alternative way of probing cosmology with the peaks of the matter density field is to consider WL peaks. Those also show promising results (Davies et al. 2024). This section compares our results with the ones obtained by other surveys and studies, from the smallest scales to the cosmological constraints. Indeed, due to the unknown nature of the potential modification of GR, it is important to test the theory at all physical scales.

Due to the nature of chameleon screening, the most stringent constraints are obtained from dwarf galaxies in cosmic voids. Desmond & Ferreira (2020) found $|f_{R0}| < -7.85$ using galaxy morphology, which practically rules out the Hu-Sawicki model of $f(R)$ gravity at astrophysical scales. Despite these results, testing the theory at larger scales is still relevant.

In the case of galaxy clusters, previous studies obtained consistent results with those presented in this work, but combining

cluster abundances with the CMB and other LSS probes (Schmidt et al. 2009; Lombriser et al. 2012). When combined with CMB data, the most recent study by Cataneo et al. (2015) obtained $\log|f_{R0}| < -4.79$. We reach a similar precision with clusters only, primarily due to the statistical power of the eRASS1 cosmology sample, the largest ICM-selected sample to date.

Other probes at the cosmological scale typically obtain weaker constraints than the ones probing the small scales. They include the WL 2pt statistic. Tröster et al. (2021), using the KIDS 3 × 2 pt measurements of cosmic shear from KiDS-1000, did not find significant constraints on $f(R)$ gravity using their data. Conversely, Liu et al. (2016) obtained $\log(f_{R0}) < -5.15$ combining measurements of the WL peaks with CMB priors. Hu et al. (2016) used supernovae data from JLA, WiggleZ, and CFHTLenS, in combination with CMB and BAO data, to find $\log(f_{R0}) < -3.2$. More recently, Kou et al. (2024) found an upper limit of $\log|f_{R0}| < -4.61$, cross-correlating CMB lensing with galaxy clustering.

Consequently, our results are consistent with the latest cosmological constraints. Astronomical probes on astrophysical and cosmological scales have found no evidence of a deviation from GR. Recent forecasts established for the Euclid survey (Casas et al. 2023) show that next-generation surveys will have the ability to distinguish GR from $f(R)$ gravity with extreme precision, potentially ruling out the model definitively.

## 6. Discussion and conclusions

In this work, we present results on potential deviations from standard GR and ΛCDM using X-ray detected eRASS1 clusters, in combination with LS DR10-South for optical confirmation and redshift measurements, and DES, KiDS, and HSC for WL mass calibration. We have derived constraints on the HS model of $f(R)$ gravity. We have considered both the case of massive and massless neutrinos.

All the models we have studied are statistically consistent with GR. Given the constraining power of eRASS1, we have chosen to present constraints from the cluster abundance only. This choice allows us to measure our parameters with a consistent physical framework for the first time. Indeed, this work shows that cluster counts alone have a great potential to constrain MG at the scale at which they are sensitive. Deeper surveys with eROSITA will significantly increase the constraints that we obtain here and might rival the results obtained at smaller scales. Overall, combined with the follow-up WL observations, the first eROSITA’s All-Sky survey has the statistical power to constrain the cosmological parameters with percent-level precision and to test the GR at large scales.

Lastly, we have investigated the consistency between our $S_8$ measurements and the one given in the standard ΛCDM analysis and found good agreement. Modifying gravitation at the scale covered by cluster abundance does not significantly affect the $S_8$ measurement. Our error bars are increased, but all posteriors are fully compatible. This emphasizes that the $S_8$ measurement provided in the ΛCDM analysis is robust against modifications and in good agreement with Planck Collaboration VI (2020).

Finally, the precision of our analysis mostly relies on the available HMF models. We have used the framework proposed by Hagstotz et al. (2019), while others like Cataneo et al. (2016) propose alternative approaches. Additional work will be necessary to increase the precision of our constraints by improving the reliability of the HMF, not only for $f(R)$ gravity but also for any additional model of interest.
eROSITA, thus far, has collected data for more than four all-sky surveys. The final eRASS will detect about 10⁵ galaxy clusters and groups (projected from the eFEDS results in Liu et al. 2022; Bulbul et al. 2022). The deeper eROSITA data and extensive cluster catalogs, in combination with state-of-the-art WL mass calibration, will allow us to tighten our constraints on the *f(R)* gravity models.

*Acknowledgement.* The full acknowledgments are available in Appendix B.

## References

Aihara, H., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S8

Alam, S., Ata, M., Bailey, S., et al. 2017, MNRAS, 470, 2617

Amon, A., Gruen, D., Troxel, M. A., et al. 2022, Phys. Rev. D, 105, 023514

Artis, E., Melin, J.-B., Bartlett, J. G., & Murray, C. 2021, A&A, 649, A47

Asgari, M., Lin, C.-A., Joachimi, B., et al. 2021, A&A, 645, A104

Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167

Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291

Bayer, J., Osendorfer, C., Diot-Girard, S., Rueckstiess, T., & Urban, S. 2015, TUM, Tech. Rep.

Bekenstein, J. D., & Sanders, R. H. 1994, ApJ, 429, 480

Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20

Bocquet, S., Saro, A., Mohr, J. J., et al. 2015, ApJ, 799, 214

Bocquet, S., Dietrich, J. P., Schrabback, T., et al. 2019, ApJ, 878, 55

Bocquet, S., Heitmann, K., Habib, S., et al. 2020, ApJ, 901, 5

Bocquet, S., Grandis, S., Bleem, L. E., et al. 2024a, Phys. Rev. D, submitted [arXiv:2310.12213v2]

Bocquet, S., Grandis, S., Bleem, L. E., et al. 2024b, Phys. Rev. D, submitted [arXiv:2401.02075]

Brunner, H., Liu, T., Lamer, G., et al. 2022, A&A, 661, A1

Buchner, J. 2021, The Journal of Open Source Software, 6, 3001

Bulbul, E., Liu, A., Pasini, T., et al. 2022, A&A, 661, A10

Bulbul, E., Liu, A., Kluge, M., et al. 2024, A&A, 685, A106

Casas, S., Cardone, V. F., Sapone, D., et al. 2023, A&A, submitted [arXiv:2306.11053]

Cataneo, M., Rapetti, D., Schmidt, F., et al. 2015, Phys. Rev. D, 92, 044009

Cataneo, M., Rapetti, D., Lombriser, L., & Li, B. 2016, J. Cosmology Astropart. Phys., 2016, 024

Chisari, N. E., Alonso, D., Krause, E., et al. 2019, CCL: Core Cosmology Library, Astrophysics Source Code Library [record ascl:1901.0003]

Clerc, N., & Finoguenov, A. 2023, in Handbook of X-ray and Gamma-ray Astrophysics, eds. C. Bambi, & A. Santangelo, 123

Clerc, N., Comparat, J., Seppi, R., et al. 2024, A&A, 687, A238

Costanzi, M., Villaescusa-Navarro, F., Viel, M., et al. 2013, J. Cosmology Astropart. Phys., 2013, 012

Costanzi, M., Rozo, E., Simet, M., et al. 2019, MNRAS, 488, 4779

Davies, C. T., Harnois-Déraps, J., Li, B., et al. 2024, MNRAS, 533, 3546

De Felice, A., & Tsujikawa, S. 2010, Liv. Rev. Relativ., 13, 3

Desmond, H., & Ferreira, P. G. 2020, Phys. Rev. D, 102, 104060

Despali, G., Giocoli, C., Angulo, R. E., et al. 2016, MNRAS, 456, 2486

Duffy, A. R., Schaye, J., Kay, S. T., & Dalla Vecchia, C. 2008, MNRAS, 390, L64

Eisenstein, D. J., Zehavi, I., Hogg, D. W., et al. 2005, ApJ, 633, 560

Euclid Collaboration (Castro, T., et al.) 2023, A&A, 671, A100

Euclid Collaboration (Giocoli, C., et al.) 2024, A&A, 681, A67

Freedman, W. L., Madore, B. F., Hatt, D., et al. 2019, ApJ, 882, 34

Fumagalli, A., Costanzi, M., Saro, A., Castro, T., & Borgani, S. 2024, A&A, 682, A148

Gannouji, R., Moraes, B., & Polarski, D. 2009, J. Cosmology Astropart. Phys., 2009, 034

Garrel, C., Pierre, M., Valageas, P., et al. 2022, A&A, 663, A3

Gatti, M., Sheldon, E., Amon, A., et al. 2021, MNRAS, 504, 4312

Ghirardini, V., Bulbul, E., Artis, E., et al. 2024, A&A, 689, A298

Giblin, B., Heymans, C., Asgari, M., et al. 2021, A&A, 645, A105

Giocoli, C., Baldi, M., & Moscardini, L. 2018, MNRAS, 481, 2813

GPy 2012, GPy: A Gaussian process framework in python, http://github.com/SheffieldML/GPy

Grandis, S., Bocquet, S., Mohr, J. J., Klein, M., & Dolag, K. 2021a, MNRAS, 507, 5671

Grandis, S., Mohr, J. J., Costanzi, M., et al. 2021b, MNRAS, 504, 1253

Grandis, S., Ghirardini, V., Bocquet, S., et al. 2024, A&A, 687, A178

Gupta, S., Hellwing, W. A., Bilicki, M., & García-Farieta, J. E. 2022, Phys. Rev. D, 105, 043538

Hagstotz, S., Costanzi, M., Baldi, M., & Weller, J. 2019, MNRAS, 486, 3927

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Hildebrandt, H., van den Busch, J. L., Wright, A. H., et al. 2021, A&A, 647, A124

Hu, W., & Kravtsov, A. V. 2003, ApJ, 584, 702

Hu, W., & Sawicki, I. 2007, Phys. Rev. D, 76, 064004

Hu, B., Raveri, M., Frusciante, N., & Silvestri, A. 2014, Phys. Rev. D, 89, 103530

Hu, B., Raveri, M., Rizzato, M., & Silvestri, A. 2016, MNRAS, 459, 3880

Huff, E., & Mandelbaum, R. 2017, arXiv e-prints [arXiv:1702.02600]

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Ider Chitham, J., Comparat, J., Finoguenov, A., et al. 2020, MNRAS, 499, 4768

Khoury, J., & Weltman, A. 2004, Phys. Rev. D, 69, 044026

Kluge, M., Comparat, J., Liu, A., et al. 2024, A&A, 688, A210

Knebe, A., Pearce, F. R., Lux, H., et al. 2013, MNRAS, 435, 1618

Kopp, M., Appleby, S. A., Achitouv, I., & Weller, J. 2013, Phys. Rev. D, 88, 084015

Kou, R., Murray, C., & Bartlett, J. G. 2024, A&A, 686, A193

Kuijken, K., Heymans, C., Dvornik, A., et al. 2019, A&A, 625, A2

Lacasa, F., & Grain, J. 2019, A&A, 624, A61

Lesci, G. F., Nanni, L., Marulli, F., et al. 2022, A&A, 665, A100

Lesgourgues, J. 2011, arXiv e-prints [arXiv:1104.2932]

Lewis, A., & Challinor, A. 2011, CAMB: Code for Anisotropies in the Microwave Background, Astrophysics Source Code Library [record ascl:1102.026]

Li, B., Zhao, G.-B., & Koyama, K. 2012, MNRAS, 421, 3481

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, X., Li, B., Zhao, G.-B., et al. 2016, Phys. Rev. Lett., 117, 051101

Liu, A., Bulbul, E., Ghirardini, V., et al. 2022, A&A, 661, A2

Lombriser, L., Slosar, A., Seljak, U., & Hu, W. 2012, Phys. Rev. D, 85, 124038

Lombriser, L., Li, B., Koyama, K., & Zhao, G.-B. 2013, Phys. Rev. D, 87, 123511

Mantz, A. B., von der Linden, A., Allen, S. W., et al. 2015, MNRAS, 446, 2205

Merloni, A., Lamer, G., Liu, T., et al. 2024, A&A, 682, A34

Mitchell, M. A., He, J.-h., Arnold, C., & Li, B. 2018, MNRAS, 477, 1133

Mitchell, M. A., Arnold, C., & Li, B. 2021, MNRAS, 508, 4157

Nakamura, T. T., & Suto, Y. 1997, Progr. Theor. Phys., 97, 49

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Oyaizu, H. 2008, Phys. Rev. D, 78, 123523

Percival, W. J., & White, M. 2009, MNRAS, 393, 297

Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565

Planck Collaboration VI. 2020, A&A, 641, A6

Predehl, P., Andritschke, R., Arefiev, V., et al. 2021, A&A, 647, A1

Preston, C., Amon, A., & Efstathiou, G. 2023, MNRAS, 525, 5554

Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009

Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., & Scolnic, D. 2019, ApJ, 876, 85

Ruan, C.-Z., Cuesta-Lazaro, C., Eggemeier, A., et al. 2023, MNRAS, 527, 2490

Salvati, L., Douspis, M., & Aghanim, N. 2020, A&A, 643, A20

Sanders, J. S., Fabian, A. C., Russell, H. R., & Walker, S. A. 2018, MNRAS, 474, 1065

Schmidt, F. 2010, Phys. Rev. D, 81, 103002

Schmidt, F., Vikhlinin, A., & Hu, W. 2009, Phys. Rev. D, 80, 083505

Seppi, R., Comparat, J., Ghirardini, V., et al. 2024, A&A, 686, A196

Sheldon, E. S., & Huff, E. M. 2017, ApJ, 841, 24

Sheth, R. K., & Tormen, G. 2002, MNRAS, 329, 61

Sunayama, T., Miyatake, H., Sugiyama, S., et al. 2023, arXiv e-prints [arXiv:2309.13025]

Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709

Tröster, T, Asgari, M., Blake, C., et al. 2021, A&A, 649, A88

Troxel, M. A., MacCrann, N., Zuntz, J., et al. 2018, Phys. Rev. D, 98, 043528

Tsujikawa, S., Gannouji, R., Moraes, B., & Polarski, D. 2009, Phys. Rev. D, 80, 084044

Vikhlinin, A., Burenin, R. A., Ebeling, H., et al. 2009, ApJ, 692, 1033

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261

Vogt, S. M. L., Bocquet, S., Davies, C. T., Mohr, J. J., & Schmidt, F. 2024, Phys. Rev. D, 109, 123503

von Braun-Bates, F., Winther, H. A., Alonso, D., & Devriendt, J. 2017, J. Cosmology Astropart. Phys., 2017, 012

Wang, Z., Mirpoorian, S. H., Pogosian, L., Silvestri, A., & Zhao, G.-B. 2023, J. Cosmology Astropart. Phys., 2023, 038

Wolf, W. J., & Ferreira, P. G. 2023, Phys. Rev. D, 108, 103519

Wright, C. O., & Brainerd, T. G. 2000, ApJ, 534, 34

Wright, A. H., Hildebrandt, H., van den Busch, J. L., & Heymans, C. 2020, A&A, 637, A100

Xu, L. 2015, arXiv e-prints [arXiv:1506.03232]

Zhao, G.-B., Li, B., & Koyama, K. 2011, Phys. Rev. D, 83, 044007

Zubeldia, I., & Challinor, A. 2019, MNRAS, 489, 401

Zucca, A., Pogosian, L., Silvestri, A., & Zhao, G. B. 2019, J. Cosmology Astropart. Phys., 2019, 001

Zumalacárregui, M., Bellini, E., Sawicki, I., Lesgourgues, J., & Ferreira, P. G. 2017, J. Cosmology Astropart. Phys., 2017, 019
# Appendix A: Impact of the errors on the halo mass function

In this work, we use correcting factors for the HMF errors. Indeed, several methods yield different results for predicting the HMF in the HS-$f(R)$ scenario. We use the prediction of Hagstotz et al. (2019), but alternative approaches exist like the one introduced by Cataneo et al. (2016). These different approaches yield up to 20% differences in the prediction of the HMF.

We use the analytic form introduced in Costanzi et al. (2019), where

$$ \frac{d\tilde{n}^{f(R)}}{d \ln M_{500c}} = \frac{dn^{f(R)}}{d \ln M_{500c}}(s \ln(M/10^{14} M_{\odot}) + q), $$

We compare the posterior distribution of our parameters s and

<figure><img src="image_8.png" /><figcaption>Fig. A.1: Relative difference of the HMF in different MG scenarios. $\Delta n = dn/d \ln M(f_{R0}) - dn/d \ln M$, where $dn/d \ln M(f_{R0})$ is the HMF described with Eq. 20 and $dn/d \ln M$ is obtained with the GR fitting function from Tinker et al. (2008). We use the posterior distribution of the correcting parameters s and q and compare their impact to one of the MG parameters. The dark gray area corresponds to the 68% confidence level area, while the light gray one corresponds to the 95% confidence area.</figcaption></figure>

q and their impact on the HMF. This is shown in figure A.1. We fix the cosmology to the results of Planck Collaboration VI (2020). Then, we compute the HMF for GR (using the fit provided by Tinker et al. (2008)) and for different values of $f_{R0}$. Assuming that the GR HMF is correct, we then use the posterior distribution of the parameters s and q to quantify their impact on the HMF. We show that including the slope and amplitude correction of the HMF corrects most of the expected difference between the different models of $f(R)$, i.e., 20% difference. In the future, a more comprehensive modeling of the HMF uncertainties will be necessary to obtain high-precision measurements on the MG parameters.

# Appendix B: Acknowledgements

The authors thank Prof. Catherine Heymans for her help and availability for the implementation of the blinding strategy and her helpful comments on the draft.

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 the 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 Arge-lander Institute for Astronomy of the 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 (Brunner et al. 2022).

V. 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). CNES financially supported N. Clerc. M. Cataneo acknowledges the financial support provided by the Alexander von Humboldt Foundation through the Humboldt Research Fellowship program, as well as support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max Planck-Humboldt Research Award endowed by the Federal Ministry of Education and Research. T. Schraback and F. Kleinebreil acknowledge support from the German Federal Ministry for Economic Affairs 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. The Innsbruck group also acknowledges support provided by the Austrian Research Promotion Agency (FFG) and the Federal Ministry of the Republic of Austria for Climate Action, Environment, Energy, Mobility, Innovation and Technology (BMK) via the Austrian Space Applications Programme with grant numbers 899537 and 900565.

The 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.
Funding 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 Facilities Council 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 program 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).

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 work made use of the following Python software packages: SciPy<fnref n="1" /> (Virtanen et al. 2020), Matplotlib<fnref n="2" /> (Hunter 2007), Astropy<fnref n="3" /> (Astropy Collaboration 2022), NumPy<fnref n="4" /> (Harris et al. 2020), CAMB (Lewis & Challinor 2011), pyCCL<fnref n="5" /> (Chisari et al. 2019), GPy<fnref n="6" /> (GPy 2012), climin<fnref n="7" /> (Bayer et al. 2015), ultranest<fnref n="8" /> (Buchner 2021)

<footnote n="1">¹ https://scipy.org/</footnote>

<footnote n="2">² https:// matplotlib.org/</footnote>

<footnote n="3">³ https://www.astropy.org/</footnote>

<footnote n="4">⁴ https:// numpy.org/</footnote>

<footnote n="5">⁵ https://github.com/LSSTDESC/CCL</footnote>

<footnote n="6">⁶ https://github.com/SheffieldML/GPy</footnote>

<footnote n="7">⁷ https://github.com/BRML/climin</footnote>

<footnote n="8">⁸ https://github.com/JohannesBuchner/UltraNest/</footnote>