The Redshift-Space Cluster-Galaxy Cross-Correlation as a Probe of Gravity on Cosmic Scales

Introduction

After three decades of intensive study, the origin of cosmic acceleration remains the most profound mystery of contemporary cosmology. Does accelerating expansion signal a breakdown of General Relativity (GR) on cosmological scales, or is it caused by an exotic form of energy that exerts repulsive gravity within GR? The example of a cosmological constant shows that the answer to this question may be ambiguous, as $\Lambda g_{\mu\nu}$ can be added to the curvature or stress-energy side of the GR field equation [1, 2]. However, many modifications of GR give rise to distinctive gravitational effects in addition to changing the expansion equation, such as growth of structure at a rate that depends on the redshift, scale, or environment (see reviews in refs. [3, 4]). Furthermore, the bending of light and acceleration of non-relativistic particles depend on distinct gravitational potentials that are equal at late times in GR but may differ in modified gravity theories, so that mass estimates based on gravitational lensing and peculiar velocities can yield incompatible results [5, 6]. Building on previous studies by Zu & Weinberg [7] and Zu et al. [8], we investigate the cluster-galaxy redshift-space cross-correlation function $\xi^{s}_{\mathrm{cg}}(r_{\mathrm{p}}, r_{\pi})$ as a probe of modified gravity and a test of GR.

Many modified gravity theories take the form of tensor-scalar models, in which a scalar field sources an acceleration acting as a “fifth force” alongside the standard gravitational interaction of the tensor field [e.g., 3, 4, 9, 10]. Since GR makes extremely accurate predictions on solar system and stellar scales, including binary pulsars and merging black holes, any scalar field contribution must be strongly suppressed in these environments. In empirically viable theories, this suppression emerges naturally from non-linear scalar field interactions, giving rise to screening mechanisms that operate in regions of high potential or high density [11, 12]. As the universe evolves, screening is effective at early times but weakens in regions near the cosmic mean density at late times, where gravitational accelerations can depart from GR predictions at order unity and produce substantial deviations in galaxy peculiar velocities. In $f(R)$ gravity, for instance, where the Ricci scalar $R$ in the Einstein-Hilbert action is replaced by $R + f(R)$, accelerations in unscreened regions are enhanced by a factor of 4/3 [13]. The same enhancement appears in the unscreened, static regime of the normal branch of the DGP model [14, 15], where modification of gravity on cosmic scales arises from the presence of an additional spatial dimension.

The infall regions around massive galaxy clusters offer a novel opportunity to test the GR-predicted velocity scale. Within the virial radius, the velocity dispersion of member galaxies produces the well-known Finger-of-God (FoG) distortion, stretching clusters along the line of sight in redshift space. At larger distances, coherent infall produces a squashing distortion, including caustics near the turnaround radius, where infall velocities cancel the Hubble flow [16, 17, 18, 19, 20]. Zu & Weinberg [7] developed an analytic model for $\xi^{s}_{\mathrm{cg}}(r_{\mathrm{p}}, r_{\pi})$, and used it to extract the characteristic radial velocity profile $v_{\mathrm{c}}(r)$ as a function of the three-dimensional radius, $r = (r_{\mathrm{p}}^{2} + r_{\pi}^{2})^{1/2}$. Applying this model to group catalogs [21] from the Sloan Digital Sky Survey (SDSS; [22, 23]), they recovered $v_{\mathrm{c}}(r)$ profiles in good agreement with predictions from the MILLENNIUM $N$-body simulation [24], with maximum infall velocities of 600–800 $\mathrm{km}\,\mathrm{s}^{-1}$.

Zu et al. [8] applied the same analysis to $N$-body simulations of $f(R)$ gravity and a Galileon model, which follow the chameleon [12] and Vainshtein [11] screening mechanisms, respectively. In both cases, the characteristic infall velocities $v_{\mathrm{c}}(r)$ were enhanced by $\sim50$–150 $\mathrm{km}\,\mathrm{s}^{-1}$ relative to the GR predictions for cluster populations of equal space density, attributable to the higher gravitational accelerations in unscreened environments at large distances from clusters. However, these enhancements of 10–30% were comparable in magnitude to the uncertainties in the SDSS analysis of Zu & Weinberg [7], preventing a definitive test of the modified gravity models.

This work revisits the analysis of Zu & Weinberg [7] and Zu et al. [8] for two reasons. First, the Dark Energy Spectroscopic Instrument [DESI; 25] Bright Galaxy Survey (BGS; [26]) provides a galaxy redshift survey ten times larger than the SDSS main galaxy survey at a similar number density [27], substantially reducing the statistical uncertainties that limited the earlier analysis. Second, computational advances now make it feasible to model $\xi^{s}_{\mathrm{cg}}(r_p,r_\pi)$ directly within a parameter inference framework, bypassing the need for the analytic model of Zu & Weinberg [7]. Recent studies of galaxy-galaxy lensing, galaxy redshift-space correlations, and cluster weak lensing have increasingly relied on simulation-based halo modeling, using emulators [28, 29, 30, 31, 32, 33], subhalo abundance matching [34, 35, 36], or direct halo occupation modeling of simulations [37, 38, 39]. We plan to construct an emulator for $\xi^{s}_{\mathrm{cg}}(r_p,r_\pi)$ in future work.

For this paper, we evaluate derivatives of observables from perturbations around a fiducial model and use these to carry out a Fisher matrix forecast, illustrating the general principles of our approach and the sensitivity achievable with DESI BGS.

In this work, we use GR simulations from the AbacusSummit $N$-body suite [40, 41], adopting a simple parameterization to capture modified gravity effects by scaling all halo velocities by a factor $F_G$, where $F_G = 1$ recovers the standard GR prediction. This choice is motivated by the linear growth rate parameter $f(\Omega_{\mathrm{m}}) \approx\Omega_{\mathrm{m}}^{0.55}$, which scales peculiar velocities relative to displacements in linear perturbation theory. Following the results of Zu et al. [8], applying this methodology to modified gravity simulations is expected to yield best-fit values of $F_G > 1$, with 1.1–1.3 typical for the models examined therein. We defer analysis of modified gravity simulations to future work. While our parameterization does not map onto a specific modified gravity theory, our forecast addresses a well-defined question: with what precision can a survey like the DESI BGS test the GR-predicted scale of peculiar velocities in cluster infall regions?

We focus on the $z = 0.3$ snapshots of the AbacusSummit simulations and present forecasts scaled to the DESI DR3 BGS volume in the redshift range $z = 0.25$–$0.35$. Redshift evolution across $z$-bins and the potential gains from analysing samples at lower and higher redshifts are left to future work. We adopt a conventional but idealized cluster selection model, in which each observed cluster is drawn from a population of massive halos whose observational mass estimates are related to their true three-dimensional masses by log-normal scatter. Anisotropic selection biases are known to affect predicted weak lensing profiles significantly [42, 43]. Analogous biases are expected to be important for $\xi^{s}_{\mathrm{cg}}(r_p,r_\pi)$, even if clusters are identified from spectroscopic galaxy catalogs, since galaxies within $\sim1000\ \mathrm{km\ s^{-1}}$ velocity separations can reside within the cluster or in the foreground or background [42, 44, 45]. Experience from cluster weak lensing analyses demonstrates that full forward modeling of cluster selection within an $N$-body simulation framework can mitigate or eliminate such biases, enabling robust cosmological inference [46, 47]. We anticipate that a similar framework will be required for accurate modeling of $\xi^{s}_{\mathrm{cg}}(r_p,r_\pi)$ with DESI BGS, but we defer its development to future work.

With the high precision potentially achievable by DESI BGS, cosmological constraints from $\xi^{s}_{\mathrm{cg}}(r_p,r_\pi)$ can complement those from cluster weak lensing in large photometric surveys (e.g., refs. [46, 47, 48, 49, 50, 51, 52]). Like stacked weak lensing profiles, redshift-space distortions of $\xi^{s}_{\mathrm{cg}}(r_p,r_\pi)$ probe the virial masses and extended mass profiles of cluster halos. In combination with cluster number counts, these measurements constrain the mean matter density $\Omega_{\mathrm{m}}$ and the amplitude of the matter power spectrum $\sigma_8$ at the typical cluster mass scale [53]. We give preliminary estimates of the sensitivity to $\sigma_8$ and $\Omega_{\mathrm{m}}$ in our Fisher forecasts, and we find only moderate correlation between the constraints on these parameters and those on $F_G$.

Ultimately, we aim to combine $\xi^{s}_{\mathrm{cg}}(r_p,r_\pi)$ with weak lensing to test the consistency of the gravitational potentials that govern light bending and galaxy motions. Light bending and galaxy motions respond differently to the two scalar potentials that describe the space and time components of the perturbed Friedmann-Robertson-Walker metric (see, e.g., ref. [5]). These potentials are equal in magnitude in GR¹ but can differ in modified gravity models. The $E_{\mathrm{G}}$ statistic, also known as the gravitational slip test, combines galaxy-galaxy lensing and galaxy redshift-space distortions to test the consistency of these potentials [5, 6, 54]. Our $F_{\mathrm{G}}$ parameterization is motivated partly by analogy to $E_{\mathrm{G}}$, and future joint analyses of $\xi_{\mathrm{cg}}^{\mathrm{s}}(r_p,r_\pi)$ and cluster weak lensing could carry out an analogous test in the physically distinct environment of cluster infall regions.

The remainder of this paper is organized as follows. § 2 describes the AbacusSummit simulations and mock galaxy and cluster catalogs. § 3 characterizes the mean radial infall velocity profile and its dependence on each model parameter. § 4 introduces $\xi_{\mathrm{cg}}^{\mathrm{s}}(r_p,r_\pi)$ and our measurement methodology. § 5 quantifies the sensitivity of the full data vector to each parameter through derivatives of $\log\xi_{\mathrm{cg}}^{\mathrm{s}}$, $\log n_{\mathrm{gal}}$, and $\log n_{\mathrm{halo}}$. § 6 presents the Fisher matrix formalism and the Full Forecast constraints. § 7 examines the information content of the data vector, the role of the infall velocity field in driving the constraints, and the steps required to apply this methodology to DESI BGS data. § 8 summarizes our main conclusions.

Simulations and Mock Catalogs

We construct our mock galaxy and cluster catalogs from the AbacusSummit $N$-body simulation suite, which provides the large-scale dark matter distribution needed to model the cluster infall regime. To populate these simulations, we model the galaxy–halo connection using a Halo Occupation Distribution (HOD) framework. We also account for the stochastic mapping between observable cluster proxies and true halo mass through a log-normal mass–observable relation. From these catalogs, we compute the clustering ($\log\xi_{\mathrm{cg}}^{\mathrm{s}}$) and number density ($\log n_{\mathrm{gal}}$ and $\log n_{\mathrm{halo}}$) observables that form our data vector.

AbacusSummit Simulations

AbacusSummit [40, 41] is a suite of large volume, high-accuracy cosmological $N$-body simulations designed to meet the strict modeling requirements of DESI [55]. The standard boxes contain $6912^3$ particles in a $(2\ h^{-1}\ \mathrm{Gpc})^3$ volume, corresponding to a particle mass resolution of $m_p \approx2 \times10^9\ h^{-1}\ M_\odot$, sufficient to resolve the $M \gtrsim10^{12}\ h^{-1}\ M_\odot$ halos that host typical bright galaxies, which serve as tracers of cluster infall regions.

We adopt the Planck cosmology [56] as our fiducial cosmology. We use the 25 independent phase realizations of the fiducial cosmology (c000) to estimate the covariance of $\log\xi_{\mathrm{scg}}^{s}$, $\log n_{\mathrm{gal}}$, and $\log n_{\mathrm{halo}}$, and their derivatives with respect to the HOD, modified gravity, and mass–observable parameters. We compute derivatives with respect to the cosmological parameters using paired simulation boxes that vary each parameter above and below the fiducial value: c102 and c103 for $\Omega_{\mathrm{m}}$, and c112 and c113 for $\sigma_{8}$. Table 1 lists the parameter configurations for each box, including $\Omega_{\mathrm{m}}$, the primordial scalar amplitude $A_{\mathrm{s}}$, $\sigma_{8}$, and the dimensionless Hubble parameter $h$.² In the AbacusSummit cosmology grid, $h$ is not an independent parameter. Instead, $H_{0}$ is chosen in each box so that the angular size of the sound horizon at last scattering, $\theta_{\ast}$, matches the value measured by Planck. Therefore, the values of $h$ listed for c102 and c103 shift along with $\Omega_{\mathrm{m}}$, rather than varying independently.

Box	Description	$\Omega_m$	$A_s$ [$10^{-9}$]	$\sigma_8$	$h$
c000	Planck	0.3152	2.0830	0.8080	0.6736
c102	high $\Omega_m$	0.3378	2.0205	0.8083	0.6597
c103	low $\Omega_m$	0.2942	2.1541	0.8081	0.6877
c112	high $\sigma_8$	0.3152	2.1672	0.8241	0.6736
c113	low $\sigma_8$	0.3152	2.0021	0.7921	0.6736

In this paper, we focus on the simulation snapshots at $z = 0.3$. To forecast constraints from DESI BGS (§ 6.3), we scale the covariance matrix to a comoving volume $V_{\mathrm{survey}} = 0.76\ h^{-3}\ \mathrm{Gpc}^{3}$, corresponding to the redshift interval $z = 0.25$–$0.35$ over a sky area of 14,000 $\mathrm{deg}^{2}$. We choose this interval because DESI BGS maintains a high galaxy number density $\bar{n} > 3 \times10^{-3}\ h^{3}\ \mathrm{Mpc}^{-3}$ out to $z = 0.35$ [26], roughly corresponding to a volume-limited sample with $M_{r} \lesssim-20.5$ in the DESI Y3 BGS-BRIGHT sample [57]. The high number density reduces shot noise in $\xi_{\mathrm{cg}}^{s}(r_{p}, r_{\pi})$ measurements, and it would allow a split into red and blue galaxy subsamples to test whether modified gravity signals and other parameter constraints are insensitive to the modeling of galaxy bias and the galaxy–halo connection.

