# Dynamics of pairwise motions in the fully non-linear regime in LCDM and
Modified Gravity cosmologies

Mariana Jaber  0000-0001-7507-9516<EMAIL_ADDRESS>Center for Theoretical
Physics, Polish Academy of Sciences, Al. Lotników 32/46, 02-668 Warsaw, Poland
Wojciech A. Hellwing  0000-0003-4634-4442 Center for Theoretical Physics,
Polish Academy of Sciences, Al. Lotników 32/46, 02-668 Warsaw, Poland  Jorge
E. García-Farieta  0000-0001-6667-5471 Instituto de Astrofísica de Canarias,
s/n, E-38205, La Laguna, Tenerife, Spain  Departamento de Astrofísica,
Universidad de La Laguna, E-38206, La Laguna, Tenerife, Spain  Suhani Gupta
0000-0003-0955-6761 Center for Theoretical Physics, Polish Academy of
Sciences, Al. Lotników 32/46, 02-668 Warsaw, Poland  Maciej Bilicki
0000-0002-3910-5809 Center for Theoretical Physics, Polish Academy of
Sciences, Al. Lotników 32/46, 02-668 Warsaw, Poland

###### Abstract

In contrast to our understanding of density field tracers, the modeling of
direct statistics pertaining to the cosmic velocity field remains open to
significant opportunities for improvement. The lack of accurate modeling for
the non-linear domain of pairwise velocities restricts our capacity to fully
exploit the information encoded in this observable. We present a robust
approach for modelling the mean infall velocities, $v_{12}(r,a)$, with broad
applicability spanning sub-megaparsec scales and cosmologies extending beyond
the standard $\mathrm{\Lambda CDM}$ paradigm. Our approach involves solving
the full pair-conservation equation using accurate non-linear power spectrum
descriptions. To assess the robustness of our model, we extend it to
cosmologies beyond the standard $\mathrm{\Lambda CDM}$, in particular, the Hu-
Sawicki $f(R)$-gravity and Dvali-Gabadadze-Porrati (DGP) modified gravity
models. Remarkably, our predictions for pairwise velocities of dark matter
particles at kilo-parsec scales exhibit excellent agreement with N-body
simulations throughout the entire dynamical range ($0.1\lesssim\xi\lesssim
1000$, or $r\geq 0.4$$h^{-1}\,\mbox{Mpc}$). Furthermore we show that different
gravity models leave distinct signatures in the shape and dynamics of the mean
pairwise velocities, providing a potent test of cosmological gravity laws.

##  I Introduction

Pairwise velocities and the large-scale velocity fields serve as crucial
statistics for large-scale structure (LSS), offering valuable insights that
complement the widely utilized density-based measurements. However, unlike our
understanding of cosmic density, the modelling and understating of direct
statistics of the cosmic velocity field still have considerable room for
improvement. This deficiency is one of the main reasons why the rich
information encoded in the velocity field and its potential as a cosmological
probe is still largely untapped.

Various approaches have been employed to model the non-linear regime of
structure formation in the density field. These include the standard
perturbation theory (SPT) (see [1, 2, 3, 4, 5]), regularized perturbation
theory [6], fitting formulae like Halofit and HMCODE (see [7, 8, 9]), and
emulator approaches (see eg [10, 11, 12]).

As we delve into the small-scale features of the Universe, it becomes
increasingly vital to harness the complementary information embedded in
peculiar velocity fields. Pairwise velocities, measuring the typical relative
velocity between objects at specific separations, contain information about
the peculiar motions induced by gravity. As such the information encoded in
the moments of the pairwise velocity distribution may be considered as crucial
for modelling and deciphering the non-linear regime.

The velocity field data has been already effectively utilized in various
applications, including redshift-space distortions (RSD) growth-rate estimates
(for instance [13, 14]), bulk-flow measures ([15]), cosmic web analysis (see
[16]), measuring the velocity power spectrum (e.g. [17]), or in field
reconstructions (see [18, 19] or [20] where the authors reconstruct the
density and velocity fields using [21], the most recent update to the Cosmic
flows catalog [22]).

In the late 1970s, a potent method for modeling mean pairwise velocities (MPV)
emerged [23, 24, 25, 26]. This approach led to the application of the
Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy of equations, which
describe the dynamic evolution of self-gravitating particle systems, for
modelling of the large-scale cosmic velocity fields. Within this hierarchy,
the first moment is known as the pair-conservation equation [24, 25],
represented as:

$\frac{\partial\xi}{\partial t}+\frac{1}{x^{2}a}\frac{\partial}{\partial
x}[x^{2}(1+\xi)v]=0.$ (1)

It expresses the conservation of particle pairs separated by a comoving
distance $x$ in terms of the two-point correlation function (2PCF) of density
fluctuations, $\xi(x,t)$, and the relative velocity of pairs, $v$. Here, $t$
denotes cosmological time, and $a$ stands for the scale factor.

The result derived from Eq. (1) was employed by [27] to propose an
interpolation-based ansatz for the mean pairwise velocity (MPV), bridging the
gap between linear and non-linear regimes, and improving the results proposed
by [28]. The physical reason behind the coefficients of [27]’s fitting formula
was later explained in [29] and extended to treat galaxies in [30].

The BBGKY formalism was used to constrain the cosmic matter density parameter
$\Omega_{\rm M}$ (as seen in [31]) and for testing cosmological models,
including those with non-zero curvature [27]. However, its performance
significantly deteriorates on highly non-linear scales.

MPVs capture the dynamics of peculiar motions driven by gravitational
interactions in an expanding background Universe. This sensitivity has
motivated research efforts aimed at using MPVs to identify deviations from GR
in large-scale structure formation (as demonstrated in [32, 33, 34, 35, 36,
37]). The absence of precise modeling for the non-linear regime of pairwise
velocities hinders our ability to fully harness the physical information
inherent in this observable.

In this paper, we introduce a precise modeling approach for pairwise
velocities. Our primary focus lies in enhancing predictions on non-linear
scales. We validate our predictions against N-body simulations, and show the
versatility and strength of this modelling, by stress-testing it for non-
standard cosmologies such as Modified Gravity (MG) scenarios.

We aim to construct our model using well-established and readily available
methodologies rather than introducing new analytical or numerical tools. By
demonstrating the effective application of these methods for precise
predictions on sub-megaparsec separations, we seek to provide the community
with a robust framework for modeling MPVs and using them to test the
underlying gravity model.

This paper is organized as follows: in Sec. II we review the pair conservation
equation, describe the models for non-linear scales included in our work and
our simulation data. In Sec. III we present our results for the MPV in
$\mathrm{\Lambda CDM}$ and for MG scenarios, which we discuss in Sec. IV and
we present our main conclusions in Sec. V.

##  II  Methods

We proceed to describe the analytical models, numerical techniques and
simulation data we used to build and test our pairwise velocity modelling.

###  II.1 Model

Our model for computing the mean pairwise velocities, $v_{12}(x,a)$, relies in
the consistent numerical solution of the following equation [27]:

$\frac{a}{3[1+\xi(x,a)]}\frac{\partial\bar{\xi}(x,a)}{\partial
a}=-\frac{v_{12}(x,a)}{H(a)r}$ (2)

where $\bar{\xi}(r,a)=3x^{-3}\int^{x}_{0}\xi(y,a)y^{2}dy$ is the volume-
averaged 2PCF, $H(a)$ is the Hubble function, and $r=ax$ is the proper
separation between pairs.

Before describing the different approaches for the non-linear clustering, let
us briefly examine the limiting values of Eq. (2). In the case of pairs
separated at $x\ll 1$$h^{-1}\,\mbox{Mpc}$ or, equivalently, high values of
$\xi\gg 1$ we recover the stable clustering regime. This regime consists of
scales where the natural tendency of a pair to collapse, as induced by the
gravitational clustering, dominates over the background Hubble expansion, thus
a regime where the term $(v_{12}(r)+Hr)$ is negative. The opposite case, i.e.
at large $x$ where $\xi(x,t)\ll 1$, is known as the linear-regime (see for
e.g. [26]). Here the growing mode of structure formation dominates and we have
$\xi(r,t)=D(t)^{2}\xi(r,t=0)=D(t)^{2}\xi_{0}(r)$, with which we get:

$v_{12}(x,a)=-\frac{2}{3}axHf\bar{\bar{\xi}}(x,a),$ (3)

where $\bar{\bar{\xi}}(x,a)\equiv\bar{\xi}(x,a)/[1+\xi(x,a)]$ and $f$ is the
logarithmic derivative of the growth function, $f\equiv d\ln D/d\ln a$.

