# Constraining $f(R)$ gravity using future galaxy cluster abundance and weak-
lensing mass calibration datasets

Sophie M. L. Vogt<EMAIL_ADDRESS>University Observatory, Faculty of
Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, 81679 Munich, Germany
Excellence Cluster ORIGINS, Boltzmannstr. 2, 85748 Garching, Germany  Max
Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, 85748 Garching,
Germany  Sebastian Bocquet  University Observatory, Faculty of Physics,
Ludwig-Maximilians-Universität, Scheinerstr. 1, 81679 Munich, Germany
Christopher T. Davies  University Observatory, Faculty of Physics, Ludwig-
Maximilians-Universität, Scheinerstr. 1, 81679 Munich, Germany  Joseph J. Mohr
University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität,
Scheinerstr. 1, 81679 Munich, Germany  Max Planck Institute for
extraterrestrial Physics, Giessenbach-Str. 2, 85748 Garching, Germany  Fabian
Schmidt  Max Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1,
85748 Garching, Germany  Excellence Cluster ORIGINS, Boltzmannstr. 2, 85748
Garching, Germany

###### Abstract

We present forecasts for constraints on the Hu & Sawicki $f(R)$ modified
gravity model using realistic mock data representative of future cluster and
weak lensing surveys. We create mock thermal Sunyaev-Zel’dovich effect
selected cluster samples for SPT-3G and CMB-S4 and the corresponding weak
gravitational lensing data from next-generation weak-lensing (ngWL) surveys
like Euclid and Rubin. We employ a state-of-the-art Bayesian likelihood
approach that includes all observational effects and systematic uncertainties
to obtain constraints on the $f(R)$ gravity parameter $\log_{10}|f_{R0}|$. In
this analysis we vary the cosmological parameters $[\Omega_{\rm
m},\Omega_{\nu}h^{2},h^{2},A_{s},n_{s},\log_{10}|f_{R0}|]$, which allows us to
account for possible degeneracies between cosmological parameters and $f(R)$
modified gravity. The analysis accounts for $f(R)$ gravity via its effect on
the halo mass function which is enhanced on cluster mass scales compared to
the expectations within general relativity (GR). Assuming a fiducial GR model,
the upcoming cluster dataset SPT-3G$\times$ngWL is expected to obtain an upper
limit of $\log_{10}|f_{R0}|<-5.95$ at $95\,\%$ credibility, which
significantly improves upon the current best bounds. The CMB-S4$\times$ngWL
dataset is expected to improve this even further to $\log_{10}|f_{R0}|<-6.23$.
Furthermore, $f(R)$ gravity models with $\log_{10}|f_{R0}|\geq-6$, which have
larger numbers of clusters, would be distinguishable from GR with both
datasets. We also report degeneracies between $\log_{10}|f_{R0}|$ and
$\Omega_{\mathrm{m}}$ as well as $\sigma_{8}$ for $\log_{10}|f_{R0}|>-6$ and
$\log_{10}|f_{R0}|>-5$ respectively. Our forecasts indicate that future
cluster abundance studies of $f(R)$ gravity will enable substantially improved
constraints that are competitive with other cosmological probes.

cosmology, modified gravity, clusters, weak-lensing

##  I Introduction

The cause of the late-time accelerated expansion of the universe is one of the
most pertinent and challenging questions in cosmology [1, 2]. The phenomenon
can be explained within the framework of general relativity (GR) if a fluid
with negative pressure is introduced, or alternatively if a cosmological
constant $\Lambda$ is added to the Einstein-Hilbert action. The latter
scenario is known as the $\Lambda$ cold dark matter model, $\Lambda$CDM. An
alternative approach to explain cosmic acceleration is through the
modification of the Einstein-Hilbert action, in a manner which seeks to avoid
the addition of a cosmological constant, referred to as modified gravity. A
key feature of modified gravity is that the clustering of matter is in general
different from that of GR [see, e.g., 3, 4, 5, for reviews]. This means that
studies of the growth of structure can be a useful test of GR on cosmic
scales, which may then shed light on the underlying cause of the late-time
acceleration.

One of the most popular modified gravity models alters the Einstein-Hilbert
action using a non-linear function $f(R)$ of the scalar curvature $R$ [6].
This leads to an additional, gravitational-strength, fifth force. This fifth
force affects structure formation, introducing a scale-dependence and
enhancing the structure growth on galaxy cluster scales. A further feature of
$f(R)$ gravity models is the chameleon screening mechanism, which suppresses
the modification in high-density regions [7], ensuring that $f(R)$ remains
consistent with tests of GR in the Solar system [8]. In this work we use the
widely studied $f(R)$ model of Hu & Sawicki [9] with $n=1$ for which self-
consistent N-body simulations [10, 11, 12, 13] as well as semi-analytical
models for the halo mass function [14, 15] exist.

The abundance of massive galaxy clusters and their evolution over cosmic time
has long been recognized as a powerful probe of the cosmic acceleration [16,
17]. Cluster surveys have been used to constrain $\Lambda$CDM [18, 19, 20, 21,
22, 23, 24] as well as modified gravity models [25, 26, 27, 28]. The abundance
of collapsed haloes is described by the halo mass function (HMF), which,
through the scale-dependent structure growth, depends on the strength of the
$f(R)$ gravity model, where the strength of the model is encoded in the single
parameter $f_{R0}$. This, together with the nonlinear screening mechanism,
results in a mass-dependent enhancement of the HMF. We incorporate this via a
prefactor which depends on the critical overdensity of halo collapse in $f(R)$
gravity [27]. The critical overdensity in $f(R)$ gravity can be calculated via
a semi-analytical model [29, 14]. In this approach, the critical overdensity
is computed for each halo mass and redshift by solving a system of coupled
differential equations. This makes the calculation computationally expensive,
and so we construct an emulator to make the computation more efficient.

In practice, the mass of a galaxy cluster, and hence the HMF is not measured
directly. Instead, observable signatures are used to infer the cluster mass.
These observables are measured and then mapped to the halo mass through
observable–mass relations, which are generally of power-law form in mass and
redshift and include a model of the scatter in the observable at a fixed halo
mass and redshift. To calibrate these relations, a robust mass measurement is
needed. In this analysis, we employ weak gravitational lensing data (WL) for
the mass calibration of the cluster ensemble. The advantage of using the WL
signal to calibrate the halo masses is that the WL data provide mass
constraints with well-characterized and controllable biases and associated
uncertainties on those biases. Moreover, no assumption of hydrostatic or
virial equilibrium is required.

One way to detect massive galaxy clusters is through the thermal Sunyaev-
Zel’dovich effect (tSZE) [30]. This phenomenon arises from the upscattering of
the cosmic microwave background (CMB) photons by energetic electrons within
the intracluster medium (ICM). Note that the resulting spectral distortion
leads to an observable signature that is approximately independent of
redshift, contributing to the fact that tSZE surveys provide a clean probe to
study the growth of structure up to the highest redshifts where massive
clusters exist (i. e. , $z\simeq 2$). Ongoing and planned tSZE surveys such as
those from SPT-3G [31], conducted with the South Pole Telescope (SPT) [32],
Simons Observatory [33], or CMB-S4 [34] will detect thousands to tens of
thousands of galaxy clusters [35].

Datasets from next-generation galaxy weak-lensing (hereafter ngWL) surveys as
from the ongoing Euclid mission [36, 37] or the Vera C. Rubin Observatory [38,
39] will allow for improved and more robust measurements of cluster halo
masses. Therefore, combining future cluster tSZE survey data with ngWL data
will yield powerful probes of cosmology and modified gravity, greatly
improving upon the existing constraints from cluster abundance analyses.

Modified gravity models such as $f(R)$ in general not only modify the halo
mass function, but also the observable–mass relation [40]. Many models,
including $f(R)$, do not modify gravitational lensing, but modify the dynamics
of non-relativistic tracers such as gas and stars. Provided with an
independent calibration of the mass proxy, this could be used to put
additional constraints on $f(R)$. Here, we conservatively assume that this
effect is completely absorbed by the weak lensing mass calibration, such that
the calibrated mass function is the same as the intrinsic halo mass function
(within statistical and systematic errors) for both GR and $f(R)$.

In this work, we make realistic forecasts for the constraining power of weak-
lensing informed galaxy cluster abundance studies on $f(R)$ gravity. We use
tSZE surveys from SPT-3G and CMB-S4 combined with cluster mass calibration
constraints from ngWL datasets similar to those expected from Euclid and
Rubin. For this, we create mock cluster and WL data for different values of
$\log_{10}|f_{R0}|$ as well as for a GR cosmology. We analyze the mock
datasets with a cluster-by-cluster likelihood approach in two steps. First, we
use the ngWL dataset for the full cluster sample to constrain the
observable–mass relation parameters, and then we adopt the posteriors from
that first calculation as priors and employ the cluster abundance likelihood
to obtain constraints on the $f(R)$ gravity parameter $\log_{10}|f_{R0}|$ and
other cosmological parameters of interest.

The analysis presented in this work employs the state-of-the-art framework
developed for the cosmological analyses of galaxy clusters selected in the
SPT-SZ and SPTpol surveys, with weak-lensing mass calibration using data from
the Dark Energy Survey data and the Hubble Space Telescope (HST) [41]. This
analysis framework and dataset was used to obtain competitive cosmological
constraints on $\Lambda$CDM and $w$CDM cosmologies [24]. Constraining $f(R)$
gravity using SPT clusters with DES and HST weak-lensing is in progress and
will presented in a future work.

The paper is organized as follows. Section II presents the $f(R)$ modified
gravity model, briefly describing the Hu & Sawicki model and discusses the
semi-analytical spherical collapse model used to calculate the critical
overdensity and subsequently quantify how the HMF changes with respect to the
GR version. Moreover, we present emulators which are used to speed up the
calculation of the critical overdensity and its derivative with respect to
$\mathrm{ln}M$. We describe in Section III the different cosmological surveys
we use for the forecasts. The observable–mass relations used in our analysis
are summarized in Section IV. The generation of the mocks is described in
Section V. We present the likelihood and the analysis method in Section VI.
The results of our analysis are presented in Section VII. Finally, Section
VIII gives a brief summary of our work.

Throughout this paper $\mathcal{U}(a,b)$ denotes a uniform distribution
between limits $a$ and $b$, and $\mathcal{N}(\mu,\sigma)$ is a Gaussian
distribution with mean $\mu$ and standard deviation $\sigma$. In this analysis
we adopt the halo mass definition $M_{200\mathrm{c}}$, which is the mass
within the cluster region where the enclosed mean density is 200 times the
critical density.

##  II $f(R)$ Modified gravity

In $f(R)$ gravity models the Einstein-Hilbert action of general relativity
(GR) is modified to include an arbitrary function $f(R)$ of the scalar
curvature $R$ [6]

$S=\int\textnormal{d}^{4}x\sqrt{-g}\left[\frac{R+f(R)}{16\pi
G}+\mathcal{L}_{m}\right]\,.$ (1)

Here $g$ is the determinant of the GR metric tensor $g_{\mu\nu}$, $G$ the
gravitational constant and $\mathcal{L}_{m}$ is the matter Lagrangian. Note,
that we use natural units $c=\hbar=1$. The field equation for $f(R)$ gravity
can be obtained by varying the action with respect to the metric tensor

$\begin{split}G_{\mu\nu}+f_{R}R_{\mu\nu}-&\left(\frac{f}{2}-\square
f_{R}\right)g_{\mu\nu}\\\ &\hskip
42.67912pt-\nabla_{\mu}\nabla_{\nu}f_{R}=8\pi GT_{\mu\nu}\,,\end{split}$ (2)

where $G_{\mu\nu}$ denotes the Einstein tensor, $R_{\mu\nu}$ represents the
Ricci tensor, $T_{\mu\nu}$ is the energy-momentum tensor and
$f_{R}=\textnormal{d}f(R)/\textnormal{d}R$ which behaves as an additional
scalar degree of freedom and is named the scalaron field.

Under the quasistatic and weak-field approximation, the trace of the field
equation gives the equation of motion for $f_{R}$

$\nabla^{2}\delta f_{R}=\frac{1}{3}(\delta R-8\pi G\delta\rho)\,,$ (3)

with $\delta\rho=\rho-\bar{\rho}$. The modified Poisson equation is obtained
from the time-time component of Eq. (2)

$\nabla^{2}\Phi=\frac{16\pi G}{3}\delta\rho-\frac{1}{6}\delta R\,,$ (4)

with $\delta R=R-\bar{R}$ and the Newtonian potential $\Phi$ is defined via
$2\Phi=\delta g_{00}/g_{00}$. Combining these two equations gives for the
modified Poisson equation

$\nabla^{2}\Phi=4\pi G\delta\rho-\frac{1}{2}\nabla^{2}\delta f_{R}\,.$ (5)

So, in $f(R)$ gravity the Poisson equation includes an extra term directly
proportional to $\nabla^{2}\delta f_{R}$.

To be effective at late times and large scales in cosmology, $f(R)$ has to be
a decreasing function of $R$, so that $f_{R}<0$. Depending on the field value
$f_{R}$ we can distinguish two different regimes: the large-field regime where
$|f_{R}|\gg|\Phi|$, which corresponds to low curvature, and the small-field
regime, i. e. $|f_{R}|\ll|\Phi|$ and thus is related to high curvature. In the
first case, $\delta R\ll 8\pi G\delta\rho$ and the Poisson Eq. (4) corresponds
to an enhancement of $\Phi$, and hence gravitational forces, by a factor of
$4/3$. On the other hand, in the small field regime one obtains that $\delta
R\approx 8\pi G\delta\rho$ and thus Eq. (4) reduces to the GR Poisson
equation. Hence, structure growth becomes scale-dependent. The recovery of the
GR Poisson equation in the high curvature regimes in $f(R)$ gravity is induced
by the so-called chameleon screening mechanism [7] which makes $f(R)$ gravity
models consistent with solar system tests [8].

In this paper we use the Hu & Sawicki form for the function $f(R)$ which is
given by [9]

$f(R)=-m^{2}\frac{c_{1}\left(\frac{R}{m^{2}}\right)^{n}}{c_{2}\left(\frac{R}{m^{2}}\right)^{n}+1}\,,$
(6)

with $m^{2}=\Omega_{\mathrm{m}}H^{2}_{0}$, the free parameters $n$, $c_{1}$,
$c_{2}$ and $H_{0}$ the Hubble constant. When $c_{2}^{1/n}R/m^{2}\gg 1$, the
above equation is approximately

$f(R)\approx-m^{2}\frac{c_{1}}{c_{2}}-\frac{f_{R0}R_{0}^{n+1}}{nR^{n}}\,,$ (7)

where $R_{0}$ is the present background curvature and $f_{R0}:=f_{R}(R_{0})$
which is a parameter that quantifies the strength of the $f(R)$ gravity model.
Note that $f_{R0}<0$. This approximation is correct up to order
$\sim(f_{R0})^{2}$; since current constraints are already at the level of
$|f_{R0}|\lesssim 10^{-4}$ or better, Eq. (7) is entirely sufficient. Since we
do not have a strong theory prior on a particular scale of $f_{R0}$, we choose
to work with $\log_{10}|f_{R0}|$ for numerical convenience. We will also
impose a uniform prior on $\log_{10}|f_{R0}|$ in our analysis.

Since we recover GR in the limit $|f_{R0}|\to 0$, we obtain from Eq. (1)