We assume that the evolution of structure and HOD parameters can be neglected over a $\Delta z = 0.1$ interval, though this assumption will require verification in future work. DESI BGS would also enable measurements in the redshift intervals $0.15 < z < 0.25$ and $0.35 < z < 0.45$, with comoving volumes approximately $0.5\times$ and $1.6\times$ that of the fiducial $0.25 < z < 0.35$ interval, respectively. In the higher-redshift bin, it may be preferable to use the faint BGS sample, which offers a higher number density at the cost of more complex selection and completeness. At still higher redshifts, applying our method to luminous red galaxy (LRG) samples may be feasible, though their substantially lower number density and higher galaxy bias will require dedicated forecasts to assess sensitivity.

Mock Cluster Catalogs

Cluster samples are selected via observable mass proxies such as X-ray luminosity, the Sunyaev–Zel’dovich signal [58], or optical richness, rather than the underlying true halo mass. Each of these proxies only approximately traces the true halo mass [59, 60]. Halos of the same true mass can differ in formation history, structure, and environment, so even a noiselessly measured proxy retains an intrinsic scatter with respect to the true mass [59, 60]. Observational measurement uncertainty adds further scatter to the mass–observable relation. Intrinsic and observational scatter shift halos across any fixed selection threshold on the mass proxy, scattering lower-mass systems into the sample and higher-mass systems out of it. Neglecting this effect would introduce a significant bias into our forecast, as the resulting cluster sample would not reflect the realistic mass distribution of an observed cluster catalog.

To construct a physically realistic cluster catalog, we begin with the halo catalogs identified by the CompaSO halo finder [61]. We assign each halo an observed mass proxy $M_{\mathrm{obs}} by drawing from a log-normal distribution centered on the true halo mass $M_{\mathrm{true}}$,

where $\sigma_{\log M_c}$ denotes the combined logarithmic scatter (intrinsic plus measurement) between the proxy and the true cluster mass, in dex. We then define the cluster sample by selecting halos with $\log M_{\mathrm{obs}} \ge\log M_{\mathrm{obs}}^{\mathrm{cut}}$. This selection allows halos with true masses below the threshold to scatter into the sample and halos above the threshold to scatter out.

We treat both $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ and $\sigma_{\log M_c}$ as free parameters in our forecast, so the resulting halo number density $\log n_{\mathrm{halo}}$ varies across the parameter space. We include $\log n_{\mathrm{halo}}$ as an observable in the data vector. Figure 1 illustrates how variations in $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ and $\sigma_{\log M_c}$ alter the mean mass and comoving number density of the cluster sample.

In each panel, the red curve shows the limiting case of zero scatter ($\sigma_{\log M_c}=0$), which gives the highest mean mass $\langle\log M_{\mathrm{true}}\rangle$ and the lowest cumulative number density $\log n_{\mathrm{halo}}$ of clusters above the observable threshold $\log M_{\mathrm{obs}}^{\mathrm{cut}}$. Increasing $\sigma_{\log M_c}$ allows lower-mass halos to enter the sample, increasing $\log n_{\mathrm{halo}}$ while reducing $\langle\log M_{\mathrm{true}}\rangle$. For our fiducial model parameters, $\log M_{\mathrm{obs}}^{\mathrm{cut}}=14.0$ and $\sigma_{\log M_c}=0.25$ dex, the number density is $n_{\mathrm{halo}}=2.20\times10^{-5}\,h^3\,\mathrm{Mpc}^{-3}$, marked by the horizontal dashed line in panel (b). Points along the curves in panel (a) show $\langle\log M_{\mathrm{true}}\rangle$ at this number density for each value of $\sigma_{\log M_c}$. At the fiducial number density, increasing the scatter from $0$ to $0.3$ dex lowers the mean true mass by $\approx0.16$ dex.

Mock Galaxy Catalogs

We construct galaxy catalogs by populating the CompaSO halos using a Halo Occupation Distribution (HOD) framework³ [e.g., 62–68], which specifies the statistical distribution of galaxies within halos as a function of halo mass. The HOD formalism allows us to marginalize over the uncertain physics governing the galaxy–halo connection when deriving constraints on $F_{\mathrm{G}}$ or on cosmological parameters ($\sigma_8$, $\Omega_{\mathrm{m}}$). We include the galaxy number density $\log n_{\mathrm{gal}}$ as an observable in the data vector alongside $\log\xi_{\mathrm{cg}}$ and $\log n_{\mathrm{halo}}$, providing additional leverage to constrain the HOD parameters. In principle, the galaxy autocorrelation function could provide further HOD constraints [e.g. 69, 70], but we do not include such constraints here.

The mean occupation numbers for central and satellite galaxies are given by

following ref. [71]. The parameter $\log M_{\mathrm{min}}$ sets the characteristic halo mass at which the central occupation reaches 0.5, and $\sigma_{\log M}$ controls the sharpness of this transition. The satellite occupation follows a power law with slope $\alpha$ above a low-mass cutoff $M_0$, normalized such that $\left\langle N_{\mathrm{sat}}(M_0+M_1)\right\rangle=1$, where $M_1$ is the characteristic satellite mass scale parametrized by $\log M_1\equiv\log_{10}(M_1/h^{-1}M_\odot)$. We fix $\log_{10}(M_0/h^{-1}M_\odot)=10.0$, well below the minimum resolved halo mass, so the satellite occupation is effectively a pure power law, $\left\langle N_{\mathrm{sat}}\right\rangle\approx\left\langle N_{\mathrm{cen}}\right\rangle(M/M_1)^\alpha$. We draw the number of centrals per halo from a Bernoulli distribution with mean $\left\langle N_{\mathrm{cen}}(M)\right\rangle$, and the number of satellites from a Poisson distribution with mean $\left\langle N_{\mathrm{sat}}(M)\right\rangle$.

We place central galaxies at the center of their host halo, and distribute satellite galaxies spatially according to a Navarro–Frenk–White profile [NFW; 72],

where $c_{\mathrm{vir}}=r_{\mathrm{h}}/r_{\mathrm{s}}$ is the halo concentration. We assign halo concentrations using the mass–concentration relation of ref. [73], rescaled from $M_{200\mathrm{c}}$ to $M_{200\mathrm{b}}$ by multiplying by $\sqrt{2}$ [74].

Galaxies in spectroscopic surveys are observed in redshift space, where peculiar velocities displace their apparent line-of-sight positions. In our calculations, we treat each of the three Cartesian axes of the simulation cube in turn as the line-of-sight direction and average the results (§ 4). Since baryons and dark matter experience the same gravitational acceleration, the mean velocity of galaxies in a halo should on average match that of the halo itself. However, the velocity dispersion of satellite galaxies within a halo can differ from that of the dark matter, and the central galaxy may have a non-zero velocity relative to the halo center of mass. Following standard practice for redshift-space galaxy clustering, we characterize these differences through velocity bias parameters [75]. For a halo with three-dimensional dark matter velocity dispersion $\sigma_{\mathrm{v}}$, we assign each Cartesian component of the central galaxy velocity a stochastic offset drawn from a Laplace distribution with scale parameter $\eta_{\mathrm{cen}}\sigma_{\mathrm{v}}/\sqrt{6}$, where a factor of $\sqrt{3}$ converts $\sigma_{\mathrm{v}}$ to the one-dimensional dispersion and a factor of $\sqrt{2}$ converts the one-dimensional dispersion to the Laplace scale parameter. Similarly, we draw each Cartesian component of the satellite velocity relative to the halo mean velocity from a Gaussian with dispersion $\eta_{\mathrm{sat}}\,\sigma_v/\sqrt{3}$, so that $\eta_{\mathrm{sat}}=1$ corresponds to satellites tracing the velocity dispersion of the dark matter. Velocity bias may depend on the halo mass, and the satellite kinematics of infalling group-mass halos may differ systematically from those of massive cluster-mass halos. We therefore introduce separate parameters $\eta_{\mathrm{sat,low}}$ and $\eta_{\mathrm{sat,high}}$ that describe the satellite velocity bias in halos below and above $10^{13.5}\,h^{-1}\,M_\odot$, respectively.

We list the fiducial values for all varied parameters in Table 2, with cosmological parameter configurations given separately in Table 1. At $z=0.3$, our fiducial HOD choices give a comoving galaxy space density $n_{\mathrm{gal}}=2.76\times10^{-3}\,h^3\,\mathrm{Mpc}^{-3}$ ($\log n_{\mathrm{gal}}=-2.559$), averaged over 25 simulation phases, comparable to the density of BGS-BRIGHT galaxies at this redshift [26]. The large-scale galaxy bias factor computed from $(\xi_{\mathrm{gg}}/\xi_{\mathrm{mm}})^{1/2}$ on scales $20$–$30\,h^{-1}\,\mathrm{Mpc}$, with $\xi_{\mathrm{mm}}$ the linear-theory matter correlation function of our fiducial cosmology at $z=0.3$, is $b_g\approx1.4$.

Parameter	Definition	Fiducial	$\delta\theta$
$\sigma_{\log M}$	Log-normal scatter width of central occupation transition	0.60	0.10
$\log M_{\min}$	$\log_{10}$ of characteristic minimum halo mass [$h^{-1}\,M_\odot$]	12.50	0.10
$\log M_1$	$\log_{10}$ of satellite mass scale $M_1$ [$h^{-1}\,M_\odot$]	13.35	0.10
$\alpha$	Power-law index of satellite occupation	1.00	0.10
$\eta_{\mathrm{cen}}$	Central-galaxy velocity bias	0.10	0.10
$\eta_{\mathrm{sat,low}}$	Satellite velocity bias for $M \leq10^{13.5}\,h^{-1}\,M_\odot$	1.00	0.10
$\eta_{\mathrm{sat,high}}$	Satellite velocity bias for $M \geq10^{13.5}\,h^{-1}\,M_\odot$	1.00	0.10
$\sigma_{\log M_c}$	Scatter in mass–observable relation [dex]	0.25	0.05
$\log M_{\mathrm{obs}}^{\mathrm{cut}}$	$\log_{10}$ of observable mass threshold [$h^{-1}\,M_\odot$]	14.00	0.25
$F_{\mathrm{G}}$	Phenomenological modified gravity parameter	1.00	0.33

Infall Velocities

The mean radial peculiar velocity profile $v_{\mathrm{rad}}(r)$ of galaxies around massive clusters is the central physical quantity of this analysis. It encodes the coherent kinematic response of matter to the cluster gravitational potential and, as discussed in § 1, is directly sensitive to gravitational accelerations in the quasi-linear infall regime. Before connecting $v_{\mathrm{rad}}(r)$ to the observable redshift-space cross-correlation function $\xi_{\mathrm{cg}}^{s}$ in § 4, we first characterize how it responds to variations in each model parameter.

Physical Motivation and Parameterization

Galaxy clusters are the most massive gravitationally bound structures in the Universe, and their deep potential wells drive coherent infall of surrounding matter out to large distances. In many modified gravity theories, a scalar field contributes to the gravitational acceleration alongside the metric, enhancing infall velocities relative to the GR prediction at fixed cluster mass. Acceleration is enhanced in unscreened regions, while GR is restored in screened regions; for some broad classes of models, the enhancement in unscreened regions is a factor of $4/3$ [13]. Since galaxy velocities accumulate over time, they may remain enhanced relative to GR even in regions that have since been screened [8].

Rather than simulating specific modified gravity models, we adopt a phenomenological approach in which departures from GR are parameterized by a single multiplicative factor $F_{\mathrm{G}}$ that uniformly rescales the mean velocities of halos relative to the GR prediction encoded in the AbacusSummit simulations. This factor consequently rescales mean infall velocities,

where $v_{\mathrm{rad}}^{\mathrm{GR}}(r)$ denotes the baseline velocity profile measured in our GR simulation suite. In modified gravity models we generally expect $1 \leq F_{\mathrm{G}} \leq4/3$, though our parameterization permits $F_{\mathrm{G}} < 1$. Zu et al. [8] find that $f(R)$ and Galileon gravity produce qualitatively similar enhancements in the cluster infall velocity profile relative to GR, despite their distinct screening mechanisms. At $r = 5\,h^{-1}\,\mathrm{Mpc}$, the infall velocity enhancement in the $f(R)$ model is approximately 30% for $f_0 = 10^{-4}$ and 15% for $f_0 = 10^{-5}$ (see their figure 5). Comparable but slightly smaller enhancements are seen in their Galileon models. A uniform velocity rescaling captures the dominant signature of modified gravity in a model-independent manner. In future work, we plan to extend this parameterization to allow for scale- or mass-dependent $F_{\mathrm{G}}$ calibrated directly against modified gravity simulations.

Dependence of the Infall Velocity Profile on Model Parameters

Figure 2 presents $v_{\mathrm{rad}}(r)$ computed from the halo and galaxy catalogs described in § 2, evaluated at the fiducial Planck cosmology using a single simulation phase (ph000). For the cosmological variants, we use a single matched phase from the corresponding boxes listed in Table 1. The profile exhibits the characteristic structure of the cluster infall regime. At large separations ($r \gtrsim10\,h^{-1}\,\mathrm{Mpc}$), the mean radial velocity approaches zero, consistent with the background Hubble flow. Moving inward, galaxies accelerate toward the cluster, and the mean infall speed peaks at $r \sim1$–$3\,r_{\mathrm{vir}}$. At smaller radii, coherent infall gives way to virialized, approximately isotropic motions within the cluster interior, and the mean radial velocity returns toward zero. Zu & Weinberg [7] decompose cluster velocity fields into infalling and orbiting components, as have several subsequent studies [76, 77, 78, 79, 80, 81, 82, 83].