An essential aspect is the transition point where the influence of the Hubble
flow on large scales gives way to the stable clustering regime at smaller
separations, and then further into the virialized regime, where the orbits
inside the collapsed haloes are essentially randomized. These transitions are
of particular significance, especially when considering cosmologies with
different growth-rate parameters, such as MG theories.

###  II.2  Modified gravity

In our study, as a way of stress-testing our model, we include
non-$\mathrm{\Lambda CDM}$ cosmologies. In particular, we explore MG models as
extensions to $\mathrm{\Lambda CDM}$. These MG theories introduce alterations
to the standard framework, leading to distinct phenomenological consequences
that are relevant for our research. Specifically, we focus on two MG theories:

  1. 1.

The normal branch of the Dvali-Gabadadze-Porrati (nDGP) model. This theory
seeks to explain the accelerated expansion of the universe by introducing
extra dimensions that influence gravity on large scales while preserving
standard General Relativity (GR) at small scales. It does so by introducing a
scalar field, called the brane-bending mode, which influences the
gravitational interaction on large scales. This modification allows for a
departure from standard gravitational laws at cosmic distances. [38].

  2. 2.

The Hu-Sawicki form of $f(R)$-gravity [39]. This widely adopted extension of
the Hilbert-Einstein action introduces a functional dependence of the Ricci
scalar ($R$) within the Einstein-Hilbert equations, which characterizes the
curvature of space-time. The function $f(R)$ allows for variations in the
gravitational force as a function of the curvature so that this modification
aims to explain cosmic acceleration without invoking dark energy (for a review
see for instance [40, 41]).

To deviate from GR on cosmological scales while respecting both high-density
and strong-field regime constraints, these MG theories introduce “screening
mechanisms”. These mechanisms are theoretical concepts designed to reconcile
the predictions of MG with experimental observations on cosmological and
astrophysical scales, effectively concealing modifications to Einstein’s field
equations in various environments (the regime of high densities or small
distances). The first family illustrates the Vainshtein screening mechanism
[42], implemented in the nDGP gravity model. The Vainshtein mechanism relies
on the non-linearity of equations of motion for the scalar field, leading to a
suppression of its effects at small scales. On the other hand, the specific
form of $f(R)$ gravity we consider employs the Chameleon screening mechanism
[43]. In the Chameleon mechanism, a scalar field possesses a variable
effective mass dependent on local matter density, allowing it to be screened
in high-density environments while remaining active in low-density regions.
For a thorough review on this type of theories we refer the reader to [44].

###  II.3 Data

For calibration and testing of our pairwise-velocity model, we use data
generated from the suite of DM-only MG $N$-body simulations: elephant(Extended
LEnsing PHysics using ANalytic ray Tracing), [45], as these simulations
provide a good test-bed to study the impact of both the Chameleon and
Vainshtein screening mechanisms on the large scale non-linear clustering of
matter.

This set of $N$-body simulations assumes a $\mathrm{\Lambda CDM}$ background,
and implements on top of this the solution of the scalar field and modified
Einstein equations in the MG models described above: normal branch of the DGP
theory, or nDGP, and the Hu-Sawicki form of $f(R)$. In the first case, the
specific values of the extra parameter are: $Hr_{c}=1,5$ (referred to as N1
and N5, respectively), and for the second model, the strength of the
modification is codified by the present-day value of the derivative of the
function $f(R)$, $|f_{R,0}|$, taking the values $10^{-5}$, and $10^{-6}$
(referred to as F5 and F6, respectively). N1 codifies a model that deviates
more strongly from GR than N5, in a similar way just as F5 deviates more from
GR than F6.

For each model and redshift, we have five independent realisations. The
different MG models have the same fiducial $\mathrm{\Lambda CDM}$ background.
For a detailed explanation of the numerical parameters used in this suite of
simulations we refer the reader to Sec. 2.2 from [46].

For the direct calculation of the pairwise velocities in the simulation data,
we use the positions and velocities of the dark-matter particles and analyse
the snapshots at $z=0,0.3$ and $0.5$.

###  II.4 Non-linear Power spectrum

To model the MPV accurately for sub-Mpc separations, we must address the non-
linear clustering component in Eq. (2). To achieve this, we employ various
proposals for the non-linear power spectrum, denoted as $P_{\text{nl}}(k,z)$.
These $P_{\text{nl}}(k,z)$ are subjected to an inverse Fourier transform (IFT)
to provide the corresponding non-linear two-point correlation function (2PCF),
denoted as $\xi_{\text{nl}}(r,a)$, as required in Eq. (2).

Non-linear clustering in $\mathrm{\Lambda CDM}$:  In our study, we test and
employ the Halofit solution which is implemented in the cosmological Boltzmann
solver CAMB [47, 48]. We fix the cosmological parameters to the values of the
elephant simulations. The Halofit formula is an accurate fitting formula for
the non-linear matter power spectrum, presented in [49] and re-calibrated in
[8]. We choose to use the Halofit model for the $\mathrm{\Lambda CDM}$ case
due to its established accuracy, tested to be up to $\sim 5\%$ for $k\leq
10$$h^{-1}\,\mbox{Mpc}$, and $z\leq 2$, for a variety of $\mathrm{\Lambda
CDM}$ and wCDM cosmologies (see for instance [8]).

Additionally, we also compare our results against the solution from the
Convolution Lagrangian Perturbation Theory (CLPT) [50], which is a non-
perturbative resummation of Lagrangian perturbation theory and provides as an
output, an estimation for $v_{12}(r,a)$.

Non-linear clustering in Modified Gravity: In the context of MG scenarios,
which deviate from the standard paradigm, there is less consensus on the
optimal approach for modeling the non-linear power spectrum. Various proposals
exist, some of which compute it directly using N-body simulations [51, 52],
perturbation theory [53], post-Friedman formalism (PPF) [54], or via the
spherical collapse model [55]. Specifically for $f(R)$ gravity, a MG version
of the standard Halofit, known as “MG-Halofit”, has been introduced [56]. It
is essential to note that these approaches are based on certain assumptions
and have limitations, which restrict their applicability and generality.
Instead, we primarily utilize the following two tools: the standard Halofit
solution and MGCAMB111https://github.com/HAWinther/FofrFittingFunction [57], a
modified version of CAMB specifically designed to handle both the $f(R)$ and
nDGP gravity models. Our choice of background cosmology remains consistent
between the $\mathrm{\Lambda CDM}$ and MG cases, aligning with the setup of
the elephant simulations. Consequently, the primary differences between GR and
the MG models manifest in the growth of perturbations. It is important to
stress that the version of Halofit that we use is agnostic to MG as it was
developed and tested for $\mathrm{\Lambda CDM}$ and wCDM cosmologies [8]. An
in-depth exploration of the limitations to model the non-linear $P(k)$ for MG
models can be found in [58]. However, and as we will demonstrate in following
sections, albeit its lower precision to model non-linear $P(k)$ in MG models,
is already enough for getting percent-level predictions for $v_{12}(r,a)$.

For a quick comparison with our MGCAMB+Halofit solutions, we also include a
recent approach proposed by [58] to compute the non-linear matter power
spectrum in the context of MG. This work introduces a non-linear matter power
spectrum for MG, $P(k)_{\text{nl, MG}}(k)$, in terms of a $\mathrm{\Lambda
CDM}$ non-liner power spectrum, $P(k)_{\text{nl},\Lambda\text{CDM}}(k)$, and a
halo model response function, $\Upsilon(k)$, which has been calibrated against
the elephant suite of simulations. The resulting $P(k)_{\text{nl, MG}}(k)$ can
be expressed as:

$P(k)_{\text{nl, MG}}(k)=\Upsilon(k)\times
P(k)_{\text{nl},\Lambda\text{CDM}}(k).$ (4)

For the $P(k)_{\text{nl},\Lambda\text{CDM}}(k)$ in Eq. (4), we also employed
the Halofit prediction from CAMB, consistent with the cosmological parameters
used in the elephant simulations. The specific components of the halo model
for MG were detailed in [58], where the authors report a 5% accuracy up to
non-linear scales of $k\lesssim 2.5-3$$h^{-1}\,\mbox{Mpc}$ in their resulting
power spectra.

Figure 1:  Pairwise velocity models for DM tracers, $v_{12}(r,a)$, in units of
$\mathrm{km\,s}^{-1}$ and for the $\mathrm{\Lambda CDM}$ cosmology. Points
correspond to velocities from elephant runs, lines correspond to the different
theoretical models: the full solution to Eq. (2) using the linear CAMB matter
power spectrum (dashed green), solution obtained via non-linear Halofit based
power spectrum (solid red), the quasilinear approximation given by eq. Eq. (3)
(dot-dashed purple), and the prediction from CLPT (dotted blue).