$\frac{c_{1}}{c_{2}}=6\frac{\Omega_{\Lambda}}{\Omega_{\mathrm{m}}}\,.$ (8)

In this paper, we further adopt $n=1$, as this is by far the most widely
studied $f(R)$ scenario; see Ref. [42] for an approach to approximately
rescale constraints from $n=1$ to other values of $n$.

###  II.1  Spherical collapse in $f(R)$ gravity

As described in the previous section, structure growth is scale-dependent in
$f(R)$ gravity, and consequently the formation of dark matter halos is too.
For the formation of clusters we use a spherical collapse model developed and
described in detail elsewhere [29, 14, 12].

We approximate a dark matter halo as a spherical top-hat overdensity, with a
radius $R_{\mathrm{TH}}$ and constant density inside $\rho_{\mathrm{in}}$ and
outside $\rho_{\mathrm{out}}$. The mass and environment dependence of the
spherical collapse model [29] with the chameleon screening effect is
implemented through the thin-shell approximation [7].

The physical radius of the spherical top-hat halo is defined as $\xi(a)$,
where $\xi(a_{\mathrm{i}})=a_{\mathrm{i}}R_{\mathrm{TH}}$ at initial scale
$a_{\mathrm{i}}$. Due to the nonlinear evolution of the overdensity, the
radius $\xi(a)$ deviates from the linear relation for larger scale factor $a$
and we define the variable $y$ to be the deviation from this linear relation,
$y(a)=\xi(a)/aR_{\mathrm{TH}}$. Note that due to mass conservation, we have
$\rho/\bar{\rho}=y^{-3}$. The equation of motion for the spherical shell is
given by [29, 43, 14]

$\frac{\ddot{\xi}}{\xi}=-\frac{4\pi
G}{3}(\bar{\rho}_{m}-2\bar{\rho}_{\lambda})-\frac{4\pi
G}{3}(1+F)\delta\rho_{\mathrm{m}}\,,$ (9)

where dots denote derivatives with respect to time. Here $F$ is the extra
force from modified gravity and is given by the thin-shell approximation [7,
14, 12]

$\small
F\left(\frac{\Delta\xi}{\xi}\right)=\frac{1}{3}\mathrm{min}\left(3\frac{\Delta\xi}{\xi}-3\left(\frac{\Delta\xi}{\xi}\right)^{2}+\left(\frac{\Delta\xi}{\xi}\right)^{3},1\right)\,.$
(10)

Here $\frac{\Delta\xi}{\xi}$ is the thickness of the thin shell and can be
expressed as [see 7, 14]

$\displaystyle\frac{\Delta\xi}{\xi}\approx$
$\displaystyle\frac{|f_{R0}|a^{3n+4}}{\Omega_{\mathrm{m}}(H_{0}R_{\mathrm{TH}})^{2}}y_{\mathrm{h}}\Biggl{[}\biggl{(}\frac{1+4\frac{\Omega_{\Lambda}}{\Omega_{\mathrm{m}}}}{y_{\mathrm{env}}^{-3}+\frac{\Omega_{\Lambda}}{\Omega_{\mathrm{m}}}a^{3}}\biggr{)}^{n+1}$
$\displaystyle\hskip
59.75095pt-\biggl{(}\frac{1+4\frac{\Omega_{\Lambda}}{\Omega_{\mathrm{m}}}}{y_{\mathrm{h}}^{-3}+\frac{\Omega_{\Lambda}}{\Omega_{\mathrm{m}}}a^{3}}\biggr{)}^{n+1}\Biggr{]}\,,$
(11)

where $y_{\mathrm{h}}^{-3}$, $y_{\mathrm{env}}^{-3}$ track the inner and outer
overdensities respectively and $n=1$ due to our choice of the Hu & Sawicki
model. The evolution for the inner overdensity $y_{\mathrm{h}}^{-3}$ can be
expressed as [see 14]:

$\begin{split}y_{\mathrm{h}}^{{}^{\prime\prime}}&+\left(2-\frac{3}{2}\Omega_{\mathrm{m}}(a)\right)y_{\mathrm{h}}^{{}^{\prime}}\\\
&\hskip
51.21504pt+\frac{1}{2}\Omega_{\mathrm{m}}(a)(1+F)(y_{\mathrm{h}}^{-3}-1)y_{\mathrm{h}}=0\,,\end{split}$
(12)

with primes denoting derivatives with respect to $\mathrm{ln}a$. For the outer
overdensity, we assume that it follows a GR evolution, thus $F\to 0$ and one
obtains [14]

$\begin{split}y_{\mathrm{env}}^{{}^{\prime\prime}}&+\left(2-\frac{3}{2}\Omega_{\mathrm{m}}(a)\right)y_{\mathrm{env}}^{{}^{\prime}}\\\
&\hskip
51.21504pt+\frac{1}{2}\Omega_{\mathrm{m}}(a)(y_{\mathrm{env}}^{-3}-1)y_{\mathrm{env}}=0\,.\end{split}$
(13)

Because the equation of motion for the inner region of the top-hat region
depends on $F$ and thus by Eq. (II.1) on the outer region, Eqs. (12) and (13)
are a system of coupled differential equations. The initial conditions at
$a_{\mathrm{i}}\ll 1$ are in the matter-dominated regime and given by

$\displaystyle y_{\mathrm{h/env,i}}$
$\displaystyle=1-\frac{\delta_{\mathrm{h/env,i}}}{3}\,,$ (14) $\displaystyle
y^{{}^{\prime}}_{\mathrm{h/env,i}}$
$\displaystyle=-\frac{\delta_{\mathrm{h/env,i}}}{3}\,.$ (15)

If we want to find the critical overdensity $\delta_{\mathrm{crit}}$ which
causes spherical collapse at scale factor $a_{\mathrm{c}}$, we have to choose
the initial conditions $y_{\mathrm{h,i}}$ and $y^{{}^{\prime}}_{\mathrm{h,i}}$
such that the solution of Eq. (12) gives $y_{\mathrm{h}}(a_{\mathrm{c}})=0$
with the requirement $y_{\mathrm{h}}(a)>0$ for $a<a_{\mathrm{c}}$. The
critical overdensity is then defined as the linearly extrapolated value of
$\delta_{\mathrm{h/env,i}}$, i. e.
$\delta_{\mathrm{crit}}=D(a)/D(a_{\mathrm{i}})\delta_{\mathrm{h/env,i}}$ with
the linear growth factor $D(a)$.

The initial overdensity $\delta_{\mathrm{env,i}}$ of the outer region is set
by the peak of the probability distribution of an Eulerian environmental
density $\delta_{\mathrm{env}}$ with an Eulerian (physical) radius
$\xi=5\,h^{-1}\,\mathrm{Mpc}$ [for more details see 14, 44] and is derived
from a physical model to be [44]

$\displaystyle P_{\xi}(\delta_{\mathrm{env}})=$
$\displaystyle\frac{\beta^{\omega/2}}{\sqrt{2\pi}}\left[1+(\omega-1)\frac{\delta_{\mathrm{env}}}{\delta_{\mathrm{crit,\mathrm{GR}}}}\right]\left(1-\frac{\delta_{\mathrm{env}}}{\delta_{\mathrm{crit,\Lambda}}}\right)^{\omega/2-1}$
$\displaystyle\times\exp\left(-\frac{\beta^{\omega}\delta_{\mathrm{env}}}{2(1-\delta_{\mathrm{env}}/\delta_{\mathrm{crit,\mathrm{GR}}})^{\omega}}\right)\,,$
(16)

with $\beta=(\xi/8)^{3/\delta_{\mathrm{crit,\mathrm{GR}}}}(\sigma_{8}^{\rm
GR})^{-2/\omega}$, $\delta_{\mathrm{crit,\mathrm{GR}}}$ the linear
extrapolated GR overdensity and
$\omega=\delta_{\mathrm{crit,\mathrm{GR}}}\gamma$ with

$\gamma=\frac{n_{\mathrm{s}}+3}{3}\,.$ (17)

In $f(R)$ gravity the critical overdensity for halo collapse
$\delta_{\mathrm{crit}}$ is a function of the mass of the spherical
overdensity $M=4/3\pi\bar{\rho}_{\mathrm{m}}R^{3}_{\mathrm{RT}}$ via Eq. (10),
(II.1) and the local overdensity of the environment $\delta_{\mathrm{env}}$ is
given by Eq. (II.1). Figure 1 shows the critical overdensity as a function of
mass for different values of $\log_{10}|f_{R0}|$ for the Planck cosmology
given in Table 1 [45, hereafter Planck 2018]. For reference, Fig. 1 also shows
the constant critical overdensity in a corresponding GR universe. Note that
the critical overdensity in $f(R)$ gravity approaches the GR value for high
halo masses, because $\frac{\Delta\xi}{\xi}\to 0$ for high masses and thus
$F\to 0$. The critical overdensity in $f(R)$ gravity and in GR as well as the
derivative of $\delta_{\mathrm{crit}}$ with respect to $\mathrm{dln}M$ are
used to calculate the halo mass function in $f(R)$ gravity models.

Figure 1: The critical overdensity $\delta_{\mathrm{crit}}$ for spherical collapse in $f(R)$ gravity (Eq. (12)) for different values of $\log_{10}|f_{R0}|$ at collapse redshift $z_{\mathrm{c}}=0$ in colored solid lines. The dashed black line represents $\delta_{\mathrm{crit}}$ in a corresponding GR cosmology (Eq. (13)).  Table 1: Fiducial cosmological parameters and their values for the flat Planck 2018 cosmology [45, table 2, column 4]. $\Omega_{\mathrm{m}}$ is the matter density parameter, $\Omega_{\mathrm{b}}h^{2}$ is the reduced baryon density, $h$ is the Hubble parameter, $N_{\textnormal{eff}}$ is the effective number of relativistic species, $\sum m_{\nu}=\Omega_{\nu}h^{2}\cdot 94\,\mathrm{eV}$ is the sum of neutrino masses, $A_{s}$ is the power spectrum amplitude parameter, $\sigma_{8}^{\rm GR}$ is the amplitude of the linear power spectrum at a scale $8\,h^{-1}\,\mathrm{Mpc}$ at redshift $z=0$ and $n_{s}$ is the scalar spectral index. Parameter | Fiducial value
---|---
$\Omega_{\mathrm{m}}$ | 0.3166
$\Omega_{\mathrm{b}}h^{2}$ | 0.02236
$h$ | 0.6727
$N_{\textnormal{eff}}$ | 3.046
$\Sigma m_{\nu}$ [eV] | 0.06
$\ln(10^{10}A_{s})$ | $2.1\times 10^{-9}$
$\sigma_{8}^{\rm GR}$ | 0.812
$n_{s}$ | 0.9649

###  II.2  Halo mass function in $f(R)$ gravity

To study and forecast cosmological constraints with cluster surveys we have to
understand the distribution of halos in mass and redshift, i. e. the
(differential) HMF. We model the HMF in an $f(R)$ gravity cosmology with two
components. The first component is the GR HMF, and the second is the
enhancement (or suppression) factor of the GR HMF due to $f(R)$ gravity. This
factor accounts for the scale-dependent clustering of matter [46, 27]. We
adopt the GR halo mass function from [47], which is given by

$\left.\frac{\textnormal{d}n}{\textnormal{d}\mathrm{ln}M}\right|_{\mathrm{T}}=-\frac{\bar{\rho}_{\mathrm{m}}}{2M}f(\sigma)_{\mathrm{T}}\frac{\textnormal{d}\mathrm{ln}\sigma^{2}}{\textnormal{d}\mathrm{ln}M}\,.$
(18)

Here $f(\sigma)_{\mathrm{T}}$ is the multiplicity function [47]

$f(\sigma)_{\mathrm{T}}=\tilde{A}\left[\left(\frac{\sigma}{\tilde{b}}\right)^{-\tilde{a}}+1\right]e^{-\frac{\tilde{c}}{\sigma^{2}}}\,,$
(19)

where $\tilde{A},\ \tilde{a},\ \tilde{b}$ and $\tilde{c}$ are parameters
calibrated using N-body simulations [see 47, table 2] and $\sigma=\sigma(M)$
is the variance of the overdensity on a mass scale $M$ in the corresponding GR
cosmology.

For the ratio of the HMF in $f(R)$ gravity to GR we use the Sheth & Tormen
halo mass function [48], which can account for the scale-dependent collapse
through $\delta_{\mathrm{crit}}$. Furthermore, Refs. [43, 14] showed that the
Sheth & Tormen prescription with modified collapse thresholds provides a good
fit to the HMF in N-body simulations of $f(R)$ gravity. The Sheth & Tormen HMF
is given by

$\left.\frac{\textnormal{d}n}{\textnormal{d}\mathrm{ln}M}\right|_{\mathrm{ST}}=\frac{\bar{\rho}_{\mathrm{m}}}{M}f(\nu)_{\mathrm{ST}}\left[\frac{\textnormal{d}\mathrm{ln}\delta_{\mathrm{crit}}}{\textnormal{d}\mathrm{ln}M}-\frac{1}{2}\frac{\textnormal{d}\mathrm{ln}\sigma^{2}}{\textnormal{d}\mathrm{ln}M}\right]\,,$
(20)

with $\nu=\delta_{\mathrm{crit}}/\sigma$ the peak height and
$f(\nu)_{\mathrm{ST}}$ the Sheth & Tormen multiplicity function, which is
parameterized as

$f(\nu)_{\mathrm{ST}}=A\sqrt{\frac{a\nu}{2\pi}}\left[1+(a\nu^{2})^{-p}\right]e^{-\frac{a\nu^{2}}{2}}\,,$
(21)

where $A,\ a$, $p$ are free parameters. We adopt fitting formulae for these
parameters [49]. This allows one to calculate the HMF for different halo
overdensities $\Delta_{\mathrm{c}}$, i. e. different halo definitions
$M_{\Delta_{\mathrm{c}}}=4/3\pi\Delta_{\mathrm{c}}\bar{\rho}_{\mathrm{crit}}R^{3}$.
The fitting formulae are given by

$\displaystyle A$ $\displaystyle=-0.1362x+0.3292\,,$ $\displaystyle a$
$\displaystyle=0.4332x^{2}+0.2263x+0.7665\,,$ (22) $\displaystyle p$
$\displaystyle=-0.1151x^{2}+0.2554x+0.2488\,.$

Here $x(z)=\log_{10}(\Delta/\Delta_{\mathrm{vir}}(z))$ and
$\Delta_{\mathrm{vir}}(z)$ is the virial overdensity. In this model we use an
effective virial overdensity from the corresponding GR cosmology which is
given by [50]

$\Delta_{\mathrm{c,vir}}(z)=18\pi^{2}-82\left[1-\Omega_{\mathrm{m}}(z)\right]-39\left[1-\Omega_{\mathrm{m}}(z)\right]^{2}\,.$
(23)

The difference between the Sheth & Tormen HMF in $f(R)$ gravity and GR is
completely encoded in the critical overdensity coming from Eq. (12) and (13).
Therefore, we use an effective variance $\sigma=\sigma(M)$ and thus employ in
both cases the GR value for the variance [14, 15]. Furthermore, if we include
massive neutrinos, the shape of the HMF is more universal if only the cold
dark matter and the baryon power spectrum is used to calculate the variance
[51]. We adopt this also in the case of an $f(R)$ gravity cosmology to account
for the effect of massive neutrinos, because there is a known degeneracy
between the $f(R)$ gravity parameter and massive neutrinos [see e. g. 52, 53,
54]. In summary, our $f(R)$ HMF has the following form