Panel (a) shows how $v_{\mathrm{rad}}(r)$ responds to variations in the mass–observable relation parameters $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ and $\sigma_{\log M_c}$. Increasing $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ selects halos of higher virial mass, deeper gravitational potential, and larger virial radius, raising the peak infall speed and shifting the peak to larger radius. The peak infall speed increases from roughly 700 km s$^{-1}$ at the fiducial $\log M_{\mathrm{obs}}^{\mathrm{cut}} = 14.0$ (navy curve) to roughly 1250 km s$^{-1}$ at $\log M_{\mathrm{obs}}^{\mathrm{cut}} = 15.0$ (orange curve), reflecting the faster infall into more massive and spatially extended potential wells, consistent with the expected mass–velocity scaling. Increasing $\sigma_{\log M_c}$ at fixed $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ (dotted navy curve, $\sigma_{\log M_c} = 0.35$ dex) introduces a larger fraction of lower-mass halos into the sample through upward scatter across the selection threshold, reducing the effective mean cluster mass and lowering the peak infall speed by roughly 50 km s$^{-1}$ relative to the fiducial.

Panel (b) shows how $v_{\mathrm{rad}}(r)$ responds to variations in $F_{\mathrm{G}}$. Since $F_{\mathrm{G}}$ rescales the mean halo velocities without altering the internal velocity dispersions, its effect on $v_{\mathrm{rad}}(r)$ is nearly identical to a direct multiplication by $F_{\mathrm{G}}$. An increase in $F_{\mathrm{G}}$ from 1.0 to 1.33 raises the peak infall speed from approximately 700 to 950 km s$^{-1}$, comparable to the effect of increasing $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ from 14.0 to 14.5. However, the two effects are distinguishable: raising $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ also shifts the peak to larger radius by increasing the halo virial radius, whereas increasing $F_{\mathrm{G}}$ uniformly rescales the velocity field at fixed spatial structure, raising the peak infall speed without appreciably shifting its location.

Panel (c) compares $v_{\mathrm{rad}}(r)$ between the fiducial cosmology (navy curve) and the high-$\Omega_{\mathrm{m}}$ (red curve) and high-$\sigma_8$ (blue curve) variants listed in Table 1. All three profiles are nearly indistinguishable, unlike the $F_{\mathrm{G}} = 1.33$ curve (orange), whose peak infall speed is about 250 km s$^{-1}$ higher.

Figure 3 shows the ratio $v_{\mathrm{rad}}/v_{\mathrm{rad}}^{\mathrm{fid}}$ for variations in $\Omega_{\mathrm{m}}$ (top row) and $\sigma_{8}$ (bottom row). These ratios explain why $v_{\mathrm{rad}}(r)$ is so insensitive to cosmology. For each parameter, we compare fixed-number-density samples (the 100,000 most massive halos; left column) and samples selected at a fixed mass threshold, $\log M_{\mathrm{obs}}^{\mathrm{cut}} = 14.0$ (right column). This decomposition isolates changes in the velocity field due to the cosmological parameter itself from those arising indirectly through shifts in the selected cluster population. Figure 2 corresponds to the fixed $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ case.

When cluster abundance is held fixed (panel (a)), a higher $\Omega_{\mathrm{m}}$ requires a larger mass threshold, since a denser universe produces more massive structures. The infall velocities around these more massive halos are enhanced by 2.9% at $r = 2$ $h^{-1}\,\mathrm{Mpc}$, decreasing to 1.6% by $r = 10$ $h^{-1}\,\mathrm{Mpc}$ (green curve, $\Omega_m = 0.338$). At fixed $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ (panel (b)), the differences are small at $r \gtrsim2$ $h^{-1}\,\mathrm{Mpc}$ because, relative to the fiducial model, the high-$\Omega_m$ model raises

the mean mass of the selected halos by only 1.0% (0.0044 dex) and the linear growth rate $f(\Omega_{\rm m}) \approx\Omega_{\rm m}^{0.55}$ at $z = 0.3$ by only 3%. Below $r \approx2\ h^{-1}\ \mathrm{Mpc}$ the curves diverge more strongly, with the high-$\Omega_{\rm m}$ model up to $\sim8\%$ above and the low-$\Omega_{\rm m}$ model up to $\sim9\%$ below the fiducial near $r = 1\ h^{-1}\ \mathrm{Mpc}$. Lower $\Omega_{\rm m}$ also shifts matter-radiation equality to larger scales, boosting large-scale power and producing slightly larger $v_{\rm rad}$ at $r \gtrsim15\ h^{-1}\ \mathrm{Mpc}$, reaching a 2.1% excess by $r = 30\ h^{-1}\ \mathrm{Mpc}$ (the high-$\Omega_{\rm m}$ model is correspondingly 2.3% lower).

Raising $\sigma_8$ increases the amplitude of the matter power spectrum without changing its shape. At fixed number density (panel (c)), a higher $\sigma_8$ (green curve) selects more massive halos, producing enhanced infall velocities of 1.9–2.4% at $r = 10$–$30\ h^{-1}\ \mathrm{Mpc}$, consistent with the 2.0% fractional difference in $\sigma_8$ itself. The ordering inverts at $r \lesssim1.6\ h^{-1}\ \mathrm{Mpc}$, where the high-$\sigma_8$ model falls below unity, by $\sim9\%$ at $r = 1\ h^{-1}\ \mathrm{Mpc}$. This inversion arises because more massive halos have larger virial radii, shifting the transition from coherent infall to virialized orbital motion to larger physical scales. Fixing $\log M_{\rm obs}^{\rm cut}$ suppresses the difference between models to $\lesssim1.4\%$ at $r \gtrsim5\ h^{-1}\ \mathrm{Mpc}$ (panel (d)), and the small-scale inversion persists at reduced amplitude. The remaining difference may be caused by the slight increase in the mean halo mass above the selection threshold, which is about 0.0078 dex (1.8%) higher in the high-$\sigma_8$ model.

Returning to Figure 2, panel (d) shows the response of $v_{\rm rad}(r)$ to variations in the HOD parameters $\sigma_{\log M}$ (green), $\log M_{\min}$ (yellow), $\alpha$ (orange), and $\eta_{\rm sat}$ (red), each evaluated at the perturbed values listed in Table 2. None of these variations produces a discernible deviation from the fiducial profile, indicating that $v_{\rm rad}(r)$ is insensitive to the details of the galaxy–halo connection captured by our HOD parameterization. This is broadly consistent with the expectation that gravitational acceleration at radius $r$ is primarily sourced by the cluster mass profile rather than by the internal structure of the galaxy population. Our HOD model does not include assembly bias, in which galaxy properties depend on halo formation history or large-scale environment beyond halo mass alone. Since such effects could preferentially select galaxies tracing different regions of the infall zone, we cannot entirely rule out that a more realistic treatment of assembly bias would introduce environmental sensitivity in $v_{\rm rad}(r)$ that our current model does not capture. However, because $v_{\rm rad}(r)$ is insensitive to the standard HOD parameters, which also change how galaxies trace the matter distribution, a significant impact of assembly bias on $v_{\rm rad}(r)$ seems unlikely.

Taken together, the four panels show that $v_{\rm rad}(r)$ responds most strongly to $F_{\rm G}$, $\log M_{\rm obs}^{\rm cut}$, and $\sigma_{\log M_c}$, while remaining largely insensitive to the HOD and cosmological parameters explored here. Unlike $\log M_{\rm obs}^{\rm cut}$ and $\sigma_{\log M_c}$, which also shift the peak location, $F_{\rm G}$ produces a spatially uniform velocity enhancement. HOD and cosmological parameters influence our observable primarily through the real-space correlation function and the number densities of galaxies and halos, rather than through $v_{\rm rad}$ itself.

The Redshift-Space Galaxy–Cluster Cross-Correlation Function

Having characterized how $v_{\rm rad}(r)$ responds to variations in the model parameters, we now turn to the observable quantity that encodes this velocity field in spectroscopic data. The real-space infall profile is not directly accessible; instead, peculiar velocities displace the apparent line-of-sight positions of galaxies, imprinting a characteristic anisotropy on the redshift-space galaxy–cluster cross-correlation function $\xi_{\rm cg}^{s}(r_p,r_\pi)$.

We estimate $\xi_{\mathrm{cg}}^{\mathrm{s}}(r_p,r_\pi)$ from galaxy–cluster pair counts as

where $DD(r_p,r_\pi)$ is the number of galaxy–cluster pairs in a given $(r_p,r_\pi)$ bin and $RR(r_p,r_\pi)$ is the expected number of pairs if both the galaxies and clusters were independently drawn from uniform random distributions with the same mean number densities. Assuming uniform densities, we compute the random counts analytically:

where $N_{\mathrm{c}}$ is the number of clusters, $n_{\mathrm{g}}$ is the galaxy number density, and the factor of 2 enters because values of $\pm r_p$ are assigned to the same bin. Since we compute the $RR$ term analytically and the simulations employ periodic boundary conditions, we do not require the Landy–Szalay estimator [84] to reduce variance.

We compute galaxy–cluster pair counts, $DD(r_p,r_\pi)$, using the Corrfunc package [85, 86] in cylindrical coordinates. The transverse separation spans $r_p\in[0.3,30]\,h^{-1}\mathrm{Mpc}$ with a nonlinear binning scheme,

where $\alpha_{\mathrm{bins}} = 2$ concentrates bins toward small scales where the signal varies most rapidly. We divide the line-of-sight separation into 30 uniform bins of width $1\,h^{-1}\mathrm{Mpc}$, covering $r_\pi\in[0,30]\,h^{-1}\mathrm{Mpc}$.

We incorporate redshift-space distortions by displacing galaxy and halo positions along the line of sight,

where $r_\parallel$ is the real-space line-of-sight separation, $v_\parallel$ is the peculiar velocity component along the line of sight, and $H(z)$ is the Hubble parameter. To increase the precision of the measurement, we treat each of the three Cartesian axes of the simulation box in turn as the line of sight and average the results, enforcing periodic boundary conditions by wrapping displaced positions back into the box. For $F_{\mathrm{G}} \ne1$, we multiply the halo peculiar velocities by $F_{\mathrm{G}}$, and galaxies inherit these enhanced bulk motions. This rescales the mean infall profile (equation (5)) without changing the internal velocity dispersions within halos.

Figure 4 shows $\log\xi_{\mathrm{cg}}^{\mathrm{s}}(r_p,r_\pi)$ at $z = 0.3$, averaged over 25 independent realizations of the fiducial AbacusSummit cosmology and over the three line-of-sight directions. In the absence of peculiar velocities, contours of constant $\xi_{\mathrm{cg}}^{\mathrm{s}}$ would be circular in the $(r_p,r_\pi)$ plane. Deviations from circular symmetry therefore directly reflect the influence of the velocity field. On small transverse scales, large virial velocities of galaxies within cluster halos elongate the contours along the line of sight, producing the FoG. On larger scales, coherent infall compresses the contours along $r_\pi$ [16, 17], a distortion commonly referred to as the Kaiser effect.

Sensitivity of Observables to Model Parameters

To quantify how each model parameter affects our observables, we compute derivatives of the full data vector ($\log\xi_{\mathrm{cg}}^{\mathrm{s}}$, $\log n_{\mathrm{gal}}$, and $\log n_{\mathrm{halo}}$) with respect to each parameter $\theta$, evaluated at the fiducial values listed in Table 2. We calculate these derivatives using the centered finite-difference approximation,

where $\mathcal{O}$ denotes any element of the data vector, and Table 2 lists the step sizes $\delta\theta$. For $\log\xi_{\mathrm{cg}}^{\mathrm{s}}$ and the number densities, we derive the HOD, modified gravity, and mass–observable parameter derivatives by averaging over the 25 Planck fiducial cosmology phases of AbacusSummit, and obtain the $\Omega_{\mathrm{m}}$ and $\sigma_8$ derivatives from the paired simulation boxes listed in Table 1. We show the real-space counterparts of the $\log\xi_{\mathrm{cg}}^{\mathrm{s}}$ derivatives in appendix § 9 for reference.

The Redshift-Space Cross-Correlation Function

Figure 5 and Figure 6 show $\partial\log\xi_{\mathrm{cg}}^{\mathrm{s}}/\partial\theta$ throughout the $(r_p,r_\pi)$ plane. Figure 5 presents derivatives with respect to the HOD and mass–observable parameters, while Figure 6 presents those with respect to the cosmological parameters, with $F_{\mathrm{G}}$ included in both for reference. In each figure,

red and blue hues indicate positive and negative derivatives, corresponding to enhancement and suppression of the clustering signal, respectively. The black dots in the $\sigma_{\log M}$ panel of Figure 5 and in all panels of Figure 6 mark the discrete $(r_p,r_\pi)$ bins at which we evaluate the derivatives. These values are the quantitative output used in the Fisher analysis of § 6. The contours interpolated between these points serve only to visualize the spatial structure of each derivative for the reader and play no role in the quantitative analysis itself. We present derivatives of the real-space correlation function $\xi^{\mathrm{r}}_{\mathrm{cg}}(r_p,r_\pi)$ in appendix § 9 and reference them throughout the discussion below.

Galaxy Occupation Parameters