##  III Results

In this section, we present the outcomes of our study, which revolves around
the modeling of MPVs in the non-linear regime. To achieve this, we utilize the
previously mentioned methodologies for computing the non-linear matter power
spectrum and solve for $v_{12}(r,a)$ using Eq. (2). Our presentation begins
with a focus on the $\mathrm{\Lambda CDM}$ model, followed by an exploration
of the applicability of our MPV model in the context of MG scenarios.

###  III.1 MPV for the $\mathrm{\Lambda CDM}$ model.

Figure 1 displays the solution for Eq. (2) for $\mathrm{\Lambda CDM}$ across
different panels. Moving from top to bottom, we present the MPV, denoted as
$v_{12}(r,a)$ in units of $\mathrm{km\,s}^{-1}$, at varying redshift values.
The panels labelled with their redshift value: $z=0,0.3,0.5$, encompass our
findings for separations ranging between $r\in 0.05-140$$h^{-1}\,\mbox{Mpc}$,
shown in linear scaling, allowing us to cover both the linear and non-linear
regimes simultaneously. Below these, we zoom in on the region $r\in
0.05-40$$h^{-1}\,\mbox{Mpc}$, now in logarithmic x-axis scaling, providing a
more detailed view of our solution’s behavior in the non-linear regime.
Finally, the panels marked with $\Delta v_{12}$, also in logarithmic x-axis
scaling, depict the ratio of our numerical solution for $v_{12}(r,a)$ with
respect to the simulation data, represented as
$\Delta$$v_{12}$$\equiv(v_{12}-v_{12,\text{sim}})/v_{12,\text{sim}}$. The
simulation MPVs, $v_{12,\text{sim}}(r)$, are determined directly calculating
the projected velocity differences of pairs of dark matter particles. These
simulation results are indicated by black dots, with error bars illustrating
the variance across five independent realizations of each snapshot.

To obtain numerical solutions for $v_{12}(r,a)$, we utilize Eq. (2) with
specific prescriptions. We employ the non-linear output from Halofit,
represented by red solid lines labeled as “BBKGY-Halofit”. We also include the
linear power spectrum from CAMB, shown as green dashed lines and labeled as
“BBKGY-linear”. To facilitate comparison with the results of [27], we present
a solution to Eq. (3) using the Halofit power spectrum as input, indicated by
purple dot-dashed lines and labeled as the “quasi-linear” solution. Lastly, we
include the solution for $v_{12}(r,a)$ obtained from the numerical
implementation of the Convolution Lagrangian Perturbation Theory
(CLPT)222https://github.com/wll745881210/CLPT-GSRSD, providing the real-space
pairwise in-fall velocity in units of $v/aH(a)f$. These results are
represented by the dotted blue lines.

It is worth noting the agreement among all the solutions for separations
$r\geq 80$$h^{-1}\,\mbox{Mpc}$. This agreement is expected since these scales
fall well within the linear regime, where the dynamics of pairwise velocities
are primarily driven by the Hubble expansion. It is important to mention that
the results from CLPT were originally optimized for galaxy clustering
analysis, specifically for extracting the baryonic acoustic oscillations (BAO)
feature from the 2PCF of biased tracers. As a result, CLPT provides a robust
solution around BAO scales ($r_{BAO}\sim 100$ $h^{-1}\,\mbox{Mpc}$). However,
we observe deviations for separations below $r\sim 80$$h^{-1}\,\mbox{Mpc}$.
Both the perturbative CLPT and the “quasi-linear” solution start to diverge
from the simulation data and the other solutions, which becomes especially
noticeable in the middle row panel where we focus on separations below
$r=40$$h^{-1}\,\mbox{Mpc}$. On these smaller scales, the CLPT solution follows
a trajectory similar to the quasi-linear solution up to $r\sim
11$$h^{-1}\,\mbox{Mpc}$, below which it deviates the most within the non-
linear regime. The “BBKGY-linear” solution performs better than the previous
two cases, tracking the evolution of the simulation data within 5$\%$ of
accuracy up to separations of $r\sim 11$$h^{-1}\,\mbox{Mpc}$ across all
redshift values, but it deviates in the intermediate and non-linear regimes.
Lastly, the fully non-linear solution derived from Eq. (2) with $P_{\text{non-
lin}}(k)$ from Halofit closely mirrors the full non-linear dynamics of dark
matter particles in the simulation, achieving an accuracy of approximately
$\sim 10\%$ up to separation of $r\sim 1.1$$h^{-1}\,\mbox{Mpc}$, for $z=0$
and, $r\sim 1$$h^{-1}\,\mbox{Mpc}$, for $z=0.3,0.5$. Our solution exhibits a
consistent trend with redshift when compared to the simulation data, showing a
gradual decrease in the $v_{12}(r,a)$ value at its minimum during earlier
snapshots. We delve deeper into the analysis of this characteristic within the
framework of various MG models.

For the remainder of the paper, unless explicitly specified otherwise, it is
presumed that the solutions presented pertain to the “BBKGY-Halofit” method.

Figure 2:  Pairwise velocities, $v_{12}(r,a)$ in units of
$\mathrm{km\,s}^{-1}$, and models in MG theories. Points correspond to DM
particles from the elephant runs, solid lines to the result from BBKGY-Halofit
prediction, and dotted lines to the solution using Eq. (4), labeled as “BBKGY-
Gupta(2023)”. The results for the GR case are represented in black, the F5
model with turquoise, F6 in blue, N1 in orange and N5 in magenta.  Figure 3:
Source (left-hand side) term of Eq. (2) with a focus on the time derivative of
the volume-averaged correlation function part, for the different gravity
models (GR, F5, F6, N1, and N5), and for $z=0$ (solid lines), $z=0.5$ (dashed
lines). As an inset plot we show the factor multiplying
$\partial\bar{\xi}/\partial a$, for the same models, $z$ values, and $r$
range.

###  III.2 MPV in MG theories.

In our treatment of MG theories, the modifications to GR are contained in the
linear power spectra from MGCAMB, which properly takes into account the
modified gravitational interactions for a Hu-Sawicky $f(R)$-gravity, and the
nDGP model. The non-linear part of the solution, however, is generated from
the standard $\mathrm{\Lambda CDM}$-Halofit fitting formulae. This setting
will be referred to as “BBGKY-Halofit” in our presentation.

The resulting $v_{12}$ models are visualized in Fig. 2 in units of
$\mathrm{km\,s}^{-1}$, which, as done previously for the $\mathrm{\Lambda
CDM}$ case, consists of three rows corresponding to three snapshots from the
simulations: from top to bottom, $z=0,0.3,0.5$. The rows labeled with the
corresponding snapshot value showcase the trend of $v_{12}(r,a)$ in these MG
models over the range of separations $r\in[0.05,40]$$h^{-1}\,\mbox{Mpc}$.
Specifically, the six panels on the left illustrate the results for $f(R)$
variants (F5 and F6), while the six panels on the right depict the results for
the nDGP variants (N1 and N5). The data points with associated errors
represent the simulation outputs obtained directly calculating the projected
velocity differences of pairs of dark matter particles. To facilitate
comparison, we also present the simulation results for GR in all the figures,
denoted by the solid points, as well as the corresponding Halofit solutions to
Eq. (2), in black lines.

In each case, the bottom sub-panel displays the ratio of the solution for
$v_{12}(r,a)$ with respect to the corresponding MG model simulation, expressed
as $\Delta$$v_{12}$$\equiv(v_{12}-v_{12,\text{sim}})/v_{12,\text{sim}}$. The
shaded regions in these subplots represent the corresponding relative errors
associated with the simulation data.

Starting with the plots which illustrate the full evolution of $v_{12}(r,a)$,
we observe that the simulation data for the F6 (blue points) and N5 (magenta
points) models closely track the evolution seen in the GR case (black points).
This trend persists as we move from top to bottom in Fig. 2, indicating that
our model captures the temporal evolution even in these non-$\mathrm{\Lambda
CDM}$ scenarios.

Regarding our model for $v_{12}(r,a)$ (solid lines), we see how our solution
from the “BBGKY-Halofit” model captures the trends of the simulation data in
all cases and for all snapshots. For the case of Hu-Sawicki $f(R)$ gravity, we
find that our model follows more closely the simulation data in the weaker
version of the theory (F6, indicated by green lines and green points). For the
F5 variant, we notice that our solutions deviate from the simulation data at
the larger separations $r>10,5,1$ $h^{-1}\,\mbox{Mpc}$, for the snapshots
$z=0.5,0.3,0$, respectively. To better appreciate how the dynamics of pairwise
velocities is captured by our model, we refer to the ratios $\Delta v_{12}$
shown in the narrow subplots at each panel. Starting from the left side, or
the $f(R)$ model, we notice a constant offset of the “BBKGY-Halofit” solution
at large scales ($r\approx 40$ $h^{-1}\,\mbox{Mpc}$). This offset reaches
$\Delta v_{12}\sim 10\%$ for the F6 case, and it is more prominent ($\Delta
v_{12}\sim 20\%$) for the stronger variant of the model, F5. In contrast, the
solutions obtained with Eq. (4) or “BBKGY-Gupta(2023)”, recover the expected
values at these large separations.