$\displaystyle\frac{\textnormal{d}n}{\textnormal{d}\mathrm{ln}M}$
$\displaystyle=\left.\frac{\textnormal{d}n}{\textnormal{d}\mathrm{ln}M}\right|_{\mathrm{T}}\times\mathcal{R}\,,$
(24) $\displaystyle\mathrm{with}\quad\mathcal{R}$
$\displaystyle=\frac{\left.\frac{\textnormal{d}n}{\textnormal{d}\mathrm{ln}M}\right|_{\mathrm{ST},\
f(R)}}{\left.\frac{\textnormal{d}n}{\textnormal{d}\mathrm{ln}M}\right|_{\mathrm{ST},\
\mathrm{GR}}}\,.$

Figure 2: The ratio $\mathcal{R}$ of the $f(R)$ HMF and the GR Sheth & Tormen
HMF in Eq. (20), for different values of the $f(R)$ parameter
$\log_{10}|f_{R0}|$ in colored lines.

The difference between the HMF in an $f(R)$ gravity model to GR for different
values of $\log_{10}|f_{R0}|$, i. e. the ratio $\mathcal{R}$ from Eq. (24), is
shown in Fig. 2. Generally, $f(R)$ gravity enhances the growth of structure
and thus $f(R)$ gravity predicts more clusters compared to a GR cosmology. As
expected, models with larger $|f_{R0}|$ show larger differences in the HMF.
Moreover, the shape of the enhancement depends on the strength of the $f(R)$
model. If the value of $|f_{R0}|$ is comparable to the cosmological potential
$\Phi$ of massive dark matter halos (see Eq. (5)), the fifth force is screened
within these halos, so that the abundance of high mass halos is not increased
substantially. This is the case for $|f_{R0}|\leq 10^{-6}$. In the case of
$\log_{10}|f_{R0}|=-5$, a high enhancement for the most massive halos is seen,
because for this the value of $f_{R0}$ even massive halos are at most
partially screened. Furthermore, Fig. 2 shows that the HMF enhancement is
weaker for smaller $f_{R0}$. Consequently, distinguishing between such weak
modified gravity models and GR using a cluster abundance analysis will be
challenging (see Sec. VI.2.2).

###  II.3  Emulating the critical overdensity quantities

The computational bottleneck in our HMF calculation is the computation of the
critical overdensity in $f(R)$ gravity and GR as well as the derivative of
$\delta_{\mathrm{crit}}$ with respect to $\mathrm{ln}M$. These quantities are
obtained by solving a system of coupled differential equations for each mass
and redshift individually in the range where we need the HMF and thus a lot of
computations must be performed to obtain the HMF. To speed up the calculations
of $\delta_{\mathrm{crit}}$,
$\textnormal{d}\delta_{\mathrm{crit}}/\textnormal{d}\mathrm{ln}M$ and
$\delta_{\mathrm{crit},\mathrm{GR}}$ we create three separate emulators for
these quantities

We choose to emulate the critical density and its derivative rather than the
HMF directly, because a direct emulation of the HMF would require us to
additionally sample baryon and neutrino density parameters, thereby making the
emulation more complex. In the end, executing our emulators and computing the
HMF is essentially as fast as emulating the HMF.

We use Gaussian processes regression (GPR), a supervised learning method, to
build the emulators. The emulators are trained on a data set that samples the
desired parameter space and then verified on an independent validation dataset
to assess the performance of the emulators. The parameter space of
$\delta_{\mathrm{crit}}$ and
$\textnormal{d}\delta_{\mathrm{crit}}/\textnormal{d}\mathrm{ln}M$ that we
sample is seven-dimensional, with five cosmological parameters together with
the halo mass $M$ and redshift $z$. This allows us to emulate the behavior of
the coupled system of differential equations (12) and (13). The ranges we
choose for the parameters are

$\displaystyle z$ $\displaystyle\in[0,2]\,,$ $\displaystyle\log_{10}M$
$\displaystyle\in[13,16]\,,$ $\displaystyle\Omega_{\mathrm{m}}$
$\displaystyle\in[0.11,0.4]\,,$ $\displaystyle h$
$\displaystyle\in[0.6,0.82]\,,$ (25) $\displaystyle n_{\mathrm{s}}$
$\displaystyle\in[0.8,1.1]\,,$ $\displaystyle\sigma_{8}^{\mathrm{GR}}$
$\displaystyle\in[0.6,0.9]\,,$ $\displaystyle\log_{10}|f_{R0}|$
$\displaystyle\in[-7,-3]\,.$

The ranges for the cosmological parameters are large enough to ensure that we
are not hitting the boundaries in our likelihood sampling, see Sec. VI.2. In
the case of the halo mass and redshift, we choose ranges such that we cover
the interval of masses and redshifts of clusters detectable with the surveys
from SPT-3G and CMB-S4.

The parameter space for $\delta_{\mathrm{crit},\mathrm{GR}}$ is only two-
dimensional with the parameters $\Omega_{\mathrm{m}}$ and $z$, because the
spherical collapse equation, Eq. (13), for GR only depends on these two
quantities. The ranges for these parameters are the same as for the two other
emulators.

We sample the points in the parameter space with a Sobol sequence algorithm.
Here we use 512 points for training and 512 points for validating the
emulators. Both datasets are evenly distributed in the parameter space and are
disjoint. Figure 3 shows the parameter space with training data in blue and
the validation set in orange; the characteristic pattern for the Sobol
sequence algorithm can be seen.

Figure 3:  The 7-dimensional parameter space for the emulation of
$\delta_{\mathrm{crit}}$ and
$\textnormal{d}\delta_{\mathrm{crit}}/\textnormal{d}\mathrm{ln}M$ drawn from a
Sobol sequence algorithm. In total 1024 points are sampled. The training set
(first half of the points) is marked in blue and the validation data are shown
in orange. The ranges of the parameters are given in Eq. (II.3).

The performance of the emulator is quantified by the normalized median
absolute deviation (nMAD), which is a robust estimator of the scatter around
the true value. The nMAD is a better quantity for the deviation from the
average than the standard deviation if the dataset shows a large scatter.

Our validation tests show emulator accuracy for $\delta_{\mathrm{crit}}$,
$\textnormal{d}\delta_{\mathrm{crit}}/\textnormal{d}\mathrm{ln}M$ and
$\delta_{\mathrm{crit,GR}}$ of 0.2 %, 22.9 % and 0.0002 %, respectively. Thus,
the critical overdensity in $f(R)$ gravity and GR computed with the emulator
is accurate at the sub-percent level. However, the logarithmic mass derivative
of the critical overdensity can only be predicted at the $\sim 25\,\%$
accuracy level by the emulator.

We test also whether increasing the accuracy of the
$\textnormal{d}\delta_{\mathrm{crit}}/\textnormal{d}\mathrm{ln}M$ emulator by
reducing the dynamical range of the quantity and emulating
$\mathrm{ln}(\textnormal{d}\delta_{\mathrm{crit}}/\textnormal{d}\mathrm{ln}M)$
instead. Unfortunately, this approach does not increase the accuracy of the
emulator. The relatively high inaccuracy can be explained by the inherent
complexity of
$\textnormal{d}\delta_{\mathrm{crit}}/\textnormal{d}\mathrm{ln}M$, which is
the halo mass as a free parameter, along with a quantity that is emulated
while involving differentiation with respect to the halo mass. We can reach a
higher accuracy by either using more points in the training dataset or by
further developing or replacing the GPR emulator used to emulate
$\textnormal{d}\delta_{\mathrm{crit}}/\textnormal{d}\mathrm{ln}M$. However,
our forecast depends on the accuracy of the HMF rather than the accuracy of
the critical overdensity and derivative. As discussed in the next section, the
error of the emulated HMF is mainly driven by the error of
$\delta_{\mathrm{crit}}$, and the error of
$\textnormal{d}\delta_{\mathrm{crit}}/\textnormal{d}\mathrm{ln}M$ does not
have a big impact.

###  II.4 HMF validation

Because we use the emulated critical overdensities in $f(R)$ gravity and GR as
well as the logarithmic derivative with respect to mass, we examine the impact
of inaccuracies in the emulation on the $f(R)$ HMF, Eq. (24). We compute the
ratio $\mathcal{R}$, Eq. (24), for the points in the validation dataset once
with the analytical result and compare them to the results using the emulated
values.

Figure 4: Ratio of the true HMF to the HMF calculated using the emulator
values of $\delta_{\mathrm{crit},\mathrm{GR}}$, $\delta_{\mathrm{crit}}$ and
$\textnormal{d}\delta_{\mathrm{crit}}/\textnormal{d}\mathrm{ln}M$. The ratio
is calculated for the validation dataset and shown as a function of halo mass;
the validation dataset spans the full range of redshift, mass and cosmological
parameters adopted for the emulation. The gray dashed and dotted lines
represent the one and two $\sigma$ (actually nMAD) median absolute differences
within six mass bins.

Figure 4 shows the ratio
$\mathcal{R}_{\mathrm{true}}/\mathcal{R}_{\mathrm{emu}}$ of the semi-
analytical and the emulator results for the validation points as a function of
mass. These validation points span the full range of redshift over which the
emulators were trained. One can see that as the mass increases the scatter in
the ratio $\mathcal{R}_{\mathrm{true}}/\mathcal{R}_{\mathrm{emu}}$ also
increases, indicating a greater characteristic uncertainty in the HMF
emulation for higher masses. Moreover, the accuracy of the HMF is worse for
locations where the amplitude of the halo mass function is lower (as indicated
by the color bar). We quantify the accuracy of the emulated HMF by the nMAD in
six mass bins evenly spaced in logspace between
$[10^{13},10^{16}]\,h^{-1}M_{\odot}$ and we obtain $\rm
nMAD=(0.7,0.5,0.7,0.8,0.9,3.0)\,\%$. Note that the nMAD is computed from the
test set which is sampled in redshift, mass and cosmological parameters and
thus averaged over the entire ranges of redshift and cosmological parameters.
As seen the nMAD in the last mass bin, i. e.
$[15.5,16]\,\log[M/(h^{-1}M_{\odot})]$, is significantly larger than that of
the other bins. However, these ultra-high mass clusters are exceedingly rare
in the Universe, and therefore the statistical uncertainties due to shot noise
are larger than the two to five percent inaccuracies of the HMF in this
regime.

Furthermore, Fig. 4 shows that the relatively large uncertainty in the
derivative of the $f(R)$ critical overdensity with respect to $\mathrm{dln}M$
does not significantly affect the accuracy of the HMF. This is expected.
First, the exponential form in the overdensity of the multiplicity function,
Eq. (21), is sensitive to the absolute difference of the emulated and semi-
analytical result of $\delta_{\mathrm{crit}}$. And second,
$\delta_{\mathrm{crit}}$ has a much lower dynamical range than $\sigma$, i. e.
$\textnormal{d}\delta_{\mathrm{crit}}/\textnormal{d}\mathrm{ln}M$ is smaller
than $\textnormal{d}\sigma/\textnormal{d}\mathrm{ln}M$ in the Sheth & Tormen
HMF, Eq. (20) and therefore makes only a small contribution to the HMF.

The semi-analytical HMF has previously been shown to be consistent with
simulations in the $20\,\%$ range [14]. New results from the FORGE HMF
emulator show agreement with the semi-analytical HMF at the $15\,\%$ level
[55]. Therefore, the inaccuracy in our HMF due to employing an emulator for
the critical overdensities and logarithmic mass derivative is not a limiting
factor in our analysis.

The calculation of the HMF with the semi-analytical model on a mass and
redshift grid of twenty times twenty points and a Planck 2018 cosmology with
$\log_{10}|f_{R0}|=-6$ needs around 1920 seconds. In comparison the HMF with
the emulated quantities needs roughly 1.4 seconds, corresponding to more than
three orders of magnitude reduction in runtime.

##  III  Cosmological surveys

In this analysis, we focus on the constraining power on $f(R)$ gravity of two
different galaxy cluster samples when combined with weak gravitational lensing
mass calibration. We describe here the two galaxy cluster samples, which are
representative of SPT-3G and CMB-S4, and the ngWL weak lensing dataset needed
for mass calibration, which could come from either the Euclid satellite or
Rubin observatory.

###  III.1 tSZE-selected cluster catalogs

This section describes the cluster survey specifics for SPT-3G and CMB-S4,
which are used to forecast constraints on $f(R)$ gravity. In modeling the
cluster samples we follow the approach adopted in the recent SPT$\times$DES
cluster cosmological analysis [24]. The massive galaxy clusters in the two
samples

are selected using their thermal SZE signature (tSZE). A matched filter tuned
to detect the cluster tSZE is applied to the multi-frequency maps over a broad
range of cluster core radii [56, 57] to identify peaks. Given the noise in the
maps and the amplitude of the tSZE peak, an associated detection significance
$\hat{\zeta}$ is assigned. The approach produces a candidate list of tSZE
selected clusters with an initial contamination fraction. The $\hat{\zeta}$
observable serves as a mass proxy for the cluster analysis as described in
Section IV.

Each cluster candidate is then confirmed by applying a matched-filter
technique to the optical and near-infrared (NIR) galaxy catalog at the
location of each cluster [e.g., 58, 59, 60]. In this process, each cluster is
assigned a redshift $z$ and a cluster richness $\hat{\lambda}$, which
corresponds to the color and position weighted number of passive galaxies
within the cluster. For follow-up tools like the multi-component matched
filter [MCMF; 61, 60], the chance that each candidate is a random
superposition of a tSZE noise fluctuation and a physically unassociated
optical/NIR system is also quantified. With this approach it is possible to
exclude the likely contaminants from the tSZE selected list, producing a
catalog with a specific targeted contamination (typically percent level).

Each confirmed cluster then has three observables used for the selection: tSZE
significance $\hat{\zeta}$, richness $\hat{\lambda}$ and redshift $z$, and the
sample for analysis is defined by the selection thresholds $\hat{\zeta}_{\rm
min}$, $\hat{\lambda}_{\rm min}(z)$ and $z_{\mathrm{min}}<z<z_{\mathrm{max}}$.
Here the redshift dependent richness selection comes from MCMF and is required
to reduce the contamination in the tSZE candidate cluster list. The values of
these thresholds depend on the cluster survey and are therefore specified in
the next sections.

In our analysis, each cluster has a fourth observable, which is the weak
lensing inferred mass, which we discuss further below.

####  III.1.1  SPT-3G cluster survey

The South Pole Telescope (SPT) is a $10$ m telescope operating in the
microwave wavelength and is located near the Amundsen-Scott South Pole Station
[32]. Since 2018 the SPT has been equipped with a new, third-generation
camera, SPT-3G [31], which has been used to detect approximately an order of
magnitude more clusters per sky solid angle than from the previous SPT
surveys: SPT-SZ [59, 62], SPTpol ECS [63] and SPTpol 500d [64]. In this work,
we consider a $4,000\,\mathrm{deg}^{2}$ survey by SPT-3G that consists of two
regions: a deep $1,500\,\mathrm{deg}^{2}$ field for which we model the
selected cluster sample to be those with detection significance
$\hat{\zeta}>4.25$, and a shallower $2,500\,\mathrm{deg}^{2}$ field for which
we adopt a selection $\hat{\zeta}>5$. Despite its smaller survey area, there
are more detected clusters in the deep field because of its greater depth.