The parameters $\log M_{\min}$ and $\sigma_{\log M}$ govern which halos host a central galaxy by setting the characteristic mass scale and the sharpness of the occupation transition (equation (2)). We interpret their effects on $\xi^{\mathrm{s}}_{\mathrm{cg}}$ most directly through the way they shift the effective bias of the galaxy sample. On large scales, where linear bias applies, the galaxy–cluster cross-correlation approaches

where $\xi_m(r)$ is the matter auto-correlation function, $b_c$ is the effective large-scale bias of the cluster sample, and $b_g$ is the effective large-scale bias of the galaxy sample, given by [67]

where $dn/dM$ is the halo mass function, $\langle N(M)\rangle= \langle N_{\mathrm{cen}}(M)\rangle+ \langle N_{\mathrm{sat}}(M)\rangle$ is the mean total occupation, $b_h(M)$ is the halo bias, and $\bar{n}_g$ normalizes the weighted average. Any HOD parameter that shifts the occupation-weighted mean of $b_h(M)$ directly alters $b_g$ and hence the large-scale amplitude of $\xi^{\mathrm{s}}_{\mathrm{cg}}$.

Raising $\log M_{\min}$ restricts central occupation to more massive, more strongly biased halos, increasing $b_g$ through equation (12) and amplifying the two-halo clustering term through equation (11). Since these halos lie far above the $M_{\min}$ threshold, raising $\log M_{\min}$ does not affect the number of galaxies within clusters. However, it lowers the overall galaxy number density $n_g$, which increases the one-halo term of $\xi^{\mathrm{s}}_{\mathrm{cg}}$ by reducing the $RR$ denominator in equation (6). This effect amplifies the derivative at small $r$ in real space (appendix § 9) and produces a strong positive response in the FoG regime (small $r_p$, wide range of $r_\pi$) in redshift space (Figure 5).

Increasing $\sigma_{\log M}$ allows a larger fraction of lower mass, less biased halos to host central galaxies, producing effects nearly opposite to those of $\log M_{\min}$ in both real and redshift space. The primary difference is that the relative response in the one-halo FoG regime is weaker for $\sigma_{\log M}$ than for $\log M_{\min}$.

Satellite Galaxy Occupation Parameters

The satellite occupation parameters $\log M_1$ and $\alpha$ govern the mean number of satellite galaxies per halo (equation (3)). Decreasing $\log M_1$ increases the satellite count across halos of all masses, while increasing $\alpha$ preferentially boosts satellites in high-mass halos relative to lower-mass halos. Both changes enhance the large-scale bias by increasing the contribution of high-mass, strongly biased halos in equation (12), and they amplify the one-halo term by increasing the satellite population within clusters. The latter effect is stronger for $\alpha$, since the clusters preferentially reside in high-mass halos where this parameter has the greatest leverage.

In both real space (appendix § 9) and redshift space (Figure 5), the derivatives for $\log M_1$ and $\alpha$ resemble those for $\log M_{\min}$, with a sign reversal for $\log M_1$. While $\log M_{\min}$ affects the one-halo term through the $RR$ denominator of equation (6), $\log M_1$ and $\alpha$ act through the $DD$ numerator. The satellite parameters, particularly $\alpha$, have a stronger influence on the one-halo term relative to the large-scale bias, producing a more pronounced FoG response. In particular, increasing $\alpha$ produces a large enhancement of $\xi_{\mathrm{cg}}^{\mathrm{s}}(r_p,r_\pi)$ at small $r_p$ and large $r_\pi$, since it increases the satellite population in the highest-mass clusters, which also have the largest velocity dispersions.

Cluster Sample Selection Parameters

The parameters $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ and $\sigma_{\log M_c}$ govern the properties of the cluster sample. Either increasing $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ or decreasing $\sigma_{\log M_c}$ raises the mean cluster mass, so the two parameters have nearly opposite effects on $\xi_{\mathrm{cg}}^{\mathrm{s}}(r_p,r_\pi)$. More massive clusters have a larger bias factor $b_c$, and in real space the resulting amplification or suppression of the cross-correlation is nearly scale-independent beyond $r \approx2 h^{-1}\,\mathrm{Mpc}$. In the one-halo regime, the effect is stronger because more massive clusters host larger satellite populations, which in redshift space produce the FoG distortion. Consequently, the largest derivatives for the cluster parameters occur at small $r_p$ over a wide range of $r_\pi$. Although more massive clusters also have larger infall velocities, the redshift-space derivatives at large $r_p$ show only marginally more structure than their real-space counterparts. The cluster selection parameters therefore affect $\xi_{\mathrm{cg}}^{\mathrm{s}}(r_p,r_\pi)$ primarily by modulating the FoG strength through the satellite population within clusters.

Velocity Parameters

The velocity parameters modulate the line-of-sight motions of galaxies and halos without altering their positions, so derivatives of the real-space log $\xi_{\mathrm{cg}}^{\mathrm{r}}(r_p,r_\pi)$ vanish with respect to these parameters, and their signatures are entirely encoded in the anisotropy of the redshift-space log $\xi_{\mathrm{cg}}^{\mathrm{s}}(r_p,r_\pi)$. Figure 5 presents derivatives with respect to $\eta_{\mathrm{sat,high}}$ and $F_G$; we omit those with respect to $\eta_{\mathrm{sat,low}}$ and $\eta_{\mathrm{cen}}$, whose amplitudes are negligible compared to the other parameters. In the one-halo regime, $\eta_{\mathrm{sat,low}}$ has no effect because it does not modulate cluster-mass halos, and $\eta_{\mathrm{cen}}$ has little effect because it governs only the central galaxy. In the two-halo regime, these parameters affect $\xi_{\mathrm{cg}}^{\mathrm{s}}(r_p,r_\pi)$ in principle by modifying the line-of-sight velocity distribution (see equation 1 of ref. [7]), but since they do not alter the mean infall velocity, the impact of a slightly broader distribution proves negligible in practice.

The parameter $\eta_{\mathrm{sat,high}}$ scales the velocity dispersions of satellite galaxies within cluster halos above $M = 10^{13.5}h^{-1}\,M_\odot$ relative to the virial velocity dispersion. Increasing $\eta_{\mathrm{sat,high}}$ broadens the line-of-sight velocity distribution of satellites orbiting within these massive halos, redistributing pairs along $r_\pi$ within the one-halo regime. This produces a positive derivative at small $r_p$ and large $r_\pi$, where pairs scatter to larger line-of-sight separations, and a negative derivative at small $r_p$ and small $r_\pi$, where those same pairs are depleted, resulting in the dipole structure visible at small $r_p$ in Figure 5. The derivative remains non-zero at large $r_p$ and large $r_\pi$ because a significant fraction of excess pairs in this regime involves the high-velocity tail of the cluster population, but it vanishes at large $r_p$ and smaller $r_\pi$. We emphasize that 10–20% departures from $\eta_{\mathrm{sat,high}} = 1$ are entirely plausible within GR as a consequence of orbital anisotropy, dynamical friction, preferential disruption of satellites on low eccentricity orbits, or departures from dynamical equilibrium.

The modified gravity parameter $F_{\mathrm{G}}$ enhances the coherent infall velocity field at fixed cluster mass and internal velocity dispersion, mimicking a gravitational fifth force that accelerates matter toward cluster overdensities without disturbing the cluster’s virial equilibrium. Its signature is therefore confined to large-scale galaxy streaming, with no corresponding change in the FoG signal. A residual response does appear at small $r_{\mathrm{p}}$, but this arises from galaxy pairs in the two-halo regime rather than from pairs within the cluster, as is the case for $\eta_{\mathrm{sat,high}}$. Increasing $F_{\mathrm{G}}$ shifts a fraction of these small-$r_{\mathrm{p}}$ pairs from small to larger $r_{\pi}$. At large $r_{\mathrm{p}}$, compression of structure along the line of sight produces a positive derivative at small to intermediate $r_{\pi}$ across all transverse scales. In contrast, at large $r_{\pi}$, where pairs are drawn inward from large line-of-sight separations, the derivative becomes negative. Appendix § 10 presents an idealized spherical model for cluster redshift-space distortions, which explains the qualitative features of the $F_{\mathrm{G}}$ derivative in Figure 5 but does not reproduce it quantitatively.

Cosmological Parameters

Figure 6 shows the derivatives with respect to $\sigma_{8}$ and $\Omega_{\mathrm{m}}$, alongside the $F_{\mathrm{G}}$ derivative included for reference. We present the corresponding real-space derivatives in appendix § 9. As noted above, we compute the cosmological derivatives from single paired simulation phases rather than the full average over 25 realizations. They are consequently much noisier than the HOD and mass–observable derivatives in Figure 5. Since the volume of a single AbacusSummit box is approximately ten times larger than the DESI BGS volume in the range $z = 0.25$–$0.35$, we expect this noise to have a small effect relative to the forecast statistical uncertainties. Nonetheless, a more robust assessment will require analysis of additional simulations. We consider our $\Omega_{\mathrm{m}}$ and $\sigma_{8}$ forecasts preliminary and subject to revision in future work.

The derivative with respect to $\sigma_{8}$ exhibits an anisotropic structure concentrated primarily in the one-halo regime. At small $r_{\mathrm{p}}$, the topology resembles the dipole-like signature of $\eta_{\mathrm{sat,high}}$, with a negative derivative at small $r_{\pi}$ and a positive derivative at larger $r_{\pi}$. This structure reflects the response of the cluster population to increasing $\sigma_{8}$ at fixed $\log M_{\mathrm{cut}}$. While a higher $\sigma_{8}$ globally enhances matter clustering, it also shifts the halo mass function toward more massive clusters above the selection threshold. These more massive clusters have deeper potential wells and larger internal velocity dispersions, driving stronger FoG stretching along $r_{\pi}$ and redistributing pairs from small to large $r_{\pi}$, locally competing with and exceeding the global amplitude boost. On large scales, the real- and redshift-space derivatives lack coherent structure and likely reflect noise from the limited simulation volume.

The derivative with respect to $\Omega_{\mathrm{m}}$ is negative across both the one-halo and two-halo regimes. Because the cosmological parameters are defined at $z = 0$, increasing $\Omega_{\mathrm{m}}$ at fixed $\sigma_{8}(z = 0)$ results in a lower growth factor at the epoch of observation ($z = 0.3$). This suppresses the amplitude of matter fluctuations, reducing $\log\xi_{\mathrm{csg}}$ throughout all scales and producing the negative derivative visible in Figure 6. The derivatives are qualitatively similar in real and redshift space, suggesting that changes in 3D clustering, rather than redshift-space distortions, primarily drive them. This inference is also supported by Figure 2c, which shows that $v_{\mathrm{rad}}(r)$ is only weakly sensitive to changes in $\Omega_{\mathrm{m}}$ and $\sigma_{8}$. However, the shapes of derivative contours for $\Omega_{\mathrm{m}}$ are somewhat different in real and redshift space.

An important feature of Figure 6 is that the $\sigma_{8}$ and $\Omega_{\mathrm{m}}$ derivatives exhibit a different structure from that of $F_{\mathrm{G}}$. Furthermore, $\sigma_{8}$ and $\Omega_{\mathrm{m}}$ influence the cluster space density (§ 5.2), while $F_{\mathrm{G}}$ does not. Together, these distinctions suggest that little degeneracy exists between $F_{\mathrm{G}}$ and the cosmological parameters $\sigma_{8}$ and $\Omega_{\mathrm{m}}$.

Galaxy and Cluster Number Densities

In addition to $\log\xi_{\mathrm{cg}}^{\mathrm{s}}$, we include $\log n_{\mathrm{gal}}$ and $\log n_{\mathrm{halo}}$ in our data vector, calculated directly from the simulation as $n_{\mathrm{gal}} = N_{\mathrm{gal}}/V_{\mathrm{box}}$ and $n_{\mathrm{halo}} = N_{\mathrm{halo}}/V_{\mathrm{box}}$. Figure 7 shows their derivatives with respect to each model parameter. We average the derivatives with respect to the HOD and mass–observable parameters over the 25 Planck fiducial cosmology phases, and compute those with respect to $\Omega_{\mathrm{m}}$ and $\sigma_8$ using the paired simulation boxes in Table 1. We omit $F_{\mathrm{G}}$ and the velocity bias parameters ($\eta_{\mathrm{sat,high}}$, $\eta_{\mathrm{sat,low}}$, $\eta_{\mathrm{cen}}$) from Figure 7, since they modify the motions of galaxies and halos without affecting their total counts.

The halo number density ($\log n_{\mathrm{halo}}$; shown in red) responds only to the cluster selection and cosmological parameters, since HOD variations only change the distribution of galaxies within halos and have no effect on the halo population itself. Increasing $\log M_{\mathrm{cut}}^{\mathrm{obs}}$ restricts the cluster sample to rarer and more massive halos, producing a large negative derivative. Increasing $\sigma_{\log M_{\mathrm{c}}}$ allows more low-mass halos to scatter above the selection threshold than high-mass halos to scatter below it, resulting in a net increase in $n_{\mathrm{halo}}$. Both $\sigma_8$ and $\Omega_{\mathrm{m}}$ produce positive derivatives, as higher values of either parameter increase the overall abundance of halos above the selection threshold (see refs. [53, 87] for a more detailed discussion). For $\sigma_8$, the enhanced clustering amplitude shifts the halo mass function toward more massive halos, directly increasing the abundance above $\log M_{\mathrm{cut}}^{\mathrm{obs}}$. For $\Omega_{\mathrm{m}}$, the increased background matter density shifts the mass function to higher mass at all space densities, which outweighs the slight suppression of the growth factor at $z = 0.3$.