Generally speaking, our model with the “BBKGY-Halofit” solution accurately
captures the dynamics from the $f(R)$ simulations with better than 10$\%$
accuracy for separations of a few $h^{-1}\,\mbox{Mpc}$ (between
$2-10$$h^{-1}\,\mbox{Mpc}$). However, as we approach the fully non-linear
regime or sub-megaparsec separations, there are noticeable deviations from the
data, particularly for F5 at $z=0.3$. The results from using the “BBKGY-
Gupta(2023)” prescription perform better in all the cases mentioned above.

We notice how both solutions perform well for the $z=0.5$ snapshot, achieving
an accuracy of of $10\%$ or better for scales down to $r\sim
0.3(0.15)$$h^{-1}\,\mbox{Mpc}$, for F5 and F6, respectively.

Turning to our nDGP solutions (right column of Fig. 2), we observe a better
agreement between our model and the simulation data, for the specific case of
the “BBKGY-Halofit” solution. In each sub-panel presenting the ratio $\Delta
v_{12}$, we see that our prediction initially shows an offset of $\sim 10\%$
at $40$$h^{-1}\,\mbox{Mpc}$ separations but then converges to the result from
direct calculation in the simulation data for separations between a few and
$10$ $h^{-1}\,\mbox{Mpc}$. At $z=0$ our model captures the non-linear
evolution of $v_{12}(r,a)$ with 10$\%$ accuracy or better for pairs separated
between $3-40$ $h^{-1}\,\mbox{Mpc}$. This level of accuracy is similarly
observed at $z=0.3$, and it improves further at $z=0.5$, where we achieve an
agreement better than 10$\%$ for scales $r\geq 1$$h^{-1}\,\mbox{Mpc}$. For our
model using the “BBKGY-Gupta(2023)” prescription (shown in dotted lines) we
notice that these solutions under-perform in comparison to the “BBKGY-Halofit”
solution, particularly for $z=0.3$, at which the pairwise velocities offset by
10$\%$ or more the direct calculation from our simulations, for separations
$r\sim 5$$h^{-1}\,\mbox{Mpc}$.

To understand better where these differences among gravity models come from,
in Fig. 3 we illustrate the various components of the source term (left-hand
side) of Eq. (2). The clustering itself, denoted by $\xi(r,a)$, exhibits
minimal variation among gravity models (depicted in the inset figure for
redshift values, $z=0$, and $0.5$, for variants F5, F6, N1, N5, and the GR
case). However, its time evolution, represented by $\partial\bar{\xi}/\partial
a$, reveals notable differences across models.

Specifically, the more pronounced modifications of GR, N1, and F5 exhibit a
maximum deviation of approximately $13\%$ from GR, at $z=0$ (shown in solid
lines), in contrast to the milder variants (N5 and F6), which deviate around
$3-4\%$ from the GR prediction at their maximum.

Notably, while at $z=0$, the values for $\partial\bar{\xi}/\partial a$ for
models F6 and N1 coincide, these models display different values at $z=0.5$.
In this case, their relative deviations from GR are approximately $11\%$ and
$14\%$, respectively, while the N5 and F6 variants show deviations of around
$2-4\%$ compared to the GR prediction.

###  III.3  Signatures of MG in the MPVs

Comparing pairwise dynamics across different gravity models reveals
interesting characteristics. As we can see in Fig. 2, the minimum value of
$v_{12}(r,a)$ undergoes variations based on the gravity model. Specifically,
for the stronger deviations from GR, $v_{12}(r,a)$ reaches a more pronounced
minimum compared to scenarios involving weaker modifications of gravity or
adhering to the GR framework.

To quantify this effect, we present the maximum value of $|$$v_{12}(r,a)$$|$,
denoted as $|\hat{v}_{12}|$, as a function of redshift, $z$, in the top row of
Fig. 4. As anticipated, during earlier cosmic times, relative velocities
attain smaller magnitudes in comparison to the present epoch. Across all cases
displayed in the top row of Fig. 4, there is a consistent trend toward smaller
$|\hat{v}_{12}|$ (indicated by shallower curves for $v_{12}(r,a)$ in earlier
snapshots, as seen in Fig. 1–2). The noteworthy observation lies in how this
value varies across distinct gravity models.

First we focus on the values from the simulation data. These are represented
by the error-bar points at each snapshot value: $|\hat{v}_{12}|(z=0,0.3,0.5)$.
The errors express the standard deviation relative to the mean calculated from
all realizations, also shown in the respective shaded regions. The
distinctions between the values of $|\hat{v}_{12}|$ for GR and MG models
become quite apparent. While the value of $|\hat{v}_{12}(z)|$ for the weaker
modifications of GR (the variants F6 and N5 denoted by blue crosses and purple
diamonds, respectively) closely tracks the trend of GR (indicated by black
circles), we see that for the more pronounced MG variants (F5 and N1,
represented by green triangles and orange squares, respectively)
$|\hat{v}_{12}(z)|$ exhibits significant deviations from the GR values. In
particular, at $z=0$, the relative difference, defined as
$\Delta|\hat{v}_{12}|\equiv(|\hat{v}_{12,\text{MG}}|/|\hat{v}_{12,\text{GR}}|-1)$,
reaches $16\%$ (or an increment of $52\mathrm{km\,s}^{-1}$), for the N1 case.
Similarly, the F5 variant has $\Delta|\hat{v}_{12}|\sim 14\%$ with respect to
the GR value (an increment of $\sim 47\mathrm{km/s}$). For the models N5 and
F6, this deviations are negligible with respect of the uncertainties
($\Delta\hat{v}_{12}\approx 3.6,1.2\%$, respectively). For earlier snapshots,
the deviation between N1 or F5 and GR is of the same order of magnitude:
$\Delta|\hat{v}_{12}|\approx 12-13\%$. Importantly, these deviations between
gravity models are larger than the uncertainties in $\Delta|\hat{v}_{12}|$
from the simulation data. Conversely, for the weaker modifications of GR, as
expressed in the F6 and N5 variants, the relative deviations are contained
within the data uncertainties. We add the solutions from our model, shown as
dotted lines.

Figure 4: Prediction of the maximum value of the infall pairwise velocity,
$|\hat{v}_{12}(z)|$, and the stable-clustering crossing scale, $R_{*}$, as
function of redshift. The expectation of a decreasing function for higher
values of $z$ is present for all models; however, we notice differences
compared to GR.  Figure 5: Stable-clustering crossing scale, $R_{*}$, at $z=0$
for the different models used in this work. In both panels we show the term
$-v_{12}(r)$ for each model, and the Hubble flow $rH(z)$, which is identical
in all models. The data points represent the values extracted from the
simulation data with error-bars displaying the variance from five independent
box realizations. In the top panel we show the $v_{12}(r,a)$ for the F5 and F6
variants of $f(R)$-gravity, and the bottom panel shows the N1 and N5 variants
of the nDGP model.

For the GR case, and the N1 and F5 variants, our prediction for
$|\hat{v}_{12}|$ agrees within the data uncertainties, except for N1 at
$z=0.3$, while our prediction for $|\hat{v}_{12}|$ in the case F6 and N5
presents an offset with respect to the simulation data. However, in all cases
we recover the two main trends shown in our data: a decreasing value of
$|\hat{v}_{12}|$ for earlier snapshots, and more importantly, the relative
increment with respect to GR in the different MG scenarios. In particular, we
recover the maximum deviation for the N1 variant, followed by the F5 case,
then the N5 variant, and lastly, the F6 model.

Another interesting aspect to consider is the scale at which each model
transitions into the stable clustering regime. We denote this specific scale
as $R_{*}$ and defined to be such that $aR_{*}H(z)-v_{12}(R_{*})=0$. To
illustrate this point, in Fig. 5, we explicitly show the values for
$-v_{12}(r,z=0)$ for all the models under study, together with the Hubble
velocity, $rH_{0}$ (dashed line). In vertical dotted lines we indicate the
value of $r$ at which both lines cross each other, in other words, the value
for $r=R_{*}$ for each model.