For the sample analyzed here, we assume that the full survey from SPT-3G is
covered by the upcoming ngWL survey datasets similar to those that will be
available from Euclid and Rubin. These datasets are crucial not only for the
weak-lensing information to calibrate cluster masses, but also for the
optical/NIR confirmation of SPT-3G cluster candidates in the range
$z\in[0.25,2]$.111We adopt the same lower redshift limit as adopted in
previous cosmological analyses of SPT clusters [21, 22, 24]. The selection
threshold associated with the MCMF exclusion of contaminants is
$\hat{\lambda}_{\mathrm{min}}(z)$, which we take to be the same as that used
in the recent SPT$\times$DES analysis [41, 24]. These thresholds are chosen to
lead to a contamination fraction in the final SPT-3G confirmed cluster sample
at the percent level, where no explicit modeling of the contaminants will be
required.

####  III.1.2  CMB-S4 cluster survey

CMB-S4 is a future cosmic microwave background survey currently in the design
and construction phase that will start operation at the end of this decade
[34]. The survey will cover roughly 50 % of the sky, and the resulting mm-wave
maps will be ideal for the detection of tens of thousands of clusters through
their tSZE signatures, as shown in a recent forecast for CMB-S4 [35].

Following this previous forecast, we adopt a cluster redshift range
$z\in[0.1,3]$ and a tSZE detection significance limit $\hat{\zeta}>5$. Note
that this sample extends to lower redshift than the SPT samples because CMB-S4
has more frequency bands, enabling greater reliance on frequency filtering to
remove the primary CMB anisotropies which affects cluster detection at low
redshift if unmitigated. We only consider clusters up to redshift $z=2$ for
the following two reasons. First, follow-up of clusters at $z>2$ is not yet
demonstrated and well understood. As previously noted, we expect the deep ngWL
imaging datasets in the optical/NIR to be sufficient for follow-up to $z=2$.
Second, we have chosen the Tinker HMF (Eq. (18)) for the mock generation, and
it is only calibrated at $z\lesssim 2$.

We adopt a survey area corresponding to the overlap between CMB-S4 and the
upcoming ngWL survey from the Euclid satellite. We use only the overlapping
region to ensure that every confirmed CMB-S4 cluster will have available weak-
lensing data. This is a conservative approach, because the HMF constraints
from cluster samples do not require that each individual cluster have
available weak-lensing data. With a Euclid-like footprint for our ngWL dataset
the overlapping region is roughly $10,100\,\mathrm{deg}^{2}$ [35, 37]. The
Rubin coverage would be even larger and would enhance the CMB-S4 sample
relative to what we adopt here. For the optical/NIR selection
$\hat{\lambda}_{\mathrm{min}}(z)$ from MCMF, we adopt also the selection
thresholds recently used in the SPT$\times$DES analysis [41, 24].

###  III.2 Next-generation weak-lensing (ngWL) data

Next-generation weak-lensing surveys will collect an order of magnitude more
lensing data than current WL surveys and will thus provide improved
constraints on the parameters of the observable–mass relation relative to what
is possible today.

In this analysis we focus on ngWL data similar to that which we expect from
the Euclid mission [36, 37]. Note that a similar analysis could be carried out
for other ngWL surveys like the Legacy Survey of Space and Time conducted with
the Vera C. Rubin Observatory [39]. Given the similarities of the weak-lensing
datasets and photometric redshifts between the two surveys, we carry out a
single analysis.

The most significant expected advantage in constraining power for Rubin
compared to Euclid comes from the survey footprint. In the case of Rubin, the
overlap with CMB-S4 is expected to be $\sim 25\,\%$ larger than for Euclid,
and thus constraints would be approximately $10\,\%$ tighter. For SPT-3G both
Euclid and Rubin fully cover the planned survey region, and so there would be
no expected significant differences in constraining power.

We assume a lensing source density of $30\,\mathrm{arcmin}^{-2}$, a shape
noise of 0.3, and a source redshift distribution with a median redshift
$z_{\mathrm{m}}=0.9$ [36]. These characteristics are similar to the adopted
goals for the ngWL surveys from Euclid and Rubin. We assume an uncertainty in
source redshifts of $\sigma_{z}=0.06$ [36], and bin the lensing source
galaxies into ten tomographic bins evenly spaced in redshift between $z=0$ and
$z=2$, and add an eleventh bin in the redshift range 2.0–2.6. Figure 5 shows
the redshift distributions of the ngWL source galaxies and of the CMB-S4
clusters.

Figure 5: Normalized redshift distribution of the ngWL sources and the CMB-S4
cluster sample. The full ngWL source redshift distribution is split up into
eleven tomographic bins. Colors are used to delineate redshifts.

##  IV  Observable-mass relations

Galaxy cluster ensembles exhibit considerable regularity, shown first in tight
observable–observable scaling relationships involving the X-ray sizes,
temperatures and ICM masses [65, 66] and then more recently in observable–mass
scaling relations calibrated using weak lensing data [e.g. 67, 68, 69, 70].
Because galaxy clusters are typically identified and selected using their
observable properties, while the HMF is expressed in terms of halo mass and
redshift, scientific analyses of cluster samples usually rely on the existence
of these observable–mass scaling relations. The analysis method we employ here
is based on an empirical calibration of the observable–mass scaling relations
using information from weak gravitational lensing based mass estimates as
previously demonstrated with tSZE selected cluster samples from SPT [22, 41,
24].

As mentioned in the introduction, while $f(R)$ gravity, in general, modifies
the dynamics and hence also the tSZE signal, it does not affect gravitational
lensing directly. To be precise, the lensing signal is rescaled by a factor of
$(1+|f_{R0}|)^{-1}$ [71, 72], which for the range of $f_{R0}$ we consider is
negligible, and we can thus assume that gravitational lensing is the same in
$f(R)$ gravity and GR. Hence, we assume in the following that all modified
gravity effects on the mass proxy are calibrated via the weak-lensing mass
calibration. In this work, we assume that $f(R)$ gravity does not strongly
affect the halo shapes and that the GR-based weak-lensing mass calibration we
will introduce in Section IV.3 is valid. A fully self-consistent approach will
require $f(R)$ numerical simulations to properly account for the impact of
halo shape changes on the inferred weak-lensing masses. We defer such an
analysis to future work.

Here we summarize the observable–mass scaling relations relevant for our
forecasts.

###  IV.1 tSZE $\zeta$–mass relation

We assume the following relation for the mean intrinsic tSZE galaxy cluster
detection significance $\zeta$

$\begin{split}\langle\ln\zeta\rangle=\ln{\zeta_{0}}&+{\zeta_{M}}\ln\left(\frac{M_{200\mathrm{c}}}{3\times
10^{14}\,h^{-1}M_{\odot}}\right)\\\
&+{\zeta_{z}}\ln\left(\frac{E(z)}{E(0.6)}\right)\,.\end{split}$ (26)

Here ${\zeta_{0}}$, ${\zeta_{M}}$ and ${\zeta_{z}}$ are the parameters
corresponding to the normalization, mass trend and redshift trend of the
relation. We assume that the intrinsic detection significance $\zeta$ scatters
around the mean relation in a log-normal fashion with a width described by
${\sigma_{\ln\zeta}}$.

The intrinsic detection significance is given for a survey of a fiducial
depth. To account for field-to-field variation in the survey noise level and
to use the cluster data from multiple fields of different depths to constrain
the underlying $\zeta$–mass relation, we rescale the normalization
${\zeta_{0}}$ of Eq. (26) using a factor $\gamma$; $\gamma=1$ corresponds to
the average depth of a specific field within the SPT-SZ survey [21, 63, 64].
Scaling up from the existing SPT-SZ and SPTpol cluster surveys, we adopt
$\gamma=3.5$ for the deep SPT-3G field and $\gamma=1.5$ for the shallower
field. For CMB-S4 we adopt a scale factor $\gamma=4$. This factor is similar
to the SPT-3G main field because at the depths of SPT-3G and CMB-S4 we expect
the noise to be dominated by the cosmic infrared background. We find that with
this scale factor, we recover similar numbers of clusters as in a previous
forecast [73].

The tSZE observed detection significance $\hat{\zeta}$ is related to the
intrinsic detection significance $\zeta$ by a normal distribution222Note that
in many past SPT analyses, the observed detection significance $\hat{\zeta}$
has been designated $\xi$. As we move to multi-observable analyses, we find it
simplest to create observed-intrinsic pairs of observables using a single
variable name with and without the hat.

$P(\hat{\zeta}|\zeta)=\mathcal{N}\left(\sqrt{\zeta^{2}+3},1\right)\,.$ (27)

The normal distribution is due to the Gaussian noise in the survey maps, and
the bias correction of 3 accounts for a noise bias introduced in the matched-
filter search for peaks within three dimensions: location on the sky (2
parameters) and effective core radius of the tSZE signature [74].

###  IV.2 Cluster richness $\lambda$–mass relation

The mean intrinsic richness $\lambda$–mass scaling relation is described by a
power law in mass and redshift

$\begin{split}\langle\ln\lambda\rangle=\ln{\lambda_{0}}&+{\lambda_{M}}\ln\left(\frac{M_{200\mathrm{c}}}{3\times
10^{14}\,h^{-1}M_{\odot}}\right)\\\
&+{\lambda_{z}}\ln\left(\frac{1+z}{1.6}\right)\,,\end{split}$ (28)

with ${\lambda_{0}}$, ${\lambda_{M}}$ and ${\lambda_{z}}$ describing the
normalization, mass trend and redshift trend. We assume that the cluster
intrinsic richness scatters around this relation in a log-normal fashion
described by its RMS variation ${\sigma_{\ln\lambda}}$.

The observed richness $\hat{\lambda}$ is related to the intrinsic richness
$\lambda$ by

$P(\hat{\lambda}|\lambda)=\mathcal{N}(\lambda,\sqrt{\lambda})\,,$ (29)

which posits the Poisson sampling noise associated with the realization of a
particular number of observed galaxies $\hat{\lambda}$ in a galaxy cluster
with an intrinsic richness $\lambda$. This expression is the Gaussian
approximation to Poisson noise that we assume to be valid for $\lambda\gtrsim
10$.

###  IV.3 Weak lensing $M_{\mathrm{WL}}$–mass relation

The six scaling-relation parameters and the two intrinsic scatter parameters
describing the $\zeta$–mass and $\lambda$–mass relations described above are
calibrated using weak-lensing data, which have well characterized and
controllable biases and uncertainties. The WL observable for each cluster is
the reduced tangential shear $\boldsymbol{g}_{\mathrm{t}}(R)$ or a collection
of shear profiles, each associated with a different tomographic bin of weak
lensing source galaxies. This observable is not employed in the cluster
selection, but it is used to extract a weak-lensing halo mass
$M_{\mathrm{WL}}$ by fitting the reduced shear profile to a Navarro-Frenk-
White profile [NFW; 75] in the radial range
$500\,h^{-1}\,\mathrm{kpc}<R<3.2/(1+z)\,h^{-1}\,\mathrm{Mpc}$. The lower bound
of the radial ranges ensures that the complex inner structure of the halo has
minimal impact on our analysis.

The derived WL mass $M_{\mathrm{WL}}$ will differ from the true halo mass,
taken here to be $M_{200\mathrm{c}}$, due to model uncertainties and
observational noise. The model uncertainties represent systematic
uncertainties that do not average down with increasing number of weak lensing
source galaxies. Therefore, we track these systematics by introducing an
additional observable–mass scaling relation relating the WL and halo masses
[76, 77, 78].

$\begin{split}&\left\langle\mathrm{ln}\left(\frac{M_{\mathrm{WL}}}{2\times
10^{14}\,h^{-1}M_{\odot}}\right)\right\rangle={\ln M_{\mathrm{WL}_{0}}}\\\
&\hskip
71.13188pt+{M_{\mathrm{WL}_{M}}}\mathrm{ln}\left(\frac{M_{200\mathrm{c}}}{2\times
10^{14}\,h^{-1}M_{\odot}}\right)\,,\end{split}$ (30)

where ${\ln M_{\mathrm{WL}_{0}}}$ is the logarithmic mass bias normalization
and ${M_{\mathrm{WL}_{M}}}$ is the mass trend in this bias. We discuss the
redshift dependence in Section V.2. In addition, the weak lensing mass
$M_{\mathrm{WL}}$ exhibits a mass dependent log-normal scatter about the mean
relation given by the variance $\sigma^{2}_{\ln\mathrm{WL}}$

$\begin{split}&\mathrm{ln}\,\sigma_{\ln\mathrm{WL}}^{2}={\ln\sigma^{2}_{\ln\mathrm{WL}_{0}}}\\\
&\hskip
56.9055pt+\sigma^{2}_{\ln\mathrm{WL}_{M}}\mathrm{ln}\left(\frac{M_{200\mathrm{c}}}{2\times
10^{14}\,h^{-1}M_{\odot}}\right)\,,\end{split}$ (31)

where ${\ln\sigma^{2}_{\ln\mathrm{WL}_{0}}}$ is the normalization and
$\sigma^{2}_{\ln\mathrm{WL}_{M}}$ is the mass trend. The determination of the
posterior distributions of the parameters of these two relations is discussed
in Section V.2.

##  V Generating mock data

The mock catalogs used in this work are created by drawing the data from our
model, and thus we consider all statistical and systematic uncertainties. The
model we have adopted is fully consistent with the recently analyzed
SPT$\times$DES dataset [41, 24], and therefore we expect it to be an excellent
baseline description of these future cluster and weak-lensing datasets. In the
following two sections, we describe in detail how we create the galaxy mock
catalogs and the follow-up next-generation weak-lensing data.

###  V.1  Mock cluster catalog

Mock tSZE catalogs in a fiducial cosmology are created by first computing the
HMF in the mass range $M_{200\mathrm{c}}\in[10^{13},10^{16}]\,h^{-1}M_{\odot}$
and in the redshift range of the given survey. We scale this HMF by the
appropriate redshift-dependent volume, creating a function containing the
expected number of halos of a given mass and redshift within mass and redshift
bins:

$\left<N(M,z)\right>\simeq\frac{\textnormal{d}n(\boldsymbol{p},M,z)}{\textnormal{d}M}\frac{\textnormal{d}V(\boldsymbol{p},z)}{\textnormal{d}z}\mathop{\textnormal{d}M}\mathop{\textnormal{d}z}\,,$
(32)

where the first factor is the HMF and the second factor is the survey solid
angle $\Omega_{\rm s}$ dependent differential volume. We then create clusters
of particular mass and redshift by drawing a Poisson realization of the
expected number of halos within each mass and redshift bin.

For each halo we assign a tSZE intrinsic detection significance $\zeta$ and an
intrinsic optical richness $\lambda$ using the observable–mass relations Eqs.
(26) and (28) together with the associated log-normal scatter. For the mock
catalog, the parameters of the observable–mass relations are fixed to the
values
$(\ln{\zeta_{0}},{\zeta_{M}},{\zeta_{z}},{\sigma_{\ln\zeta}})=(0.96,1.5,0.5,0.2)$
for the $\zeta$–mass relation and
$(\ln{\lambda_{0}},{\lambda_{M}},{\lambda_{z}},{\sigma_{\ln\lambda}})=(4.25,1.0,0.0,0.2)$
for the $\lambda$–mass relation. The observed tSZE detection significance
$\hat{\zeta}$ and the observed richness $\hat{\lambda}$ are then drawn using
the intrinsic values and the measurement and sampling noise described in Eqs.
(27) and (29), respectively.

