We present constraints on the fR gravity model using a sample of 1005 galaxy clusters in the redshift range 0.25  1.78 that have been selected through the thermal SunyaevZeldovich effect tSZE from South Pole Telescope SPT data and subjected to optical and nearinfrared confirmation with the Multicomponent Matched Filter MCMF algorithm. We employ weak gravitational lensing mass calibration from the Dark Energy Survey DES Year 3 data for 688 clusters at z lt 0.95 and from the Hubble Space Telescope HST for 39 clusters with 0.6 lt z lt 1.7. Our cluster sample is a powerful probe of fR gravity because this model predicts a scaledependent enhancement in the growth of structure which impacts the halo mass function HMF at cluster mass scales. To account for these modified gravity effects on the HMF our analysis employs a semianalytical approach calibrated with numerical simulations. Combining calibrated cluster counts with primary cosmic microwave background CMB temperature and polarization anisotropy measurements from the Planck2018 release we derive robust constraints on the fR parameter fR0. Our results log10 fR0 lt 5.32 at the 95  credible level are the tightest current constraints on fR gravity from cosmological scales. This upper limit rules out fRlike deviations from general relativity that result in more than a sim20  enhancement of the cluster population on mass scales Mmathrm200cgt3times1014Modot.
['Astrophysics - Cosmology and Nongalactic Astrophysics']
Constraints on fR gravity from tSZEselected SPT galaxy clusters and weak lensing mass calibration from DES and HST
{'Constraints on f ( R ) gravity from tSZE-selected SPT galaxy clusters and weak lensing mass calibration from DES and HST': "S. M. L. Vogt, 1, 2, 3, ∗ S. Bocquet, 1 C. T. Davies, 1 J. J. Mohr, 1, 4 F. Schmidt, 3 C.-Z. Ruan, 5 B. Li, 6 C. Hern'andez-Aguayo, 3 S. Grandis, 7, 1 L. E. Bleem, 8, 9 M. Klein, 1, 4 T. Schrabback, 7 M. Aguena, 10 D. Brooks, 11 D. L. Burke, 12, 13 A. Campos, 14, 15 A. Carnero Rosell, 16, 10 J. Carretero, 17 M. Costanzi, 18, 19, 20 L. N. da Costa, 10 M. E. S. Pereira, 21 J. De Vicente, 22 P. Doel, 11 S. Everett, 23 I. Ferrero, 24 J. Frieman, 25, 26 J. Garc'ıa-Bellido, 27 M. Gatti, 28 G. Giannini, 17, 26 D. Gruen, 1 R. A. Gruendl, 29, 30 S. R. Hinton, 31 D. L. Hollowood, 32 S. Lee, 33 M. Lima, 34, 10 J. L. Marshall, 35 J. Mena-Fern'andez, 36 R. Miquel, 37, 17 J. Myles, 38 M. Paterno, 25 A. Pieres, 10, 39 A. A. Plazas Malag'on, 12, 13 C. L. Reichardt, 40 A. K. Romer, 41 S. Samuroff, 42 A. Sarkar, 43 E. Sanchez, 22 I. Sevilla-Noarbe, 22 M. Smith, 44 E. Suchyta, 45 M. E. C. Swanson, 29 G. Tarle, 46 V. Vikram, 47 N. Weaverdyck, 48, 49 and J. Weller 1, 4  \n(the SPT and DES Collaborations) 1 University Observatory, Faculty of Physics, Ludwig-Maximilians-Universitat, Scheinerstr. 1, 81679 Munich, Germany 2 Excellence Cluster Origins, Boltzmannstr. 2, 85748 Garching, Germany 3 Max Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany 4 Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 2, 85748 Garching, Germany 5 Institute of Theoretical Astrophysics, University of Oslo, 0315 Oslo, Norway 6 Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK 7 Universitat Innsbruck, Institut fur Astro- und Teilchenphysik, Technikerstr. 25/8, 6020 Innsbruck, Austria 8 High-Energy Physics Division, Argonne National Laboratory, 9700 South Cass Avenue, Lemont, IL 60439, USA 9 Kavli Institute for Cosmological Physics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA 10 Laborat'orio Interinstitucional de e-Astronomia - LIneA, Rua Gal. Jos'e Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil 11 Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK 12 Kavli Institute for Particle Astrophysics & Cosmology, P. O. Box 2450, Stanford University, Stanford, CA 94305, USA 13 SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA 14 Department of Physics, Carnegie Mellon University, Pittsburgh, Pennsylvania 15312, USA 15 NSF AI Planning Institute for Physics of the Future, Carnegie Mellon University, Pittsburgh, PA 15213, USA 16 Instituto de Astrofisica de Canarias, E-38205 La Laguna, Tenerife, Spain 17 Institut de F'ısica d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona) Spain 18 Astronomy Unit, Department of Physics, University of Trieste, via Tiepolo 11, I-34131 Trieste, Italy 19 INAF-Osservatorio Astronomico di Trieste, via G. B. Tiepolo 11, I-34143 Trieste, Italy 20 Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014 Trieste, Italy 21 Hamburger Sternwarte, Universitat Hamburg, Gojenbergsweg 112, 21029 Hamburg, Germany 22 Centro de Investigaciones Energ'eticas, Medioambientales y Tecnol'ogicas (CIEMAT), Madrid, Spain 23 California Institute of Technology, 1200 East California Blvd, MC 249-17, Pasadena, CA 91125, USA 24 Institute of Theoretical Astrophysics, University of Oslo. P.O. Box 1029 Blindern, NO-0315 Oslo, Norway 25 Fermi National Accelerator Laboratory, P. O. Box 500, Batavia, IL 60510, USA 26 Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA 27 Instituto de Fisica Teorica UAM/CSIC, Universidad Autonoma de Madrid, 28049 Madrid, Spain 28 Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA 29 Center for Astrophysical Surveys, National Center for Supercomputing Applications, 1205 West Clark St., Urbana, IL 61801, USA 30 Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA 31 School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia 32 Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA 33 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA 34 Departamento de F'ısica Matem'atica, Instituto de F'ısica, Universidade de S˜ao Paulo, CP 66318, S˜ao Paulo, SP, 05314-970, Brazil 35 George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, and Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA 36 LPSC Grenoble - 53, Avenue des Martyrs 38026 Grenoble, France 37  \nInstituci'o Catalana de Recerca i Estudis Avan¸cats, E-08010 Barcelona, Spain  \n38 Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA 39 Observat'orio Nacional, Rua Gal. Jos'e Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil 40 School of Physics, University of Melbourne, Parkville, VIC 3010, Australia 41 Department of Physics and Astronomy, Pevensey Building, University of Sussex, Brighton, BN1 9QH, UK 42 Department of Physics, Northeastern University, Boston, MA 02115, USA 43 Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, 70 Vassar St, Cambridge, MA 02139 44 Physics Department, Lancaster University, Lancaster, LA1 4YB, UK 45 Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831 46 Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA 47 Argone National Laboratory, 9700 S. Cass Avenue, Lemont, IL 60439, USA 48 Department of Astronomy, University of California, Berkeley, 501 Campbell Hall, Berkeley, CA 94720, USA 49 Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA  \nWe present constraints on the f ( R ) gravity model using a sample of 1,005 galaxy clusters in the redshift range 0 . 25 -1 . 78 that have been selected through the thermal Sunyaev-Zel'dovich effect (tSZE) from South Pole Telescope (SPT) data and subjected to optical and near-infrared confirmation with the Multi-component Matched Filter (MCMF) algorithm. We employ weak gravitational lensing mass calibration from the Dark Energy Survey (DES) Year 3 data for 688 clusters at z < 0 . 95 and from the Hubble Space Telescope (HST) for 39 clusters with 0 . 6 < z < 1 . 7. Our cluster sample is a powerful probe of f ( R ) gravity, because this model predicts a scale-dependent enhancement in the growth of structure, which impacts the halo mass function (HMF) at cluster mass scales. To account for these modified gravity effects on the HMF, our analysis employs a semianalytical approach calibrated with numerical simulations. Combining calibrated cluster counts with primary cosmic microwave background (CMB) temperature and polarization anisotropy measurements from the Planck 2018 release, we derive robust constraints on the f ( R ) parameter f R 0 . Our results, log 10 | f R 0 | < -5 . 32 at the 95 % credible level, are the tightest current constraints on f ( R ) gravity from cosmological scales. This upper limit rules out f ( R )-like deviations from general relativity that result in more than a ∼ 20% enhancement of the cluster population on mass scales M 200c > 3 × 10 14 M ⊙ .", 'I. INTRODUCTION': "One of the most challenging questions in modern cosmology is understanding the nature of the accelerating expansion of the Universe [1, 2]. Various cosmological theories have been proposed to explain this phenomenon. Within the framework of general relativity (GR), this acceleration can be explained by introducing a cosmological constant Λ to the Einstein-Hilbert action, leading to the well-known Λ cold dark matter (ΛCDM) model. However, adding a cosmological constant to the EinsteinHilbert action offers little physical insight into the nature of Dark Energy. Therefore, there is strong motivation to consider modifications to the Einstein-Hilbert action that give rise to modified gravity models (see e.g., reviews [35]).  \nThese modifications impact the growth of cosmic structures. Consequently, the abundance of massive galaxy clusters, as the end products of the hierarchical growth of cosmic structures, is sensitive to the different matter clustering and therefore, serves as an excellent probe for constraining modified gravity models and offering an independent test of GR.  \nIn this work, we focus on a specific modified gravity model that introduces a non-linear function f ( R ) of the  \nscalar curvature R into the Einstein-Hilbert action [6]. We employ the widely studied Hu & Sawicki model for the function f ( R ) [7]. Physically, this model introduces an additional gravitational-strength fifth force, altering structure formation in a scale-dependent manner and enhancing structure formation on galaxy cluster scales. The extent to which the f ( R ) model deviates from GR is encoded in the single parameter f R 0 , which has been constrained using various observations on cosmological scales.  \nBecause the effects of f ( R ) gravity persist to very small scales, constraints on galactic and solar-system scales are very stringent for this model. Studies using galaxy rotation curves and morphology report log 10 | f R 0 | < -6 . 1 and log 10 | f R 0 | < -7 . 55 at 95% credible level, respectively [8, 9]. However, because these studies probe small scales, systematics such as uncertainties in galaxy formation play a key role. This makes constraints from larger scales that are less sensitive to galaxy formation a complementary test.  \nThe current tightest constraint from clusters comes from a combination of ROSAT clusters, primary cosmic microwave background (CMB) anisotropy data, Supernovae (SNe), and baryonic acoustic oscillations (BAO), with an upper bound of log 10 | f R 0 | < -4 . 79 at the 95% credible level [10]. The recent analysis of eROSITA clusters reports log 10 | f R 0 | < -4 . 12 at the 95 % credible level using clusters alone and marginalizing over the neutrino  \nmass [11].  \nStronger constraints from large scales probes are obtained from a weak-lensing peak analysis, which used data from the Canada-France-Hawaii-Telescope Lensing Survey (CFHTLenS), and are given by log 10 | f R 0 | < -5 . 16 [12]. Similar constraints come from the cross-correlation of galaxies with CMB lensing and galaxy weak-lensing, CMB, SNe, and BAO, reporting log 10 | f R 0 | < -4 . 61 and log 10 | f R 0 | < -4 . 5 at the 95 % credible level respectively [13, 14].  \nIn the analysis presented here, we focus on the weak lensing informed galaxy cluster abundance to constrain f ( R ) gravity. We are motivated to pursue this study partly because it has long been recognized that galaxy cluster surveys would be powerful probes of cosmic growth and therefore the action of gravity [15, 16] and also because in recent years cluster surveys have been successfully employed to study the standard ΛCDM model [17-25] as well as modified gravity models [11, 26-30].  \nOne of the most promising cluster samples currently available has been constructed using South Pole Telescope (SPT) [31] survey data and the thermal SunyaevZel'dovich effect (tSZE) [32]. The tSZE is caused by high-energy electrons in the intra-cluster medium (ICM) scattering off CMB photons, resulting in a spectral distortion of the CMB at the cluster position. Because the tSZE is a direct tracer of the hot ICM, it enables the detection of massive galaxy clusters. Moreover, the cluster tSZE signature is strongly mass dependent and approximately redshift-independent and thus can be employed to identify galaxy clusters up to the highest redshifts where clusters of sufficient mass exist. To constrain cosmological parameters with galaxy clusters, one has to relate observables such as the tSZE detection significance to the underlying halo mass at all relevant redshifts. These observable-mass relations can be empirically calibrated using weak gravitational lensing data and are typically modeled as power laws in cluster mass and redshift.  \nIn this study, we employ the sample of 1,005 galaxy clusters detected using SPT data and confirmed using the MCMF algorithm [33-35] with optical and near-infrared data from the Dark Energy Survey (DES) [36-38] and the Wide-field Infrared Survey Explorer (WISE) [39]. To obtain mass estimates for the cluster sample, we use weak-lensing measurements from DES and targeted observations from the Hubble Space Telescope (HST). The f ( R ) analysis framework employed for the SPT clusters with mass calibration from DES and HST is based on the state-of-the-art method developed for the recent ΛCDM analyses of this same sample [24, 40] (hereafter SB24a and SB24b). In a recent paper, this same framework was modified and employed to carry out validation tests and forecasts for f ( R ) gravity constraints from upcoming Stage-III and Stage-IV surveys [41] (hereafter SV24a).  \nFollowing SV24a, we incorporate f ( R ) gravity into our analysis by modifying the halo mass function (HMF), which is enhanced relative to GR, in a mass- and redshiftdependent way. We implement this modification by  \nintroducing a multiplicative factor to the GR HMF [28, 42], dependent on the spherical collapse threshold for halo collapse in f ( R ) gravity for which we use a semianalytical model [43, 44]. In the present analysis we calibrate this semi-analytical HMF model against the f ( R ) FORGE numerical simulations to obtain a more accurate halo mass function [45]. The simulations are not used directly to predict the HMF, because the available mass range from the simulations is limited, whereas the HMF from the semi-analytical model can be calculated for the wide mass range needed in this analysis.  \nThis paper is organized as follows. Section II presents a summary of the SPT cluster dataset and the DES and HST weak-lensing data. We review in Sec. III f ( R ) gravity and the f ( R ) HMF model used in this work as well as the calibration to the FORGE simulations. In Sec. IV, we discuss our analysis method, including the observablemass relations, weak-lensing model, likelihood approach, and priors. The results are presented in Sec. V, and we conclude with a summary in Sec. VI.  \nThroughout this paper, U ( a, b ) denotes a uniform distribution between limits a and b , and N ( µ, σ 2 ) is a Gaussian distribution with mean µ and variance σ 2 . We adopt the halo mass definition M 200c , which is the mass within a radius where the mean density is 200 times the critical density.", 'II. DATA': 'This section gives a brief summary of the cluster and weak-lensing data we use in this work. A detailed description of the data products is presented in SB24a.', 'A. SPT cluster catalog': 'The tSZE selected cluster catalogs from the SPT-SZ, SPTpol ECS and SPTpol 500d surveys employed here cover a total solid angle of 5 , 270 deg 2 of the southern sky [34, 35, 46, 47]. Note that the whole SPTpol 500d survey lies within the SPT-SZ footprint, and we use only the data from the deeper SPTpol 500d survey in the overlapping region. Cluster candidates of these surveys are selected in tSZE detection significance ˆ ζ and confirmed using optical and infrared data, which also add redshift information. Over the whole SPT survey region, only clusters with z > 0 . 25 are included in the sample, because the impact of atmospheric and primary CMB noise on the selection function at low redshift becomes dominant.  \nIn the SPT survey region that is not covered by DES (1 , 327 deg 2 , ca. 27 % of the total solid angle), the cluster candidates are confirmed by targeted observations that also provide redshift measurements. Candidates in this  \nregion are selected by  \nˆ ζ > 5 z > 0 . 25 , (1)  \nwhich results in a sample with 110 clusters and a purity ≳ 95 % [46, 47].  \nFor the region covered by DES (3 , 567 deg 2 , 75 % of the SPT area) cluster candidates are optically confirmed using the Multi-Component Matched Filter cluster confirmation tool [MCMF; 33, 34], and we follow the work of Refs. [34, 35]. Moreover, measurements for cluster redshift z , optical richness ˆ λ and optical center position are obtained using MCMF. Because DES data are only reliable for redshifts z ≤ 1 . 1, WISE data are used to compute richnesses and redshifts for clusters with z > 1 . 1. There is the chance that an overdensity of galaxies is a random superposition along the line of sight of a tSZE noise fluctuation. To exclude chance associations, we apply a redshift dependent richness cut, ˆ λ min ( z ), which ensures a > 98 % purity over the entire redshift range probed and a cluster candidate is then confirmed if the measured richness is larger than this threshold. Using MCMF for cluster confirmation allows us to validate clusters with lower tSZE detection significance while achieving a high sample purity, and results in a cluster sample, which is 30 % larger using the MCMF tool.  \nDue to the different depths of the individual tSZE surveys, different selection thresholds in ˆ ζ are applied to obtain an approximately constant purity for the combined SPT sample. The selection criteria are  \nˆ ζ > 4 . 25 / 4 . 5 / 5 (500d / SZ / ECS) , ˆ λ > ˆ λ min ( z ) , z > 0 . 25 . (2)  \nThese selections result in a sample of 895 confirmed clusters over this region.  \nTo summarize, the total cluster sample consists of 1,005 clusters, each characterized by the observables tSZE detection significance ˆ ζ and redshift z . In the SPT area covered by DES, these clusters also have additional measurements of richness ˆ λ and cluster center position from either DES ( z ≤ 1 . 1) or WISE ( z > 1 . 1).  \nFigure 1 shows the 1,005 confirmed cluster sample in the space of tSZE detection significance ˆ ζ and redshift (left) as well as in optical richness ˆ λ and redshift (middle). The clusters are color coded according to which of the three SPT surveys they originate from. In addition, the corresponding selection thresholds are shown as lines in each observable. One can see in the ˆ ζ -z distribution that while the sample extends to z ∼ 1 . 8 the bulk of the clusters lie at z < 1. In addition, this distribution makes clear that the majority of the cluster lies close to the tSZE detection thresholds of the surveys. The ˆ λ -z distribution displayed in the middle panel shows that the bulk of the clusters have relatively high richnesses  \nˆ λ > 30, and that the optical selection thresholds are not significantly impacting the completeness of the confirmed cluster sample. Indeed, the MCMF selection thresholds largely impact the high purity of the sample, because the noise fluctuations present in the tSZE selected candidate clusters are efficiently removed due to their low optical richnesses ( ˆ λ < 10).', 'B. DES Y3 weak-lensing data': 'The DES was conducted in the griz Y bands and covers a sky area of 5 , 000 deg 2 . The DES Y3 weak lensing shape catalog [48] utilized data from the first three years of observations and covers approximately 4 , 143 deg 2 of the sky after masking. 3 , 567 deg 2 of the DES region overlaps with the SPT surveys, corresponding to 75 % of the whole SPT survey area. The weak-lensing shape catalog of DES Y3 is built with the Metacalibration pipeline from the r , i and z bands [49, 50]. Lensing source galaxies are selected in four tomographic redshift bins as employed in the 3 × 2pt analysis of DES [51].  \nFor each cluster in the overlapping region of SPT and DES, we use the weak-lensing shear profiles within the radial range 0 . 5 < r/ ( h -1 Mpc) < 3 . 2 (1 + z cluster ) -11 around the optical cluster center. The lower limit on the radial range excludes the inner region of the cluster, which is largely affected by feedback from active galactic nuclei, miscentering, blending, cluster member contamination and non-linear shear. The upper limit on the radial range guarantees that only the one-halo term region is used for mass calibration [52]. We only use weaklensing data from DES for clusters with z < 0 . 95, corresponding to the median redshift of the highest redshift tomographic bin [SB24b, SB24a].  \nIn our analysis, we account for systematic and statistical uncertainties such as cluster member contamination, miscentering of the shear profile, shear and photo-z calibration, halo mass modeling and the impact of largescale structure. A detailed description of the modeling of these uncertainties can be found in detail in SB24a, Section V. Note that the calibration of these uncertainties was performed within the ΛCDM paradigm, but we expect these to not change significantly in f ( R ) gravity [53]. In total, our analysis includes 688 cluster shear profiles from 555 , 912 source galaxies with an average of 808 shear source galaxies per cluster [SB24a]. We show for illustration the averaged matter density profile of the DES Y3 data in the right panel of Fig. 1. The profiles are broken into four redshift ranges of comparable signalto-noise. The combined dataset corresponds to a 31 . 2 σ detection of the matter profiles of these clusters.  \nFIG. 1. Left : tSZE detection significance ˆ ζ and redshift distribution for the three SPT surveys. Dashed-colored lines show the detection threshold for the corresponding survey in the region that overlaps with DES. The black dashed line shows the ˆ ζ threshold for the clusters outside of the DES region (same as the cut for the SPTpol ECS survey). Middle : optical richness ˆ λ and redshift distribution of the cluster sample color coded by the three surveys. Colored lines correspond to the ˆ λ detection threshold of the given survey. The SPTpol 500d is significantly deeper than the other two surveys (yellow-green dashed line) and thus a lower ˆ λ ( z ) threshold is applied (solid blue line). Right : averaged weak-lensing inferred projected matter profiles for DES Y3 data shown for four redshift bins in purple to yellow. The redshift bins are chosen such that the signal-to-noise is approximately equal in each bin. The HST data are shown in blue. Over the radial range used for weak-lensing mass calibration, the DES Y3 data have a signal-to-noise of 31 . 2 and the HST data have 9 . 7.  \n<!-- image -->', 'C. HST weak-lensing data': 'DES lensing data are only reliable for z ≲ 0 . 95 and therefore, we complement the weak-lensing dataset with HST data to obtain weak-lensing information for the mass calibration at high redshift. We use the HST39 dataset [54-56] to obtain weak-lensing shear profiles. This dataset contains 39 clusters of our tSZE selected sample in the redshift range 0 . 6 -1 . 7. More details about the dataset and the analysis can be found in Refs. [54-59]. The averaged matter density profile from the 39 clusters from HST is shown in Fig. 1 in the right panel in blue. These data correspond to a 9 . 7 σ detection of the matter profiles of these halos.', 'III. f ( R ) MODIFIED GRAVITY': 'f ( R ) gravity modifies GR by introducing an arbitrary function f ( R ) of the Ricci scalar R into the EinsteinHilbert action [6]  \nS = ∫ d 4 x √ -g [ R + f ( R ) 16 πG + L m ] , (3)  \nwhere g denotes the determinant of the metric tensor, G is the gravitational constant and L m is the matter Lagrangian density. Note that we use natural units where c = ℏ = 1 and if f ( R ) = -2Λ we recover GR plus a cosmological constant, hence a ΛCDM cosmology. The  \nfield equation obtained from this modified action takes the form  \nG µν + f R R µν -( f 2 -□ f R ) g µν -∇ µ ∇ ν f R = 8 πGT µν . (4)  \nHere G µν is the Einstein tensor, R µν is the Ricci tensor, T µν represents the energy-momentum tensor and f R = d f ( R ) / d R is an additional scalar degree of freedom, which indicates the strength of the modifications to GR.  \nThe equation of motion for f R is derived from the trace of the field equation in the quasistatic and weak-field limit  \n∇ 2 δf R = 1 3 ( δR -8 πGδρ ) , (5)  \nwhere δx = x -¯ x is the perturbation of the quantity x with respect to the cosmic mean. Additionally, from the time-time component of the field equation, Eq (4), we derive the modified Poisson equation in f ( R ) gravity:  \n∇ 2 Φ = 16 πG 3 δρ -1 6 δR, (6)  \nwith Φ the Newtonian potential which is defined via 2Φ = δg 00 /g 00 . Combining Eqs. (5) and (6) shows how the Poisson equation and thus the structure growth depends on the strength of the f ( R ) gravity model  \n∇ 2 Φ = 4 πGδρ -1 2 ∇ 2 δf R . (7)  \nCompared to GR the Poisson equation includes a term proportional to the Laplacian of δf R and thus depends directly on the strength of the model.  \nThe modified Poisson equation, Eq. (6), shows that the strength of the gravitational force depends on the environment. In low curvature regions, δR ≪ 8 πGδρ , and consequently the f ( R ) Poisson equation reduces to a modified Poisson equation with gravitational forces enhanced by a factor of 4 / 3 compared to GR. In high curvature environments, the Ricci scalar R is approximately 8 πGδρ and thus Eq. (6) corresponds to the unmodified GR Poisson equation, i. e. in high-density regions f ( R ) gravity falls back to GR. This property of driving f R → 0 in high-density regions is the so-called chameleon screening mechanism [60], which makes f ( R ) gravity consistent with Solar systems tests [7, 61]. The two limits show that structure growth in f ( R ) gravity is environmentdependent in contrast to GR.  \nIn this paper, we adopt the widely used and studied Hu & Sawicki model [7]  \nf ( R ) = -m 2 c 1 ( R m 2 ) n c 2 ( R m 2 ) n +1 , (8)  \nwith m 2 = Ω m H 2 0 , H 0 the Hubble constant and the free parameters n , c 1 , c 2 . In the high curvature regime we have c 1 /n 2 R/m 2 ≫ 1 and Eq. (8) is approximately given by  \nf ( R ) ≈ -m 2 c 1 c 2 -f R 0 R n +1 0 nR n . (9)  \nHere R 0 is the present background curvature and f R 0 := f R ( R 0 ), which is the parameter that quantifies the strength of this f ( R ) gravity model.  \nTo have a modified gravity scenario that is active at late times and large scales, when the acceleration of the Universe happens, and that does not spoil the successful description of early Universe observables such as BBN and CMB the function f ( R ) has to satisfy f R < 0, i. e. it is a decreasing function of R . Therefore, the approximation above is correct up to order ∼ ( f R 0 ) 2 , and since the constraints from our cluster sample are better than | f R 0 | ∼ 10 -4 , the approximation in Eq. (9) is entirely sufficient. We work with log 10 | f R 0 | for numerical convenience, and, given that the theory gives no strong prior on the scale of f R 0 , we impose a uniform prior on log 10 | f R 0 | .  \nTo obtain an expansion history consistent with ΛCDM in the limit | f R 0 | → 0, the parameters c 1 and c 2 are given by  \nc 1 c 2 = 6 Ω Λ Ω m . (10)  \nWe further adopt n = 1 for the Hu & Sawicki model; see Ref. [62] for an approach to approximately rescale constraints from n = 1 to other values of n .', 'A. The halo mass function in f ( R ) gravity': 'As we constrain f ( R ) gravity with the help of cluster abundance datasets we need a model for the distribution of halo mass and redshift within this theory, in other words the (differential) halo mass function. In this work we adopt a two-component model [28, 42] for which the first component is the HMF in GR, and the second factor models the enhancement of the f ( R ) gravity, which accounts for the scale-dependent clustering due to the scale-dependent structure growth described in Sec. III:  \nd n dln M = d n dln M ∣ ∣ ∣ ∣ GR ×R . (11)  \nFor the GR HMF, we use the halo mass function from Ref. [63]  \nd n dln M ∣ ∣ ∣ ∣ GR = -¯ ρ m 2 M f ( σ ) T dln σ 2 dln M , (12)  \nwith the multiplicity function of the form  \nf ( σ ) T = ˜ A [ ( σ ˜ b ) -˜ a +1 ] e -˜ c σ 2 , (13)  \nwhere ˜ A, ˜ a, ˜ b and ˜ c are parameters calibrated using Nbody simulations [see 63, table 2] and σ = σ ( M ) is the variance of the overdensity field on a mass scale M in the corresponding GR cosmology.  \nThe enhancement factor R is calculated from the ratio of the Sheth-Tormen HMF [64] in f ( R ) gravity to GR, i. e.  \nR = d n dln M ∣ ∣ ST , f ( R ) d n dln M ∣ ∣ ST , GR . (14)  \nThis HMF can account for the scale-dependent clustering through the spherical collapse threshold δ crit . In this work, we calculate the spherical collapse threshold δ crit from the spherical collapse model of Refs. [43, 44]. With this, the Sheth-Tormen HMF is given by  \nd n dln M ∣ ∣ ∣ ∣ ST = ¯ ρ m M f ( ν ) ST [ dln δ crit dln M -1 2 dln σ 2 dln M ] . (15)  \nHere ν = δ crit /σ is the peak height and f ( ν ) ST refers to the Sheth-Tormen multiplicity function [64],  \nf ( ν ) ST = A √ aν 2 2 π [ 1 + ( aν 2 ) -p ] e -aν 2 2 , (16)  \nwith A, a , p are free parameters for which we adopt the parameters of Ref. [65].  \nThe change in the clustering and thus in the HMF is completely captured by the scale-dependent spherical collapse threshold δ crit and therefore, we use an effective variance σ derived from the corresponding GR cosmology [44, 66]. Moreover, when including massive neutrinos, the HMF shape is closer to universal if the halo  \nmass and variance are calculated neglecting the neutrino component [67, 68], i. e. using the baryon and cold dark matter power spectrum instead of the total matter power spectrum for the variance. Since there is a degeneracy between massive neutrinos and log 10 | f R 0 | (at least for high values of Σ m ν and log 10 | f R 0 | ) [see e. g. 69-71] we use the approach of Ref. [67, 68] to account for the non-zero total neutrino mass.  \nThe derivation of δ crit with the spherical collapse model requires solving a computationally expensive system of coupled differential equations for each value of mass, redshift and cosmology [43, 44]. To speed up the calculations we use emulators to predict the values of δ crit and d δ crit / dln M . The emulators adopted in this analysis are presented in SV24a.  \nFigure 2 shows the f ( R ) HMF on top and the ratio R on the bottom for different values of log 10 | f R 0 | and redshift. In f ( R ) gravity, the HMF is larger compared to GR as expected from the enhanced structure growth. The shape and strength of the enhancement depend on the value of log 10 | f R 0 | with a larger enhancement for stronger f ( R ) models. Furthermore, the enhancement becomes smaller at higher redshifts, because modified gravity falls back to GR in the early universe. A large enhancement of the most massive halos is only seen for the strongest model with log 10 | f R 0 | = -4 . 5 and at lower redshifts, because even the most massive halos are only partially screened in this scenario. In the other cases, the potential Ψ of the high mass halos becomes compatible with f R and the fifth force is screened in these massive halos, reducing the modified gravity effect on the halo abundance.', 'B. HMF Calibration using FORGE simulations': 'The model for the HMF presented above follows a semi-analytical approach where the f ( R ) gravity is incorporated via the spherical collapse threshold δ crit calculated from spherical collapse theory in f ( R ) gravity. Besides semi-analytical models, the HMF can also be derived from simulations. Therefore, we make a comparison with f ( R ) simulations to validate the semi-analytical model and to check for any kind of discrepancy.  \nIn this work, we use the state-of-the-art FORGE N -body simulations which encompass 49 f ( R ) gravity cosmologies (nodes) [45]. The simulations sample the cosmological parameters Ω m , h , S GR 8 = σ GR 8 √ Ω m / 0 . 3, and log 10 | f R 0 | (hereafter FORGE parameters) with a latin hypercube, while the other cosmological parameters are fixed to n s = 0 . 9652 Ω b = 0 . 049199, Ω ν = 0 and Ω Λ = 1 -Ω m [45]. Each FORGE f ( R ) gravity node has a ΛCDM counterpart. The parameter ranges explored in  \nFIG. 2. The f ( R ) HMF, Eq. (11), (top) and the ratio R of the f ( R ) and the GR Sheth-Tormen HMF, Eq. (15), (bottom) for different values of log 10 | f R 0 | (different colors) and three redshifts (different lifestyles) corresponding to the minimum, mean and maximum redshift of the SPT cluster sample. The deviation from GR depends on the strength of the f ( R ) gravity and redshift, whereby weaker f ( R ) models and higher redshifts show less enhancement in the HMF relative to GR.  \n<!-- image -->  \n/circledot  \nthe FORGE simulations are  \n0 . 11 < Ω m < 0 . 54 , 0 . 61 < h < 0 . 81 , 0 . 6 < S GR 8 < 0 . 9 , (17) -6 . 17 < log 10 | f R 0 | < -4 . 51 .  \nThe ranges in Ω m and S GR 8 translate to a range in σ GR 8 of 0 . 49 < σ GR 8 < 1 . 31.  \nRef. [53] presents an emulator for the f ( R ) HMF that is based upon the FORGE simulations. In this work, the authors trained a neural network to directly predict the enhancement factor of the HMF due to f ( R ) gravity for each mass, redshift and set of FORGE parameters. However, by design, the resulting halo mass function is only valid within the mass, redshift and cosmological parameter range used in the neural network training set. To make predictions outside the mass range available within the FORGE simulations, we calibrate the semi-analytical HMF model, which is valid up to halo masses of 10 16 h -1 M ⊙ , with the HMF retrieved from the FORGE simulations.  \nTo calibrate the semi-analytical HMF model with the  \nFORGE simulations we use the high-resolution simulations with 1024 3 dark matter particles in a box with length L = 500 h -1 Mpc and a mass resolution of 9 . 5 × 10 9 (Ω m / 0 . 3) h -1 M ⊙ . We extract halo catalogs from all f ( R ) simulations as well as from their corresponding ΛCDM nodes at redshifts z = 0 . 00 , 0 . 25 , 0 . 5 , 0 . 75 , 1 . 00 , 1 . 25 , 1 . 50 , 1 . 75, and 2 . 00, using a bin width of 0.1 in log 10 M within the mass range 10 13 h -1 M ⊙ ≤ M 200c ≤ 5 × 10 15 h -1 M ⊙ . To ensure no empty bins at the high-mass end, we combine the last high-mass bins into one (large) bin such that it contains at least 20 halos 2 . We compute covariance matrices for each halo catalog using 5 3 jackknife samples to account for noise due to sample variance and shot noise.  \nWe characterize the difference between the semianalytical model and the simulation by comparing the enhancement in the HMF, R , from the simulations and the semi-analytical HMF model. This has the advantage that the cosmic variance in the simulations partially cancels in this ratio. With the halo catalogs from f ( R ) gravity and ΛCDM FORGE simulations, we can calculate the enhancement in the cluster counts, R FORGE , for each node and redshift. The enhancement from the semi-analytical model, R SAM , is also calculated for each FORGE cosmology and redshift.  \nThe top panels of Fig. 3 show the ratios between the enhancements in each mass bin, i. e. R FORGE / R SAM , for four of the nine different redshifts with error bars derived from the jackknife covariance of the simulations. There is a bias between the semi-analytical HMF and the FORGE simulations, which varies with mass, redshift and the FORGE parameters. In Fig. 3 the ratio is color-coded based on the value of log 10 | f R 0 | and the discrepancy is larger with higher values of log 10 | f R 0 | . Note that also the other cosmological parameters can drive the difference between the semi-analytical HMF and the simulations. Therefore, we assume that the bias depends on all FORGE parameters. Overall the semi-analytical HMFmodel predicts more clusters than the FORGE simulations, and the agreement is better for weaker f ( R ) models, i. e. smaller values of log 10 | f R 0 | .  \nBased on the top panels of Fig. 3 we model the ratio and thus the correction to the semi-analytical HMF with a broken linear function in log 10 M with a pivot scale log 10 M piv and smooth transition with strength k between the two linear functions at given logarithmic mass log 10 M 1 . To be precise, for log 10 M smaller than log 10 M 1 , the ratio is modeled with a linear function in log 10 M with slope a and intercept b . For log 10 M larger than log 10 M 1 , the same ratio is modeled by a different linear function with slope a 2 = a + ∆ a and intercept fixed by continuity at the transition mass M 1 . We then smooth the transition between the two linear relations by a power-law interpolation controlled by a parameter k . For fixed pivot mass, transition mass and smoothing  \nstrength, the fitting function for the correction has three free parameters and is given by  \nc ( M, p , z ) = a ( p , z )(log 10 M -log 10 M piv ) + b ( p , z ) +∆ a ( p , z ) ln(1 + e -k (log 10 M -log 10 M 1 ) ) k +∆ a ( p , z )(log 10 M -log 10 M 1 ) . (18)  \nHere p represents the vector containing the FORGE parameters, a ( p , z ), ∆ a ( p , z ) , b ( p , z ) are the fitting parameters of the correction for which we assume a dependence on the FORGE parameters and redshift. We decided to fit for the parameter ∆ a ( p , z ) instead of a 2 ( p , z ) = a ( p , z ) + ∆ a ( p , z ) because it showed a better behavior during fitting.  \nThe parameters of the broken linear function are obtained by fitting the correction function, Eq. (18), to the ratio between the enhancements, R FORGE /R SAM , where we fix the parameters k = 2, log 10 M 1 = 14 . 5 and log 10 M piv = 13. With this correction, we achieve good agreement with the simulations within the error bars as shown in the bottom panels of Figure 3.  \nBecause we assume that the fitting parameters depend on cosmology and redshift we have to predict these parameters for an arbitrary cosmology and redshift to use the calibrated HMF in our analysis. To do so we build emulators based on Gaussian process regression to predict the three fitting parameters as a function of the FORGE parameters and redshift.  \nWith the emulators, we can calculate the correction to the HMF for an arbitrary cosmology (within the parameter ranges Eq. (17)), and our calibrated semi-analytical HMF is given by  \nd n dln M = d n dln M ∣ ∣ ∣ ∣ GR × c ×R , (19)  \nwhere R is the enhancement factor in the semi-analytical halo mass function from Eq. (14) and c is the correction function, Eq. (18), calibrated with the FORGE simulations. As we model the correction function c ( M, p , z ) with a broken linear function, the enhancement of the HMF, c R , may be reduced to values below one, i. e. the f ( R ) HMF is suppressed rather than enhanced. Because the semi-analytical model always predicts an enhancement and also the FORGE simulations show no sign of suppression apart from noise at the high-mass end, we enforce c R ≥ 1 for the mass range considered in this work. Specifically, we set c R = 1 for all masses above the mass where the corrected HMF factor drops below unity.  \nUsing the 49 FORGE simulations at 9 different redshifts we calibrate the semi-analytical HMF model presented in Sec. III A; however, there are limitations to this approach that we have to address. First, due to the relatively small box size of 500 h -1 Mpc the halo catalogs contain only halos with M ≲ 4 × 10 15 h -1 M ⊙ depending on cosmology and redshift. Thus the correction function is only calibrated in the mass range  \nFIG. 3. Comparison of the enhancement in the HMF from the FORGE simulations, R FORGE , and the semi-analytical HMF model, R SAM for different redshifts, color-coded by the log 10 | f R 0 | values. Grey dashed lines are plotted at 10 % deviation to guide the eye. The upper four plots show the comparison of the FORGE simulations to the semi-analytical model, Eq. (11), without a correction. The four plots at the bottom show the comparison to the corrected semi-analytical HMF, Eq. (19), leading to a better agreement between the two HMF enhancements. Error bars are derived from the jackknife covariance of the FORGE simulations.  \n<!-- image -->  \n/circledot  \n/circledot  \n/circledot  \n/circledot  \nthat is available from the simulations. However, the semi-analytical prediction gives an overall shape of the f ( R ) HMF. Moreover, the halo mass function is exponentially suppressed at high masses and therefore a correction in this regime has negligible impact on our cosmological analysis. The second limitation is the available range for log 10 | f R 0 | values. FORGE samples log 10 | f R 0 | in the range -6 . 17 < log 10 f R 0 < -4 . 51, so we can only predict the parameters of the correction function in this range. Therefore, we adopt a hard upper limit prior of log 10 | f R 0 | = -4 . 51 in our analysis. We extend our analysis to log 10 | f R 0 | = -7 based on the fact that the f ( R ) HMF approaches the GR HMF for log 10 | f R 0 | → -∞ and thus a, ∆ a → 0 and b → 1. Thus, we interpolate the parameters between log 10 f R 0 = -6 . 17 and log 10 f R 0 = -7 under the assumption that a model with log 10 f R 0 = -7 is indistinguishable from GR with the SPT dataset. Moreover, the correction from the simulation is smaller as we approach GR. We will show in our analysis that the upper bound has no impact on the f R 0 posterior, as the data strongly disfavor values of log 10 f R 0 > -4 . 51. Given these limitations on the calibration, we present both the results derived when using the semi-analytical HMF model and when using the FORGE-informed calibrated semi-analytical HMF model.  \nIn the analysis presented below, we account for the impact of remaining uncertainties in the HMF by following the approach of Ref. [72] and introducing uncertainties  \nin the amplitude and the logarithmic mass trend of the HMF, i. e. the HMF with uncertainties given by  \nd n dln M = d n dln M ( q + s ln ( M 200c 10 14 h -1 M ⊙ )) , (20)  \nwhere q is the uncertainty in the amplitude and s is the uncertainty in the trend with logarithmic mass and we marginalize over the q and s in our analysis (see Sec. IV D)', 'IV. ANALYSIS METHOD': 'The method we employ in this work is based on the state-of-the-art weak lensing informed cluster cosmological analysis of the SPT sample [SB24a]. The method was also used and validated in the recent f ( R ) gravity forecast [SV24a] for SPT-3G [73] and CMB-S4 [74] cluster samples with next-generation weak-lensing data like those expected to come from the Euclid mission [75, 76] or the Vera C. Rubin Observatory [77, 78].', 'A. Observable-mass relations': 'In tSZE cluster surveys, galaxy clusters are identified and selected by observables such as the tSZE detection significance and richness and observable-mass relations link these observables to the halo mass [e.g., 79, 80].  \nThrough gravitational weak lensing calibration of these relations, we can relate the observed cluster sample to the HMF, which describes the abundance of halos depending on cosmology, mass, and redshift. In this analysis, we employ observable-mass relations that are empirically calibrated with weak-lensing data [e.g. 24, 81-86]. This section outlines the observable-mass relation for tSZE detection significance and optical/NIR richness.', '1. tSZE ζ -mass relation': 'As in previous SPT studies, we first relate the observed tSZE detection significance ˆ ζ to the intrinsic detection significance ζ to account for noise in the data. The relation between ˆ ζ and ζ is given by [87]  \nP ( ˆ ζ | ζ ) = N ( √ ζ 2 +3 , 1 ) . (21)  \nThe distribution accounts for the Gaussian noise present in the survey maps, with a correction factor of 3 due to the noise resulting from the matched-filter search for peaks in three dimensions. The mean intrinsic tSZE detection significance ζ is then modeled by  \n⟨ ln ζ ⟩ = ln ζ 0 + ζ M ln ( M 200c 3 × 10 14 h -1 M ⊙ ) + ζ z ln ( E ( z ) E (0 . 6) ) , (22)  \nwhere ζ 0 , ζ M and ζ z are the parameters corresponding to the normalization, mass and redshift trend of the scaling relation and E ( z ) = H ( z ) /H 0 . Adiitionally, we assume a lognormal intrinsic scatter in ζ with width σ ln ζ .  \nBecause the SPT surveys vary in depth, and we want to employ one ζ -mass relation for all surveys, the normalization ζ 0 and the redshift trend ζ z are rescaled for each field [35, 46, 47], i. e. ζ 0 , field = γ field ζ 0 and ζ z, field = ζ z +constant. In the case of the SPTpol ECS survey fields the normalization is difficult to calibrate and thus the parameter γ ECS is allowed to vary in the analysis [SB24a]. For the redshift trend ζ z the variation of the rescaling parameter c across fields is negligible in the SPT-SZ and SPTpol surveys. Therefore, we rescale ζ z for each survey where the SPT-SZ survey is taken as the reference [35]:  \nζ z, SPT-SZ = ζ z , ζ z, SPTpol ECS = ζ z -0 . 09 , ζ z, SPTpol 500d = ζ z +0 . 26 . (23)', '2. Cluster richness λ -mass relation': 'As for the tSZE detection significance, the observed richness ˆ λ is related to the intrinsic richness λ by a Gaussian distribution of the form  \nP ( ˆ λ | λ ) = N ( λ, √ λ ) . (24)  \nThis relationship accounts for Poisson sampling noise. Note that the Gaussian approximation of a Poisson distribution approximately holds for λ ≳ 10, which is below the richness cut we apply to our sample. Similar to the ζ -mass scaling relation we assume for the mean intrinsic richness a power law, in mass and (1 + redshift):  \n⟨ ln λ ⟩ = ln λ 0 + λ M ln ( M 200c 3 × 10 14 h -1 M ⊙ ) + λ z ln ( 1 + z 1 . 6 ) . (25)  \nThe parameters λ 0 , λ M and λ z govern the normalization, mass and redshift trend, respectively. The cluster intrinsic richness varies around this relation by a log-normal distribution with a width σ ln λ .  \nAs mentioned in Sec. II A we use richness measurements from DES for clusters with z ≤ 1 . 1 and data from WISE for high-redshift clusters. As matching two distinct types of richness measurements is challenging, we use two separate λ -mass relations for the DES and WISE data [SB24b], hereafter denoted by subscripts DES and WISE respectively.', 'B. Weak-lensing model in f ( R ) gravity': 'With the above described observable-mass relations, we can model the cluster sample in the ζ -λ -z space by transforming the halo mass function into the halo observable function and using it to predict cosmological parameters. However, there are no informative priors on the parameters of the ζ and λ scaling relations and their scatters. To empirically calibrate these relations and the corresponding scatters we rely on weak-lensing data. It has been shown that weak-lensing measurements are a robust way to measure halo masses with well-characterized and controllable biases [SB24b, SB24a].  \nIn this analysis, we assume that any f ( R )-gravity modification of the mapping from cluster potential to lensing signal can be neglected, as in previous works [11, 41]. First, the f ( R ) effect on the lensing signal for a given fixed mass distribution is given by a rescaling of the GR signal by a factor of (1 + | f R 0 | ) -1 [88, 89], which is negligible for the values of | f R 0 | we consider in this work. Second, while the cluster observables and halo profiles do undergo modifications in f ( R ) gravity, these effects are small [53, 90, 91]. Any changes to the cluster observables ˆ ζ and ˆ λ will be accounted by the empirical calibration of the observable mass relations. Changes to the halo profiles are more concerning, because they could impact the weak lensing inferred cluster masses. Accounting for these effects self-consistently within the weak lensing model described below would require a study of the halo shapes using f ( R ) numerical simulations and measurement of any changes in the inferred weak-lensing masses. Because we know these effects are smaller than the current uncertainties on the weak lensing model, which are  \ndominated by uncertainties in the hydrodynamical effects and on photometric redshift systematics, we adopt the GR-based calibration of the weak lensing model presented in SB24a and described below.', '1. DES weak-lensing model': 'The model we adopt for DES weak-lensing data was studied and described in detail in SB24a and works referenced therein. Here we provide a summary of the method. The weak-lensing observable is the reduced tangential shear profile, which is related to the underlying projected halo mass distribution Σ by  \ng t ( r, M WL ) = ∆Σ( r, M WL ) Σ -1 crit 1 -Σ( r, M WL ) Σ -1 crit . (26)  \nHere ∆Σ( r ) ≡ ⟨ Σ( < r ) ⟩-Σ( r ) is the surface density contrast and Σ -1 crit is the lensing efficiency or inverse critical surface mass density, given by  \nΣ -1 crit = 4 πG c 2 D l D s × max[0 , D ls ] , (27)  \nwhere c is the speed of light and D s , D l , D ls are the angular diameter distances between the observer and the source, the observer and the lens, and the source and the lens, respectively. We model Σ by the line of sight integral of a Navarro-Frenk-White profile (NFW) [92, 93] and we refer to the associated mass as the weak-lensing mass, M WL . We account for possible miscentering of the selected cluster center by assuming a constant density within the cluster miscentering radius R min , i. e. Σ( R ) = Σ( R min ) for R ≤ R min (see SB24a section IVC for more details). Cluster member contamination f cl ( r ) is corrected by a factor (1 -f cl ( r )) to the reduced tangential shear profile [SB24a].  \nGiven that the cluster profiles are not perfectly described by an NFW profile, the computed weak-lensing mass M WL is a biased and noisy estimator of the true cluster mass M 200c [94, 95]. To account for the bias we use a scaling relation between M WL and M 200c with a mean relation of [52]  \n〈 ln ( M WL M 0 )〉 = ln M WL 0 ( z ) + M WL M ln ( M 200c M 0 ) . (28)  \nHere ln M WL 0 is the logarithmic mass bias normalization and M WL M is the mass trend in this bias at a pivot mass M 0 = 2 × 10 14 h -1 M ⊙ . We assume a log-normal scatter of the true relations with a width described by  \nln σ ln WL 2 = ln σ 2 ln WL 0 ( z ) + σ 2 ln WL M ln ( M 200c M 0 ) , (29)  \nwhere ln σ 2 ln WL 0 is the normalization and σ 2 ln WL M is the mass trend of the scatter.  \nThe parameters of the above mean scaling relation and scatter are calibrated from simulations by extracting the  \nTABLE I. Normalization and uncertainties of the amplitude and scatter of the weak-lensing-mass-to-halo-mass relation derived from the simulations at redshifts z ∈ { 0 . 252 , 0 . 470 , 0 . 783 , 0 . 963 } .  \nweak-lensing inferred mass from hydrodynamical simulations and calculating the corresponding cluster mass from the matched N -body simulation at different redshifts [52]. The calibration results in a mean value and uncertainty obtained from the posterior for each of the above parameters.  \nIn this model, the logarithmic mass bias normalization, ln M WL 0 , and the normalization of the scatter, ln σ 2 ln WL 0 , are functions of redshift and calibrated from the simulations at four redshift values: z ∈ { 0 . 252 , 0 . 470 , 0 . 783 , 0 . 963 } . Therefore, we model the two parameters in the analysis as  \np = N (¯ p, (∆ p ) 2 ) = ¯ p ( z ) + ∆ p ( z ) N (0 , 1) , (30)  \nwhere ¯ p ( z ) is the mean value and ∆ p ( z ) is the uncertainty of the corresponding parameter p ( z ) at redshift z . We interpolate linearly to obtain the values for these parameters at any intermediate redshift. To accurately describe the uncertainty of the logarithmic mass bias ln M WL 0 ( z ), the uncertainty in this parameter, ∆ln M WL 0 ( z ), is modeled as a linear combination of two redshift-dependent functions [SB24a]:  \n∆ln M WL 0 ( z ) = ∆ 1 ln M WL 0 ( z ) + ∆ 2 ln M WL 0 ( z ) . (31)  \nThe values of the bias and scatter normalization parameters, Eqs. (28) and (29), as well as their uncertainties at the simulation redshifts used in this work are summarized in Tab. I. The uncertainties of these parameters include various elements such as uncertainties from baryonic effects, photoz calibration, miscentering and shear calibration [52]. The total uncertainty is primarily influenced by uncertainties in baryonic effects at low redshifts, while at high redshifts, the uncertainty in photoz calibration becomes dominant. Overall uncertainty from the weak-lensing model remains small across the calibrated redshift range, contributing to approximately 1 % of the total uncertainty (see Figure 10 in SB24a).', '2. HST weak-lensing model': 'A similar model is applied to the HST-39 dataset. The shear profiles from HST are modeled by the line of sight integral of an NFW profile with a concentration from  \nRef. [96]. From the NFW a weak-lensing mass M WL is calculated and related to the true halo mass with a mean relation [54].  \n⟨ ln M WL ⟩ = ln M WL 0 +ln M 200c . (32)  \nThe true relation scatters around the mean by a Gaussian distribution with width σ ln WL . The scatter σ ln WL accounts for all sources of uncertainties in the M WL -M 200c relation. Here each cluster has its own bias and scatter and associated uncertainties by calibrating Eq. (32) for each cluster individually. We refer the reader to the original works for a more detailed explanation of the cluster lensing model employed in the HST dataset [54-56, 59].', 'C. Multivariate observable-mass relation': 'To account for possible correlation among the three observables, unbiased tSZE detection significance ζ , intrinsic richness λ and weak-lensing mass M WL we employ the multivariant observable-mass relation from the work of SB24a. For this, the lognormal scatters of the observables, σ ln ζ , σ ln λ and σ ln WL are combined into a covariance matrix of the form  \nΣ =   σ 2 ln ζ ρ SZ , WL σ ln ζ σ ln WL ρ SZ ,λ σ ln ζ σ ln λ ρ SZ , WL σ ln ζ σ ln WL σ ln WL 2 ρ WL ,λ σ ln WL σ ln λ ρ SZ ,λ σ ln ζ σ ln λ ρ WL ,λ σ ln WL σ ln λ σ 2 ln λ   (33)  \nwhere ρ SZ ,λ , ρ SZ , WL and ρ WL ,λ are the correlation coefficients between ζ and M WL , ζ and λ , and λ and M WL respectively. The joint multi-observable-mass relation is then given by a multivariate Gaussian with correlation matrix Σ  \nP (   ln ζ ln M WL ln λ   ∣ ∣ M,z, p ) = N (   ⟨ ln ζ ⟩ ( M,z, p ) ⟨ ln M WL ⟩ ( M,z, p ) ⟨ ln λ ⟩ ( M,z, p )   , Σ ) . (34)', 'D. Likelihood and priors': 'The analysis relies on Bayesian statistics and we obtain cosmological and scaling relation parameters p using a cluster population model. The likelihood model employed in this analysis is based on the recent ΛCDM SPT × DES+HST analysis of SB24b, SB24a and was verified for an f ( R ) cosmology in our forecast work [SV24a]. Following SB24a and SB24b, we approximate the multiobservable cluster abundance likelihood with a Poisso-  \n,  \nnian distribution:  \nln L ( p ) = ∑ i ln ∫ ∞ ˆ λ min d ˆ λ d 3 N ( p ) d ˆ ζ d ˆ λ d z ∣ ∣ ∣ ˆ ζ i ,z i -∫ z max z min d z ∫ ∞ ˆ ζ min d ˆ ζ ∫ ∞ ˆ λ min d ˆ λ d 3 N ( p ) d ˆ ζ d ˆ λ d z + ∑ i ln d 4 N ( p ) d ˆ ζ d ˆ λ d g t d z ∣ ∣ ∣ ˆ ζ i , ˆ λ i , g t ,i ,z i ∫ ∞ ˆ λ min d ˆ λ d 3 N ( p ) d ˆ ζ d ˆ λ d z ∣ ∣ ∣ ˆ ζ i ,z i +const . , (35)  \nwhere both sums run over all clusters i . The differential cluster numbers d 3 N dobs in the above likelihood are the differential halo observable function (HOF) and is given in the ˆ ζ -ˆ λ -z space by  \nd 3 N ( p ) d ˆ ζ d ˆ λ d z = ∫ dΩ s ∫∫∫ d M d λ d ζ P ( ˆ ζ | ζ ) P ( ˆ λ | λ ) P ( ζ, λ | M,z, p ) d 2 N ( M,z, p ) d M d V d 2 V ( z, p ) d z dΩ s , (36)  \nand in the ˆ ζ -ˆ λ -g t -z by  \nd 4 N ( p ) d ˆ ζ d ˆ λ d g t d z = ∫ dΩ s ∫∫∫∫ d M d ζ d λ d M WL P ( g t | M WL , p ) P ( ˆ ζ | ζ ) P ( ˆ λ | λ ) P ( ζ, λ, M WL | M,z, p ) d 2 N ( M,z, p ) d M d V d 2 V ( z, p ) d z dΩ s . (37)  \nHere Ω s is the survey solid angle, the factors d 2 N ( M,z, p ) d M d z and d 2 V ( z, p ) d z dΩ s are the HMF and the differential volume element for the corresponding cosmology whereas P ( ˆ ζ | ζ ) and P ( ˆ λ | λ ) relate the observed quantity to the intrinsic one given by Eqs. (21) and (24) respectively. Moreover, P ( ζ, λ | M,z, p ) and P ( ζ, λ, M WL | M,z, p ) are obtained from the multivariant observable-mass relation Eq. (34), and P ( g t | M WL , p ) is given by the product of Gaussian probabilities in each radial bin i of the tangential shear profiles g t with the shape noise ∆ g t ,i  \nP ( g t | M WL , p ) = ∏ i ( √ 2 π ∆ g t ,i ) -1 exp [ -1 2 ( g t ,i -g t ,i ( M WL , p ) ∆ g t ,i ) 2 ] . (38)  \nNote that the first two terms of the total Poisson likelihood, Eq. (35), are independent of the weak-lensing data and thus are associated with the Poisson likelihood of the cluster sample in the ˆ ζ -z space with the condition ˆ λ > ˆ λ min ( z ). The last term in Eq. (35) includes the information of the mass calibration from the weak-lensing data and is therefore often called the mass-calibration likelihood. The likelihood is implemented in the CosmoSIS framework as a Python module [97] and we use the  \nnested sampling algorithm nautilus to run the MCMC chains [98].  \nGiven the large number of cosmological and nuisance parameters considered in the analysis, we combine the SPT cluster data set with the primary CMB data from Planck (Planck 2018 TT,TE,EE+lowE) [99] to break parameter degeneracies and recover meaningful constraints on f ( R ) gravity. This combination of data is sound, because the standard cosmological analysis from SB24b showed no statistically significant tension between the SPT-clusters × WL dataset and the Planck data.  \nWe emphasize that primary CMB data like those from Planck 2018 place only weak constraints on f ( R ) gravity of the order of log 10 | f R 0 | ≲ -3 [13, 100], and thus the constraints in log 10 | f R 0 | are primarily coming from the SPT cluster sample. However, CMB data are essential to constrain the remaining cosmological parameters such as Ω m h 2 . The Planck 2018 likelihood is implemented in CosmoSIS , and because one can assume that the cluster likelihood and the Planck 2018 likelihood are independent, we multiply the two likelihoods in a joint analysis. We account for the effect of f ( R ) gravity in the Planck 2018 likelihood by using the f ( R ) power spectrum computed with MGCAMB [101-104].  \nIn this analysis, we vary 23 nuisance parameters and eight cosmological parameters. All parameters with their priors are listed in Table II. For the standard cosmological parameters, we adopt uniform priors with ranges that are based on the Planck 2018 posteriors, since these parameters will be best constrained by the Planck 2018 dataset. For the f ( R ) gravity parameter log 10 | f R 0 | we apply a uniform prior U ( -7 , -3). Note that the GR limit at | f R 0 | = 0 cannot be reached when using a logarithmic prior and thus would introduce an infinitely large parameter volume below our lower bound when we calculate the upper limit of the f ( R ) parameter. To avoid the sensitivity of the upper limit of the f ( R ) parameter on the choice of the lower prior boundary, we transform the parameter space from logarithmic to linear. In this linear space, the parameter volume between 0 and 10 -7 is negligible. We account for the uncertainty in the weak-lensing model as described in Sec. IV B 1 by Gaussian priors on these parameters. To ensure a positive-definite covariance matrix of the multivariate observable-mass relation we assume priors on the correlation coefficients of U ( -0 . 5 , 0 . 5). No informative priors are applied for the tSZE and richness observable-mass relation parameters, and we adopt sufficiently wide uniform priors for these parameters. To account for systematic uncertainties in the HMF, we use Gaussian priors on the amplitude and slope of the HMF as defined in Eq. (20).', 'V. RESULTS': 'In this section, we present our constraints on the f ( R ) gravity model derived from the analysis of DES and HST weak lensing calibrated SPT clusters combined with  \nTABLE II. Parameters and priors of our f ( R ) gravity analysis of the DES and HST weak-lensing informed SPT cluster sample and Planck 2018. The prior on Ω ν h 2 corresponds to a prior on the sum of neutrino masses Σ m ν ∼ U (0 , 0 . 6) eV.  \nPlanck 2018 data. For our baseline analysis, we employ the FORGE-calibrated HMF model, because it is considered to be more accurate due to its empirical calibration with the FORGE numerical simulations. For comparison, we also include results obtained using the uncorrected semi-analytical HMF model. Note that all reported upper limits are provided at the 95 % credible level and uncertainties at 68 % credibility.', 'A. Comparison to Λ CDM results': 'We first compare our results to those from the ΛCDM analysis of SB24b. Because we report a tight upper limit on log 10 | f R 0 | (see next section) and thus the deviation from a ΛCDM cosmology is relatively small, we expect consistent results between the two analyses for all standard cosmological parameters. Figure 4 shows the posterior distribution of the cosmological parameters from both analyses. Overall, the standard cosmological parameters are in good agreement with the ΛCDM results. The only discernible parameter shifts are observed in Ω b h 2 and n s , which are both primarily constrained by the Planck 2018 data. However, the shift of the two parameters is within the 68 % credible contours and within the statistical uncertainties. An explanation for the shift can be given by the fact that Ω b h 2 and n s as well as log 10 | f R 0 | change the power spectrum on small scales, and in opposite directions.  \nA complete comparison of all parameters from the f ( R ) and ΛCDM analysis can be found in Appendix A, Fig. 6. We summarize in Tab. III the constraints on the ζ -mass and λ -mass relation parameters, as well as on the correlation coefficients for the f ( R ) and ΛCDM analyses. The results are in good agreement with ΛCDM, indicating no significant deviations in the observable-mass relations under the f ( R ) gravity model. Compared to the ΛCDM constraints, we obtain slightly weaker constraints for most of the observable-mass relation parameters (see the third column of Tab. III). The largest increase in the uncertainties is observed in the ζ -mass relation parameters, the amplitude mass trend of the λ -mass relation for the DES richness, and the scatter of the λ -mass relation for the WISE richness. An overall increase in the uncertainties is expected due to the additional degrees of freedom introduced by the f ( R ) gravity model.  \nWe observe a 1 . 5 σ shift in the correlation parameter between the tSZE detection significance and richness, ρ SZ,λ , compared to the ΛCDM results, which is attributed to different models for the scatter in richness assumed in the two analyses: in the ΛCDM analysis, a lognormal richness scatter with width λ -1 / 2 was used, which resulted in a negative correlation between the tSZE detection significance and richness. In this f ( R ) analysis, we use a Gaussian approximation of Poisson noise (see Eq (24)) and report a vanishing correlation between the two quantities. When using the same richness scatter model as in the ΛCDM analysis, we obtain the same results for this parameter. The same is true for the scatter parameter of the WISE richness σ ln λ, WISE .', 'B. f ( R ) gravity constraints': "Our f ( R ) analysis of the combination of weak-lensing mass calibrated SPT clusters and the Planck 2018 dataset results in the current tightest constraints available from clusters and CMB data and other probes of the large-  \nTABLE III. Constraints on the observable-mass relation parameters, the scatter correlation coefficients and cosmological parameters for the f ( R ) and ΛCDM analyses in the second and third column respectively. The last column shows the relative increase of the parameter uncertainties in the f ( R ) analysis with respect to ΛCDM. Note that all reported upper limits are provided at the 95 % credible level and uncertainties at 68 % credibility.  \nscale structure. With our baseline analysis using the FORGE-calibrated HMF, we obtain an upper bound on the f ( R ) parameter:  \nlog 10 | f R 0 | < -5 . 32 (95 % credible level) . (39)  \nThis result is consistent with a ΛCDM cosmology and excludes all f ( R ) parameter space that does not lead to substantial screening of halos. When applying the semianalytical HMF model, we achieve a 3 % tighter upper bound on log 10 | f R 0 | :  \nlog 10 | f R 0 | < -5 . 46 (95 % credible level) . (40)  \nThe slightly stronger constraint from the semi-analytical HMF model can be explained by the larger enhancement in the halo mass function compared to the FORGE simulations, see Sec. III B. Therefore, one expects a greater sensitivity to log 10 | f R 0 | , which leads to tighter constraints on the f ( R ) parameter. The remaining cosmological parameter posteriors are found to be consistent  \nFIG. 4. Posterior distribution of the cosmological parameters for our f ( R ) analysis (gray and red) and the ΛCDM analysis from SB24b (blue) of the DES and HST weak-lensing mass calibrated SPT clusters combined with Planck 2018 CMB data. Because the two f ( R ) HMF models give similar results, the gray contours from the semi-analytical HMF (see Sec. III A) are hidden by the red ones from the FORGE calibrated HMF (see Sec. III B). As expected, the results from the analysis show consistent results in the nonf ( R ) cosmological parameters.  \n<!-- image -->  \nΩ  \nb  \nh  \nwith those in the ΛCDM analysis. This is because we use Planck 2018 data in our analysis, which tightly constrains the ΛCDM cosmological parameters and thus eliminates potential degeneracies. Moreover, degeneracies between log 10 | f R 0 | and other cosmological parameters that have previously been found are more pronounced for higher  \nvalues of log 10 | f R 0 | , which are excluded by our dataset [30, 70, 105].  \nFor comparison, the best previous constraints on f ( R ) gravity from clusters are presented in Ref. [28]. The authors obtained log 10 | f R 0 | < -4 . 79 using clusters from ROSAT and the Massive Cluster Survey, combined with  \nprimary CMB, SN, and BAO data. Our analysis improves upon this result by a factor of 3.4 in f R 0 without using any information from SNe or BAO. The improved constraints are due to the large cluster sample and the weak lensing mass calibration dataset that we use in this analysis, with 1 , 005 clusters compared to 224 clusters in the analysis of Ref. [28]. Recent results from the eROSITA cluster analysis reported an upper limit of log 10 | f R 0 | < -4 . 12, considering the neutrino mass as a free parameter [11]. Although this constraint is significantly weaker than ours, it is important to note that their analysis was based on clusters alone. eROSITA can place meaningful constraints from the clusters alone, because of the larger weak-lensing calibrated cluster sample of 5 , 259 clusters and lower redshift range for which the effect of f ( R ) gravity is larger. On the other hand, they employed a different HMF based on the model of Ref. [30], which predicts a slightly smaller enhancement of the f ( R ) HMF. For a fairer comparison with the eROSITA results, we also perform our analysis using the HMF of Ref. [30] but still including Planck 2018 data. We obtain log 10 | f R 0 | < -5 . 11 (95 % credible level), which is an order of magnitude better than the results from Ref. [11] and 62 % weaker in f R 0 than our baseline result. These different results show that a reliable f ( R ) HMF model is needed to obtain accurate constraints on this modified gravity model. Note that different work apply different prior on log 10 | f R 0 | which affects this parameter's upper limit due to the infinite volume below the lower prior. Therefore, a comparison is always affected by the prior choice when using a uniform prior in log 10 | f R 0 | .  \nThe strongest constraints from large-scale probes of f ( R ) gravity are derived from weak-lensing peak abundance using weak-lensing data from CFHTLenS [12]. In this study, the authors found an upper bound of log 10 | f R 0 | < -5 . 16 (95 % credible level) with priors from the Planck 15 analysis on Ω m and A s . This result is comparable, but slightly weaker than the constraints presented here.  \nOther recent cosmological constraints on f ( R ) gravity have been derived from the cross-correlation of galaxies from BOSS combined with primary CMB and lensing data, yielding log 10 | f R 0 | ≤ -4 . 61 [13] and a combination of galaxy weak-lensing shear from the Canada-France-Hawaii Telescope Lensing Survey, CMB, SN and BAO was used to obtain log 10 | f R 0 | ≤ -4 . 50 [12]. Figure 5 shows the comparison of the constraints from the different cosmological probes discussed in this section.  \nModified gravity models such as f ( R ) gravity generally enhance structure formation and thus lead to larger values of σ 8 . Therefore, they do not provide a solution to the S 8 tension. Additionally, galaxy clusters are not sensitive to H 0 , so our analysis does not offer any insights into the H 0 tension.  \nFIG. 5. Comparison of log 10 | f R 0 | constraints from SPTclusters + CMB with the FORGE calibrated HMF (top bar) and the semi-analytical HMF model (second bar) with other recent results on cosmological scales. All limits are given at the 95 % credible level.  \n<!-- image -->  \n|  \n|", 'VI. SUMMARY': 'This work presents constraints on f ( R ) gravity derived from the DES and HST weak lensing informed SPT cluster abundance combined with primary CMB data from Planck 2018. We use a sample of 1 , 005 galaxy clusters selected from the SPT-SZ, SPTpolECS, and SPT500d [34, 35, 46, 47] surveys with redshifts z > 0 . 25. 688 of these clusters have weak-lensing information from DES [SB24b, SB24a] and 39 from HST [54, 55, 59]. Our analysis framework is based on the methodology established by SB24a and in the recent f ( R ) gravity forecast for upcoming Stage-III and -IV surveys [SV24a].  \nf ( R ) gravity alters gravity and leads, compared to GR, to a scale-dependent enhancement of structure formation and thus modifies the HMF. This enhancement in the abundance of massive galaxy clusters makes cluster samples, such as those from the SPT surveys, powerful probes for testing modified gravitational models like f ( R ) gravity.  \nTo capture the effects of f ( R ) gravity on the HMF, we employ a semi-analytical approach for calculating the mass-dependent spherical collapse threshold δ crit [44]. The f ( R ) HMF is then given by the GR HMF scaled by an enhancement factor, which includes the massdependent spherical collapse threshold. We emphasize that this model is designed to also capture the nontrivial screening effects in this modified gravity model, which play a major role at small modified gravity amplitudes that our cluster sample is able to probe. We compare the predictions of the HMF from this semi-analytical approach with those from the FORGE simulations [45], and we find a discrepancy between the two that depends on cosmology, redshift and cluster mass. We use the simulations to calibrate the semi-analytical model to obtain a more robust HMF, while still allowing for an analysis within a broader parameter range than the simulations allow.  \nIn this analysis, we neglect the f ( R ) gravity effects on the gravitational lensing potential, which are subdomi-  \nnant compared to the weak-lensing mass calibration uncertainties derived from GR simulations [52]. Furthermore, modifications to the observable-mass relations and the halo profiles are minimal, keeping the weak-lensing parameters similar to those in GR within their uncertainties [53, 90, 91].  \nWe achieve consistent results in all parameters compared to the ΛCDM analysis presented in SB24b. This is reassuring and expected, since our constraints on log 10 | f R 0 | do not indicate a preference for a deviation from the ΛCDM cosmology.  \nWe report an upper bound of log 10 | f R 0 | < -5 . 32 at the 95% credible level with the HMF calibrated by the FORGE simulations. A slightly tighter constraint of log 10 | f R 0 | < -5 . 46 is obtained when using the semianalytical HMF. The difference in the constraints shows the necessity for a reliable f ( R ) HMF to place accurate constraints on f ( R ) gravity. The constraints reported here are the tightest constraints on f ( R ) gravity from clusters and on cosmological scales published to date.  \nUpcoming Stage-III and Stage-IV surveys, such as SPT-3G [73] and CMB-S4 [74] or the Simons Observatory [106], will provide significantly larger cluster samples, which cover a broader redshift range [107]. In combination with next-gravitational lensing data like from the Euclid [75, 76] satellite or the Vera C. Rubin observatory [77, 78], these cluster data sets will lead to improved constraints on f ( R ) gravity [SV24a]', 'ACKNOWLEDGMENTS': "This research was supported by 1) the Excellence Cluster ORIGINS, which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC2094-390783311, by 2) the Max Planck Society Faculty Fellowship program at MPE, and by 3) the LudwigMaximilians-Universitat in Munich.  \nThe FORGE simulations and analyses of this project made use of the DiRAC@Durham facility managed by the Institute for Computational Cosmology (ICC) on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National eInfrastructure.  \nBL is supported by STFC via Consolidated Grant ST/X001075/1  \nThe Innsbruck authors acknowledge support provided by the Austrian Research Promotion Agency (FFG) and the Federal Ministry of the Republic of Austria for Climate Action, Environment, Mobility, Innovation and Technology (BMK) via the Austrian Space Applications Programme with grant numbers 899537, 900565, and 911971.  \nThe South Pole Telescope program is supported by the National Science Foundation (NSF) through the Grant No. OPP-1852617 and 2332483. Partial support is also provided by the Kavli Institute of Cosmological Physics at the University of Chicago. Work at Argonne National Lab is supported by UChicago Argonne LLC, Operator of Argonne National Laboratory (Argonne). Argonne, a U.S. Department of Energy Office of Science Laboratory, is operated under contract no. DE-AC02-06CH11357  \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 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¸c˜ao Carlos Chagas Filho de Amparo 'a Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Cient'ıfico e Tecnol'ogico and the Minist'erio da Ciˆencia, Tecnologia e Inova¸c˜ao, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey.  \nThe Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energ'eticas, Medioambientales y Tecnol'ogicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenossische Technische Hochschule (ETH) Zurich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ci'encies de l'Espai (IEEC/CSIC), the Institut de F'ısica d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians-Universitat Munchen and the associated Excellence Cluster Origins, the University of Michigan, NSF's NOIRLab, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.  \nBased in part on observations at Cerro Tololo InterAmerican Observatory at NSF's NOIRLab (NOIRLab Prop. ID 2012B-0001; PI: J. Frieman), which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.  \nThe DES data management system is supported by the National Science Foundation under Grant Numbers AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MICINN under grants ESP2017-89838, PGC2018094773, PGC2018-102021, SEV-2016-0588, SEV-2016-  \n0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. Research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Program (FP7/20072013) including ERC grant agreements 240672, 291329, and 306478. We acknowledge support from the Brazilian Instituto Nacional de Ciˆencia e Tecnologia (INCT) do e-Universo (CNPq grant 465376/2014-2).  \nThis work is based on observations made with the NASA/ESA Hubble Space Telescope , using imaging data from the SPT follow-up GO programs 12246 (PI: C. Stubbs), 12477 (PI: F. W. High), 13412 (PI: T. Schrabback), 14252 (PI: V. Strazzullo), 14352 (PI: J. Hlavacek-Larrondo), and 14677 (PI: T. Schrabback). STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. It is also based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programs 086.A-0741 (PI: Bazin), 088.A0796 (PI: Bazin), 088.A-0889 (PI: Mohr), 089.A-0824 (PI: Mohr), 0100.A-0204 (PI: Schrabback), 0100.A-0217 (PI: Hern'andez-Mart'ın), 0101.A-0694 (PI: Zohren), and  \n- [1] S. Perlmutter, G. Aldering, G. Goldhaber, R. A. Knop, P. Nugent, P. G. Castro, S. Deustua, S. Fabbro, A. Goobar, D. E. Groom, et al., Astrophys. J. 517 , 565 (1999), arXiv:astro-ph/9812133.\n- [2] A. G. Riess, A. V. Filippenko, P. Challis, A. Clocchiatti, A. Diercks, P. M. Garnavich, R. L. Gilliland, C. J. Hogan, S. Jha, R. P. Kirshner, et al., ApJ 116 , 1009 (1998), arXiv:astro-ph/9805201.\n- [3] A. Joyce, L. Lombriser, and F. Schmidt, Annual Review of Nuclear and Particle Science 66 , 95 (2016), arXiv:1601.06133.\n- [4] K. Koyama, International Journal of Modern Physics D 27 , 1848001 (2018).\n- [5] T. Baker, A. Barreira, H. Desmond, P. Ferreira, B. Jain, K. Koyama, B. Li, L. Lombriser, A. Nicola, J. Sakstein, et al., arXiv e-prints arXiv:1908.03430 (2019), arXiv:1908.03430.\n- [6] H. A. Buchdahl, MNRAS 150 , 1 (1970).\n- [7] W. Hu and I. Sawicki, Phys. Rev. D 76 , 064004 (2007), URL https://link.aps.org/doi/10.1103/PhysRevD. 76.064004 .\n- [8] A. P. Naik, E. Puchwein, A.-C. Davis, D. Sijacki, and H. Desmond, MNRAS 489 , 771 (2019), arXiv:1905.13330.\n- [9] H. Desmond and P. G. Ferreira, Phys. Rev. D 102 , 104060 (2020), URL https://link.aps.org/doi/10. 1103/PhysRevD.102.104060 .\n- [10] M. Cataneo, D. Rapetti, L. Lombriser, and B. Li, JCAP 2016 , 024 (2016), arXiv:1607.08788.\n- [11] E. Artis, V. Ghirardini, E. Bulbul, S. Grandis, C. Garrel, N. Clerc, R. Seppi, J. Comparat, M. Cataneo, Y. E. Bahar, et al., arXiv e-prints arXiv:2402.08459 (2024), arXiv:2402.08459.  \n0102.A-0189 (PI: Zohren). It is also based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnolog'ıa e Innovaci'on Productiva (Argentina), Minist'erio da Ciˆencia, Tecnologia e Inova¸c˜ao (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea), under programs 2014B-0338 and 2016B-0176 (PI: B. Benson).  \nThis manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics.", 'Appendix A: Full triangle plot': "Figure 6 shows the posterior distribution for all parameters varied in the analysis, both for the f ( R ) and ΛCDM cases of the same data SB24b.  \n- [12] X. Liu, B. Li, G.-B. Zhao, M.-C. Chiu, W. Fang, C. Pan, Q. Wang, W. Du, S. Yuan, L. Fu, et al., Phys. Rev. Lett. 117 , 051101 (2016), arXiv:1607.00184.\n- [13] R. Kou, C. Murray, and J. G. Bartlett, arXiv e-prints arXiv:2311.09936 (2023), arXiv:2311.09936.\n- [14] B. Hu, M. Raveri, M. Rizzato, and A. Silvestri, Monthly Notices of the Royal Astronomical Society 459 , 3880 (2016), ISSN 0035-8711, arXiv:https://academic.oup.com/mnras/articlepdf/459/4/3880/8191276/stw775.pdf, URL https://doi.org/10.1093/mnras/stw775 .\n- [15] L. Wang and P. J. Steinhardt, Astrophys. J. 508 , 483 (1998), arXiv:astro-ph/9804015.\n- [16] Z. Haiman, J. J. Mohr, and G. P. Holder, Astrophys. J. 553 , 545 (2001), arXiv:astro-ph/0002336.\n- [17] A. Vikhlinin, A. V. Kravtsov, R. A. Burenin, H. Ebeling, W. R. Forman, A. Hornstrup, C. Jones, S. S. Murray, D. Nagai, H. Quintana, et al., Astrophys. J. 692 , 1060 (2009), arXiv:0812.2720.\n- [18] B. A. Benson, T. de Haan, J. P. Dudley, C. L. Reichardt, K. A. Aird, K. Andersson, R. Armstrong, M. L. N. Ashby, M. Bautz, M. Bayliss, et al., Astrophys. J. 763 , 147 (2013), arXiv:1112.5435.\n- [19] S. Bocquet, A. Saro, J. J. Mohr, K. A. Aird, M. L. N. Ashby, M. Bautz, M. Bayliss, G. Bazin, B. A. Benson, L. E. Bleem, et al., Astrophys. J. 799 , 214 (2015), arXiv:1407.2942.\n- [20] T. de Haan, B. A. Benson, L. E. Bleem, S. W. Allen, D. E. Applegate, M. L. N. Ashby, M. Bautz, M. Bayliss, S. Bocquet, M. Brodwin, et al., Astrophys. J. 832 , 95 (2016), arXiv:1603.06522.\n- [21] S. Bocquet, J. P. Dietrich, T. Schrabback, L. E. Bleem, M. Klein, S. W. Allen, D. E. Applegate, M. L. N. Ashby,  \nFIG. 6. Posterior distribution for all parameters in the f ( R ) and ΛCDM analyses.  \n<!-- image -->  \n- M. Bautz, M. Bayliss, et al., Astrophys. J. 878 , 55 (2019), arXiv:1812.01679.\n- [22] T. M. C. Abbott, M. Aguena, A. Alarcon, S. Allam, S. Allen, J. Annis, S. Avila, D. Bacon, K. Bechtol, A. Bermeo, et al., Phys. Rev. D 102 , 023509 (2020), arXiv:2002.11124.\n- [23] I. N. Chiu, M. Klein, J. Mohr, and S. Bocquet, MNRAS 522 , 1601 (2023), arXiv:2207.12429.\n- [24] S. Bocquet, S. Grandis, L. E. Bleem, M. Klein, J. J. Mohr, T. Schrabback, T. M. C. Abbott, P. A. R.\n- Ade, M. Aguena, A. Alarcon, et al., PRD submitted arXiv:2401.02075 (2024), arXiv:2401.02075.\n- [25] V. Ghirardini, E. Bulbul, E. Artis, N. Clerc, C. Garrel, S. Grandis, M. Kluge, A. Liu, Y. E. Bahar, F. Balzer, et al., arXiv e-prints arXiv:2402.08458 (2024), arXiv:2402.08458.\n- [26] F. Schmidt, A. Vikhlinin, and W. Hu, Phys. Rev. D 80 , 083505 (2009), arXiv:0908.2457.\n- [27] L. Lombriser, A. Slosar, U. Seljak, and W. Hu, Phys. Rev. D 85 , 124038 (2012), arXiv:1003.3009.\n- [28] M. Cataneo, D. Rapetti, F. Schmidt, A. B. Mantz, S. W. Allen, D. E. Applegate, P. L. Kelly, A. von der Linden, and R. G. Morris, Phys. Rev. D 92 , 044009 (2015), arXiv:1412.0133.\n- [29] S. Peirone, M. Raveri, M. Viel, S. Borgani, and S. Ansoldi, Phys. Rev. D 95 , 023521 (2017), URL https: //link.aps.org/doi/10.1103/PhysRevD.95.023521 .\n- [30] S. Hagstotz, M. Costanzi, M. Baldi, and J. Weller, MNRAS 486 , 3927 (2019), arXiv:1806.07400.\n- [31] J. E. Carlstrom, P. A. R. Ade, K. A. Aird, B. A. Benson, L. E. Bleem, S. Busetti, C. L. Chang, E. Chauvin, H. M. Cho, T. M. Crawford, et al., PASP 123 , 568 (2011), arXiv:0907.4445.\n- [32] R. A. Sunyaev and Y. B. Zeldovich, Comments on Astrophysics and Space Physics 4 , 173 (1972).\n- [33] M. Klein, J. J. Mohr, S. Desai, H. Israel, S. Allam, A. Benoit-L'evy, D. Brooks, E. Buckley-Geer, A. Carnero Rosell, M. Carrasco Kind, et al., MNRAS 474 , 3324 (2018), arXiv:1706.06577.\n- [34] M. Klein, J. J. Mohr, S. Bocquet, M. Aguena, S. W. Allen, O. Alves, B. Ansarinejad, M. L. N. Ashby, D. Bacon, M. Bayliss, et al., MNRAS 531 , 3973 (2024).\n- [35] L. E. Bleem, M. Klein, T. M. C. Abbot, P. A. R. Ade, M. Aguena, O. Alves, A. J. Anderson, F. AndradeOliveira, B. Ansarinejad, M. Archipley, et al., The Open Journal of Astrophysics 7 , 13 (2024), arXiv:2311.07512.\n- [36] B. Flaugher, H. T. Diehl, K. Honscheid, T. M. C. Abbott, O. Alvarez, R. Angstadt, J. T. Annis, M. Antonik, O. Ballester, L. Beaufore, et al., ApJ 150 , 150 (2015), arXiv:1504.02900.\n- [37] Dark Energy Survey Collaboration, T. Abbott, F. B. Abdalla, J. Aleksi'c, S. Allam, A. Amara, D. Bacon, E. Balbinot, M. Banerji, K. Bechtol, et al., MNRAS 460 , 1270 (2016), arXiv:1601.00329.\n- [38] T. M. C. Abbott and et al., ApJ 239 , 18 (2018), arXiv:1801.03181.\n- [39] E. L. Wright, P. R. M. Eisenhardt, A. K. Mainzer, M. E. Ressler, R. M. Cutri, T. Jarrett, J. D. Kirkpatrick, D. Padgett, R. S. McMillan, M. Skrutskie, et al., ApJ 140 , 1868 (2010), arXiv:1008.0031.\n- [40] S. Bocquet, S. Grandis, L. E. Bleem, M. Klein, J. J. Mohr, M. Aguena, A. Alarcon, S. Allam, S. W. Allen, O. Alves, et al., PRD submitted (2023), arXiv:2310.12213.\n- [41] S. M. L. Vogt, S. Bocquet, C. T. Davies, J. J. Mohr, and F. Schmidt, Phys. Rev. D 109 , 123503 (2024), arXiv:2401.09959.\n- [42] S. Shandera, A. Mantz, D. Rapetti, and S. W. Allen, JCAP 2013 , 004 (2013), arXiv:1304.1216.\n- [43] B. Li and G. Efstathiou, MNRAS 421 , 1431 (2012), arXiv:1110.6440.\n- [44] L. Lombriser, B. Li, K. Koyama, and G.-B. Zhao, Phys. Rev. D 87 , 123511 (2013), arXiv:1304.6395.\n- [45] C. Arnold, B. Li, B. Giblin, J. Harnois-D'eraps, and Y.C. Cai, MNRAS 515 , 4161 (2022), arXiv:2109.04984.\n- [46] L. E. Bleem, B. Stalder, T. de Haan, K. A. Aird, S. W. Allen, D. E. Applegate, M. L. N. Ashby, M. Bautz, M. Bayliss, B. A. Benson, et al., ApJ 216 , 27 (2015), arXiv:1409.0850.\n- [47] L. E. Bleem, S. Bocquet, B. Stalder, M. D. Gladders, P. A. R. Ade, S. W. Allen, A. J. Anderson, J. Annis, M. L. N. Ashby, J. E. Austermann, et al., ApJ 247 , 25 (2020), arXiv:1910.04121.  \n- [48] M. Gatti, G. Giannini, G. M. Bernstein, A. Alarcon, J. Myles, A. Amon, R. Cawthon, M. Troxel, J. DeRose, S. Everett, et al., MNRAS 510 , 1223 (2022), arXiv:2012.08569.\n- [49] E. Huff and R. Mandelbaum (2017), arXiv:1702.02600.\n- [50] E. S. Sheldon and E. M. Huff, Astrophys. J. 841 , 24 (2017), arXiv:1702.02601.\n- [51] DES Collaboration, T. M. C. Abbott, M. Aguena, A. Alarcon, S. Allam, O. Alves, A. Amon, F. AndradeOliveira, J. Annis, S. Avila, et al., Phys. Rev. D 105 , 023520 (2022), arXiv:2105.13549.\n- [52] S. Grandis, S. Bocquet, J. J. Mohr, M. Klein, and K. Dolag, MNRAS 507 , 5671 (2021), arXiv:2103.16212.\n- [53] C.-Z. Ruan, C. Cuesta-Lazaro, A. Eggemeier, B. Li, C. M. Baugh, C. Arnold, S. Bose, C. Hern'andezAguayo, P. Zarrouk, and C. T. Davies, MNRAS 527 , 2490 (2024), arXiv:2301.02970.\n- [54] T. Schrabback, D. Applegate, J. P. Dietrich, H. Hoekstra, S. Bocquet, A. H. Gonzalez, A. von der Linden, M. McDonald, C. B. Morrison, S. F. Raihan, et al., MNRAS 474 , 2635 (2018), arXiv:1611.03866.\n- [55] T. Schrabback, S. Bocquet, M. Sommer, H. Zohren, J. L. van den Busch, B. Hern'andez-Mart'ın, H. Hoekstra, S. F. Raihan, M. Schirmer, D. Applegate, et al., MNRAS 505 , 3923 (2021), arXiv:2009.07591.\n- [56] H. Zohren, T. Schrabback, S. Bocquet, M. Sommer, F. Raihan, B. Hern'andez-Mart'ın, O. Marggraf, B. Ansarinejad, M. B. Bayliss, L. E. Bleem, et al., AnA 668 , A18 (2022), arXiv:2208.10232.\n- [57] S. F. Raihan, T. Schrabback, H. Hildebrandt, D. Applegate, and G. Mahler, MNRAS 497 , 1404 (2020), arXiv:2007.01211.\n- [58] B. Hern'andez-Mart'ın, T. Schrabback, H. Hoekstra, N. Martinet, J. Hlavacek-Larrondo, L. E. Bleem, M. D. Gladders, B. Stalder, A. A. Stark, and M. Bayliss, AnA 640 , A117 (2020), arXiv:2007.00386.\n- [59] M. W. Sommer, T. Schrabback, D. E. Applegate, S. Hilbert, B. Ansarinejad, B. Floyd, and S. Grandis, MNRAS 509 , 1127 (2022), arXiv:2105.08027.\n- [60] J. Khoury and A. Weltman, Phys. Rev. D 69 , 044026 (2004), arXiv:astro-ph/0309411.\n- [61] C. Burrage and J. Sakstein, Living Reviews in Relativity 21 , 1 (2018), arXiv:1709.09071.\n- [62] S. Ferraro, F. Schmidt, and W. Hu, Phys. Rev. D 83 , 063503 (2011), arXiv:1011.0992.\n- [63] J. Tinker, A. V. Kravtsov, A. Klypin, K. Abazajian, M. Warren, G. Yepes, S. Gottlober, and D. E. Holz, Astrophys. J. 688 , 709 (2008), arXiv:0803.2706.\n- [64] R. K. Sheth and G. Tormen, MNRAS 308 , 119 (1999), arXiv:astro-ph/9901122.\n- [65] G. Despali, C. Giocoli, R. E. Angulo, G. Tormen, R. K. Sheth, G. Baso, and L. Moscardini, MNRAS 456 , 2486 (2016), arXiv:1507.05627.\n- [66] L. Lombriser, K. Koyama, and B. Li, JCAP 2014 , 021 (2014), arXiv:1312.1292.\n- [67] K. Ichiki and M. Takada, Phys. Rev. D 85 , 063521 (2012), URL https://link.aps.org/doi/10. 1103/PhysRevD.85.063521 .\n- [68] M. Costanzi, F. Villaescusa-Navarro, M. Viel, J.-Q. Xia, S. Borgani, E. Castorina, and E. Sefusatti, JCAP 2013 , 012 (2013), arXiv:1311.1514.\n- [69] H. Motohashi, A. A. Starobinsky, and J. Yokoyama, Phys. Rev. Lett. 110 , 121302 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.  \n110.121302 .  \n- [70] M. Baldi, F. Villaescusa-Navarro, M. Viel, E. Puchwein, V. Springel, and L. Moscardini, MNRAS 440 , 75 (2014), arXiv:1311.2588.\n- [71] B. S. Wright, K. Koyama, H. A. Winther, and G.-B. Zhao, JCAP 2019 , 040 (2019), arXiv:1902.10692.\n- [72] M. Costanzi, E. Rozo, M. Simet, Y. Zhang, A. E. Evrard, A. Mantz, E. S. Rykoff, T. Jeltema, D. Gruen, S. Allen, et al., MNRAS 488 , 4779 (2019), arXiv:1810.09456.\n- [73] B. A. Benson, P. A. R. Ade, Z. Ahmed, S. W. Allen, K. Arnold, J. E. Austermann, A. N. Bender, L. E. Bleem, J. E. Carlstrom, C. L. Chang, et al., 9153 , 91531P (2014), arXiv:1407.2973.\n- [74] K. Abazajian, G. Addison, P. Adshead, Z. Ahmed, S. W. Allen, D. Alonso, M. Alvarez, A. Anderson, K. S. Arnold, C. Baccigalupi, et al., arXiv e-prints arXiv:1907.04473 (2019), arXiv:1907.04473.\n- [75] R. Laureijs, J. Amiaux, S. Arduini, J. L. Augu'eres, J. Brinchmann, R. Cole, M. Cropper, C. Dabin, L. Duvet, A. Ealet, et al., arXiv e-prints arXiv:1110.3193 (2011), arXiv:1110.3193.\n- [76] Euclid Collaboration, R. Scaramella, J. Amiaux, Y. Mellier, C. Burigana, C. S. Carvalho, J. C. Cuillandre, A. Da Silva, A. Derosa, J. Dinis, et al., AnA 662 , A112 (2022), arXiv:2108.01201.\n- [77] Z. Ivezic, T. Axelrod, W. N. Brandt, D. L. Burke, C. F. Claver, A. Connolly, K. H. Cook, P. Gee, D. K. Gilmore, S. H. Jacoby, et al., Serbian Astronomical Journal 176 , 1 (2008).\n- [78] The LSST Dark Energy Science Collaboration, R. Mandelbaum, T. Eifler, R. Hloˇzek, T. Collett, E. Gawiser, D. Scolnic, D. Alonso, H. Awan, R. Biswas, et al., arXiv e-prints arXiv:1809.01669 (2018), arXiv:1809.01669.\n- [79] N. Kaiser, MNRAS 222 , 323 (1986).\n- [80] R. E. Angulo, V. Springel, S. D. M. White, A. Jenkins, C. M. Baugh, and C. S. Frenk, MNRAS 426 , 2046 (2012), arXiv:1203.3216.\n- [81] H. Hoekstra, R. Herbonnet, A. Muzzin, A. Babul, A. Mahdavi, M. Viola, and M. Cacciato, MNRAS 449 , 685 (2015), arXiv:1502.01883.\n- [82] A. B. Mantz, S. W. Allen, R. G. Morris, A. von der Linden, D. E. Applegate, P. L. Kelly, D. L. Burke, D. Donovan, and H. Ebeling, MNRAS 463 , 3582 (2016), arXiv:1606.03407.\n- [83] H. Miyatake, N. Battaglia, M. Hilton, E. Medezinski, A. J. Nishizawa, S. More, S. Aiola, N. Bahcall, J. R. Bond, E. Calabrese, et al., Astrophys. J. 875 , 63 (2019), arXiv:1804.05873.\n- [84] F. Bellagamba, M. Sereno, M. Roncarelli, M. Maturi, M. Radovich, S. Bardelli, E. Puddu, L. Moscardini, F. Getman, H. Hildebrandt, et al., MNRAS 484 , 1598 (2019), arXiv:1810.02827.\n- [85] I. N. Chiu, V. Ghirardini, A. Liu, S. Grandis, E. Bulbul, Y. E. Bahar, J. Comparat, S. Bocquet, N. Clerc, M. Klein, et al., AnA 661 , A11 (2022), arXiv:2107.05652.\n- [86] S. Grandis, V. Ghirardini, S. Bocquet, C. Garrel, J. J. Mohr, A. Liu, M. Kluge, L. Kimmig, T. H. Reiprich,  \n- A. Alarcon, et al., arXiv e-prints arXiv:2402.08455 (2024), arXiv:2402.08455.\n- [87] K. Vanderlinde, T. M. Crawford, T. de Haan, J. P. Dudley, L. Shaw, P. A. R. Ade, K. A. Aird, B. A. Benson, L. E. Bleem, M. Brodwin, et al., Astrophys. J. 722 , 1180 (2010), arXiv:1003.0003.\n- [88] T. P. Sotiriou and V. Faraoni, Reviews of Modern Physics 82 , 451 (2010), arXiv:0805.1726.\n- [89] P. Zhang, Phys. Rev. D 76 , 024007 (2007), arXiv:astroph/0701662.\n- [90] F. Schmidt, Phys. Rev. D 81 , 103002 (2010), arXiv:1003.0409.\n- [91] M. A. Mitchell, J.-h. He, C. Arnold, and B. Li, MNRAS 477 , 1133 (2018), arXiv:1802.02165.\n- [92] J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 490 , 493 (1997), arXiv:astro-ph/9611107.\n- [93] M. Bartelmann and M. Maturi, Scholarpedia 12 , 32440 (2017), arXiv:1612.06535.\n- [94] M. R. Becker and A. V. Kravtsov, Astrophys. J. 740 , 25 (2011), arXiv:1011.1681.\n- [95] M. Oguri and T. Hamana, MNRAS 414 , 1851 (2011), arXiv:1101.0650.\n- [96] B. Diemer and M. Joyce, Astrophys. J. 871 , 168 (2019), arXiv:1809.07326.\n- [97] J. Zuntz, M. Paterno, E. Jennings, D. Rudd, A. Manzotti, S. Dodelson, S. Bridle, S. Sehrish, and J. Kowalkowski, Astronomy and Computing 12 , 45 (2015), arXiv:1409.3409.\n- [98] J. U. Lange, MNRAS 525 , 3181 (2023), arXiv:2306.16923.\n- [99] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, et al., AnA 641 , A6 (2020), arXiv:1807.06209.\n- [100] Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, N. Bartolo, et al., AnA 594 , A14 (2016), arXiv:1502.01590.\n- [101] G.-B. Zhao, L. Pogosian, A. Silvestri, and J. Zylberberg, Phys. Rev. D 79 , 083513 (2009), arXiv:0809.3791.\n- [102] A. Hojjati, L. Pogosian, and G.-B. Zhao, JCAP 2011 , 005 (2011), arXiv:1106.4543.\n- [103] A. Zucca, L. Pogosian, A. Silvestri, and G. B. Zhao, JCAP 2019 , 001 (2019), arXiv:1901.05956.\n- [104] Z. Wang, S. H. Mirpoorian, L. Pogosian, A. Silvestri, and G.-B. Zhao, JCAP 2023 , 038 (2023), arXiv:2305.05667.\n- [105] J. Harnois-D'eraps, C. Hernandez-Aguayo, C. CuestaLazaro, C. Arnold, B. Li, C. T. Davies, and Y.-C. Cai, MNRAS 525 , 6336 (2023), arXiv:2211.05779.\n- [106] P. Ade, J. Aguirre, Z. Ahmed, S. Aiola, A. Ali, D. Alonso, M. A. Alvarez, K. Arnold, P. Ashton, J. Austermann, et al., JCAP 2019 , 056 (2019), arXiv:1808.07445.\n- [107] S. Raghunathan, N. Whitehorn, M. A. Alvarez, H. Aung, N. Battaglia, G. P. Holder, D. Nagai, E. Pierpaoli, C. L. Reichardt, and J. D. Vieira, Astrophys. J. 926 , 172 (2022), arXiv:2107.10250."}