As previously explained, $v_{12}(r,a)$ is dominated by the Hubble flow at
large separations, gradually transitioning into the stable clustering regime
as pairs of objects draw closer ([28, 29]). However, the scale $R_{*}$ at
which this transition takes place varies depending on the underlying gravity
theory. The outcomes of this analysis are presented in the bottom panel of
Fig. 4.

In the bottom panel of Fig. 4 we show the stable-clustering crossing scale
$R_{*}(z)$ for each gravity model. The values from our simulation data are
shown in solid lines with error bars and the corresponding model prediction,
in dotted lines. An interesting result is how this scale increases for the MG
models that deviate more prominently from the GR case. We recover a similar
behaviour as the variation of $|\hat{v}_{12}|$ with the gravity model, namely,
that the strongest deviation in the value of $R_{*}$ is found for the N1
variant, followed by the F5 case, the N5 and lastly, the F6 case. Our model
predictions lie outside the data uncertainties. Nevertheless, we recover the
data trends for each one of the MG cases.

##  IV Discussion

In this paper we have presented the modeling for the evolution of pairwise
velocities, for the linear, quasilinear, and fully non-linear regime. Our
model relies on the proper solution of the clustering in said regimes, and at
its core, it sits on the pair-conservation equation as derived from the BBGKY
hierarchy. The elegance of this approach allowed us to test its validity
beyond the standard framework of GR and $\mathrm{\Lambda CDM}$.

In Fig. 1 we have shown the solution from Eq. (2) and compared it against
known approximations, and also against the direct calculation of the projected
relative velocity differences of pairs in our N-body simulations.

In particular we see that the fully consistent solution of Eq. (2), using the
Halofit template for the the clustering (“BBKGY-Halofit”), provides excellent
agreement with the $\mathrm{\Lambda CDM}$ simulation data throughout all the
range ($r\in[0.05-140]$$h^{-1}\,\mbox{Mpc}$), and importantly so, in the sub-
megaparsec regime.

It is essential to acknowledge the limitations of our direct estimation from
simulation data, with the smallest resolved scales being at
$r=0.05$$h^{-1}\,\mbox{Mpc}$. A recent study [59] delved into the accuracy of
pairwise velocities in N-body simulations, albeit in a scale-invariant
Einstein-de-Sitter cosmology ($\Omega_{tot}=\Omega_{m}=1$). Despite this
limitation, their findings revealed that the direct calculation of projected
velocities between particle pairs in the simulation and the estimation using
the 2PCF in the conservation equation are equivalent. Importantly, their
simulations reached smaller scales than those covered in our work, which
encompasses the range of scales presented in our analysis.

While perturbative solutions, such as CLPT, were proposed to provide accurate
solutions in the BAO regime, they understandably do not match the evolution of
the $v_{12}(r,a)$ below several tens of $h^{-1}\,\mbox{Mpc}$. The ansatz
proposed by the authors of [27] (labeled in our solutions as “Quasilinear”)
shows a better performance in comparison to the perturbative treatment of
CLPT. While the behaviour of $v_{12}(r,a)$ for CLPT shows an nonphysical turn
towards positive values of $v_{12}(r,a)$ in the range of $r\leq
10$$h^{-1}\,\mbox{Mpc}$, the “Quasilinear” approximation correctly shows the
tendency to a decreasing value in the $v_{12}(r,a)$ in such regime.
Importantly too, we can notice how the “Quasilinear” converges to the solution
of the full BBKGY hierarchy using as an input a linear model for the 2PCF (a
solution we named “BBKGY-linear”). This result shows the power of consistently
solving Eq. (2), and its ability to provide better results than with other
approximations.

In the bottom panel of Fig. 1 we show that our numerical solution is within
$10\%$ accuracy for the three snapshots, all the way down to $r\approx
1$$\mbox{Mpc}/h$. In the particular case of $z=0.5$, however, our numerical
solution is consistent with numerical errors (which are below the $10\%$
lines) even at separations of $r\sim 300-400$$\mbox{kpc}/h$. For the
intermediate snapshot, $z=0.3$, this is achieved up to $r\sim
1$$\mbox{Mpc}/h$, and for $z=0$, up to $r\sim 2$$\mbox{Mpc}/h$.

The underlying assumptions for the validity of Eq. (2) are tied to the
conservation of pairs of dark matter particles in an expanding Universe (see
chapter IV of [26] and the discussion in [30]). We check whether the solutions
can be applied to cosmological models beyond $\mathrm{\Lambda CDM}$. In
particular, we have investigated the case when we relax the assumption of
General Relativity as the underlying law of gravity at all scales.

In Fig. 2 we can see the remarkable agreement between the solution of Eq. (2)
and the simulation data obtained for the MG scenarios described in Sec. II. In
this case we present the results of fully solving the BBKGY hierarchy Eq. (2)
with the Halofit prescription to obtain the non-linear power spectrum,
$P_{\text{nl}}(k)$, from the linear power spectrum in the MG scenarios under
consideration, $P(k)_{\text{lin, MG}}(z)$. While, as mentioned previously, the
Halofit ingredients were calibrated to match the non-linear clustering in
$\mathrm{\Lambda CDM}$, our $P(k)_{\text{lin, MG}}$ properly takes into
account the modified gravitational interactions for a Hu-Sawicky
$f(R)$-gravity, and the nDGP model.