The galaxy number density ($\log n_{\mathrm{gal}}$; shown in blue) is sensitive to both the HOD and cosmological parameters, since both govern the abundance of halos available to host galaxies. It is insensitive to $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ and $\sigma_{\log M_c}$, as these parameters determine which halos enter the cluster sample but do not affect galaxy assignment. Among the HOD parameters, $\log M_{\min}$ and $\sigma_{\log M}$ produce the dominant responses, reflecting the fact that central galaxies constitute approximately 80% of the total galaxy population. Raising $\log M_{\min}$ restricts central occupation to more massive halos, leaving a large population of lower-mass halos without central galaxies and producing a large negative derivative. Increasing $\sigma_{\log M}$ broadens the occupation transition, allowing more lower-mass halos to host central galaxies and yielding a positive derivative of comparable magnitude. The satellite parameters $\log M_1$ and $\alpha$ have a smaller impact on $\log n_{\mathrm{gal}}$, as satellites constitute only $\sim20\%$ of the total galaxy population. Raising $\log M_1$ suppresses satellite occupation in intermediate-mass halos, reducing $\log n_{\mathrm{gal}}$ and producing a negative derivative. Increasing $\alpha$ steepens the satellite occupation–mass relation, redistributing satellites from abundant intermediate-mass halos to rare high-mass halos. The gain in the high-mass tail does not compensate for the loss at lower masses, yielding a small negative derivative. As with $\log n_{\mathrm{halo}}$, both $\sigma_8$ and $\Omega_m$ yield positive derivatives, driven by the increased abundance of halos massive enough to host central and satellite galaxies.

The derivatives in Figure 5, Figure 6, and Figure 7 demonstrate that each parameter imprints a distinct signature across the full data vector. In particular, while the cosmological and HOD parameters have minimal impact on $v_{\mathrm{rad}}(r)$ as shown in Figure 2, their derivatives with respect to $\log\xi_{\mathrm{cg}}^{s}$ exhibit structured anisotropic patterns across the $(r_p,r_\pi)$ plane, indicating that the redshift-space clustering signal encodes additional parameter sensitivity beyond the mean infall velocity profile.

Of particular relevance for constraining $F_{\mathrm{G}}$ is the degree to which other parameters produce similar signatures in $\log\xi_{\mathrm{cg}}^{s}$. Parameters such as $\log M_{\min}$, $\log M_{\mathrm{obs}}^{\mathrm{cut}}$, and $\alpha$ all produce broadly positive derivatives throughout the $(r_p,r_\pi)$ plane, superficially similar to $F_{\mathrm{G}}$. However, these parameters also enhance the FoG elongation at small $r_p$ through their impact on the halo mass and satellite population, whereas $F_{\mathrm{G}}$ leaves the internal velocity dispersions unchanged and produces no corresponding FoG enhancement. Furthermore, $F_{\mathrm{G}}$ produces a negative derivative at large $r_\pi$ that is absent in the HOD and mass–observable derivatives. These distinctions in the anisotropic structure of $\log\xi_{\mathrm{cg}}^{s}$ allow us to constrain $F_{\mathrm{G}}$ independently of the mass–observable, HOD, and cosmological parameters, as we quantify in the following section.

Fisher Forecast

Formalism and Construction

We define the full data vector as

where the first set of elements corresponds to the two-dimensional correlation function evaluated over all $(r_p,r_\pi)$ bins, followed by the scalar observables $\log n_{\mathrm{gal}}$ and $\log n_{\mathrm{halo}}$. As described in § 5, we calculate the derivatives of each observable with respect to the model parameters using the centered finite-difference approximation in equation (10).

Given the derivatives $\partial d/\partial\theta_i$ evaluated at the fiducial parameter values (see Table 1 and Table 2) and the covariance matrix $C$ of the full data vector (see § 6.2 below), we construct the Fisher information matrix following the standard formalism [88],

where $d_k$ runs over all elements of the combined data vector. The covariance matrix takes the block form

where $C_{\xi\xi}$ is the covariance matrix of $\log\xi_{\mathrm{cg}}^{s}$, $\sigma_{n_g}^{2}$ and $\sigma_{n_h}^{2}$ are the variances of $\log n_{\mathrm{gal}}$ and $\log n_{\mathrm{halo}}$, the off-diagonal blocks $C_{\xi n_g}$ and $C_{\xi n_h}$ describe the cross-covariance between the correlation function and each number density, and $C_{n_g n_h}$ captures the covariance between the two number density observables. We obtain marginalized $1\sigma$ constraints on each parameter $\theta_i$ from the inverse Fisher matrix,

Unless otherwise stated, all quoted constraints are marginalized over the full parameter set.

Covariance matrix

We estimate $C_{\xi\xi}$ from subvolumes of the 25 Planck fiducial cosmology phases of AbacusSummit. For each phase, we divide the simulation volume into 64 equal subvolumes, each with side length $500\,h^{-1}\,\mathrm{Mpc}$, and measure $\log\xi_{\mathrm{cg}}^{s}(r_p,r_\pi)$ along the three Cartesian axes of the simulation box. We confirm that this choice of 25 phases and 64 subvolumes per phase, for a total of $N_s = 1600$ subvolumes, is sufficient by testing the convergence of the derived parameter constraints as a function of sample size, finding that all marginalized parameter errors at $N_s \approx1000$ lie within $\sim3\%$ of their $N_s = 1600$ values, and within $1.3\%$ at $N_s = 1280$. We compute a covariance matrix for each projection separately using the scatter across all subvolumes and phases,

where $d_i$ is $\log\xi_{\mathrm{cg}}^{s}(r_p,r_\pi)$ in bin $i$ for a given subvolume and $\bar{d}_i$ is its mean over all subvolumes and phases. The combined covariance matrix is the average of the three projection-specific covariances,

The DESI BGS sample covers approximately 14,000 deg$^2$ over $0.25 < z < 0.35$. In a Planck cosmology, this corresponds to an effective cubic volume with side length $L_{\mathrm{survey}} = 912.55\,h^{-1}\,\mathrm{Mpc}$. Because each subvolume is smaller than the survey volume, we rescale the covariance by the ratio of subvolume to survey volume.

The values of $\sigma_{n_g}$ and $\sigma_{n_h}$ computed from the subvolumes are 0.0030 and 0.0065 dex, respectively, after scaling to the DESI BGS volume. Although the cosmic variance scatter in $\log n_{\mathrm{gal}}$ and $\log n_{\mathrm{halo}}$ is very small, the number densities in realistic surveys have systematic uncertainties from selection effects (e.g., photometric variation or spectroscopic incompleteness) and even from imperfect definition of the survey footprint (e.g., area masked by bright galaxies). We elect to be conservative and set both to 5%, adopting $\sigma_{n_g} = \sigma_{n_h} = 0.022$ dex in our fiducial forecast. We demonstrate in § 7.1 that our forecasts for $F_G$ and the cosmological parameters are insensitive to the values of $\sigma_{n_g}$ and $\sigma_{n_h}$. Varying both from 1% to 10% changes the error on $F_G$ by only 3% and those on $\Omega_m$ and $\sigma_8$ by $\lesssim1.5\%$ (see Figure 12). The constraints on $\log M_{\mathrm{cut}}$, $\sigma_{\log M}$, and $\log M_{\mathrm{min}}$, which are set primarily by the number densities, scale more directly with the assumed uncertainty. Because we are implicitly assuming that $\sigma_{n_g}$ and $\sigma_{n_h}$ are inflated by factors not represented in our simulations, we set the cross-covariance terms $C_{\xi n_g}$ and $C_{\xi n_h}$ to zero. We verify that this choice changes the forecast errors on $F_G$, $\Omega_m$, and $\sigma_8$ by $\lesssim1\%$ and on most HOD parameters by $\lesssim2\%$; the largest effects are on $\sigma_{\log M}$ and $\log M_{\mathrm{min}}$, whose errors increase by 6% and 9%, so neglecting these terms is conservative. We retain $C_{n_g n_h}$, the covariance between the two number densities, rescaled from its measured subvolume value to preserve their measured correlation coefficient ($\rho\approx0.87$) under the assumed 5% variances. Neglecting $C_{n_g n_h}$ would change the errors on $F_G$, $\Omega_m$, and $\sigma_8$ by less than 1% but shift those on $\log M_{\mathrm{min}}$ and $\log M_1$ by $\sim7\%$ in opposite directions. In contrast, the off-diagonal structure of $C_{\xi\xi}$ is essential. A diagonal covariance would underestimate the error on $F_G$ by 24%, shrink that on $\Omega_m$ by a factor of 2.6, and inflate that on $\sigma_8$ by 42%.

With $C_{\xi n_g}$ and $C_{\xi n_h}$ set to zero, the covariance matrix in equation (15) becomes block diagonal, with $C_{\xi\xi}$ in one block and the two number density variances and their covariance in the other. We apply the Hartlap correction [89] to debias the inverse of the $C_{\xi\xi}$ block only,

where $N_s = 1600$ (25 phases $\times$ 64 subvolumes) and $N_d = 300$ (10 bins in $r_p \times$ 30 bins in $r_\pi$), giving a correction factor of 1298/1599 $\approx0.812$. We do not apply the Hartlap correction to the number density block, since its values are set to an assumed external uncertainty rather than estimated from a finite number of realizations. We invert the two blocks independently and combine them to form the full inverse covariance used in equation (14).

In a small number of bins at large $r_\pi$, $\xi_{\mathrm{cg}}^{s} \leq0$ in some subvolumes. We set $\log\xi$ to a floor value in these cases, which inflates the variance of these bins and effectively removes them from the Fisher analysis.

Figure 8 shows the covariance structure of $\log\xi_{\mathrm{cg}}^{s}$. The left panel shows $\sigma_{\log\xi_{\mathrm{cg}}^{s}}$ as a function of $r_p$, where each curve corresponds to a fixed $r_\pi$ slice and is colored by the value of $\log\xi_{\mathrm{cg}}^{s}$ in that bin. At fixed $r_p$, $\sigma_{\log\xi_{\mathrm{cg}}^{s}}$ increases with $r_\pi$, reflecting the decrease in the signal with line-of-sight separation. At small $r_\pi$, where the signal is strong, the uncertainty is below 0.03 dex and rises only gradually towards large $r_p$. Once $r_\pi$ exceeds $\sim20h^{-1}\,\mathrm{Mpc}$, the error grows rapidly in all $r_p$, most significantly in the smallest $r_p$, where $\xi_{\mathrm{cg}}^{s}$ beyond the extent of the FoG is close to zero and becomes negative in some subvolumes. The derivatives in Figure 5 exhibit distinctive behavior in this regime, but these bins carry little weight in the constraints because of their large statistical errors.

The right panel of Figure 8 shows the correlation matrix of $\log\xi_{\mathrm{cg}}^{s}$. The repeating block structure reflects the organization of the data vector into $r_\pi$ slices. Bins at the same $r_p$ are strongly correlated across different $r_\pi$ slices out to $r_\pi\sim20h^{-1}\,\mathrm{Mpc}$, since they sample the same halos and their surrounding structure. These correlations weaken at larger $r_\pi$, where the bins are dominated by noise.

Parameter constraints

We present the marginalized parameter constraints from our Fisher forecast for the DESI BGS at $z \sim0.3$, based on a joint analysis of $\log\xi^{s}_{\mathrm{cg}}(r_p,r_\pi)$, $\log n_{\mathrm{gal}}$, and $\log n_{\mathrm{halo}}$. We refer to this configuration as the Full Forecast. It simultaneously constrains the modified gravity parameter $F_{\mathrm{G}}$, the cosmological parameters $\Omega_{\mathrm{m}}$ and $\sigma_8$, the HOD parameters governing galaxy occupation and velocity bias, and the mass–observable relation parameters that describe cluster selection.

Figure 9 presents the marginalized $1\sigma$ and $2\sigma$ posterior ellipses for the Full Forecast, where the ellipses correspond to the 68.3% and 95.4% confidence regions of the joint two-parameter distributions and the diagonal panels show the marginalized one-dimensional posteriors for each parameter. We list the $1\sigma$ marginalized constraint on each parameter, obtained from the square root of the corresponding diagonal element of the inverse Fisher matrix, at the top of each diagonal panel and in the right column of Table 3.

Parameter	Fiducial value	$1\sigma$ constraint
$\log M_{\mathrm{obs}}^{\mathrm{cut}}$	14.00	0.041
$\sigma_{\log M_c}$	0.25	0.039
$\sigma_{\log M}$	0.60	0.152
$\log M_{\min}$	12.50	0.104
$\log M_1$	13.35	0.068
$\alpha$	1.00	0.085
$\eta_{\mathrm{cen}}$	0.10	0.197
$\eta_{\mathrm{sat,low}}$	1.00	0.286
$\eta_{\mathrm{sat,high}}$	1.00	0.026
$F_{\mathrm{G}}$	1.00	0.033
$\Omega_{\mathrm{m}}$	0.3152	0.010
$\sigma_8$	0.8080	0.010

Our joint data vector yields tight constraints on the cosmological parameters, with $\sigma(\Omega_{\mathrm{m}}) = 0.010$ and $\sigma(\sigma_8) = 0.010$. As previously noted, we consider these forecast constraints preliminary because the derivatives with respect to $\sigma_8$ and $\Omega_{\mathrm{m}}$ are noisy when evaluated using a single simulation volume. We constrain the phenomenological modified gravity parameter $F_{\mathrm{G}}$ to $\sigma(F_{\mathrm{G}}) = 0.033$, demonstrating that the DESI BGS can constrain deviations from GR in the cluster infall velocity field at the 3.3% level.

We now examine the most prominent parameter degeneracies visible in Figure 9, focusing on the physical mechanisms that drive each correlation. In general, the correlations are strongest between pairs of parameters that have similar structure in their derivatives in Figure 5–Figure 7.