Finally, we produce the cluster sample, modeling the selection by applying the
appropriate lower thresholds in the tSZE detection significance
$\hat{\zeta}_{\mathrm{min}}$ and observed richness
$\hat{\lambda}_{\mathrm{min}}(z)$ within the redshift range adopted for each
survey. The threshold value in tSZE detection significance and the redshift
range for each survey are presented in Secs. III.1.1 and III.1.2.

The minimum richness threshold as a function of redshift adopted here is
modeled on the MCMF follow-up method [60] and follows from the recent
SPT$\times$DES analysis [41, 24]. Because that analysis only contains the
minimum richness threshold up to $z=1.79$, we extrapolate using
$\hat{\lambda}_{\mathrm{min}}(z)$ equal to the value of
$\hat{\lambda}_{\mathrm{min}}$ at $z=1.79$ for all higher redshifts considered
in these forecasts.

We check the sensitivity of the cluster sample to different richness cuts by
assuming a constant $\hat{\lambda}_{\mathrm{min}}$ for all redshifts with
values of 10, 5 and 1. The number of clusters only varies by a few percent
among these samples, similar to the variation we see when creating mocks with
different random seeds. This is an indication that the tSZE significance
selection threshold $\hat{\zeta}_{\mathrm{min}}$ is dominating the selection,
and that the richness threshold has only a weak impact on the mock sampling
and therefore the exact values we adopt are not important for our forecasts.

To create $f(R)$ gravity mocks we use the HMF described in Eq. (24). Here the
number of clusters depends also on the strength of the gravity modification,
with more clusters being obtained for stronger $f(R)$ models, as demonstrated
in the enhancement of the HMF visible in Fig. 2. Figure 6 shows the
distribution and abundance of clusters in mass for different values of
$\log_{10}|f_{R0}|$ for our two different mock surveys, which we refer to as
SPT-3G$\times$ngWL and CMB-S4$\times$ngWL.

Figure 6: Mass distribution of the SPT-3G$\times$ngWL (above) and
CMB-S4$\times$ngWL (below) mock catalogs for different $f(R)$ gravity models
in yellow, light green and dark green as well as for the GR Planck 2018
cosmology in dark violet.

As we can see from Fig. 6, cluster catalogs with a GR cosmology and $f(R)$
gravity with $\log_{10}|f_{R0}|=-7$ have close to the same number of clusters.
This is in agreement with the difference in the HMF of these two models, which
is only a few percent in the lower mass range (see Fig. 2). Therefore, we
assume in our analysis that an $f(R)$ gravity model with
$\log_{10}|f_{R0}|=-7$ is indistinguishable from GR, and we choose the lower
bound of $\log_{10}|f_{R0}|$ in the emulators described in Sec. II.3 to be
$\log_{10}|f_{R0}|=-7$.

###  V.2  Mock ngWL data

Table 2: Priors on the weak-lensing mass to halo mass relation parameters and parameter uncertainties (see Eqs. (30) and (31)). The mean values are adopted from a recent SPT$\times$DES analysis [41], and the uncertainties are assumed to be two times smaller.  Parameter | Prior | Uncertainty | Prior
---|---|---|---
$\ln M_{\mathrm{WL},0}$ | -0.050 | $\Delta{\ln M_{\mathrm{WL}_{0}}}$ | 0.010
${M_{\mathrm{WL}_{M}}}$ | 1.029 | $\Delta{M_{\mathrm{WL}_{M}}}$ | 0.009
${\ln\sigma^{2}_{\ln\mathrm{WL}_{0}}}$ | -3.100 | $\Delta{\ln\sigma^{2}_{\ln\mathrm{WL}_{0}}}$ | 0.120
$\sigma^{2}_{\ln\mathrm{WL}_{M}}$ | -0.226 | $\Delta\sigma^{2}_{\ln\mathrm{WL}_{M}}$ | 0.300

We create mock cluster lensing data following the approach taken in the recent
SPT$\times$DES analysis [41]. The first step is to use Eqs. (30) and (31) to
assign a weak-lensing mass $M_{\mathrm{WL}}$ to each tSZE selected cluster in
the mock sample. Doing so ensures that we include the fact that the weak-
lensing mass is not the same as the halo mass. We then use this weak-lensing
mass and the corresponding NFW $\Delta\Sigma(R)$ to create a tangential shear
profile $\boldsymbol{g}_{\mathrm{t}}(R)$ for each background tomographic bin.
Source galaxy number densities and their redshift distribution are presented
in Section III.2.

In addition, we add cluster member contamination, consistent with the recent
measurement in DES data [41]. This process produces realistically noisy and
biased tangential shear profiles for each cluster. Specifically, these shear
profiles include all the known systematic and stochastic effects needed to
model cluster shear profiles in DES data.

The parameters describing the $M_{\mathrm{WL}}$-$M_{200\mathrm{c}}$ relations
(Eqs. (30) and (31)) are derived through the application of our mass
measurement technique to hydrodynamical simulations of clusters. Following a
previous work [78], the weak-lensing mass is extracted from hydrodynamical
simulations of clusters, and associated N-body simulations are used to
determine the corresponding true halo masses. Given both masses, the model for
the $M_{\rm WL}$–$M_{200\mathrm{c}}$ relation and its scatter is calibrated.
The true halo masses adopted need to match those that are used to calibrate
the HMF and its dependence on cosmology. In general, systematic uncertainties
in the shear and photometric redshift estimates of the source galaxies can
also play an important role in this weak-lensing to halo mass relation [41].

Crucial in this process of parameter estimation for the weak-lensing to halo
mass relation is that there are remaining uncertainties represented by the
parameter posteriors. These represent an effective systematic error floor in
our ability to estimate cluster halo masses using weak-lensing data. The
uncertainties on the bias and scatter parameters, $\Delta{\ln
M_{\mathrm{WL}_{0}}}$, $\Delta{M_{\mathrm{WL}_{M}}}$,
$\Delta{\ln\sigma^{2}_{\ln\mathrm{WL}_{0}}}$ and
$\Delta\sigma^{2}_{\ln\mathrm{WL}_{M}}$, are then marginalized over during the
mass calibration analysis to properly include the impact of weak lensing
systematic uncertainties on the cluster observable–mass relations.

Table 2 contains the parameter values of the weak-lensing mass to halo mass
relation that are adopted for the creation of the mock shear profiles and for
the analysis results presented below. We use as a baseline the results from
the DES weak-lensing calibration adopted in the recent SPT$\times$DES analysis
[41]. We assume no redshift dependence in any of the
$M_{\mathrm{WL}}$–$M_{200\mathrm{c}}$ scaling relation parameters in our
analysis here because that dependence is largely driven by uncertainties in
the photometric redshift biases in the DES shear catalog. The photometric
redshift requirements for both Euclid and Rubin are so tight that related
systematic uncertainties will be subdominant in an ngWL mass calibration
analyses.

Compared to the current state-of-the-art, the ngWL datasets will include
improvements in photometric redshift biases, shear measurement biases, and the
understanding of the baryonic physics impact on the cluster mass profiles.
Therefore, we assume that the uncertainties on the
$M_{\mathrm{WL}}$–$M_{200\mathrm{c}}$ parameters will improve by a factor of
two compared to DES. We view this as a relatively conservative choice for the
improvement of the WL uncertainties because in addition to better
characterization of the cluster population through simulations, we can alter
our analysis approach to reduce sensitivity to certain systematics (e.g.,
avoiding more of the cluster core where the baryonic effects are most
important and selecting background source galaxies conservatively to
dramatically reduce the cluster member contamination).

In Appendix A, we discuss the impact of adopting more optimistic parameter
uncertainties for ngWL surveys that are a factor of 10 smaller than those in
the recent DES analysis; interestingly, moving from 2 times to 10 times
reduction in systematic uncertainties has a minimal impact on the ngWL mass
calibration. This suggests that at the level of weak lensing systematics
adopted in our ngWL forecast, our observable–mass relation parameter
posteriors are shape noise dominated.

##  VI Likelihood and Analysis

The analysis method we employ in this work is based on the state-of-the-art
SPT$\times$DES cosmological analysis [41, 24]. We summarize in this section
the abundance and mass calibration likelihood and how we analyze our mock data
for the SPT-3G$\times$ngWL and CMB-S4$\times$ngWL sample.

###  VI.1  Abundance and mass calibration likelihood

The multi-observable likelihood used in this forecast is given by

$\ln\mathcal{L}(\boldsymbol{p})=\sum_{i}\ln\int_{\hat{\lambda}_{\mathrm{min}}}^{\infty}\textnormal{d}\hat{\lambda}\frac{\textnormal{d}^{3}N(\boldsymbol{p})}{\textnormal{d}\hat{\zeta}\,\textnormal{d}\hat{\lambda}\,\textnormal{d}z}\Big{|}_{\hat{\zeta}_{i},z_{i}}-\int_{z_{\mathrm{min}}}^{z_{\mathrm{max}}}\textnormal{d}z\int_{\hat{\zeta}_{\mathrm{min}}}^{\infty}\textnormal{d}\hat{\zeta}\int_{\hat{\lambda}_{\mathrm{min}}}^{\infty}\textnormal{d}\hat{\lambda}\frac{\textnormal{d}^{3}N(\boldsymbol{p})}{\textnormal{d}\hat{\zeta}\,\textnormal{d}\hat{\lambda}\,\textnormal{d}z}+\sum_{i}\ln\frac{\frac{\textnormal{d}^{4}N(\boldsymbol{p})}{\textnormal{d}\hat{\zeta}\,\textnormal{d}\hat{\lambda}\,\textnormal{d}\boldsymbol{g}_{\mathrm{t}}\textnormal{d}z}\Big{|}_{\hat{\zeta}_{i},\hat{\lambda}_{i},\boldsymbol{g}_{\mathrm{t},i},z_{i}}}{\int_{\hat{\lambda}_{\mathrm{min}}}^{\infty}\textnormal{d}\hat{\lambda}\frac{\textnormal{d}^{3}N(\boldsymbol{p})}{\textnormal{d}\hat{\zeta}\,\textnormal{d}\hat{\lambda}\,\textnormal{d}z}\Big{|}_{\hat{\zeta}_{i},z_{i}}}+\mathrm{const.}$
(33)

where both sums run over all clusters $i$.333More details and motivation of
the likelihood approach can be found in Ref. [41]. The vector $\boldsymbol{p}$
contains the cosmological and scaling relation parameters, and the observables
are tSZE detection significance $\hat{\zeta}$, richness $\hat{\lambda}$,
tangential shear profile $\boldsymbol{g}_{\mathrm{t}}$ and redshift $z$. The
first two terms of the above equation represent the Poisson likelihood
associated with the cluster abundance, which is independent of the weak-
lensing data, while the last term represents the information from the mass
calibration with the WL data.

The differential cluster number
$\frac{\textnormal{d}^{3}N}{\textnormal{d}\,\mathrm{obs}}$ that appears in the
first two terms is the differential halo observable function HOF in the
observable space $\hat{\zeta}-\hat{\lambda}-z$

$\begin{split}\frac{\textnormal{d}^{3}N(\boldsymbol{p})}{\textnormal{d}\hat{\zeta}\,\textnormal{d}\hat{\lambda}\,\textnormal{d}z}=\int\textnormal{d}\Omega_{\mathrm{s}}\iiint\textnormal{d}M\,\textnormal{d}\lambda\,\textnormal{d}\zeta\,P(\hat{\zeta}|\zeta)P(\hat{\lambda}|\lambda)P(\zeta,\lambda|M,z,\boldsymbol{p})\frac{\textnormal{d}^{2}N(M,z,\boldsymbol{p})}{\textnormal{d}M\,\textnormal{d}V}\frac{\textnormal{d}^{2}V(z,\boldsymbol{p})}{\textnormal{d}z\,\textnormal{d}\Omega_{\mathrm{s}}}\,.\end{split}$
(34)

Here $P(\hat{\zeta}|\zeta)$ and $P(\hat{\lambda}|\lambda)$ follow from Eqs.
(27) and (29), respectively, whereas $P(\zeta,\lambda|M,z,\boldsymbol{p})$ is
obtained from Eqs. (26) and (28), $\Omega_{\mathrm{s}}$ is the survey solid
angle, and the factors
$\frac{\textnormal{d}^{2}N(M,z,\boldsymbol{p})}{\textnormal{d}M\textnormal{d}z}$
and
$\frac{\textnormal{d}^{2}V(z,\boldsymbol{p})}{\textnormal{d}z\textnormal{d}\Omega_{\mathrm{s}}}$
are the HMF and the differential volume element for the corresponding
cosmology. The other differential cluster number in the observable space
$\hat{\zeta}-\hat{\lambda}-\boldsymbol{g}_{\mathrm{t}}-z$ is given by

$\begin{split}\frac{\textnormal{d}^{4}N(\boldsymbol{p})}{\textnormal{d}\hat{\zeta}\,\textnormal{d}\hat{\lambda}\,\textnormal{d}\boldsymbol{g}_{\mathrm{t}}\textnormal{d}z}=\int\textnormal{d}\Omega_{\mathrm{s}}\iiiint\textnormal{d}M\,\textnormal{d}\zeta\,\textnormal{d}\lambda\,\textnormal{d}M_{\mathrm{WL}}\,&P(\boldsymbol{g}_{\mathrm{t}}|M_{\mathrm{WL}},\boldsymbol{p})P(\hat{\zeta}|\zeta)P(\hat{\lambda}|\lambda)P(\zeta,\lambda,M_{\mathrm{WL}}|M,z,\boldsymbol{p})\frac{\textnormal{d}^{2}N(M,z,\boldsymbol{p})}{\textnormal{d}M\,\textnormal{d}V}\frac{\textnormal{d}^{2}V(z,\boldsymbol{p})}{\textnormal{d}z\,\textnormal{d}\Omega_{\mathrm{s}}}\,,\end{split}$
(35)

Where $P(\zeta,\lambda,M_{\mathrm{WL}}|M,z,\boldsymbol{p})$ follows from Eqs.
(26), (28) and (30) and
$P(\boldsymbol{g}_{\mathrm{t}}|M_{\mathrm{WL}},\boldsymbol{p})$ is given by
the product of Gaussian probabilities in each radial bin $i$ of the tangential
shear profiles

$\begin{split}P(\boldsymbol{g}_{\mathrm{t}}|M_{\mathrm{WL}},\boldsymbol{p})=\prod_{i}\left(\sqrt{2\pi}\Delta
g_{\mathrm{t},i}\right)^{-1}\exp\left[-\frac{1}{2}\left(\frac{g_{\mathrm{t},i}-g_{\mathrm{t},i}(M_{\mathrm{WL}},\boldsymbol{p})}{\Delta
g_{\mathrm{t},i}}\right)^{2}\right]\,,\end{split}$ (36)

with the shape noise $\Delta g_{\mathrm{t},i}$.

###  VI.2  Analysis

The analysis presented in this work is done with
CosmoSIS444https://cosmosis.readthedocs.io/[79] using the Multinest and
Nautilus samplers [80, 81]. In our analysis we separate the cluster abundance
and mass calibration elements into independent MCMC chains. The mass
calibration likelihood is first used to quantify posterior distributions on
the parameters of the observable–mass relations at a fixed cosmology.
Thereafter, we adopt these parameter constraints as Gaussian priors on the
parameters of the observable–mass relation parameters when evaluating
cosmological constraints with the cluster abundance likelihood. The advantage
of this approach is that it simplifies and dramatically speeds up the
likelihood calculation. The downside is that the cosmological sensitivity of
the mass calibration element (primarily due to $\Omega_{\mathrm{m}}$
sensitivity of the distance-redshift relation) is not correctly captured. As
discussed in Appendix B, this separation of the analysis has little impact on
our forecast cosmological or scaling relation parameter posteriors.