In the first MG scenario, the Hu-Sawicki $f(R)$ theory, we found a constant
shift between our solution and the simulation data, for pairs separated a few
tens of $\mbox{Mpc}/h$. While one plausible explanation for this deviation
might stem from the absence of a well-defined large-scale limit for pairwise
velocities, as discussed in [1], an alternative hypothesis implicates the
potential influence of the absent effective screening mechanism in this
intermediate regime (at approximately ($r=30$$h^{-1}\,\mbox{Mpc}$). The
observed constant shift prompts speculation about the nuanced interplay
between the characteristics of the $f(R)$ theory and the larger-scale
dynamics, necessitating further investigation to discern the precise origins
of this discrepancy.

Note however that this offset is not present when we use the solution from
[58] for the $f(R)$ model. Instead, the regime of applicability of our model
is extended using the “BBGKY-Gupta” solution in our prediction for
$v_{12}(r,a)$. We recall that the response function in Eq. (4) was calibrated
against the elephant suite of simulations, to provide a proper modeling for
the non-linear clustering in the MG scenarios under consideration.

On the other hand, for the nDGP case (bottom panel of the second row of Fig.
2) we show that our solution properly models the dynamics for pairs separated
$r\gtrapprox 2$$h^{-1}\,\mbox{Mpc}$ for the snapshots $z=0,0.3$, and for pairs
separated at $r\gtrapprox 1$$h^{-1}\,\mbox{Mpc}$, for $z=0.5$. We compare our
model against the numerical results obtained from the simulation data, and
shown as shaded regions around the value
$\Delta$$v_{12}(r,a)$$\equiv\left[v_{12}-v_{12,\text{sim}}\right]/v_{12,\text{sim}}$.
This level of agreement is achieved with our “BBKGY-Halofit” prediction. In
the case of the “BBKGY-Gupta(2023)” solution the results deviate beyond the
propagated error range of the simulation data across all values of $r$, in the
specific case of N1, at the snapshots $z=0,0.3$, but we obtained a better
match for the value $z=0.5$. The solution for the N5 variant from the “BBKGY-
Gupta(2023)” was found more consistent with the uncertainties from the
simulation particles in the three $z$ values. We speculate that this can be an
effect of the offset in their response functions (figure 2 of [58]) for the
specific case of N1.

In all the cases that we discussed, both in GR and MG cosmologies, we find
better agreement between our solutions and the simulation data for the higher
redshift snapshot ($z=0.5$), in comparison to the subsequent values,
$z=0.3,0$. This can be attributed to the fact that the more recent snapshots
represent highly non-linear stages of the evolution of the Universe,
therefore, the scale at which our prediction is accurate decreases with the
redshift.

From examining the left hand side term of Eq. (2), we conclude that it is the
time evolution of $\bar{\xi}$ what determines the dependence of $v_{12}(r,a)$
on the gravity model.

A perhaps more interesting result can be seen in the results shown in Fig. 4,
where we study the differences in the dynamics of $v_{12}(r,a)$ between the
various gravity prescriptions. We focus on two distinct features, the value of
$v_{12}(r,a)$ at its minimum as a function of redshift, $\hat{v}_{12}(z)$, and
the scale at which the pairs enter the stable clustering regime, $R_{*}(z)$,
instead of being dominated by the Hubble flow. This scale is defined by the
ratio of the streaming velocity to the Hubble expansion. It is important to
remember that the modifications of the underlying gravity model are taken into
consideration only at perturbation level in the elephant suite of simulations,
keeping the same background as in GR for all the models under scrutiny. In
other words, the effect of changes in the expansion rates that such models
have (see for instance [60, 61, 62] where the specific case of $f(R)$-gravity
theories in the cosmological context is discussed) are not considered in our
simulations. Therefore, the differences found in the value of $R_{*}(z)$
express changes arising only from gravitational clustering in these MG
scenarios, which helps to pinpointing changes attributed to the modified force
law rather than by the modified expansion dynamics. The same claim is valid
for the maximum value of $|v_{12}(z)|$, denoted as $|\hat{v}_{12}(z)|$, which
encodes the enhanced infall velocities of pairs of galaxies in this
alternative gravity models, with respect to GR.

Interestingly, we see a connection between the effect of the strength of the
modification of gravity and the changes in the values of $|\hat{v}_{12}(z)|$
for a given $z$:
$\Delta|\hat{v}_{12}|\equiv|\hat{v}_{12,\text{MG}}|/|\hat{v}_{12,\text{GR}}|$.
The stronger modification of GR in the nDGP model, the case N1, displays the
larger difference $\Delta|\hat{v}_{12}|$, followed closely by the F5 case of
the $f(R)$ model.

These findings point to a distinct and recognizable signature from
modifications of gravity on cosmological scales, which can potentially help
detect the clear effect predicted by the implementation of these models,
supporting the results presented in the letter [32], where different
statistics for the pairwise velocities were analyzed. In that work, the
authors focus on the amplitude of $\sigma_{12}$, the line-of-sight centered
pairwise dispersion, also derived under the BBKGY formalism. A direct
comparison to their results is not possible, as their analysis was based on
the HOD (Halo Occupation Distribution) mock galaxies, while we have kept our
analysis on the dark matter particles, for which we are guaranteed that the
conservation of pairs is fulfilled throughout cosmic history in the different
snapshots we analysed. However, some indirect comparison can be made as a part
of the signal in $\sigma_{12}$ that originates, in fact, from
$\Delta$$v_{12}(r,a)$. Their results showed an increase in the amplitude of
$\sigma_{12}$, for F5 with respect to GR, of approximately $25\%$ at
separations of $r=1$$h^{-1}\,\mbox{Mpc}$ and $r=5$$h^{-1}\,\mbox{Mpc}$, the
two cases probed in their analysis. This reinforces the primary finding that
pairwise velocities serve as a potent tool for assessing the validity of
General Relativity (GR) on cosmological scales.

##  V Conclusions

In summary, this study has outlined a robust methodology for accurately
computing pairwise velocities across a wide range of regimes, encompassing
linear, mildly non-linear, and fully non-linear stages. We have demonstrated
that by adequately considering the clustering aspects and their temporal
evolution, our approach can effectively predict the mean pairwise velocities
(MPVs) within the range covered by our clustering model. Our analysis relies
on the fundamental equation derived from the BBGKY hierarchy and employs non-
linear power spectrum models. Through Fourier transformation, we obtain the
non-linear two-point correlation function $\xi_{\text{non-lin}}(r,a)$, which
serves as the input to our core equation.

It is important to emphasize that this approach does not necessitate
simplifications or approximations to generate reliable predictions for
$v_{12}(r,a)$.

Furthermore, we have shown that this model is applicable to gravity models
other than GR. Specifically, we have established that by appropriately
accounting for non-linear clustering in these alternative gravity models, our
equation can seamlessly provide predictions for the infall velocities of pairs
in these modified gravity (MG) scenarios. This highlights the versatility and
robustness of our methodology, making it suitable for a broader range of
cosmological investigations.

Even more, the different MG scenarios have a particular imprint in the
dynamics of in-fall velocities. We have shown that the physical scale at which
the dynamics between pairs is dominated by their gravitational attraction, as
opposed to being dominated by the Hubble flow (crossing scale in the stable-
clustering regime, $R_{*}$), is distinctively affected by the MG families of
theories we analysed. As screening mechanisms are a relatively generic
prediction of viable MG theories, detecting deviations in the pairwise
dynamics as we have described them would be a signature of physics beyond GR.

Furthermore, the pairwise velocity dispersion (PVD), $\sigma_{12}$, appears in
the first moment of the second BBKGY equation [23, 32]. The analysis of PVD is
a vital component of redshift space distortions (RSD) models. To harness the
full potential of ongoing and future data collection cosmological surveys like
DESI [63], Euclid [64], and 4MOST surveys (CRS [65], and 4HS [66]), we must
push the boundaries of our current RSD models. An extension of our model could
result in physically-motivated model for RSD that could be applicable to a
variety of gravity theories.

However, as we delve into smaller separations, the interplay between baryonic
physics and dark matter distribution becomes increasingly significant.
Specifically, the effect of baryons on RSDs, cosmic density and velocity
fields has been thoroughly investigated in [67] by employing comprehensive
hydrodynamical simulations. Their findings reveal that the impact of baryonic
matter on halo and galaxy velocities becomes notable only at very small
separations, typically in the order of kiloparsecs. In fact, for separation
scales smaller than $r\sim 0.6$$h^{-1}\,\mbox{Mpc}$, the deviation in the
total matter power spectrum is anticipated to be less than $20\%$.
Furthermore, over a broad range of distances spanning from $r\sim 1.6$ to $62$
$h^{-1}\,\mbox{Mpc}$, the overall amplitude of the power spectrum deviates in
less than $1\%$ when compared to DM-only simulations. Therefore, the impact of
baryons on the velocities of the DM samples employed in our work remains
negligible over the broad range of scales under consideration.

However, for a realistic detection of such signatures we need to consider a
number of systematic effects. The most immediate one is that we have modelled
the signal directly on the pairs of dark matter particles, which can be
guaranteed to be a closed system therefore ensuring the conservation of pairs.
Before attempting any direct comparison with biased density tracers, we would
need to understand how the history mergers affect the temporal evolution of
the clustering signal, which is the cornerstone of our model. Even more, a
non-linear bias model is indispensable to properly account for the dynamics of
pairs at sub-megaparsec scales. These two directions are immediate steps we
plan to tackle in a follow-up paper.

###### Acknowledgements.

The authors would like to thank Adi Nusser for useful discussions, and to Hans
A. Winther for kindly providing us with mgcamb version for specific forms of
$\mu(a,k)$ and $\gamma(a,k)$ functions implementing our nDGP models. The
authors acknowledge the support of the Polish Ministry of Science and Higher
Education MNiSW grant DIR/WK/2018/12, as well as the research project
“VErTIGO” funded by the National Science Center, Poland, under agreement
number 2018/30/E/ST9/00698. JEGF is supported by the Spanish Ministry of
Universities through a María Zambrano grant with reference UP2021-022, funded
within the European Union-Next Generation EU. JEGF acknowledge the IAC
facilities and the personnel of the Servicios Informáticos Comunes (SIC) of
the IAC. MB is supported by the Polish National Science Center through grants
no. 2020/38/E/ST9/00395, 2018/30/E/ST9/00698, 2018/31/G/ST9/03388 and
2020/39/B/ST9/03494

## References

  * Scoccimarro [2004] R. Scoccimarro, “Redshift-space distortions, pairwise velocities, and nonlinearities,” Phys. Rev. D 70, 083007 (2004), arXiv:astro-ph/0407214 [astro-ph].
  * Taruya et al. [2010] A. Taruya, T. Nishimichi, and S. Saito, “Baryon acoustic oscillations in 2D: Modeling redshift-space power spectrum from perturbation theory,” Phys. Rev. D 82, 063522 (2010), arXiv:1006.0699 [astro-ph.CO].
  * Bernardeau et al. [2002] F. Bernardeau, S. Colombi, E. Gaztañaga, and R. Scoccimarro, “Large-scale structure of the Universe and cosmological perturbation theory,” Phys. Rep. 367, 1 (2002), arXiv:astro-ph/0112551 [astro-ph].
  * Gil-Marín et al. [2012] H. Gil-Marín, C. Wagner, L. Verde, C. Porciani, and R. Jimenez, “Perturbation theory approach for the power spectrum: from dark matter in real space to massive haloes in redshift space,” J. Cosmology Astropart. Phys 2012, 029 (2012), arXiv:1209.3771 [astro-ph.CO].
  * Crocce and Scoccimarro [2006] M. Crocce and R. Scoccimarro, “Renormalized cosmological perturbation theory,” Phys. Rev. D 73, 063519 (2006), arXiv:astro-ph/0509418 [astro-ph].
  * Taruya [2014] A. Taruya, “RegPT: Regularized cosmological power spectrum,” Astrophysics Source Code Library, record ascl:1404.012 (2014), ascl:1404.012.
  * Smith and Angulo [2019] R. E. Smith and R. E. Angulo, “Precision modelling of the matter power spectrum in a Planck-like Universe,” MNRAS 486, 1448 (2019), arXiv:1807.00040 [astro-ph.CO].
  * Takahashi et al. [2012] R. Takahashi, M. Sato, T. Nishimichi, A. Taruya, and M. Oguri, “Revising the Halofit Model for the Nonlinear Matter Power Spectrum,” ApJ 761, 152 (2012), arXiv:1208.2701 [astro-ph.CO].
  * Mead et al. [2015] A. J. Mead, J. A. Peacock, C. Heymans, S. Joudaki, and A. F. Heavens, “An accurate halo model for fitting non-linear cosmological power spectra and baryonic feedback models,” MNRAS 454, 1958 (2015), arXiv:1505.07833 [astro-ph.CO].
  * Heitmann et al. [2014] K. Heitmann, E. Lawrence, J. Kwan, S. Habib, and D. Higdon, “The Coyote Universe Extended: Precision Emulation of the Matter Power Spectrum,” ApJ 780, 111 (2014), arXiv:1304.7849 [astro-ph.CO].
  * Winther et al. [2019] H. A. Winther, S. Casas, M. Baldi, K. Koyama, B. Li, L. Lombriser, and G.-B. Zhao, “Emulators for the nonlinear matter power spectrum beyond $\Lambda$ CDM,” Phys. Rev. D 100, 123540 (2019), arXiv:1903.08798 [astro-ph.CO].
  * Euclid Collaboration and et. al. [2021] Euclid Collaboration and et. al., “Euclid preparation: IX. EuclidEmulator2 - power spectrum emulation with massive neutrinos and self-consistent dark energy perturbations,” MNRAS 505, 2840 (2021), arXiv:2010.11288 [astro-ph.CO].
  * Beutler et al. [2012] F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, G. B. Poole, L. Campbell, Q. Parker, W. Saunders, and F. Watson, “The 6dF Galaxy Survey: z$\approx$ 0 measurements of the growth rate and $\sigma$8,” MNRAS 423, 3430 (2012), arXiv:1204.4725 [astro-ph.CO].
  * Howlett et al. [2015] C. Howlett, A. J. Ross, L. Samushia, W. J. Percival, and M. Manera, “The clustering of the SDSS main galaxy sample - II. Mock galaxy catalogues and a measurement of the growth of structure from redshift space distortions at z = 0.15,” MNRAS 449, 848 (2015), arXiv:1409.3238 [astro-ph.CO].
  * Watkins et al. [2009] R. Watkins, H. A. Feldman, and M. J. Hudson, “Consistently large cosmic flows on scales of 100h-1Mpc: a challenge for the standard $\Lambda$CDM cosmology,” MNRAS 392, 743 (2009), arXiv:0809.4041 [astro-ph].
  * Hoffman et al. [2012] Y. Hoffman, O. Metuki, G. Yepes, S. Gottlöber, J. E. Forero-Romero, N. I. Libeskind, and A. Knebe, “A kinematic classification of the cosmic web,” MNRAS 425, 2049 (2012), arXiv:1201.3367 [astro-ph.CO].
  * Howlett et al. [2017] C. Howlett, L. Staveley-Smith, P. J. Elahi, T. Hong, T. H. Jarrett, D. H. Jones, B. S. Koribalski, L. M. Macri, K. L. Masters, and C. M. Springob, “2MTF - VI. Measuring the velocity power spectrum,” MNRAS 471, 3135 (2017), arXiv:1706.05130 [astro-ph.CO].
  * Ganeshaiah Veena et al. [2023] P. Ganeshaiah Veena, R. Lilow, and A. Nusser, “Large-scale density and velocity field reconstructions with neural networks,” MNRAS 522, 5291 (2023), arXiv:2212.06439 [astro-ph.CO].
  * Lilow and Nusser [2021] R. Lilow and A. Nusser, “Constrained realizations of 2MRS density and peculiar velocity fields: growth rate and local flow,” MNRAS 507, 1557 (2021), arXiv:2102.07291 [astro-ph.CO].
  * Courtois et al. [2023] H. M. Courtois, A. Dupuy, D. Guinet, G. Baulieu, F. Ruppin, and P. Brenas, “Gravity in the local Universe: Density and velocity fields using CosmicFlows-4,” A&A 670, L15 (2023), arXiv:2211.16390 [astro-ph.CO].
  * Tully et al. [2023] R. B. Tully, E. Kourkchi, H. M. Courtois, G. S. Anand, J. P. Blakeslee, D. Brout, T. d. Jaeger, A. Dupuy, D. Guinet, C. Howlett, J. B. Jensen, D. Pomarède, L. Rizzi, D. Rubin, K. Said, D. Scolnic, and B. E. Stahl, “Cosmicflows-4,” ApJ 944, 94 (2023), arXiv:2209.11238 [astro-ph.CO].
  * Turnbull et al. [2012] S. J. Turnbull, M. J. Hudson, H. A. Feldman, M. Hicken, R. P. Kirshner, and R. Watkins, “Cosmic flows in the nearby universe from Type Ia supernovae,” MNRAS 420, 447 (2012), arXiv:1111.0631 [astro-ph.CO].
  * Peebles [2019] P. Peebles, Principles of Physical Cosmology, Vol. 99 (Princeton University Press, 2019).
  * Peebles [1976] P. J. E. Peebles, “The Peculiar Velocity Field in the Local Supercluster,” ApJ 205, 318 (1976).
  * Davis and Peebles [1977] M. Davis and P. J. E. Peebles, “On the integration of the BBGKY equations for the development of strongly nonlinear clustering in an expanding universe.” ApJS 34, 425 (1977).
  * Peebles [1980] P. J. E. Peebles, The Large-Scale Structure of the Universe (Princeton University Press, 1980).
  * Juszkiewicz et al. [1999] R. Juszkiewicz, V. Springel, and R. Durrer, “Dynamics of pairwise motions,” The Astrophysical Journal 518, L25 (1999).
  * Hamilton et al. [1991] A. J. S. Hamilton, P. Kumar, E. Lu, and A. Matthews, “Reconstructing the Primordial Spectrum of Fluctuations of the Universe from the Observed Nonlinear Clustering of Galaxies,” ApJ 374, L1 (1991).
  * Sheth et al. [2001a] R. K. Sheth, L. Hui, A. Diaferio, and R. Scoccimarro, “Linear and non-linear contributions to pairwise peculiar velocities,” MNRAS 325, 1288 (2001a), arXiv:astro-ph/0009167 [astro-ph].
  * Sheth et al. [2001b] R. K. Sheth, A. Diaferio, L. Hui, and R. Scoccimarro, “On the streaming motions of haloes and galaxies,” MNRAS 326, 463 (2001b), arXiv:astro-ph/0010137 [astro-ph].
  * Juszkiewicz et al. [2000] R. Juszkiewicz, P. G. Ferreira, H. A. Feldman, A. H. Jaffe, and M. Davis, “Evidence for a Low-Density Universe from the Relative Velocities of Galaxies,” Science 287, 109 (2000), arXiv:astro-ph/0001041 [astro-ph].
  * Hellwing et al. [2014] W. A. Hellwing, A. Barreira, C. S. Frenk, B. Li, and S. Cole, “Clear and Measurable Signature of Modified Gravity in the Galaxy Velocity Field,” Phys. Rev. Lett. 112, 221102 (2014), arXiv:1401.0706 [astro-ph.CO].
  * Hearin [2015] A. P. Hearin, “Assembly bias &redshift-space distortions: impact on cluster dynamics tests of general relativity.” MNRAS 451, L45 (2015), arXiv:1501.02798 [astro-ph.CO].
  * Gronke et al. [2015] M. Gronke, C. Llinares, D. F. Mota, and H. A. Winther, “Halo velocity profiles in screened modified gravity theories,” MNRAS 449, 2837 (2015), arXiv:1412.0066 [astro-ph.CO].
  * Ivarsen et al. [2016] M. F. Ivarsen, P. Bull, C. Llinares, and D. Mota, “Distinguishing screening mechanisms with environment-dependent velocity statistics,” A&A 595, A40 (2016), arXiv:1603.03072 [astro-ph.CO].
  * Bibiano and Croton [2017] A. Bibiano and D. J. Croton, “Pairwise velocities in the “Running FLRW” cosmological model,” MNRAS 467, 1386 (2017), arXiv:1701.04453 [astro-ph.CO].
  * Valogiannis et al. [2020] G. Valogiannis, R. Bean, and A. Aviles, “An accurate perturbative approach to redshift space clustering of biased tracers in modified gravity,” J. Cosmology Astropart. Phys 2020, 055 (2020), arXiv:1909.05261 [astro-ph.CO].
  * Dvali et al. [2000] G. Dvali, G. Gabadadze, and M. Porrati, “4D gravity on a brane in 5D Minkowski space,” Physics Letters B 485, 208 (2000), arXiv:hep-th/0005016 [hep-th].
  * Hu and Sawicki [2007] W. Hu and I. Sawicki, “Models of f(R) cosmic acceleration that evade solar system tests,” Phys. Rev. D 76, 064004 (2007), arXiv:0705.1158 [astro-ph].
  * Amendola et al. [2007] L. Amendola, R. Gannouji, D. Polarski, and S. Tsujikawa, “Conditions for the cosmological viability of f(R) dark energy models,” Phys. Rev. D 75, 083504 (2007), arXiv:gr-qc/0612180 [gr-qc].
  * Jaime et al. [2012] L. G. Jaime, L. Patino, and M. Salgado, “f(R) Cosmology revisited,” arXiv e-prints , arXiv:1206.1642 (2012), arXiv:1206.1642 [gr-qc].
  * Vainshtein [1972] A. I. Vainshtein, “To the problem of nonvanishing gravitation mass,” Physics Letters B 39, 393 (1972).
  * Khoury and Weltman [2004] J. Khoury and A. Weltman, “Chameleon cosmology,” Phys. Rev. D 69, 044026 (2004), arXiv:astro-ph/0309411 [astro-ph].
  * Clifton et al. [2012] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, “Modified gravity and cosmology,” Phys. Rep. 513, 1 (2012), arXiv:1106.2476 [astro-ph.CO].
  * Alam et al. [2021] S. Alam, A. Aviles, R. Bean, Y.-C. Cai, M. Cautun, et al., “Towards testing the theory of gravity with DESI: summary statistics, model predictions and future simulation requirements,” J. Cosmology Astropart. Phys 2021, 050 (2021), arXiv:2011.05771 [astro-ph.CO].
  * Cautun et al. [2018] M. Cautun, E. Paillas, Y.-C. Cai, S. Bose, J. Armijo, B. Li, and N. Padilla, “The Santiago-Harvard-Edinburgh-Durham void comparison - I. SHEDding light on chameleon gravity tests,” MNRAS 476, 3195 (2018), arXiv:1710.01730 [astro-ph.CO].
  * Lewis et al. [2000] A. Lewis, A. Challinor, and A. Lasenby, “Efficient computation of cosmic microwave background anisotropies in closed friedmann-robertson-walker models,” The Astrophysical Journal 538, 473 (2000).
  * Howlett et al. [2012] C. Howlett, A. Lewis, A. Hall, and A. Challinor, “CMB power spectrum parameter degeneracies in the era of precision cosmology,” Journal of Cosmology and Astroparticle Physics 2012, 027 (2012).
  * Smith et al. [2003] R. E. Smith, J. A. Peacock, A. Jenkins, S. D. M. White, C. S. Frenk, F. R. Pearce, P. A. Thomas, G. Efstathiou, and H. M. P. Couchman, “Stable clustering, the halo model and non-linear cosmological power spectra,” MNRAS 341, 1311 (2003), arXiv:astro-ph/0207664 [astro-ph].
  * Carlson et al. [2013] J. Carlson, B. Reid, and M. White, “Convolution Lagrangian perturbation theory for biased tracers,” MNRAS 429, 1674 (2013), arXiv:1209.0780 [astro-ph.CO].
  * Oyaizu et al. [2008] H. Oyaizu, M. Lima, and W. Hu, “Nonlinear evolution of $f(r)$ cosmologies. ii. power spectrum,” Phys. Rev. D 78, 123524 (2008).
  * Alam and et. al. [2021] S. Alam and et. al., “Towards testing the theory of gravity with desi: summary statistics, model predictions and future simulation requirements,” Journal of Cosmology and Astroparticle Physics 2021, 050 (2021).
  * Koyama et al. [2009] K. Koyama, A. Taruya, and T. Hiramatsu, “Nonlinear evolution of the matter power spectrum in modified theories of gravity,” Phys. Rev. D 79, 123512 (2009).
  * Hu and Sawicki [2007] W. Hu and I. Sawicki, “Parametrized post-friedmann framework for modified gravity,” Phys. Rev. D 76, 104043 (2007).
  * Mead et al. [2015] A. J. Mead, J. A. Peacock, C. Heymans, S. Joudaki, and A. F. Heavens, “An accurate halo model for fitting non-linear cosmological power spectra and baryonic feedback models,” Monthly Notices of the Royal Astronomical Society 454, 1958 (2015), https://academic.oup.com/mnras/article-pdf/454/2/1958/9381114/stv2036.pdf.
  * Zhao [2014] G.-B. Zhao, “Modeling the nonlinear clustering in modified gravity models. i. a fitting formula for the matter power spectrum of f(r) gravity,” The Astrophysical Journal Supplement Series 211, 23 (2014).
  * Hojjati et al. [2011] A. Hojjati, L. Pogosian, and G.-B. Zhao, “Testing gravity with CAMB and CosmoMC,” J. Cosmology Astropart. Phys 2011, 005 (2011), arXiv:1106.4543 [astro-ph.CO].
  * Gupta et al. [2023] S. Gupta, W. A. Hellwing, and M. Bilicki, “Improved analytical modeling of the nonlinear power spectrum in modified gravity cosmologies,” Phys. Rev. D 107, 083525 (2023), arXiv:2301.12016 [astro-ph.CO].
  * Maleubre et al. [2023] S. Maleubre, D. J. Eisenstein, L. H. Garrison, and M. Joyce, “Constraining accuracy of the pairwise velocities in N-body simulations using scale-free models,” MNRAS 525, 1039 (2023), arXiv:2211.07607 [astro-ph.CO].
  * Jaime et al. [2014] L. G. Jaime, L. Patiño, and M. Salgado, “Note on the equation of state of geometric dark energy in $f(r)$ gravity,” Phys. Rev. D 89, 084010 (2014).
  * Jaime et al. [2018] L. G. Jaime, M. Jaber, and C. Escamilla-Rivera, “New parametrized equation of state for dark energy surveys,” Phys. Rev. D 98, 083530 (2018), arXiv:1804.04284 [astro-ph.CO].
  * Jaber et al. [2022] M. Jaber, G. Arciniega, L. G. Jaime, and O. A. Rodríguez-López, “A single parameterization for dark energy and modified gravity models,” Physics of the Dark Universe 37, 101069 (2022), arXiv:2102.08561 [astro-ph.CO].
  * DESI Collaboration [2016] DESI Collaboration, “The DESI Experiment Part I: Science,Targeting, and Survey Design,” arXiv e-prints , arXiv:1611.00036 (2016), arXiv:1611.00036 [astro-ph.IM].
  * Laureijs and et. al. [2011] R. Laureijs and et. al., “Euclid definition study report,” (2011), arXiv:1110.3193 [astro-ph.CO].
  * Richard et al. [2019] J. Richard, J. P. Kneib, C. Blake, A. Raichoor, J. Comparat, T. Shanks, J. Sorce, M. Sahlén, C. Howlett, E. Tempel, R. McMahon, M. Bilicki, B. Roukema, J. Loveday, D. Pryer, T. Buchert, C. Zhao, and CRS Team, “4MOST Consortium Survey 8: Cosmology Redshift Survey (CRS),” The Messenger 175, 50 (2019), arXiv:1903.02474 [astro-ph.CO].
  * Taylor et al. [2023] E. N. Taylor, M. Cluver, E. Bell, J. Brinchmann, M. Colless, H. Courtois, H. Hoekstra, S. Kannappan, C. Lagos, J. Liske, E. Tempel, C. Howlett, S. McGee, K. Said, R. Skelton, M. Gunawardhana, S. Bellstedt, L. Hunt, T. Jarrett, C. Lidman, J. Lucey, S. Alam, M. Bilicki, A. de Graaff, W. Hellwing, S. Leslie, I. Loubser, L. Marchetti, M. Maseda, M. Mogotsi, P. Norberg, A. Sonnenfeld, J. G. Sorce, and 4HS Team, “The 4MOST Hemisphere Survey of the Nearby Universe (4HS),” The Messenger 190, 46 (2023).
  * Hellwing et al. [2016] W. A. Hellwing, M. Schaller, C. S. Frenk, T. Theuns, J. Schaye, R. G. Bower, and R. A. Crain, “The effect of baryons on redshift space distortions and cosmic density and velocity fields in the EAGLE simulation,” MNRAS 461, L11 (2016), arXiv:1603.03328 [astro-ph.CO].