As expected, the observable mass threshold $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ and the mass–observable scatter $\sigma_{\log M_{\mathrm{c}}}$ are positively correlated. Increasing the mass threshold increases the mean cluster

mass and bias and decreases the cluster number density, and increasing the log-normal scatter can compensate for all three effects. For similar reasons, the minimum mass and log-normal scatter for central galaxies, $\log M_{\min}$ and $\sigma_{\log M}$, are positively correlated. The satellite occupation parameters, $\log M_1$ and $\alpha$, are positively correlated because of their compensating impact on the number of satellite galaxies in cluster mass halos (recall that raising $\log M_1$ lowers the amplitude of $\langle N_{\mathrm{sat}}(M)\rangle$). The values of $\log M_{\min}$ and $\alpha$ are anti-correlated, as one would expect from their similar derivatives in Figure 5. From these derivatives, one might expect a positive correlation between $\log M_{\min}$ and $\log M_1$, but in this case the similar impact of both parameters on $\log n_{\mathrm{gal}}$ (Figure 7) leads to anti-correlation. The cluster selection parameters exhibit only weak correlations with the galaxy occupation parameters because their impact on $\xi_{\mathrm{cg}}^{s}(r_p,r_\pi)$ is very different outside the FoG regime (Figure 5).

By the same logic, we might expect an anti-correlation between $F_{\mathrm{G}}$ and $\log M_{\mathrm{obs}}^{\mathrm{cut}}$, while Figure 9 shows a mild positive correlation. We have confirmed that the expected negative correlation appears when $\sigma_{\log M_c}$ is held fixed. However, when $\sigma_{\log M_c}$ is free, the model raises $\sigma_{\log M_c}$ alongside $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ to preserve $\log n_{\mathrm{halo}}$, introducing lower-mass halos into the sample and reducing velocities. The positive $F_{\mathrm{G}}$–$\log M_{\mathrm{obs}}^{\mathrm{cut}}$ correlation is therefore largely a consequence of the stronger $F_{\mathrm{G}}$–$\sigma_{\log M_c}$ degeneracy.

Our forecast shows a positive correlation between $F_{\mathrm{G}}$ and $\eta_{\mathrm{sat,high}}$, which multiplies the velocity dispersions of satellite galaxies in halos above $M = 10^{13.5}\,h^{-1}\,M_\odot$. Although both parameters have a strong impact on the redshift-space distortion of $\xi_{\mathrm{cg}}^{s}(r_p,r_\pi)$, their effects differ enough in shape that each remains individually well bounded. The derivatives of $\log\xi_{\mathrm{cg}}^{s}(r_p,r_\pi)$ with respect to $\eta_{\mathrm{cen}}$ and $\eta_{\mathrm{sat,low}}$ are much smaller than those for $\eta_{\mathrm{sat,high}}$ (by roughly an order of magnitude), with the result that both parameters are poorly constrained by our observables. Our forecast shows moderate anti-correlation of $F_{\mathrm{G}}$ with both of these parameters. The central and high-mass satellite velocity bias parameters, $\eta_{\mathrm{cen}}$ and $\eta_{\mathrm{sat,high}}$, are also anti-correlated with each other, since both alter the small-scale line-of-sight velocity dispersion of galaxies in cluster halos. The absence of a strong degeneracy between $F_{\mathrm{G}}$ and any of the velocity bias parameters demonstrates that measuring clustering in redshift space does not introduce a fundamental degeneracy between modified gravity and the dynamics within halos, despite both effects altering the velocity field.

The cosmological parameters $\sigma_8$ and $\Omega_{\mathrm{m}}$ are both well constrained in our forecast (at the $\approx0.01$ level) with little correlation between them. This lack of correlation shows that the information content of our data vector is quite different from that of a conventional analysis of cluster number counts versus mass, where a combination $S_8 \approx\sigma_8\Omega_{\mathrm{m}}^{0.5}$ is constrained better than either parameter individually (e.g., [52, 53]). The value of $F_{\mathrm{G}}$ is moderately anti-correlated with those of $\sigma_8$ and $\Omega_{\mathrm{m}}$, which is as expected since increasing any of these parameters leads to larger velocities in cluster infall regions, though $\sigma_8$ and $\Omega_{\mathrm{m}}$ also affect the real-space clustering $\xi_{\mathrm{cg}}^{\mathrm{s}}(r_{\mathrm{p}},r_{\pi})$ (see appendix § 9). We show below (Figure 12) that fixing $\sigma_8$ and $\Omega_{\mathrm{m}}$ only leads to a small improvement in the marginalized error on $F_{\mathrm{G}}$, and fixing $F_{\mathrm{G}}$ leads to only a small improvement in the marginalized errors on $\sigma_8$ and $\Omega_{\mathrm{m}}$, confirming the impression from Figure 9 that these correlations are only moderate. We again caution that our $\sigma_8$ and $\Omega_{\mathrm{m}}$ derivatives are noisy because they come from a single set of AbacusSummit phases, and the forecast errors and error correlations could change with more accurate derivatives.

Discussion

Forecast variations

To better understand what drives the constraints in the Full Forecast, we consider a suite of variations in which we modify the composition of the data vector or hold specific parameters fixed. We first examine how the constraints change as we remove observables from the data vector, isolating the individual contributions of $\log\xi_{\mathrm{cg}}^{\mathrm{s}}$, $\log n_{\mathrm{gal}}$, and $\log n_{\mathrm{halo}}$. We then fix subsets of parameters to identify which degeneracies most limit the constraint on $F_{\mathrm{G}}$.

Figure 10 presents the marginalized $1\sigma$ contours for the cluster selection, modified gravity, and cosmological parameters as we progressively add observables to the data vector, with all parameters free in each configuration. The constraint on $F_{\mathrm{G}}$ degrades from $\sigma(F_{\mathrm{G}}) = 0.033$ in the Full Forecast to $\sigma(F_{\mathrm{G}}) = 0.043$ when we use only $\log\xi_{\mathrm{cg}}^{\mathrm{s}}$. The dominant gain from including the number densities comes from $\log n_{\mathrm{halo}}$, which alone reduces $\sigma(F_{\mathrm{G}})$ to $0.034$, while adding only $\log n_{\mathrm{gal}}$ reduces it to $0.039$. This hierarchy is consistent with the parameter sensitivities established in § 5.2. As shown there, $\log n_{\mathrm{halo}}$ is strongly sensitive to the cluster selection parameters $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ and $\sigma_{\log M_{\mathrm{c}}}$. Providing a direct constraint on the cluster abundance tightens the constraint on $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ by nearly a factor of two (from $0.077$ to $0.042$) and mitigates the $F_{\mathrm{G}}$–$\sigma_{\log M_{\mathrm{c}}}$ degeneracy identified in § 6.3 as the largest single source of uncertainty on $F_{\mathrm{G}}$. The galaxy number density, $\log n_{\mathrm{gal}}$, by contrast, is primarily sensitive to the HOD parameters, and it provides only moderate additional leverage on $F_{\mathrm{G}}$ beyond what is already encoded in $\log\xi_{\mathrm{cg}}^{\mathrm{s}}$.

A notable consequence of removing $\log n_{\mathrm{halo}}$ is the reversal of the $F_{\mathrm{G}}$–$\log M_{\mathrm{obs}}^{\mathrm{cut}}$ correlation visible in the upper-left off-diagonal panel of Figure 10. In the Full Forecast, this correlation is positive, as discussed in § 6.3. When $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ is raised, $\sigma_{\log M_{\mathrm{c}}}$ must simultaneously increase to preserve the observed $\log n_{\mathrm{halo}}$, introducing lower-mass halos into the cluster sample and suppressing the effective velocity field, which requires $F_{\mathrm{G}}$ to compensate upward. Without $\log n_{\mathrm{halo}}$ in the data vector, this compensating mechanism is no longer enforced and the correlation reverts to the physically intuitive negative sign. Raising $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ selects more massive clusters with deeper potential wells and stronger infall signatures, which requires $F_{\mathrm{G}}$ to decrease to reproduce the observed distortion of $\xi_{\mathrm{cg}}^{\mathrm{s}}(r_{\mathrm{p}},r_{\pi})$.

The cosmological parameters $\Omega_{\mathrm{m}}$ and $\sigma_8$ show little sensitivity to the inclusion of either number density, with $\sigma(\Omega_{\mathrm{m}}) = 0.010$ and $\sigma(\sigma_8) = 0.010$–$0.011$ across all configurations. The

constraining power on $\Omega_m$ and $\sigma_8$ derives almost entirely from the distinct signatures they imprint on $\log\xi^{\mathrm{cg}}$.

To identify which parameter degeneracies most limit the constraint on $F_{\mathrm{G}}$, Figure 11 shows how the contours change when specific parameters are held fixed relative to the Full Forecast. We consider three cases: fixing the cosmological parameters $\Omega_{\mathrm{m}}$ and $\sigma_{8}$, fixing $F_{\mathrm{G}}$, and fixing the mass–observable scatter $\sigma_{\log M_{\mathrm{c}}}$. Fixing $\sigma_{\log M_{\mathrm{c}}}$ produces the largest improvement on $F_{\mathrm{G}}$, reducing $\sigma(F_{\mathrm{G}})$ from 0.033 to 0.025, a $\sim23\%$ improvement. This confirms that the $F_{\mathrm{G}}$–$\sigma_{\log M_{\mathrm{c}}}$ degeneracy identified in § 6.3 is the largest limiting factor in the modified gravity constraint in the Full Forecast. Physically, $\sigma_{\log M_{\mathrm{c}}}$ and $F_{\mathrm{G}}$ are degenerate because both modulate the effective infall velocity signal. As Wu et al. [90] discuss, uncertainty in $\sigma_{\log M_{\mathrm{c}}}$ also limits inferences of $\sigma_{8}$ and $\Omega_{\mathrm{m}}$ from cluster weak lensing. However, with direct forward modeling of cluster selection, $\sigma_{\log M_{\mathrm{c}}}$ is no longer a free parameter. Thus, while modeling the cluster selection process in DESI will introduce complications, it may improve constraints on $F_{\mathrm{G}}$. Fixing the cosmological parameters $\Omega_{\mathrm{m}}$ and $\sigma_{8}$ provides a more modest improvement in the $F_{\mathrm{G}}$ constraint, reducing $\sigma(F_{\mathrm{G}})$ from 0.033 to 0.028 ($\sim14\%$). Conversely, the cosmological constraints are unchanged when $\sigma_{\log M_{\mathrm{c}}}$ is fixed and improve by only 6–9% when $F_{\mathrm{G}}$ is fixed.

Figure 12 shows how the $F_{\mathrm{G}}$, $\Omega_{\mathrm{m}}$, and $\sigma_{8}$ errors change in many variations of our standard Full Forecast. Rows 1–7 summarize the results of Figure 10 and Figure 11. Rows 8–9 will be discussed in the next section. The next four rows show the impact of fixing various combinations of HOD parameters. Fixing $\eta_{\mathrm{sat,high}}$ alone improves the $F_{\mathrm{G}}$ error by 8%. While $\xi^{s}_{\mathrm{cg}}(r_{\mathrm{p}},r_{\pi})$ is much less sensitive to $\eta_{\mathrm{cen}}$ and $\eta_{\mathrm{sat,low}}$, these parameters are also poorly constrained, so fixing them yields a noticeable 14% improvement in the $F_{\mathrm{G}}$ error. Fixing the HOD occupation parameters $\sigma_{\log M}$, $\log M_{\mathrm{min}}$, $\log M_{1}$, and $\alpha$ leads to a substantial improvement in the errors on $F_{\mathrm{G}}$, $\Omega_{\mathrm{m}}$, and $\sigma_{8}$, by factors of 1.9, 1.8, and 1.3, respectively, and fixing all the HOD parameters improves the errors by factors of 2.2 ($F_{\mathrm{G}}$), 1.9 ($\Omega_{\mathrm{m}}$), and 1.8 ($\sigma_{8}$). In a real analysis, HOD parameters cannot be fixed a priori, but the occupation parameters can usually be well constrained with the projected galaxy correlation function $w_{\mathrm{p,gg}}(r_{\mathrm{p}})$ (e.g., [70]), and it may be possible to constrain $\eta_{\mathrm{cen}}$ and $\eta_{\mathrm{sat}}$ effectively with the redshift-space galaxy auto-correlation $\xi^{s}_{\mathrm{gg}}(r_{\mathrm{p}},r_{\pi})$. Thus, an analysis that includes the galaxy auto-correlation in the data vector might approach the fixed-HOD errors on $F_{\mathrm{G}}$, $\Omega_{\mathrm{m}}$, and $\sigma_{8}$. However, we have not included assembly bias parameters in our HOD model, and these could degrade the constraints. We leave the investigation of extended data vectors and assembly bias to future work.

The next two rows of Figure 12 show that reducing the assumed fractional errors on $\log n_{\mathrm{gal}}$ and $\log n_{\mathrm{halo}}$ from 5% to 1%, or increasing them from 5% to 10%, makes almost no difference to the $F_{\mathrm{G}}$, $\Omega_{\mathrm{m}}$, and $\sigma_{8}$ errors, which change by about 2% at most. Omitting these number densities from the data vector entirely does degrade the $F_{\mathrm{G}}$ constraint, as shown in Figure 10, but the precision on $F_{\mathrm{G}}$ is not sensitive to the $\log n_{\mathrm{gal}}$ and $\log n_{\mathrm{halo}}$ uncertainties over this fairly wide range. The last row of Figure 12 shows a forecast that uses only the diagonal elements of the covariance matrix. This forecast shows a moderate reduction (24%) in the $F_{\mathrm{G}}$ error, a substantial factor of 2.6 reduction in the $\Omega_{\mathrm{m}}$ error, and a 42% increase in the $\sigma_{8}$ error. These differences show that including off-diagonal terms in the covariance is important, but also that moderate inaccuracies in the off-diagonal correlations will not strongly affect the $F_{\mathrm{G}}$ error.