####  VI.2.1 Mass calibration

The main bottleneck in our mass calibration likelihood is the four-dimensional
convolution integral from Eq. (35), which has to be evaluated for each cluster
in our current implementation. In the calculation of this integral, we use an
efficient Monte-Carlo integration method [41]; however, the total number of
clusters is so large that it is not possible to complete a mass calibration
chain for either mock cluster sample in a reasonable amount of time. One
approach to speed up the calculation would be to adopt a stacked analysis of
cluster shear profiles [82], but for this analysis, we use an approximate
approach employing the likelihood in Eq. (33) to infer the observable–mass
relation parameter posteriors. In this approach, we use 1000 randomly selected
clusters from each tSZE survey to calculate the mass calibration likelihood
and then scale up the weights of each selected cluster by a factor $\alpha$
that accounts for the number of missing clusters from the full sample. To up-
weight the clusters we reduce the shape noise $\Delta g_{\mathrm{t}}$ by a
factor of $1/\sqrt{\alpha}$ where $\alpha$ is the rescaling factor to the
total number of clusters $N$, i. e. $N=\alpha 1000$. This impacts the
probability $P(\boldsymbol{g}_{\mathrm{t}}|M_{\mathrm{WL}},\boldsymbol{p})$
(defined in Eq. (36)) of observing a particular reduced shear profile given
the WL mass, appropriately rescaling the third term in the likelihood in Eq.
(33) to account for the full sample of clusters in the survey.

Because we assume that the mass calibration analysis is independent of the
underlying cosmology, we adopt fixed cosmological parameters that equal those
used in generating the mocks and uniform priors on the eight observable–mass
scaling relation parameters. The ranges for our flat priors are chosen such
that they are larger than the $5\sigma$ results of the observable mass
relations from previous analyses [22, 83]. Importantly, we account for the
uncertainties in the weak-lensing mass to halo mass relation (see Eqs. (30)
and (31)) by adopting mean parameter values and parameter uncertainties as
listed in Table 2 and then marginalizing over those uncertainties.

Table 3: Priors on the parameters of our cluster abundance analysis for the SPT-3G$\times$ngWL (second column) and CMB-S4$\times$ngWL (third column) cluster samples. For the scaling relation parameters, the mean values are the same as the mock inputs, and the uncertainties are sampled from the posteriors of the corresponding mass calibration MCMC chain (see Sec. VI.2.1). The prior on $\Omega_{\nu}h^{2}$ corresponds to a prior on the sum of neutrino masses $\sum m_{\nu}\sim\mathcal{U}(0,0.6)~{}\mathrm{eV}$ Parameter | SPT-3G$\times$ngWL | CMB-S4$\times$ngWL
---|---|---
Cosmology
$\Omega_{\mathrm{m}}$ | $\mathcal{U}(0.232,0.4)$ | $\mathcal{U}(0.232,0.4)$
$\Omega_{\nu}h^{2}$ | $\mathcal{U}(0,0.006)$ | $\mathcal{U}(0,0.00644)$
$h$ |  $\mathcal{U}(0.6,0.8)$ | $\mathcal{U}(0.6,0.8)$
$\ln(10^{10}A_{s})$ | $\mathcal{U}(1,4)$ | $\mathcal{U}(1,4)$
$n_{s}$ | $\mathcal{U}(0.94,1.)$ | $\mathcal{U}(0.94,1.)$
$\log_{10}|f_{R0}|$ | $\mathcal{U}(-7,-3)$ | $\mathcal{U}(-7,-3)$
tSZE $\zeta$-mass relation (Eqs. (26) and (27))
$\ln{\zeta_{0}}$ | $\mathcal{N}(0.96,0.03)$ | $\mathcal{N}(0.960,0.021)$
${\zeta_{M}}$ | $\mathcal{N}(1.50,0.04)$ | $\mathcal{N}(1.50,0.03)$
${\zeta_{z}}$ | $\mathcal{N}(0.50,0.17)$ | $\mathcal{N}(0.50,0.09)$
${\sigma_{\ln\zeta}}$ | $\mathcal{N}(0.200,0.026)$ | $\mathcal{N}(0.200,0.018)$
Richness $\lambda$-mass relation (Eqs. (28) and (29))
$\ln{\lambda_{0}}$ | $\mathcal{N}(4.250,0.019)$ | $\mathcal{N}(4.250,0.014)$
${\lambda_{M}}$ | $\mathcal{N}(1.00,0.03)$ | $\mathcal{N}(1.000,0.025)$
${\lambda_{z}}$ | $\mathcal{N}(0.00,0.16)$ | $\mathcal{N}(0.00,0.06)$
${\sigma_{\ln\lambda}}$ | $\mathcal{N}(0.200,0.012)$ | $\mathcal{N}(0.200,0.008)$

####  VI.2.2 Cluster abundance

We adopt the posterior parameter distributions for the observable–mass scaling
relation parameters that are listed in Table 3 as priors for the abundance
analysis. In the cluster abundance likelihood analysis (see Eq. (33)) we vary
the cosmological parameters $\Omega_{\mathrm{m}},\ \Omega_{\nu}h^{2},\ h,\
\mathrm{ln}(10^{10}A_{s}),\ n_{s}$, and $\log_{10}|f_{R0}|$. Note that
$\sigma_{8}^{\rm GR}$ is a derived parameter. It is important to note here
that the $\sigma_{8}^{\rm GR}$ quantity is the GR value as we use the linear
power spectrum in the corresponding GR cosmology to calculate the HMF (see
Sec. II.2).

The goal of this work is to test the constraining power from cluster data with
weak-lensing mass calibration alone, and thus we adopt flat priors on all
cosmological parameters. However, as mentioned in Sec. II.2, it is difficult
to distinguish between very weak modified gravity models and GR. Thus,
efficiently sampling the likelihood for these weak models is challenging.
Therefore, we combine the cluster abundance dataset with the primary CMB
Planck 2018 (TT, TE, EE) data for the analysis of GR and the
$\log_{10}|f_{R0}|=-7$ mock catalogs. This allows us to achieve convergence on
a significantly more reasonable timescale. Note that we do not include future
primary CMB results in this work. We account for the modified gravity
sensitivity of the Planck 2018 data by using the CMB power spectrum from
$f(R)$ gravity in this part of the likelihood. The cluster abundance data and
Planck 2018 data can be combined, because the mock catalogs are generated for
a Planck 2018 cosmology. Thus, there is no tension between the datasets.

Table 3 summarizes the priors on the cosmological and scaling relation
parameters for the SPT-3G$\times$ngWL and the CMB-S4$\times$ngWL forecasts.
Note that we do not apply the lower bound from oscillation experiments for the
massive neutrinos. The neutrino mass can only be constrained when combining
with primary CMB data, and thus the upper bound for the neutrino prior is
chosen based on the Planck 2018 results.

Note that when the mock input value of the $f(R)$ parameter lies on the lower
boundary of the $\log_{10}|f_{R0}|$ prior, the credibility limits are computed
from the lower boundary, and we present only upper bounds. We apply this to
the mock catalog with $\log_{10}|f_{R0}|=-7$ and the GR mock. Furthermore, the
GR limit with $|f_{R0}|=0$ is not reachable in a log prior and thus would add
an infinitely large volume below our lower bound. To avoid dependence of the
parameter upper limits on the choice of the lower boundary, we calculate in
the case the upper bounds in the linear $|f_{R0}|$ space by transforming the
parameter space from logarithmic to linear.

##  VII  Results

In this section, we summarize and discuss our results from the mass
calibration and cluster abundance analyses for the future SPT-3G$\times$ngWL
and CMB-S4$\times$ngWL datasets.

###  VII.1  ngWL mass calibration

We expect that SPT will detect around 6,000 clusters in a Planck 2018
cosmology with the SPT-3G camera. Applying the up-weighting approach described
in Sec. VI.2.1 with $\alpha=6$ for the mass calibration, we obtain constraints
for the uncertainties of the observable–mass scaling relation parameters as
shown in the second column of Table 3. For the $\sim 32,000$ CMB-S4 clusters,
we use the up-weighting approach with $\alpha=32$. The resulting parameter
constraints for CMB-S4$\times$ngWL are shown in the third column of Table 3.
We note that the mean recovered parameters are in agreement with the mock
inputs within the uncertainties, confirming the validity of our analysis
pipeline.

The tighter posteriors for the CMB-S4$\times$ngWL analysis are expected
because SPT-3G$\times$ngWL has roughly five times fewer clusters (and
therefore much less ngWL information). Furthermore, the CMB-S4 dataset
includes clusters at $z>0.1$, while for the SPT-3G dataset we adopt $z>0.25$
(see Sec. III.1). We see the largest improvement in the uncertainty on the
${\zeta_{z}}$ and ${\lambda_{z}}$ parameters, because these two parameters
model the redshift dependence of the observables at fixed mass.

We use the uncertainties from the mass calibration analyses presented here in
the cluster abundance analyses presented in the next section. Note that in the
cluster abundance analyses we set the mean parameter values to be equal to the
input values of the mocks. In addition, note that we have ignored the impact
of modified gravity on the halo profile in our mass calibration analysis for
this forecast, but it could be included in the analysis of the real datasets.

###  VII.2 Cosmology constraints

We now present the main results of the paper: the constraints on $f(R)$
modified gravity models for the future datasets SPT-3G$\times$ngWL and
CMB-S4$\times$ngWL. We analyze four different models with our cluster
abundance pipeline: values of $\log_{10}|f_{R0}|=\\{-7,-6,-5\\}$, denoted as
F7, F6, and F5, respectively, as well as GR. The F7 model is assumed to be
equivalent to GR in our framework (see Sec. V), and the comparison of the F7
and GR datasets allows us to test this assumption.

In the sections that follow we compare our results with the $95\,\%$ upper
bound constraint of $\log_{10}|f_{R0}|<-4.79$ found from clusters of galaxies
combined with CMB data, CMB lensing, baryon acoustic oscillations (BAO) and
type Ia supernova data [27]. This constraint is also comparable to the upper
end of the range of $f_{R0}$ constraints obtained from the Planck tSZE cluster
sample [28], where it is emphasized that uncertainties in the HMF prediction
they adopted contribute significantly to their final constraint on $f_{R0}$.
As discussed in previous sections, we have attempted to fold HMF uncertainties
into our constraints.

Figure 7: $\log_{10}|f_{R0}|$ posterior distributions from the cluster
abundance analysis for SPT-3G$\times$ngWL (red) and CMB-S4$\times$ngWL (blue)
datasets. For each survey one of the three independent mocks was chosen
randomly. For the analyses of the GR and F7 models, the Planck 2018 anisotropy
data were included.

####  VII.2.1 SPT-3G$\times$ngWL forecast

In this section we present our SPT-3G$\times$ngWL forecasts for the four
different models. Figure 7 shows the posterior distribution of the
$\log_{10}|f_{R0}|$ parameter of the four models tested. The full posteriors
of all cosmological parameters from the SPT-3G$\times$ngWL dataset are
presented in Appendix C. In fact, we have created three statistically
independent mocks for each model to test the sensitivity of our conclusions to
the statistical fluctuations that arise through the Poisson sampling of the
HMF and sources of scatter in the observable–mass relations. As expected, no
statistically significant difference is found between the three realizations
per model, and we hence only show the posterior for one of the mock
realizations here.

Table LABEL:tab:constraints_fR0 shows the constraints of the
$\log_{10}|f_{R0}|$ parameter for the four analyzed models. The results of the
$95\,\%$ upper bounds show very similar constraining power for the GR and F7
models and thus validate our choice of $\log_{10}|f_{R0}|=-7$ as a lower limit
within the context of the SPT-3G dataset. Our analysis shows that an
SPT-3G$\times$ngWL dataset can be used to distinguish GR and $f(R)$ modified
gravity models down to $\log_{10}|f_{R0}|<-5.97$ ($95\,\%$ upper bound).

Figure 8: Constraints on $\log_{10}|f_{R0}|$ and $\Omega_{\mathrm{m}}$ for the
F6 model (above) and F5 model (below) from the cluster abundance analysis for
SPT-3G$\times$ngWL (red) and CMB-S4$\times$ngWL (blue) datasets. The parameter
degeneracy that is clearly visible in the F6 panel arises due to enhanced
cluster numbers in stronger $f(R)$ models being offset partially by lower
$\Omega_{\mathrm{m}}$. In the F5 panel the degeneracy is still apparent at
some level in the SPT-3G$\times$ngWL case and only mild but with orthogonal
direction for the more constraining CMB-S4$\times$ngWL dataset. We tested the
convergence of MCMC chain by running multiple chains with different settings
and no shifts were seen.

The $95\,\%$ credibility constraints from the F6 mock show that this cluster
abundance and weak-lensing dataset without Planck 2018 data can differentiate
between GR and F6 models. Note, however, that models with a higher
$\log_{10}|f_{R0}|$ value predict more clusters when leaving all other
cosmological parameters unchanged (see Fig. 6) and these increased numbers
translate into tighter constraints, as expected. Moreover, the
SPT-3G$\times$ngWL dataset with Planck 2018 data improves over the constraint
from Ref. [27] by $25\,\%$. Finally, the analysis of the F5 model shows that
such modified gravity models could be distinguished from GR and from $f(R)$
gravity models with $\log_{10}|f_{R0}|<-5.97$ at $95\,\%$ credibility.

As seen in Fig. 7, our analysis shows an asymmetric posterior distribution for
$\log_{10}|f_{R0}|$ in the F5 model. An explanation for this asymmetric shape
is given by the fact that the HMFs of weak modified gravity models are harder
to distinguish because the differences are smaller. There is effectively less
information on the $f(R)$ gravity parameter encoded in the HMF as one moves to
weaker $f(R)$ models. Therefore, one expects a wider distribution or weaker
constraints when going to smaller values of $\mathrm{log}_{10}|f_{R0}|$.

Figure 8 shows the degeneracy between $\log_{10}|f_{R0}|$ and
$\Omega_{\mathrm{m}}$ for the F6 and F5 models. In the F6 model, the
degeneracy is mild and vanishes when approaching smaller values of
$\log_{10}|f_{R0}|$. The degeneracy between the two parameters is much
stronger in the F5 model.

In addition, the analysis of the F5 model shows a degeneracy between
$\log_{10}|f_{R0}|$ and $\sigma_{8}^{\rm GR}$ as seen in Fig. 9. This
degeneracy vanishes for low $\log_{10}|f_{R0}|$ values for the
SPT-3G$\times$ngWL dataset and the degeneracies are not seen in the analysis
of the GR, F7 and F6 models.

The degeneracies between these parameters have been previously noted in, e.g.,
studies of $f(R)$ gravity using the weak-lensing power spectrum [84]. These
degeneracies arise from the fact that the three quantities $\log_{10}|f_{R0}|$
$\Omega_{\mathrm{m}}$ and $\sigma_{8}^{\rm GR}$ all change the amplitude of
the HMF, enhancing the number of clusters for higher parameter values. The
degeneracies are broken for small $\log_{10}|f_{R0}|$ values because the
enhancement of the HMF in such $f(R)$ gravity models is very small to the
point that the HMF is difficult to distinguish from GR.