Constraining Power of the Infall Velocity Field

Zu & Weinberg [7] presented a method to estimate infall velocity profiles $v_{\mathrm{c}}(r)$ from $\xi^{s}_{\mathrm{cg}}(r_{\mathrm{p}},r_{\pi})$, and Zu et al. [8] showed that the extracted profiles could provide a test of modified gravity

models. In this paper, by contrast, we directly model $\xi_{\mathrm{cg}}^{s}(r_p,r_\pi)$ as the primary observable.

We now wish to determine whether the constraining power in $F_{\mathrm{G}}$ actually derives from the infall velocity profile and whether the constraints are degraded relative to those that would apply if we could measure $v_{\mathrm{rad}}(r)$ directly (which in practice we cannot).⁴ To address this question, we compare three data vector configurations: the Full Forecast using $\log\xi_{\mathrm{cg}}^{s}$, a forecast using the mean radial infall velocity profile $v_{\mathrm{rad}}(r)$ directly, and a forecast using the real-space cross-correlation function $\log\xi_{\mathrm{cg}}^{r}$. In all cases we include $\log n_{\mathrm{gal}}$ and $\log n_{\mathrm{halo}}$ as observables with 5% (0.022 dex) uncertainties.

Figure 13 presents the marginalized $1\sigma$ contours for the cluster selection, modified gravity, and cosmological parameters under each configuration. We marginalize over the HOD parameters but do not show them here. The constraints on cluster selection parameters are substantially better when using $\log\xi_{\mathrm{cg}}^{s}$ instead of $v_{\mathrm{rad}}(r)$, though $v_{\mathrm{rad}}(r)$ outperforms the real space $\log\xi_{\mathrm{cg}}^{r}$, showing that infall dynamics provide a key part of the information in $\log\xi_{\mathrm{cg}}^{s}$. However, for our purposes $\log M_{\mathrm{obs}}^{\mathrm{cut}}$ and $\sigma_{\log M_{\mathrm{c}}}$ are nuisance parameters, so the constraints on these parameters alone are not especially important.

The key result of figure 13 is that the constraint on $F_{\mathrm{G}}$ from $\log\xi_{\mathrm{cg}}^{s}$ is only modestly weaker than the constraint from a direct measurement of $v_{\mathrm{rad}}(r)$, $\sigma(F_{\mathrm{G}})=0.033$ versus 0.029. This similarity implies that our analysis can extract most of the relevant information about $v_{\mathrm{rad}}(r)$ from $\log\xi_{\mathrm{cg}}^{s}$, even after marginalizing over HOD parameters. The real space $\log\xi_{\mathrm{cg}}^{r}$ contains no information about $F_{\mathrm{G}}$ because it is unaffected by peculiar velocities. In principle, $\log\xi_{\mathrm{cg}}^{s}$ could contain additional information about $F_{\mathrm{G}}$ from the full velocity distribution in infall regions. Still, we expect that most of the useful velocity information resides in $v_{\mathrm{rad}}(r)$, especially after marginalizing over velocity bias parameters. The 13% larger $\sigma(F_{\mathrm{G}})$ from $\log\xi_{\mathrm{cg}}^{s}$ relative to $v_{\mathrm{rad}}(r)$ likely reflects degeneracies between $F_{\mathrm{G}}$ and the HOD and velocity bias parameters, which affect $\log\xi_{\mathrm{cg}}^{s}$ but have little effect on $v_{\mathrm{rad}}(r)$ (Figure 2), even though $\log\xi_{\mathrm{cg}}^{s}$ yields considerably tighter constraints on the cluster selection parameters. The real-space correlation function $\log\xi_{\mathrm{cg}}^{r}$ already provides interestingly tight constraints on $\Omega_{\mathrm{m}}$ and $\sigma_8$, though these improve when adding the dynamical information available in $\log\xi_{\mathrm{cg}}^{s}$, and $v_{\mathrm{rad}}(r)$ plus number densities give constraints similar to or better than those from $\log\xi_{\mathrm{cg}}^{r}$. The summary comparison of constraints from these data vectors is shown in rows 1, 8, and 9 of Figure 12.

The real-space clustering constraints should be regarded with some caution because our HOD model for spatial clustering is somewhat restrictive; in particular, we do not include assembly bias parameters, and we assume that satellite occupation is a power law in halo mass. A more flexible model might show more degeneracy between HOD parameters and $(\sigma_8,\Omega_{\mathrm{m}})$, especially for $\log\xi_{\mathrm{cg}}^{r}$, and could also degrade the $F_{\mathrm{G}}$ constraint. Although $v_{\mathrm{rad}}(r)$ itself is insensitive to the HOD (Figure 2), marginalizing over the HOD parameters already increases $\sigma(F_{\mathrm{G}})$ from $\log\xi_{\mathrm{cg}}^{s}$ by a factor of 2.2 relative to the fixed-HOD case (Figure 12). Conversely, external HOD constraints from galaxy auto-correlations could substantially strengthen it.

Future Work

The results presented here establish the statistical sensitivity achievable with the DESI BGS for testing GR-predicted infall velocities within a Fisher matrix framework. Applying this methodology to real data requires several extensions, spanning improvements to the theoretical parameterization, the forward modeling pipeline, and the observational setup.

The most immediate theoretical extension is to place the $F_{\mathrm{G}}$ constraint in the context of specific modified gravity theories. The $f(R)$ and Galileon models investigated by Zu et al. [8] predict infall velocity enhancements of 15–30% relative to GR at $r = 5\,h^{-1}\,\mathrm{Mpc}$, corresponding to $F_{\mathrm{G}} \approx1.15$–$1.30$. Our Full Forecast constraint of $\sigma(F_{\mathrm{G}}) = 0.033$ would therefore distinguish these models from GR at $\gtrsim4\sigma$. A natural next step is to run dedicated modified gravity simulations and fit for the best-fitting $F_{\mathrm{G}}$ in each case, establishing a direct mapping between our phenomenological parameter and model-specific parameters such as $f_{0}$ in Hu–Sawicki $f(R)$ gravity [10].

More generally, a uniform velocity rescaling is a simplification, as modified gravity theories generically predict enhancements that depend on scale, halo mass, and environment through their screening mechanisms. Chameleon screening (e.g., $f(R)$ gravity) and Vainshtein screening (e.g., Galileon and DGP gravity) screen via qualitatively different physical mechanisms, and these differences are expected to imprint distinct radial shapes on the infall velocity profile beyond a simple change in amplitude. We therefore plan to develop a parameterization that captures these shape differences, rather than a single scale-independent rescaling, and to characterize how this shape information propagates into the anisotropic structure of $\log\xi_{\mathrm{cg}}^{s}(r_{\mathrm{p}}, r_{\pi})$. This would extend our forecast beyond detecting a deviation from GR to distinguishing between specific screening mechanisms. Allowing $F_{\mathrm{G}}$ to vary with scale or halo mass, calibrated directly against modified gravity simulations, is a first step toward this more flexible parameterization.

Applying this methodology to real DESI BGS data will require the ability to predict $\log\xi_{\mathrm{cg}}^{s}(r_{\mathrm{p}}, r_{\pi})$ at all observationally viable points in our parameter space, so that we can find best-fit parameters and their uncertainties. Following the approach successfully applied to galaxy–galaxy lensing and redshift-space galaxy clustering [28, 29], we plan to construct an emulator that interpolates predictions for $\xi_{\mathrm{cg}}^{s}(r_{\mathrm{p}}, r_{\pi})$ across a grid of cosmological and HOD parameters. The derivatives computed in § 5 guide us in prioritizing the parameter dimensions that most strongly affect the observable. We will need to validate the emulator against simulations and mock catalogs before applying it to the data.

A realistic application to DESI BGS will also require a more accurate model of the cluster selection function. Our current model assumes a simple log-normal scatter between observed and true halo mass with isotropic selection, but realistic cluster catalogs derived from optical richness, X-ray luminosity, or the SZ signal exhibit more complex selection functions, including projection effects and mass-dependent scatter that can introduce anisotropic biases in $\xi_{\mathrm{cg}}^{s}(r_{\mathrm{p}}, r_{\pi})$ [42, 44, 47, 91, 92]. Developing a forward model that accurately reproduces the selection function of a DESI BGS group or cluster catalog and validating it against mock catalogs that replicate the photometric and spectroscopic selection of the BGS will be essential for unbiased inference. Our group’s cluster weak lensing studies again provide a useful template, and they show that direct forward modeling of cluster selection from the galaxy population can improve cosmological constraints by effectively removing $\sigma_{\log M_{\mathrm{c}}}$ as a nuisance parameter [47, 93].

Our forecast targets a single redshift slice at $z = 0.3$, but the DESI BGS spans $0.1 < z < 0.5$, offering the possibility of measurements across multiple redshift bins. We expect modified gravity signatures to evolve with redshift, since the growth of structure and the efficiency of screening mechanisms depend on the cosmic epoch, so splitting the BGS into multiple bins would provide an additional discriminant between models and a consistency check against GR.

As noted in § 5.1.5, the derivatives with respect to $\Omega_{\mathrm{m}}$ and $\sigma_{8}$ come from singlepaired simulation phases and are consequently noisier than the HOD and mass–observable derivatives, introducing uncertainty into our cosmological forecasts. We hope that constructing an emulator that interpolates across a range of cosmological parameters will reduce this noise and improve the reliability of the $\Omega_{\mathrm{m}}$ and $\sigma_8$ constraints. We also need to test the possible impact of a more flexible HOD parameterization, especially the impact of assembly bias or a more general $\langle N_{\mathrm{sat}}(M)\rangle$. As discussed in section 7.2, a more flexible HOD could weaken our constraints on $\sigma_8$ and $\Omega_{\mathrm{m}}$, and it could also degrade the $F_{\mathrm{G}}$ constraint, which is already a factor of 2.2 weaker than in the fixed-HOD case. Conversely, as discussed in section 7.1, adding the galaxy projected and redshift-space autocorrelation functions to the data vector may tighten HOD constraints and, in turn, tighten constraints on $F_{\mathrm{G}}$, $\Omega_{\mathrm{m}}$, and $\sigma_8$.

Conclusions

We have presented a Fisher matrix forecast for the precision with which the DESI BGS can constrain departures from GR in the cluster infall velocity field, using the redshift-space galaxy–cluster cross-correlation function $\xi_{\mathrm{cg}}^{\mathrm{s}}(r_{\mathrm{p}}, r_{\pi})$ as the primary observable. Our analysis simultaneously marginalizes over the HOD parameters governing the galaxy–halo connection, the mass–observable relation parameters describing the cluster selection, and the cosmological parameters $\Omega_{\mathrm{m}}$ and $\sigma_8$, providing a realistic assessment of the statistical sensitivity achievable in the presence of these astrophysical uncertainties. The main conclusions of this work are as follows.

• The DESI BGS can constrain the phenomenological modified gravity parameter $F_{\mathrm{G}}$ to $\sigma(F_{\mathrm{G}}) = 0.033$ at $z \sim0.3$, enabling the detection of deviations from GR in the cluster infall velocity field at the 3.3% level. Based on the results of Zu et al. [8], this sensitivity is sufficient to distinguish $f(R)$ and Galileon gravity models from GR at $\gtrsim4\sigma$, representing a substantial improvement over the $\sim10$–30% uncertainties that limited earlier SDSS analyses.

• The modified gravity parameter $F_{\mathrm{G}}$ is only moderately correlated with the HOD parameters governing galaxy occupation and velocity bias, demonstrating that we can measure the cluster infall velocity signal without introducing a fundamental degeneracy with the galaxy–halo connection. The strongest limitation on the $F_{\mathrm{G}}$ constraint is its degeneracy with the mass–observable scatter $\sigma_{\log M_{\mathrm{c}}}$, which is partially broken by the independent constraint on the cluster abundance. Fixing $\sigma_{\log M_{\mathrm{c}}}$ improves the constraint on $F_{\mathrm{G}}$ by $\sim23\%$.

• The redshift-space cross-correlation function $\log\xi_{\mathrm{cg}}^{\mathrm{s}}$ encodes both the infall velocity information that drives the $F_{\mathrm{G}}$ constraint and the spatial clustering information that constrains the nuisance parameters. Replacing $\log\xi_{\mathrm{cg}}^{\mathrm{s}}$ with a direct measurement of the mean radial infall velocity profile $v_{\mathrm{rad}}(r)$ gives only a modestly tighter constraint on $F_{\mathrm{G}}$ (0.029 versus 0.033) but substantially weaker constraints on the cluster selection parameters, so modeling $\log\xi_{\mathrm{cg}}^{\mathrm{s}}$ directly recovers most of the infall information.

• We forecast constraints on the cosmological parameters $\Omega_{\mathrm{m}}$ and $\sigma_8$ of $\sigma(\Omega_{\mathrm{m}}) = 0.010$ and $\sigma(\sigma_8) = 0.010$, respectively, with their constraining power deriving almost entirely from the distinct radial and anisotropic signatures they imprint on $\log\xi_{\mathrm{cg}}^{\mathrm{s}}$ rather than

from the number density observables. The $F_{\mathrm{G}}$ parameter is only moderately anti-correlated with $\Omega_{\mathrm{m}}$ and $\sigma_{8}$, establishing that we can separate the modified gravity and cosmological constraints without significant mutual degradation. Fixing $\Omega_{\mathrm{m}}$ and $\sigma_{8}$ allows a slightly improved constraint on $F_{\mathrm{G}}$, and fixing $F_{\mathrm{G}}=1$ allows slightly improved constraints on $\Omega_{\mathrm{m}}$ and $\sigma_{8}$.