Figure 9: Joint constraints on $\log_{10}|f_{R0}|$ and $\sigma_{8}^{\rm GR}$
of the F5 model from the cluster abundance analysis for SPT-3G$\times$ngWL
(red) and CMB-S4$\times$ngWL (blue) datasets. The parameter degeneracy
indicates that the more numerous clusters in stronger $f(R)$ models can be
offset to some degree by a reduction in $\sigma_{8}^{\rm GR}$. This degeneracy
is not seen in the F6, F7 and GR models for both datasets. We tested the
convergence of MCMC chain by running multiple chains with different settings
and no shifts were seen.

As stated in Sec. II.2 there is a known degeneracy between the $f(R)$ gravity
parameter and massive neutrinos [see e. g. 52, 53, 54]. Our analysis of the
SPT-3G$\times$ngWL dataset on the other hand shows no degeneracy between these
parameters for all four examined models. This is owing to the fact that in our
analysis, the sum of neutrino masses is much smaller than the values
considered in Refs. [52, 53, 54], as we include the primary CMB constraint via
our prior on $\Omega_{\nu}h^{2}$, and the impact of such low neutrino masses
on the HMF is negligible.

####  VII.2.2 CMB-S4$\times$ngWL forecast

This section summarizes the results of our CMB-S4$\times$ngWL mock catalog
analysis. As with the SPT-3G$\times$ngWL analysis, we analyze three
statistically independent mocks and find no statistically significant
differences. Figure 7 shows the posterior distribution of the
$\log_{10}|f_{R0}|$ parameter of the four models for one mock realization in
blue. The full constraints in the cosmological parameters from the
CMB-S4$\times$ngWL dataset are presented in Appendix C

The $\log_{10}|_{R0}|$ posteriors for the four analyzed models are shown in
Table LABEL:tab:constraints_fR0. One can see that the $\log_{10}|f_{R0}|$
upper limit for the GR model is slightly tighter (roughly $5\,\%$) than for
the F7 model. This indicates that with the larger cluster sample of
CMB-S4$\times$ngWL the few percent difference in the HMF as seen in Fig. 2
starts to play a role. In an actual analysis of the CMB-S4 cluster sample, one
should thus extend the sampling range to lower values of $\log_{10}|f_{R0}|$.

Our results show that CMB-S4$\times$ngWL combined with Planck 2018 data can
rule out $f(R)$ gravity models down to $\log_{10}|f_{R0}|<-6.23$ at the
$95\,\%$ upper limit. These models are currently still viable based on current
published constraints. In fact, this forecast constraint represents an
improvement of more than an order of magnitude compared to previously
published cluster constraints [27]. Moreover, the analysis of the F6 model
shows that a CMB-S4$\times$ngWL dataset can distinguish $\log_{10}|f_{R0}|=-6$
cosmology from a GR cosmology at the 95% credible interval. The tightest
constraints from the CMB-S4$\times$ngWL mock dataset are for the F5 model
because stronger modified gravity models predict more clusters.

Table 4:  Parameter posteriors for $\log_{10}|f_{R0}|$ in the four analyzed models in the two different clusters surveys SPT-3G$\times$ngWL and CMB-S4$\times$ngWL. For models F7 and GR we present $95\,\%$ upper limits whereas mean and $68\,\%$ credible intervals are shown for F6 and F5 models. Model |  SPT-3G$\times$ngWL |  CMB-S4$\times$ngWL
---|---|---
GR | $<-5.97$ | $<-6.23$
F7 | $<-5.98$ | $<-5.93$
F6 | $-6.12\pm 0.44$ | $-6.18\pm 0.31$
F5 | $-5.09\pm 0.43$ | $-5.13\pm 0.17$

Figure 8 indicates a degeneracy between $\log_{10}|f_{R0}|$ and
$\Omega_{\mathrm{m}}$ in the F6 and F5 models. In the case of the F5 model
analysis, the degeneracy is mild and the direction is orthogonal to that from
the analysis of the F6 model. The change in the degeneracy direction compared
to the one in the F6 model can be explained by the behavior of the HMF at the
high mass end for $\log_{10}|f_{R0}|=-5$ $f(R)$ gravity models. The change in
the degeneracy is not seen in the F5 model for SPT-3G$\times$ngWL mock data
because the contours are wider and leak into the low $\log_{10}|f_{R0}|$
region where the degeneracy changes direction. Furthermore the F5 model
exhibits a degeneracy between $\log_{10}|f_{R0}|$ and $\sigma_{8}^{\rm GR}$ as
seen in Fig. 9. As in the SPT-3G$\times$ngWL analysis, the F7 and GR models do
not show these two degeneracies, because the enhancement of the HMF vanishes
for low values of $\log_{10}|f_{R0}|$ as discussed in Sec. VII.2.1. In
addition, the analysis of the CMB-S4$\times$ngWL datasets results in no
significant degeneracy between $\log_{10}|f_{R0}|$ and $\Omega_{\nu}h^{2}$ for
all models.

As we show in the last two sections, cluster data with weak-lensing informed
mass calibration and with additional Planck 2018 data in the weakest two
$f(R)$ models will improve upon the best current constraints from cluster plus
Planck data.

##  VIII Conclusions

In this work, we present a forecast of cluster abundance constraints with
weak-lensing informed mass calibration for a Hu & Sawicki $f(R)$ model in two
future surveys. In $f(R)$ gravity models, the clustering of matter is enhanced
compared to the GR model due to a fifth force mediated by an extra scalar
degree of freedom. Thus, the growth of structure is scale-dependent, and the
collapse threshold for halos is correspondingly mass-dependent. We account for
this in our analysis with a different HMF, which is given by the GR HMF scaled
by an enhancement factor $\mathcal{R}$ involving the critical overdensity in
$f(R)$ gravity $\delta_{\mathrm{crit}}$, and GR
$\delta_{\mathrm{crit},\mathrm{GR}}$ as well as its derivative with respect to
$\mathrm{ln}M$
$\textnormal{d}\delta_{\mathrm{crit}}/\textnormal{d}\mathrm{ln}M$ calculated
with a semi-analytical model.

The calculation of the critical overdensities and derivatives is
computationally expensive and thus is not feasible for a cosmological
analysis. Therefore, we build emulators for the three quantities needed to
calculate the enhancement factor for the HMF. We show that the emulation error
of the HMF is at the percent level and hence subdominant compared to other
limitations. Using the emulators speeds up our calculation by more than three
orders of magnitude.

Using the $f(R)$ gravity HMF we analyze galaxy cluster datasets, which are a
powerful cosmological probe able to distinguish between GR and $f(R)$ gravity
due to the enhanced structure formation in $f(R)$ gravity. We create mock
catalogs for the future tSZE cluster surveys representative of SPT-3G and
CMB-S4, and we include weak-lensing mass calibration information from next-
generation weak-lensing data from Euclid or Rubin. For each survey, we create
three $f(R)$ gravity models with $\log_{10}|f_{R0}|\in[-7,-6,-5]$ (F7, F6 and
F5 respectively) and a fourth GR model. We assume in our analysis that the F7
and GR models are indistinguishable because the HMF of the two models only
differs by a few percent in the low mass regime, which is not easily
distinguishable by either of these two future cluster surveys. In the case of
the F7 and GR models, we add the Planck 2018 CMB anisotropy data to improve
the constraints.

Our analysis consists of two steps. First, we carry out a weak-lensing mass
calibration analysis of the observable–mass relations at a fixed, fiducial
cosmology. Then we adopt the parameter posteriors from that analysis as priors
on a cluster abundance analysis with the goal of deriving cosmological
parameter posteriors. We test the robustness of our pipeline by analyzing
three statistically independent mock datasets for each of the models and find
no statistically significant difference. The analysis of the four models of an
SPT-3$\times$ngWL dataset gives the following results:

  * •

Cluster data from SPT-3G$\times$ngWL combined with primary CMB Planck 2018
data improves the current best constraints from cluster data to
$\log_{10}|f_{R0}|<-5.95$ at $95\,\%$ upper limit.

  * •

The analysis of the F7 and GR model verifies that the two models are
indistinguishable within the statistical uncertainties in the data.

  * •

Our analysis shows that the F6 model is distinguishable from GR at the
$2\sigma$ level.

  * •

Modified gravity models with $\log_{10}|f_{R0}|=-5$ are distinguishable from
GR and from $f(R)$ gravity models below $\log_{10}|f_{R0}|=-5.97$ at the
$2\sigma$ level.

  * •

We observe degeneracies between $\log_{10}|f_{R0}|$ and $\Omega_{\mathrm{m}}$
as well as $\sigma_{8}^{\rm GR}$. These are broken for values of $|f_{R0}|$
below $10^{-6}$ and $10^{-5}$ respectively.

The results of the CMB-S4$\times$ngWL mock data are the following:

  * •

CMB-S4$\times$ngWL data plus Planck 2018 data can rule out $f(R)$ gravity
models above $\log_{10}|f_{R0}|=-6.23$ at $95\,\%$ upper bound. They not only
improve upon the current best constraints from cluster data, but yield
constraints on $|f_{R0}|$ that are lower by a factor of $0.38$ than the one
from the SPT-3G$\times$ngWL dataset.

  * •

Modified gravity models with $\log_{10}|f_{R0}|=-6$ and $\log_{10}|f_{R0}|=-5$
are distinguishable from GR at the $95\,\%$ credible interval.

  * •

The analysis of the F6 and F5 models show degeneracy between
$\log_{10}|f_{R0}|$ and $\Omega_{\mathrm{m}}$ as well as
$\sigma_{8}^{\mathrm{GR}}$.

  * •

The degeneracy between $\log_{10}|f_{R0}|$ and $\sigma_{8}^{\rm GR}$ in the
analysis of the F5 model changes direction compared to the one in the F6
model. This is due to the different behavior of the HMF at the high mass end
for the F5 model.

Overall our analysis shows that upcoming tSZE-selected cluster samples of
thousands to tens of thousands of systems, combined with next-generation weak-
lensing survey data such as that from Euclid or Rubin, will enable
substantially improved constraints on modified gravity models. Furthermore,
the analysis shows that constraints on $f(R)$ gravity from cluster and WL
surveys are strongly competitive with other cosmological probes [27, 85, 86,
87, 88, 89].

###### Acknowledgements.

We thank Baojiu Li, César Hernández-Aguayo and Cheng-Zong Ruan for sharing
their FORGE HMF emulator with us. 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 - EXC-2094-390783311, by 2) the Max Planck Society Faculty
Fellowship program at MPE, and by 3) the Ludwig-Maximilians-Universität in
Munich.

##  Appendix A The impact of improved WL systematic uncertainties on the mass
calibration

Figure 10: Mass calibration analysis with two different choices on the weak-
lensing uncertainties (WL systematic floor) for the ngWL survey. The
conservative choice in red has a factor of two improvement in WL systematics
compared to DES, while the optimistic case is a factor of ten better than DES.
The posteriors on the observable mass scaling relation parameters are hardly
impacted, indicating that in neither case are the posteriors dominated by the
systematic uncertainties in the weak-lensing analysis. Figure 11: Mass
calibration analysis with two different cosmologies: Planck 2018 (green and
red) and WMAP7 (grey and blue). For each cosmology, we execute a chain with
$\Omega_{\mathrm{m}}$ fixed (red and blue) and free (green and grey). Dashed
lines indicate input values, for $\Omega_{\mathrm{m}}$ the input value is
represented in green for the Planck 2018 mock and in grey for the WMAP7 mock.
The consistency of the observable–mass scaling relation parameter posteriors
indicates that neither the exact input value of $\Omega_{\mathrm{m}}$ used in
the mocks nor fixing it or fitting $\Omega_{\mathrm{m}}$ have an important
impact.

In this appendix, we investigate the dependence of our results on the
assumptions made about systematics in the weak lensing mass estimates. Next-
generation weak-lensing surveys as expected from the Euclid satellite and the
Vera Rubin Observatory will collect a much higher amount of data compared to
current surveys like DES, KiDS and HSC SSP [90]. Therefore, much work will be
done to improve the systematic uncertainties to take advantage of the
increasing amount of data for the mass calibration analysis.

For photometric redshift and shear systematics, we simply adopt the
requirements for the ngWL surveys, which are dramatically tighter than what
has been achieved in the ongoing surveys. But, as discussed in Section V,
there are also systematics associated with how we use the observed source
galaxy shear profiles to determine the weak-lensing mass $M_{\mathrm{WL}}$,
and we model these in the so-called weak-lensing to halo mass relation (Eqs.
(30) and (31)). The exact uncertainties in the parameters of this relation are
estimated for the ngWL surveys based on the current DES results [41].

Figure 12: Posterior distributions of the cluster abundance analysis for
SPT-3G$\times$ngWL (red) and CMB-S4$\times$ngWL (blue) datasets. For the
analyses of the GR and F7 models, the Planck 2018 anisotropy data were
included.

We assume a baseline (conservative) approach where the uncertainties on the
parameters in Eqs. (30) and (31) would be improved by a factor of two as
compared to DES, and an optimistic scenario where the uncertainties would be
reduced by a factor of ten compared to DES.

The comparison of the conservative baseline and our optimistic scenario is
shown in Fig. 10. One can see that the reduction of the WL uncertainties
(corresponding to a kind of systematic floor) by a factor of five between our
conservative and optimistic scenarios does not significantly impact the
scaling relation parameter posteriors. We do see an improvement of around
$5\,\%$ in the constraints on the parameters $\ln{\zeta_{0}}$ and
$\ln{\lambda_{0}}$, which indicates that the systematic floor in our
conservative baseline is having some impact. Qualitatively, we expect this by
looking at the contours of Fig. 10 between the scaling relation parameters and
the WL uncertainties. Besides the contours of ${\ln M_{\mathrm{WL}_{0}}}$ and
$\ln{\zeta_{0}}$ or $\ln{\lambda_{0}}$ respectively, the contours are
horizontally orientated. This is an indication that an improvement of the WL
uncertainties, which only shrinks the contours in the horizontal directions,
will not lead to better constraints on the scaling parameters. In the case of
the two amplitude parameters, the contour plot shows a small degeneracy
between ${\ln M_{\mathrm{WL}_{0}}}$ and $\ln{\zeta_{0}}$ or
$\ln{\lambda_{0}}$, respectively, which helps explain the small improvement of
these two amplitude parameters when reducing the systematic floor by a factor
of five.

##  Appendix B Cosmology-dependence of mass calibration

As mentioned in Section VII.1 we assume that the ngWL based mass calibration
analysis can be separated from the SPT-3G and CMB-S4 abundance analysis
without introducing any important biases. This assumption is driven by our
need to avoid running the mass calibration, which is the bottleneck in the
current version of our code, many different times. We verify this statement by
running the mass calibration analysis with both fixed and free matter density
parameter $\Omega_{\mathrm{m}}$ with our standard Planck 2018 priors. Freedom
in $\Omega_{\mathrm{m}}$ has the most important cosmological impact on our
mass calibration analysis because the mapping from shear profiles to halo mass
depends on the distance-redshift relation. In addition, we analyze a GR mock
with a WMAP7 cosmology [91], i. e. $\Omega_{\mathrm{m}}=0.27$,
$\Omega_{\mathrm{b}}=0.0469$, $h=0.7$, $\mathrm{ln}(10^{10}A_{s})=3.155$ and
$n_{s}=0.95$ with $\Omega_{\mathrm{m}}$ fixed or free. This tests whether
there is sensitivity to the precise input value of the cosmological
parameters.