• The sensitivity demonstrated here motivates a program of follow-up work, including calibration of $F_{\mathrm{G}}$ against dedicated modified gravity simulations to connect our phenomenological parameter to specific theories such as $f(R)$ and Galileon gravity, construction of an emulator for fast likelihood evaluation across the full parameter space, development of a more realistic cluster selection forward model, and measurements across multiple DESI BGS redshift bins to probe the redshift evolution of gravitational growth.

These results establish the DESI BGS as a powerful testbed for gravitational physics on quasi-linear scales, and position $\xi^{s}_{\mathrm{cg}}(r_{\mathrm{p}},r_{\pi})$ as a compelling complement to existing probes of modified gravity. The cluster infall regime probed by $\xi^{s}_{\mathrm{cg}}(r_{\mathrm{p}},r_{\pi})$ occupies a physically distinct region in the landscape of gravity tests, as it lies beyond the screened, virialized cluster interior where fifth forces are suppressed, yet within scales where coherent infall flows are strong and well resolved by spectroscopic surveys. This is the environment where departures from GR are predicted to be most pronounced, and where the gravitational potentials governing light deflection and matter acceleration may diverge in ways that cause lensing and kinematic mass estimates to be mutually inconsistent [6, 8]. Combining $\xi^{s}_{\mathrm{cg}}(r_{\mathrm{p}},r_{\pi})$ measurements with weak gravitational lensing around the same cluster sample would enable a novel gravitational slip test in this regime, probing the ratio of the two scalar metric potentials at few-percent precision in a physically distinct environment from existing $E_{\mathrm{G}}$ measurements [54]. Looking further ahead, the Nancy Grace Roman Space Telescope will provide weak lensing measurements of unprecedented depth and source density over wide fields, enabling a direct cross-correlation between the kinematic infall signal measured in DESI BGS and the lensing mass profiles from Roman. Such a combination would constitute a powerful simultaneous test of whether the gravitational potentials governing light and matter are equal, as GR demands, offering one of the most direct probes of modified gravity on cosmological scales achievable in the near future.

Real-Space Cross-Correlation Function

We compute the real-space derivatives $\partial\log\xi^{r}_{\mathrm{cg}}(r_{\mathrm{p}},r_{\pi})/\partial\theta$ following the same procedure used for the redshift-space derivatives in § 5, applying the centered finite-difference approximation of equation (10) to $\log\xi^{r}_{\mathrm{cg}}(r_{\mathrm{p}},r_{\pi})$ in place of $\log\xi^{s}_{\mathrm{cg}}(r_{\mathrm{p}},r_{\pi})$. As in the redshift-space case, we average the HOD and mass-observable parameter derivatives over the 25 Planck fiducial cosmology phases of AbacusSummit, and compute the $\Omega_{\mathrm{m}}$ and $\sigma_{8}$ derivatives from the paired simulation boxes listed in Table 1. Figure 14 shows the resulting real-space derivatives. The derivatives with respect to the velocity bias parameters and $F_{\mathrm{G}}$ vanish in real space and are not shown. We discuss the interpretation of these derivatives and their relationship to the corresponding redshift-space derivatives in § 5.

Spherical Model for the $F_{\rm G}$ Derivative

To build intuition for the redshift-space distortion of $\xi_{\rm cg}^{s}(r_p,r_\pi)$, we consider an approximate description based on the infall region of a single, spherically symmetric cluster, extending the model presented by Kaiser [17]. As a spherical perturbation evolves, a shell of matter with Lagrangian radius $q$ in the uniform initial conditions moves inward to a radius $r$, and the cumulative overdensity profile is

where $\rho_{\rm int}(r)$ is the mean interior overdensity in units of the cosmic mean density $\bar{\rho}$. For a spherical perturbation in an $\Omega_{\rm m}=1$ universe, the overdensity $1+\delta$ can be written in terms of a timelike variable $\eta_{\rm P}$ as

with collapse occurring at $\eta_{\rm P}=2\pi$ [94]; for a textbook treatment see, e.g., ref. [95]). One can model a spherical overdensity with profile $\rho_{\rm int}(r)=1+\delta(r)$ as a nested set of perturbations with $\eta_{\rm P}$ decreasing outward as implied by equation (B.2). The overdensity $\rho_{\rm vir}=180$, commonly taken to represent the overdensity within the virial radius, corresponds to $\eta_{\rm P}\approx4.8$. The radial peculiar velocity can be written

where negative values denote radial infall. For an $\Omega_m < 1$ universe, we make the approximation of multiplying the non-linear velocity $v_{\mathrm{rad,sph}}$ by the linear growth rate factor $f(\Omega_m)$, and we further multiply by our modified gravity parameter $F_{\mathrm{G}}$ to obtain

We find the dimensionless factor $Q(\rho_{\mathrm{int}})$ by numerically solving for $\eta_{\mathrm{P}}(\rho_{\mathrm{int}})$ using equation (B.2). The value of $Q(\rho_{\mathrm{int}})$ rises from 1.0 at $\rho_{\mathrm{int}}=1$ to about 2.4 at $\rho_{\mathrm{int}}=180$, with a roughly linear dependence on $\log\rho_{\mathrm{int}}$. For $Q(\rho_{\mathrm{int}})=F_{\mathrm{G}}=1$, equation (B.5) is equivalent to the peculiar velocity implied by the Zeldovich approximation [96],

where $\mathbf{r}$ is the final position of a particle with Lagrangian position $\mathbf{q}$.

For a particle with real-space position $(r_p,r_\pi)$ in the projected and line-of-sight directions, the redshift-space position is $s_p=r_p$ and

where $\mu=r_\pi/r$ is the cosine of the angle between $\mathbf{r}$ and the line of sight. Figure 15 shows the real-space (left) and redshift-space (right) positions of particles in a spherical overdensity with profile

with $\rho_{\mathrm{vir}}=180$. The value of $r_{\mathrm{vir}}$ defines the characteristic scale of the model, and a comoving value of $r_{\mathrm{vir}}=1\,h^{-1}\mathrm{Mpc}$ corresponds to a mass $M_{\mathrm{vir}}\approx10^{14}\,h^{-1}M_\odot$. Particles are initially distributed in a grid of Lagrangian cells of size $0.333\,r_{\mathrm{vir}}$, with random positions within each cell to suppress artifacts of a regular grid. They are displaced radially inward to produce the desired density profile, and radial velocities are computed from equation (B.5) with $F_{\mathrm{G}}=1$ and $f(\Omega_m)=(0.485)^{0.55}=0.67$, corresponding to the value of $\Omega_m$ at $z=0.3$ in a flat $\Lambda$CDM model with $\Omega_m(z=0)=0.3$. We set the maximum value of $Q(\rho_{\mathrm{int}})$ to 2.4, the value at $r=r_{\mathrm{vir}}$. We plot particles with $r/r_{\mathrm{vir}}<11.5$, corresponding to $\rho_{\mathrm{int}}>2.36$.

Particles colored red, purple, and cyan are those with radius $r\approx8$, 6, or $2.5\,r_{\mathrm{vir}}$, respectively. In redshift space, the two outer shells are compressed along the line of sight, with compression at $r_p=0$ of about a factor of two for $r=8\,r_{\mathrm{vir}}$ and about a factor of three for $r=6\,r_{\mathrm{vir}}$. For $r=2.5\,r_{\mathrm{vir}}$, peculiar velocities are large enough that particles on the far side of the cluster are shifted to appear on the near side and vice versa. The full particle distribution at $r_p<5\,r_{\mathrm{vir}}$ shows the characteristic “trumpet” pattern identified by Kaiser [17] (see his figure 5), bounded by caustics at which $ds_\pi/dr_\pi=0$. These caustics formally have infinite density, but also infinitesimal volume, and thus have finite mass. At $r_p=0$ the trumpet extends to $r_\pi/r_{\mathrm{vir}}\approx\pm6.2$. For our adopted $\rho_{\mathrm{vir}}$ and $\Omega_m$, a 1D velocity dispersion $\sigma_{\mathrm{1D}}=(GM_{\mathrm{vir}}/3r_{\mathrm{vir}})^{1/2}$ corresponds to $\sigma_{\mathrm{1D}}/H\approx3.8\,r_{\mathrm{vir}}$, so the velocity elongation from infall at $r_p\approx0$ is similar in length to the FoG from the virial velocity dispersion, though it truncates more sharply.

To create a model for the $F_{\mathrm{G}}$ derivative, we average the overdensity profiles for $r_{\mathrm{vir}}=1.0$, 1.2, 1.4, 1.6, 1.8, and $2.0\,h^{-1}\mathrm{Mpc}$, since the density caustics create unrealistically sharp features if we consider a single value of $r_{\mathrm{vir}}$ (equivalent to a single cluster mass). We adopt a finer Lagrangian grid spacing of $0.16\,h^{-1}\mathrm{Mpc}$ to suppress particle noise. In Figure 16, the top row compares the contours of $\xi_{\mathrm{cg}}^{\mathrm{s}}(r_p,r_{\pi})$ from Figure 4 to the redshift-space contours of $\rho(r_p,r_{\pi})/\bar{\rho}-1$ from the spherical model for $F_G=1.0$ and $1.2$, where $\rho(r_p,r_{\pi})/\bar{\rho}$ is the local overdensity (not the mean interior overdensity), and we subtract one to get the quantity most directly comparable to $\xi_{\mathrm{cg}}^{\mathrm{s}}(r_p,r_{\pi})$. We compute densities on the same grid of points used for our $\log\xi_{\mathrm{cg}}^{\mathrm{s}}$ contours and derivative calculations (see the $\sigma_{\log M}$ panel of Figure 14). The density contours show large-scale flattening and a high-density ridge along the caustic locus, as seen in the particle slice, with a transition between these features at $r_p\approx5$–$10\,h^{-1}\mathrm{Mpc}$ and $r_{\pi}<5\,h^{-1}\mathrm{Mpc}$. At $r_p<2\,h^{-1}\mathrm{Mpc}$ the density maximum is displaced from $r_{\pi}=0$ because

the infall velocity of particles within $r_{\rm vir}$ displaces them to large $r_\pi$. Virialization within $r_{\rm vir}$ would change the line-of-sight velocity distribution to be roughly Gaussian, altering the appearance of the FoG feature. The density contours of the spherical model qualitatively resemble the contours of $\xi_{\rm cg}^{s}(r_p,r_\pi)$, but the $\xi_{\rm cg}^{s}$ contours show a more elongated FoG feature, less large-scale flattening, and less crowding of contours in the region of density caustics. Smooth radial infall is the only source of redshift-space distortion in the spherical model, and it is unsurprising that its features are diluted in the more complex 3D structure of the $N$-body simulations.

The bottom row of Figure 16 compares the $F_G$ derivatives of $\log[\rho(r_p,r_\pi)/\bar{\rho}-1]$ to those of $\log\xi_{\rm cg}^{s}(r_p,r_\pi)$. Adding a virial FoG to the spherical model would not change these derivatives because the velocity dispersions within halos are not affected by $F_G$.⁵ In the spherical model, the strongest effect of increasing $F_G$ is to shift caustics outward in $r_\pi$ at a given $r_p$, as one can see by comparing the $F_G=1.0$ and $F_G=1.2$ contours in the top row. As a result, the density derivatives are large and positive along the diagonal ridge that traces the caustics. The derivatives are negative underneath this ridge, in the triangle at small $(r_\pi,r_p)$, because larger peculiar velocities move particles out of this region towards larger $r_\pi$. The increased flattening of large-scale density contours with larger $F_G$ leads to positive derivatives at small $r_\pi$ and $r_p=5$–$15h^{-1}\mathrm{Mpc}$, but the $F_G$ dependence is small at larger $r_\pi$ because the peculiar velocities are small at these distances.

The spherical model explains several of the key qualitative features of the $\xi_{\rm cg}^{s}$ derivatives, though it does not reproduce them quantitatively. Most importantly, it explains why the derivative is largest and positive at small $r_p$ and $r_\pi= 10$–$15\,h^{-1}\mathrm{Mpc}$, because larger $F_G$ stretches the infall features (the “Kaiser trumpet”) along the line of sight. In similar fashion, it explains the triangle of negative derivatives at small $(r_p,r_\pi)$, and large-scale flattening explains the positive derivatives at large $r_p$ and small $r_\pi$. The spherical model is idealized, and it predicts sharper features in the derivatives than the full $N$-body distributions. It also does not explain the structure at large $(r_p,r_\pi)$, where our simple density profile may provide a poor match to the correlation function. These regions carry little weight in our forecasts because the measurement errors are largest there.

Acknowledgments

We thank Tessa Baker, Shaun Cole, and Peder Norberg for helpful discussions that shaped several aspects of this work. ANS was supported by Department of Energy (DOE) grants DE-SC0009913 and DE-SC0025993. CL was supported by an NSF Postdoctoral Fellowship under award 2502789. This work used the AbacusSummit suite of large, high-accuracy cosmological $N$-body simulations, produced with the Abacus $N$-body code [40, 41]. AbacusSummit was run on the Summit supercomputer at the Oak Ridge Leadership Computing Facility through projects AST135 and AST145. Data analysis and visualization were performed with Corrfunc [85, 86], NumPy [97], SciPy [98], Matplotlib [99], and SM⁶. The simulations were analyzed in part using computational resources from the Ohio Supercomputer Center [100]. This research has used NASA’s Astrophysics Data System. This work is dedicated to the loving memory of Czeslaw Kroenke. Niech spoczywa w pokoju.

References