The results for the mass calibration for all four runs give the same
constraints on the observable–mass scaling relation parameters. Figure 11
includes the posteriors from these four runs (color coded as marked in the
figure). There is no tension among the four sets of posteriors, aside from the
$\Omega_{\mathrm{m}}$ posterior where the input Planck 2018 and WMAP7 values
are indeed different. This behavior motivates our assumption that the mass
calibration is independent of the cosmology and therefore validates our
approach of separating the mass calibration and cosmology chains in the
analyses we present here.

##  Appendix C Full results for the SPT-3G$\times$ngWL and CMB-S4$\times$ngWL
datasets

In Fig. 12 we show the constraints on all cosmological parameters for all four
examined models for the SPT-3G$\times$ngWL and CMB-S4$\times$ngWL datasets in
red and blue, respectively.

## References

  * Perlmutter et al. [1999] 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.
  * Riess et al. [1998] 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.
  * Joyce et al. [2016] A. Joyce, L. Lombriser, and F. Schmidt, Annual Review of Nuclear and Particle Science 66, 95 (2016), arXiv:1601.06133.
  * Koyama [2018] K. Koyama, International Journal of Modern Physics D 27, 1848001 (2018).
  * Baker et al. [2019] 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.
  * Buchdahl [1970] H. A. Buchdahl, MNRAS 150, 1 (1970).
  * Khoury and Weltman [2004] J. Khoury and A. Weltman, Phys. Rev. D 69, 044026 (2004), arXiv:astro-ph/0309411.
  * Burrage and Sakstein [2018] C. Burrage and J. Sakstein, Living Reviews in Relativity 21, 1 (2018), arXiv:1709.09071.
  * Hu and Sawicki [2007] W. Hu and I. Sawicki, Phys. Rev. D 76, 064004 (2007), URL https://link.aps.org/doi/10.1103/PhysRevD.76.064004.
  * Oyaizu [2008] H. Oyaizu, Phys. Rev. D 78, 123523 (2008), arXiv:0807.2449.
  * Li et al. [2012] B. Li, G.-B. Zhao, R. Teyssier, and K. Koyama, JCAP 2012, 051 (2012), arXiv:1110.1379.
  * Cataneo et al. [2016] M. Cataneo, D. Rapetti, L. Lombriser, and B. Li, JCAP 2016, 024 (2016), arXiv:1607.08788.
  * Arnold et al. [2022] C. Arnold, B. Li, B. Giblin, J. Harnois-Déraps, and Y.-C. Cai, MNRAS 515, 4161 (2022), arXiv:2109.04984.
  * Lombriser et al. [2013] L. Lombriser, B. Li, K. Koyama, and G.-B. Zhao, Phys. Rev. D 87, 123511 (2013), arXiv:1304.6395.
  * Lombriser et al. [2014] L. Lombriser, K. Koyama, and B. Li, JCAP 2014, 021 (2014), arXiv:1312.1292.
  * Wang and Steinhardt [1998] L. Wang and P. J. Steinhardt, Astrophys. J. 508, 483 (1998), arXiv:astro-ph/9804015.
  * Haiman et al. [2001] Z. Haiman, J. J. Mohr, and G. P. Holder, Astrophys. J. 553, 545 (2001), arXiv:astro-ph/0002336.
  * Vikhlinin et al. [2009] 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.
  * Benson et al. [2013] 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.
  * Bocquet et al. [2015] 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.
  * de Haan et al. [2016] 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.
  * Bocquet et al. [2019] S. Bocquet, J. P. Dietrich, T. Schrabback, L. E. Bleem, M. Klein, S. W. Allen, D. E. Applegate, M. L. N. Ashby, M. Bautz, M. Bayliss, et al., Astrophys. J. 878, 55 (2019), arXiv:1812.01679.
  * Abbott et al. [2020] 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.
  * Bocquet et al. [2024] S. Bocquet, S. Grandis, L. E. Bleem, M. Klein, J. J. Mohr, T. Schrabback, T. M. C. Abbott, P. A. R. Ade, M. Aguena, A. Alarcon, et al., arXiv e-prints arXiv:2401.02075 (2024), arXiv:2401.02075.
  * Schmidt et al. [2009a] F. Schmidt, A. Vikhlinin, and W. Hu, Phys. Rev. D 80, 083505 (2009a), arXiv:0908.2457.
  * Lombriser et al. [2012] L. Lombriser, A. Slosar, U. Seljak, and W. Hu, Phys. Rev. D 85, 124038 (2012), arXiv:1003.3009.
  * Cataneo et al. [2015] 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.
  * Peirone et al. [2017] 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.
  * Li and Efstathiou [2012] B. Li and G. Efstathiou, MNRAS 421, 1431 (2012), arXiv:1110.6440.
  * Sunyaev and Zeldovich [1972] R. A. Sunyaev and Y. B. Zeldovich, Comments on Astrophysics and Space Physics 4, 173 (1972).
  * Benson et al. [2014] 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.
  * Carlstrom et al. [2011] 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.
  * Ade et al. [2019] 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.
  * Abazajian et al. [2019] 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.
  * Raghunathan et al. [2022] 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.
  * Laureijs et al. [2011] R. Laureijs, J. Amiaux, S. Arduini, J. L. Auguères, J. Brinchmann, R. Cole, M. Cropper, C. Dabin, L. Duvet, A. Ealet, et al., arXiv e-prints arXiv:1110.3193 (2011), arXiv:1110.3193.
  * Euclid Collaboration et al. [2022] 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.
  * Ivezic et al. [2008] 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).
  * The LSST Dark Energy Science Collaboration et al. [2018] The LSST Dark Energy Science Collaboration, R. Mandelbaum, T. Eifler, R. Hložek, T. Collett, E. Gawiser, D. Scolnic, D. Alonso, H. Awan, R. Biswas, et al., arXiv e-prints arXiv:1809.01669 (2018), arXiv:1809.01669.
  * Schmidt [2010] F. Schmidt, Phys. Rev. D 81, 103002 (2010), arXiv:1003.0409.
  * Bocquet et al. [2023] 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., arXiv e-prints arXiv:2310.12213 (2023), arXiv:2310.12213.
  * Ferraro et al. [2011] S. Ferraro, F. Schmidt, and W. Hu, Phys. Rev. D 83, 063503 (2011), arXiv:1011.0992.
  * Schmidt et al. [2009b] F. Schmidt, M. Lima, H. Oyaizu, and W. Hu, Phys. Rev. D 79, 083518 (2009b), arXiv:0812.0545.
  * Li and Lam [2012] B. Li and T. Y. Lam, MNRAS 425, 730 (2012), arXiv:1205.0058.
  * Planck Collaboration et al. [2020] 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.
  * Shandera et al. [2013] S. Shandera, A. Mantz, D. Rapetti, and S. W. Allen, JCAP 2013, 004 (2013), arXiv:1304.1216.
  * Tinker et al. [2008] J. Tinker, A. V. Kravtsov, A. Klypin, K. Abazajian, M. Warren, G. Yepes, S. Gottlöber, and D. E. Holz, Astrophys. J. 688, 709 (2008), arXiv:0803.2706.
  * Sheth and Tormen [1999] R. K. Sheth and G. Tormen, MNRAS 308, 119 (1999), arXiv:astro-ph/9901122.
  * Despali et al. [2016] G. Despali, C. Giocoli, R. E. Angulo, G. Tormen, R. K. Sheth, G. Baso, and L. Moscardini, MNRAS 456, 2486 (2016), arXiv:1507.05627.
  * Bryan and Norman [1998] G. L. Bryan and M. L. Norman, Astrophys. J. 495, 80 (1998), arXiv:astro-ph/9710107.
  * Ichiki and Takada [2012] K. Ichiki and M. Takada, Phys. Rev. D 85, 063521 (2012), URL https://link.aps.org/doi/10.1103/PhysRevD.85.063521.
  * Motohashi et al. [2013] H. Motohashi, A. A. Starobinsky, and J. Yokoyama, Phys. Rev. Lett. 110, 121302 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.110.121302.
  * Baldi et al. [2014] M. Baldi, F. Villaescusa-Navarro, M. Viel, E. Puchwein, V. Springel, and L. Moscardini, MNRAS 440, 75 (2014), arXiv:1311.2588.
  * Wright et al. [2019] B. S. Wright, K. Koyama, H. A. Winther, and G.-B. Zhao, JCAP 2019, 040 (2019), arXiv:1902.10692.
  * Ruan et al. [2024] C.-Z. Ruan, C. Cuesta-Lazaro, A. Eggemeier, B. Li, C. M. Baugh, C. Arnold, S. Bose, C. Hernández-Aguayo, P. Zarrouk, and C. T. Davies, MNRAS 527, 2490 (2024), arXiv:2301.02970.
  * Melin et al. [2006] J. B. Melin, J. G. Bartlett, and J. Delabrouille, AnA 459, 341 (2006), arXiv:astro-ph/0602424.
  * Staniszewski et al. [2009] Z. Staniszewski, P. A. R. Ade, K. A. Aird, B. A. Benson, L. E. Bleem, J. E. Carlstrom, C. L. Chang, H. M. Cho, T. M. Crawford, A. T. Crites, et al., Astrophys. J. 701, 32 (2009), arXiv:0810.1578.
  * Song et al. [2012] J. Song, A. Zenteno, B. Stalder, S. Desai, L. E. Bleem, K. A. Aird, R. Armstrong, M. L. N. Ashby, M. Bayliss, G. Bazin, et al., Astrophys. J. 761, 22 (2012), arXiv:1207.4369.
  * Bleem et al. [2015] 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.
  * Klein et al. [2023a] 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., arXiv e-prints arXiv:2309.09908 (2023a), arXiv:2309.09908.
  * Klein et al. [2018] M. Klein, J. J. Mohr, S. Desai, H. Israel, S. Allam, A. Benoit-Lévy, D. Brooks, E. Buckley-Geer, A. Carnero Rosell, M. Carrasco Kind, et al., MNRAS 474, 3324 (2018), arXiv:1706.06577.
  * Klein et al. [2023b] 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., arXiv e-prints arXiv:2309.09908 (2023b), arXiv:2309.09908.
  * Bleem et al. [2020] 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.
  * Bleem et al. [2023] L. E. Bleem, M. Klein, T. M. C. Abbott, P. A. R. Ade, M. Aguena, O. Alves, A. J. Anderson, F. Andrade-Oliveira, B. Ansarinejad, M. Archipley, et al., arXiv e-prints arXiv:2311.07512 (2023), arXiv:2311.07512.
  * Mohr and Evrard [1997] J. J. Mohr and A. E. Evrard, Astrophys. J. 491, 38 (1997), arXiv:astro-ph/9707184.
  * Mohr et al. [1999] J. J. Mohr, B. Mathiesen, and A. E. Evrard, Astrophys. J. 517, 627 (1999), arXiv:astro-ph/9901281.
  * Hoekstra et al. [2015] H. Hoekstra, R. Herbonnet, A. Muzzin, A. Babul, A. Mahdavi, M. Viola, and M. Cacciato, MNRAS 449, 685 (2015), arXiv:1502.01883.
  * Mantz et al. [2016] 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.
  * Chiu et al. [2018] I. Chiu, J. J. Mohr, M. McDonald, S. Bocquet, S. Desai, M. Klein, H. Israel, M. L. N. Ashby, A. Stanford, B. A. Benson, et al., MNRAS 478, 3072 (2018), arXiv:1711.00917.
  * Chiu et al. [2022] 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.
  * Sotiriou and Faraoni [2010] T. P. Sotiriou and V. Faraoni, Reviews of Modern Physics 82, 451 (2010), arXiv:0805.1726.
  * Zhang [2007] P. Zhang, Phys. Rev. D 76, 024007 (2007), arXiv:astro-ph/0701662.
  * Raghunathan [2022] S. Raghunathan, Astrophys. J. 928, 16 (2022), arXiv:2112.07656.
  * Vanderlinde et al. [2010] 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.
  * Navarro et al. [1997] J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 490, 493 (1997), arXiv:astro-ph/9611107.
  * Becker and Kravtsov [2011] M. R. Becker and A. V. Kravtsov, Astrophys. J. 740, 25 (2011), arXiv:1011.1681.
  * Dietrich et al. [2019] J. P. Dietrich, S. Bocquet, T. Schrabback, D. Applegate, H. Hoekstra, S. Grandis, J. J. Mohr, S. W. Allen, M. B. Bayliss, B. A. Benson, et al., MNRAS 483, 2871 (2019), arXiv:1711.05344.
  * Grandis et al. [2021] S. Grandis, S. Bocquet, J. J. Mohr, M. Klein, and K. Dolag, MNRAS 507, 5671 (2021), arXiv:2103.16212.
  * Zuntz et al. [2015] 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.
  * Feroz et al. [2009] F. Feroz, M. P. Hobson, and M. Bridges, MNRAS 398, 1601 (2009), arXiv:0809.3437.
  * Lange [2023] J. U. Lange, MNRAS 525, 3181 (2023), arXiv:2306.16923.
  * Singh et al. [2024] A. Singh, C. Davies, J. Mohr, S. Bocquet, M. Klein, and D. Hernandez-Lang, in preparation (2024).
  * Saro et al. [2015] A. Saro, S. Bocquet, E. Rozo, B. A. Benson, J. Mohr, E. S. Rykoff, M. Soares-Santos, L. Bleem, S. Dodelson, P. Melchior, et al., MNRAS 454, 2305 (2015), arXiv:1506.07814.
  * Harnois-Déraps et al. [2023] J. Harnois-Déraps, C. Hernandez-Aguayo, C. Cuesta-Lazaro, C. Arnold, B. Li, C. T. Davies, and Y.-C. Cai, MNRAS 525, 6336 (2023), arXiv:2211.05779.
  * Hojjati et al. [2016] A. Hojjati, A. Plahn, A. Zucca, L. Pogosian, P. Brax, A.-C. Davis, and G.-B. Zhao, Phys. Rev. D 93, 043531 (2016), arXiv:1511.05962.
  * Hu et al. [2016] B. Hu, M. Raveri, M. Rizzato, and A. Silvestri, MNRAS 459, 3880 (2016), arXiv:1601.07536.
  * Pratten et al. [2016] G. Pratten, D. Munshi, P. Valageas, and P. Brax, Phys. Rev. D 93, 103524 (2016), arXiv:1602.06711.
  * Yamamoto et al. [2010] K. Yamamoto, G. Nakamura, G. Hütsi, T. Narikawa, and T. Sato, Phys. Rev. D 81, 103517 (2010), arXiv:1004.3231.
  * Stark et al. [2016] A. Stark, C. J. Miller, N. Kern, D. Gifford, G.-B. Zhao, B. Li, K. Koyama, and R. C. Nichol, Phys. Rev. D 93, 084036 (2016), arXiv:1603.00056.
  * Aihara et al. [2018] H. Aihara, N. Arimoto, R. Armstrong, S. Arnouts, N. A. Bahcall, S. Bickerton, J. Bosch, K. Bundy, P. L. Capak, J. H. H. Chan, et al., PASJ 70, S4 (2018), arXiv:1704.05858.
  * Komatsu et al. [2011] E. Komatsu, K. M. Smith, J. Dunkley, C. L. Bennett, B. Gold, G. Hinshaw, N. Jarosik, D. Larson, M. R. Nolta, L. Page, et al., ApJ 192, 18 (2011), arXiv:1001.4538.