# Fast generation of mock galaxy catalogues in modified gravity models with
COLA

Bartolomeo Fiorini,11footnotetext: Corresponding author. Kazuya Koyama  Albert
Izard  Hans A. Winther  Bill S. Wright  Baojiu Li

###### Abstract

We investigate the viability of producing galaxy mock catalogues with COmoving
Lagrangian Acceleration (COLA) simulations in Modified Gravity (MG) models
employing the Halo Occupation Distribution (HOD) formalism. In this work, we
focus on two theories of MG: $f(R)$ gravity with the chameleon mechanism, and
a braneworld model (nDGP) that incorporates the Vainshtein mechanism. We use a
suite of full N-body simulations in MG as a benchmark to test the accuracy of
COLA simulations. At the level of Dark Matter (DM), we show that COLA
accurately reproduces the matter power spectrum up to $k\sim 1\,h\,{\rm
Mpc}^{-1}$, while it is less accurate in reproducing the velocity field. To
produce halo catalogues, we find that the rockstar halo-finder does not
perform well with COLA simulations. On the other hand, using a simple Friends-
of-Friends (FoF) finder and an empirical mass conversion from FoF to spherical
over-density masses, we are able to produce halo catalogues in COLA that are
in good agreement with those in N-body simulations. To consider the effects of
the MG fifth force on the halo profile, we derive simple fitting formulae for
the concentration-mass and the velocity dispersion-mass relations that we
calibrate using rockstar halo catalogues in N-body simulations. We then use
these results to extend the HOD formalism to modified gravity simulations in
COLA. We use an HOD model with five parameters that we tune to obtain galaxy
catalogues in redshift space. We find that despite the great freedom of the
HOD model, MG leaves characteristic imprints in the redshift space power
spectrum multipoles and these features are well captured by the COLA galaxy
catalogues.

##  1 Introduction

Stage IV dark energy surveys like Euclid222https://www.euclid-ec.org/ and the
Vera C. Rubin Observatory’s LSST333https://www.lsst.org/ are devoted to
improving our understanding of the cosmological model: the investigation of
the Large-Scale Structure (LSS) of the Universe can give insight on the nature
of the dark sector components. Because structure formation is driven by
gravitational forces, possible deviations of gravity from General Relativity
(GR) on cosmological scales may affect our interpretation of Dark Matter (DM)
and dark energy.

On the other hand, a unified description of gravity with the other fundamental
forces of nature is an active area of intense research, with many possible
theoretical candidates but with a lack of experimental probes. Despite that
this is a concern mainly for high energy physics, some of the alternative
theories to GR known as modified gravity (MG) theories manifest a non-trivial
behaviour on cosmological scales, where the theory of gravity has not been
directly tested yet (for a review see [1]).

It is crucial to make accurate theoretical predictions of the properties of
the LSS on non-linear scales, both in terms of comparison to measurements and
to construct realistic covariance matrices, to successfully constrain the
gravity model (see e.g. [2]). This can be achieved by means of cosmological
simulations, but in excess of ${\cal O}(10^{3})$ realisations must be produced
to match the volume of Stage IV surveys and compute an accurate estimate of
the covariance matrices. This poses a serious challenge for full N-body
simulations in MG models due to their high computational cost. Such models
usually have screening mechanisms that hide modified gravity effects on small
scales. To describe the screening mechanism, an additional non-linear equation
needs to be solved in N-body simulations, which significantly slows down the
MG simulations [3].

An alternative to full N-body simulations is to exploit approximate methods
(see [4] for a comprehensive review and [5] for a comparison project in GR),
which lower the computational cost at the expenses of accuracy on non-linear
scales. Amongst these, the COLA method [6, 7, 8, 9] with its extension to MG
[10, 11, 12] offers an interesting compromise between speed-up and accuracy
without introducing any additional free parameter and is therefore ideal to
access the MG information on (mildly) non-linear scales without losing
predictability (see [13] for an alternative approach based on pinocchio).

Having a bridge between observed galaxies and simulated dark matter
distribution makes it possible to extract the LSS information encoded in the
datasets from galaxy surveys. Several prescriptions have been studied in the
literature to do so. A popular and simple empirical prescription is the Halo
Occupation Distribution (HOD), which starts from the DM halos identified in
the density field by a halo-finder algorithm [14] and populates them with
galaxies with a probability distribution conditioned on some halo properties.
In particular, the HOD model proposed in [15] relies on the halo mass and the
density: using the analytical Navarro-Frenk-White (NFW) model [16] for the
halo profile, it is possible to apply this formalism to COLA simulations that
do not resolve the internal halo properties [7]. Due to the non-trivial
dynamics of the screening mechanisms, the internal halo structure in MG
theories can differ from the one in GR and this needs to be considered to
produce realistic galaxy mocks [17, 18].

In this context, we investigate the feasibility of producing mock galaxy
catalogues from COLA simulations in MG producing a simulation suite with mg-
picola [3, 12] and using a suite of full N-body simulations performed by the
ecosmog code [19] to validate the COLA results and to estimate the accuracy in
reproducing galaxy clustering statistics.

This paper is organised as follows. After introducing the MG theories of
interest to this work in Section 2, we discuss the simulation techniques and
suites in Section 3. In Section 4 we investigate the production of halo
catalogues. We then apply the HOD formalism in Section 5 to create mock galaxy
catalogues and study the multipole moments of the galaxy spectrum in redshift
space. We summarise our results in Section 6.

##  2 Modified Gravity

Einstein’s theory of General Relativity (GR) has revolutionised our
understanding of the Universe over the course of the last century and the
$\Lambda$CDM model based on GR has been very successful in reproducing various
cosmological observations. However, GR still lacks a high energy completion
and the $\Lambda$CDM model requires the existence of a highly fine-tuned
cosmological constant to explain the current accelerated expansion of the
Universe. Furthermore, the so-called $H_{0}$ and $S_{8}$ tensions [20, 21, 22]
may be hints of new physics beyond the $\Lambda$CDM model that may be due to
the underlying gravity theory [23].

These yet to be solved problems have motivated theorists to formulate
alternatives to GR often referred to as modified gravity (MG) theories. The MG
theories of interest to cosmology are the ones showing an infrared behavior,
which is often associated with an additional force (referred as the fifth
force), mediated by a scalar field. The fifth force dynamics can be
interpreted as the result of two factors:

  * •

The effective mass of the scalar field, which limits the range of action of
the fifth force.

  * •

The screening mechanism, which may suppress the fifth force [24, 1].

Solar System tests put very stringent constraints on gravity [25]. Because of
this, theories that feature deviations from GR relevant to cosmology commonly
incorporate a screening mechanism that let the theory evade these stringent
constraints. Amongst the possible screening mechanisms we focus on:

  * •

Chameleon screening [26, 27], occurring when the effective action for the
scalar field accounts a potential term.

  * •

Vainshtein screening [28], present when the scalar field action contains a
term with second derivatives of the scalar field.

In this work we consider the Hu-Sawicki $f(R)$ and the Dvali-Gabadadze-Porrati
gravity models, which implement the Chameleon and the Vainshtein screening
mechanisms respectively, and study their impact on the large-scale structure
compared to GR.

###  2.1 Hu-Sawicki $f(R)$

The class of theories where the modifications of gravity can be described by
an additional term to the Einstein-Hilbert action in the form of a function of
the Ricci scalar is known as $f(R)$ theories [29]:

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

This action evades the Solar System constraints by means of the so-called
Chameleon mechanism [26, 27], where the density controls the shape of the
scalar field potential, determining the background value of the scalar field
and its effective mass $m$ (large masses in high density environments and
vice-versa) and the fifth force is consequently screened.

In this paper we are interested in the Hu-Sawicki model with $n=1$ [30]

$f(R)=-2\Lambda+\left|f_{R0}\right|\frac{\bar{R}_{0}^{2}}{R}\,,$ (2.2)

where $\bar{R}_{0}$ is the background curvature today and $f_{R0}$ is the
value of $f_{R}\equiv df/dR$ today. For a small $|f_{R0}|$, the background
cosmology can be approximated as the one given by the $\Lambda$CDM model. In
the linear regime, it is described in Fourier space by the modified Poisson
equation for the gravitational potential $\Psi$ [31]:

$k^{2}\Psi=-4\pi G_{\rm eff}a^{2}\delta\rho,\quad G_{\rm
eff}=G\left(\frac{4+3a^{2}m^{2}/k^{2}}{3+3a^{2}m^{2}/k^{2}}\right),$ (2.3)

where $\delta\rho$ is the fluctuation around the mean energy density and
$m^{2}=\frac{1}{6\left|f_{R0}\right|}\frac{\bar{R}^{3}}{\bar{R}_{0}^{2}}$.
Depending on scale and cosmic time, the effective gravitational constant
$G_{\rm eff}$ in Eq. (2.3) has two limits:

$\displaystyle ma\gg k:G_{\rm eff}\rightarrow G\,,$ (2.4) $\displaystyle ma\ll
k:G_{\rm eff}\rightarrow\frac{4}{3}G\,,$ (2.5)

where in the former case the deviation from GR are suppressed while in the
latter gravity is enhanced by the fifth force.

In the non-linear regime, eq. (2.3) for the Poisson equation does not hold and
it is necessary to solve the non-linear scalar field equation with a
potential. The chameleon screening is determined by the gravitational
potential of the object as well as the environment. The effect of Chameleon
screening can be approximated as [32]

$G_{\rm
eff}=\left(1+\frac{1}{3}\frac{k^{2}}{k^{2}+a^{2}m^{2}}\text{Min}\left[1,f_{s}\,\left|\frac{\Phi_{\rm
crit}(a)}{\Phi_{N}}\right|\right]\right)G,$ (2.6)

where $\Phi_{N}$ is the gravitational potential, $f_{s}$ is a parameter to
tune the strength of the screening that we shall call screening efficiency,
and

$\Phi_{\rm
crit}(a)=\frac{3f_{R0}}{2}\left(\frac{\Omega_{m,0}a^{-3}+4\Omega_{\Lambda,0}}{\Omega_{m,0}+4\Omega_{\Lambda,0}}\right)^{2},$
(2.7)

with $\Omega_{m,0}$ and $\Omega_{\Lambda,0}$ the present-day density
parameters for matter and dark energy.

###  2.2 Dvali-Gabadadze-Porrati

Another way to modify gravity is by defining a theory in a higher-dimensional
space under the assumption that we are living in a 4D brane of the high
dimensional space-time. This class of model is often known as braneworld
gravity and the simplest of these is the Dvali-Gabadadze-Porrati (DGP) model
[33], which is characterised by the action:

$S=\frac{1}{16\pi G_{(5)}}\int
d^{5}x\sqrt{-\prescript{(5)}{}{g}}\prescript{(5)}{}{R}+\int
d^{4}x\sqrt{-g}\left(\frac{1}{16\pi G}R+\mathcal{L}_{m}\right)\,,$ (2.8)

where the index (5) indicates the quantity is the 5-D generalisation of a 4-D
quantity. The normal branch of this model (nDGP) still needs an additional
dark energy component to explain the late-time acceleration and the tuning of
the cross over radius $r_{c}\equiv G_{(5)}/2G$ to recover the standard
cosmology at early times.

Although this model is not theoretically well motivated, it is useful as a toy
model to study the effects of its screening mechanism. Indeed focusing on the
normal branch, under quasi-static conditions, the theory incorporates the
Vainsthein mechanism which screens the fifth force in high curvature
environments, but allows for deviation from GR in low curvature environments.
We assume that the background cosmology is the same as the $\Lambda$CDM model
by introducing the appropriate dark energy in the model [34, 35]. Linearising
the equation of motion, the resulting modified Poisson equation is [36]:

$\nabla^{2}\Psi=4\pi G_{\rm eff}a^{2}\rho\delta,\quad G_{\rm
eff}=G\left(1+\frac{1}{3\beta}\right),$ (2.9)

where $\beta=1-2Hr_{\mathrm{c}}\left(1+\frac{\dot{H}}{3H^{2}}\right)$ and the
overdot denotes the derivative with respect to the physical time. In the case
of the Vainshtein mechanism, the scalar field equation satisfies a non-linear
equation where non-linearity appears in the second derivative of the scalar
field. The Vainshtein screening may be approximated using a density dependent
effective Newton constant. For the DGP model, the approximation is given by
[32]

$G_{\rm
eff}=\left[1+\frac{1}{3\beta}\left(\frac{2(\sqrt{1+x}-1)}{x}\right)\right]G,\quad
x=\frac{8(r_{c}H_{0})^{2}\Omega_{m,0}}{9\beta^{2}}\frac{\rho}{\overline{\rho}},$
(2.10)

where $\rho$ is the average density within a given radius and $\bar{\rho}$ is
the background density. Note that this approximation holds only for a
spherical mass distribution and can violate the condition that the theory
recovers the linear prediction on large scales [37]. Also the screening
depends on the smoothing radius of the density field $\rho$, which needs to be
tuned to match the exact result.

##  3 Simulations

To produce mock catalogues for the modified gravity models previously
described, we start by running dark matter simulations of the large-scale
structure. In particular, we explore the following modified gravity models:
the normal branch DGP model with $H_{0}r_{c}=1$ (N1) and Hu-Sawicki $f(R)$
model with $|f_{R0}|=10^{-5}$ (F5), plus the vanilla GR model.

Cosmological simulations can be accurate up to non-linear scales, but they are
computationally very expensive, as forces between particles need to be
calculated for thousands of time steps typically. The cost is even higher for
MG models, which in addition also need to solve the non-linear equation for
the fifth force. In recent years, the so-called approximate methods have
become a popular alternative to full $N$-body simulations, and they are aimed
at speeding up the creation of a density field at the expense of not resolving
accurately sub-halo scales (see [5, 38, 39, 40] for studies on their accuracy,
and [4] for a review). The simulations in this paper make use of the COmoving
Lagrangian Acceleration (COLA) method [6], in particular, we use mg-picola
[11], that includes gravity models other than GR, and is an ideal tool to
efficiently run simulations for MG.

###  3.1 COLA method

The COLA method [6] uses the Particle-Mesh algorithm, a simple N-body
technique, to evolve the displacement of particles with respect to their
second-order Lagrangian Perturbation Theory (2LPT) positions (which are
obtained analytically). Thanks to this, a small number of time-steps
(typically of $\mathcal{O}(10)$) are enough to accurately recover the DM
density field up to mildly non-linear scales, as well as the halo positions
and masses. This results in a speed-up of a factor $100-1000$ with respect to
full N-body simulations, at the expense of not resolving internal halo
properties. For the production of accurate mock halo catalogues, [8] proposed
an optimal configuration of the parameters in COLA that control the trade-off
between accuracy and computational cost, such as the grid resolution used to
compute forces and the distribution of time steps, which we adopt for this
work.

###  3.2 MG extension to COLA

When MG is incorporated in cosmological simulation, the complexity increases
because the motion of particles is governed by both the Newtonian force and
the fifth force. This can lead to a slow-down of a factor of $\mathcal{O}(10)$
in theories where computing the fifth force requires the solution of non-
linear equations. One way to avoid this significant slow-down is to use
screening approximations [32]: the fifth force dynamics can be described by
means of an effective mass of the scalar field and a coupling determined by
either the value of gravitational potential or its derivatives, depending on
the type of screening mechanism. To extend the COLA method to MG, it is also
necessary to reformulate the Lagrangian Perturbation Theory to take into
account the effect of the fifth force on the dynamics [10, 11, 41].

mg-picola [11] is a publicly available code for cosmological simulations with
the COLA method in modified gravity and it implements the scale-dependent 2LPT
and some screening approximations. Its accuracy has already been tested
against full N-body simulations at the level of the Dark Matter (DM) density
field. In this work, we study the accuracy of mg-picola also for the DM
velocity field, the halo abundance and halo clustering, and we show how it can
be used in combination with a HOD algorithm to generate mock galaxy catalogues
(with some tweaks relative to applications in GR).

###  3.3 Simulation suites

We use a full N-body simulation suite called elephant, which was introduced
and validated in [42], to benchmark the results of our runs based on COLA. The
elephant suite was produced using ecosmog [19], a MG extension of the adaptive
mesh refinement code RAMSES [43] that solves the exact equations for the fifth
force; these simulations were recently used to create mock galaxy catalogues
in [44, 45, 2]444see [46] and [47] for other attempts to create mock galaxy
catalogues in MG theories. Table 2 describes the parameters of the suite such
as the box size, the mass resolution and the gravity models implemented, and
the cosmological parameters employed are:

$\begin{array}[]{ccc}\Omega_{m,0}=0.281\,,&\Omega_{\Lambda,0}=0.719\,,&\Omega_{b,0}=0.046\,,\\\
n_{s}=0.971\,,&\sigma_{8}=0.842\,,&h=0.697\,.\end{array}$ (3.1)

The initial conditions were generated using the Zel’dovich approximation at
$z=49$. For a given realisation, the same initial seed was used for GR, F5 and
N1555Due to an anomaly with one snapshot, we use only 4 realisations for N1 in
this work..

We develop a new suite of simulations with mg-picola that matches the
cosmology and parameters (such as the mass resolution, the box size and the
number of realisations) of the fiducial elephant set. We call this new COLA
suite PIpeline TEsting Run (piter). We stop the simulations at redshift
$z=0.5057$ after 30 timesteps to match the redshift of the elephant snapshot
closest to redshift $z=0.5$. As discussed in the previous section, with the
screening approximation it is possible to tune the behaviour of the screening
mechanism: in Eq. (2.6) we use a screening efficiency $f_{s}=4.0$ for F5 and
in Eq. (2.10) we use a Gaussian smoothing with scale 1 $\,h^{-1}\,{\rm Mpc}$
for N1. We have found these values to provide the best agreement for the dark
matter power spectrum with N-body simulations amongst the values tested. A
summary of the piter simulation details is given in Table 2. Initial
conditions are generated using the second-order Lagrangian perturbation theory
(2LPT) at $z=49$. As in elephant, for a given realisation, the same initial
seed is used for GR, F5 and N1. We also ran a GR COLA simulation using the
same initial condition as one of the realisations in elephant and checked that
our conclusions were not affected by the cosmic variance. For simplicity of
notation, in the following we will just use N-body and COLA to refer to
elephant and piter respectively.

Table 1: elephant simulations Models | GR, F5, N1
---|---
Realisations | 5
Box size | 1024 $\,h^{-1}\,{\rm Mpc}$
$N_{\mathrm{part}}$ | $1024^{3}$
$M_{part}$ | $7.7985\cdot 10^{10}\,h^{-1}\,M_{\odot}$
Domain grid | $1024^{3}$
Refinement criterion | $8$
Initial conditions | Zel’dovich, $z=49$

Table 2: piter simulations Models | GR, F5, N1
---|---
Realisations | 5
Box size | 1024 $\,h^{-1}\,{\rm Mpc}$
$N_{\mathrm{part}}$ | $1024^{3}$
$M_{part}$ | $7.7985\cdot 10^{10}\,h^{-1}\,M_{\odot}$
Force grid | $3072^{3}$
Timesteps | $30$
Initial conditions | 2LPT, $z=49$

We use the public codes PowerI4666https://github.com/sefusatti/PowerI4 and
DTFE [48] to compute and compare the density and velocity-divergence power
spectra in the elephant and piter simulations. In Figures 1 and 2 we compare
the clustering signals in GR (top left panel), the ratio between COLA and
N-body (top right panel) and the enhancement in MG with respect to GR, which
we refer as the boost-factor (bottom panels, F5 on the left and N1 on the
right). The measurements correspond to the mean over 5 realisations and the
shaded regions display the standard deviation (just to have a rough estimate
of the errors to guide the eye). The errors are smaller in the boost-factor
plots because the noise cancels out: the MG and the GR simulations were
started from the same initial conditions and therefore have a very similar
cosmic variance at late redshifts that vanishes in ratios. We find better than
$3\%$ agreement between COLA and N-body for the density power spectra up to
modes of $\sim 1\,h\,{\rm Mpc}^{-1}$ in GR while the boost-factors in F5 and
N1 show agreement within the variance up to $k\sim 1\,h\,{\rm Mpc}^{-1}$. This
shows that COLA simulations capture the MG effects very accurately even with
the use of approximations for screening. The agreement in the velocity
divergence power spectra is within $3\%$ up to $k\sim 0.4\,h\,{\rm Mpc}^{-1}$
in GR and for the boost factor in F5 while the boost factor in N1 agrees
within $2\%$ up to $k\sim 2\,h\,{\rm Mpc}^{-1}$. In GR, after losing power at
$k\sim 1\,h\,{\rm Mpc}^{-1}$, the velocity divergence power spectrum in COLA
shows an enhancement compared with N-body at very high $k$. This is because
COLA has less random motion at small scales that washes out the signal in
N-body more than in COLA. The better agreement in the density power spectra is
expected because the velocity divergence power spectrum is more affected by
non-linear physics.

From the boost-factors, we can already appreciate the intrinsic difference
between F5 and N1: while F5 has a scale-dependent enhancement, N1 is almost
scale-invariant across most scales. This behavior can be explained by
comparing the respective Poisson equations: $G_{\rm eff}$ is a function of $k$
in Eq. (2.3) for $f(R)$ gravity while it is scale-independent in Eq. (2.9) for
DGP.

(a)  (b)  (c)  (d)

Figure 1:  Dark Matter power spectrum comparison between COLA (orange dashed
lines) and N-body (blue solid lines) obtained by taking the average over 5
realisations. The shaded regions represent the standard deviation over the 5
realisations. The power spectrum in GR (top panels) and the boost-factors in
F5 (bottom left) and N1 (bottom right) show agreement within the variance up
to $k\sim 1\,h\,{\rm Mpc}^{-1}$.

(a)  (b)  (c)  (d)

Figure 2:  Dark Matter velocity divergence power spectrum comparison between
COLA (orange dashed lines) and N-body (blue solid lines) obtained by taking
the average over 5 realisations. The shaded regions represent the standard
deviation over the 5 realisations. The power spectrum in GR (top panels) is in
agreement within $3\%$ up to $k\sim 0.4\,h\,{\rm Mpc}^{-1}$ where COLA starts
to lose power. The boost-factor in F5 (bottom left) and N1 (bottom right) show
better than $5\%$ accuracy even deeper in the non-linear regime, up to $k\sim
2\,h\,{\rm Mpc}^{-1}$.

##  4 Halos

Halos can be found by means of specific algorithms often referred to as halo-
finders. Amongst the several halo-finders proposed in literature we focus on
two popular codes: rockstar [49] and the Friends-of-Friends (FoF) finder [50]
(as implemented in nbodykit777https://nbodykit.readthedocs.io). rockstar uses
a 6D phase-space FoF finder and is capable of measuring many halo properties
like the spherical over-density mass, the concentration parameter and the
velocity dispersion just to name a few. The FoF finder instead is purely based
on particles positions.

###  4.1 rockstar halos

We first produce rockstar halo catalogues in COLA and N-body and compare the
cumulative halo mass function (hmf) and the halo power spectrum. We use the
setting to remove unbound particles and keep only parent halos888We found
that, if we do not remove unbound particles, the agreement between N-body and
COLA is substantially worse for the halo power spectrum. The hmf is simply
defined as the number density of halos with mass greater than $M$. Before
measuring the halo power spectrum, we draw a halo sample defined by an
abundance matching procedure, in which the target density is given by a cut at
$10^{13}\,h^{-1}\,M_{\odot}$ in the N-body catalogues for GR. Figures 3 and 4
show the comparisons of the hmf and the halo power spectrum for COLA and
N-body rockstar catalogues respectively, both in GR. The hmf is $25\%$ off for
masses greater than $10^{13}\,h^{-1}\,M_{\odot}$, and the halo power spectrum
ratio is well predicted in the large-scale limit but shows a scale-dependent
deviation, up to $10\%$ at $k\sim 0.4\,h\,{\rm Mpc}^{-1}$. These differences
are likely due to the lower force resolution and accuracy in recovering the
velocity field in COLA simulations.

(a)  (b)

Figure 3: rockstar halo mass function comparison in GR. The hmf in COLA and
N-body (left panel) and their ratio (right panel) show $\sim 25\%$ discrepancy
for masses greater than $10^{13}\,h^{-1}\,M_{\odot}$.

(a)  (b)

Figure 4: rockstar power spectrum comparison between the power spectra of COLA
and N-body in GR. The discrepancy, which is already clear in the left panel,
is shown to reach the $10\%$ with a scale-dependent behaviour in the right
panel.

Given this result, in the rest of the paper, we will use only the FoF finder
to generate halo catalogues from COLA and N-body simulations. However, we
still use rockstar in the N-body simulations to quantify the effects of MG on
the concentration and velocity dispersion of halos in F5 simulations (see
Section 4.3), which we then incorporate in the prescriptions that we use to
populate the COLA simulations with galaxies.

###  4.2 Friends-of-Friends halos

The FoF algorithm links together particles that are separated by a distance
smaller than a certain linking length $b$, expressed in units of the mean
inter-particle distance. We adopt the standard value $b=0.2$. This value,
which is commonly used in N-body simulations, has been shown to also be valid
in COLA simulations (see [7, 51, 8]), although, in some other approximate
methods that use very low-resolution density fields, $b$ has been tweaked, see
[52]. At this point, we are interested in converting FoF masses, $M_{\rm
FoF}$, to spherical over-density (SO) masses, $M_{\rm 200c}$, which we adopt
in the rest of the paper when expressing halo masses. $M_{\rm 200c}$ denotes
the mass enclosed in a spherical region with density 200 times the critical
density. We use the fit in [53] to convert mass definitions with an accuracy
of $5\%$ for most halos in the mass range
$10^{12.5}-10^{15.5}\,h^{-1}\,M_{\odot}$:

$\frac{M_{\rm
FoF}}{M_{200\mathrm{c}}}=\frac{a_{1}}{c_{200}^{2}}+\frac{a_{2}}{c_{200}}+a_{3}\,,$
(4.1)

where the coefficients $a_{1}$, $a_{2}$ and $a_{3}$ depend on the number of
particles (see their Table 1), and $c_{200}$ is the concentration parameter,
which we compute using the empirical formula from [54]:

$\displaystyle\log_{10}c_{200}$ $\displaystyle=a+b\log_{10}\left(\frac{M_{\rm
200c}}{10^{12}\,h^{-1}\,M_{\odot}}\right),$ (4.2) $\displaystyle a$
$\displaystyle=0.520+(0.905-0.520)\exp\left(-0.617z^{1.21}\right),$
$\displaystyle b$ $\displaystyle=-0.101+0.026z.$

The left panel in Figure 5 shows the conversion factor of Eq. (4.1) (for a set
of values corresponding to SO masses of 100, 600, 1000, 3000, 6000 and 10000
particles). There is very little dependence on halo mass, with all the values
within $2\%$ of 0.75. Hence we use the constant value 0.75 to convert the FoF
masses to SO masses. To validate this mass conversion, we compare the mass
function for the FoF halo catalogues $n_{\rm FoF}$ with rescaled
$M_{200\mathrm{c}}=0.75M_{\rm FoF}$ to the mass function as measured by
rockstar $n_{\rm Rock}$ in terms of $M_{\rm 200c}$. The right panel of Figure
5 shows this comparison for N-body GR halo catalogues, and the agreement of
these two mass functions is within $3\%$ for masses grater than
$10^{13}\,h^{-1}\,M_{\odot}$. We thus produce halo catalogues for both the
N-body and COLA simulations using the FoF method and assign the over-density
mass $M_{200\mathrm{c}}=0.75M_{\rm FoF}$.

(a)  (b)

Figure 5: Mass conversion between FoF and $M_{200c}$ spherical over-density
masses. The conversion factor is evaluated for 6 mass values (left panel),
corresponding to the SO mass of halos composed of 100, 600, 1000, 3000, 6000
and 10000 DM particles. All the six points lie within $2\%$ from the value
$0.75$ which we select as a constant to convert the masses of the FoF halos to
SO masses. In the right panel, we use the $N$-body simulations in GR to
compare the FoF hmf resulting after this conversion with the hmf of rockstar
catalogues. The ratio between the two (right panel) shows a better than $3\%$
agreement for masses above $10^{13}\,h^{-1}\,M_{\odot}$.

To validate the COLA halo catalogues we compare against the N-body halo
catalogues for the following two summary statistics: the halo mass function
and the halo power spectrum, which are shown in Figure 6 and Figure 7
respectively. For the halo power spectrum, we draw a halo sample defined by an
abundance matching procedure, in which the target density is given by a cut at
$10^{13}\,h^{-1}\,M_{\odot}$ in the N-body catalogues for GR. We see that COLA
reproduces the N-body mass function with an accuracy better than $5\%$ for
halos with more than 130 particles (or $10^{13}\,h^{-1}\,M_{\odot}$) in GR,
while the hmf boost-factors agree to better than $10\%$ for all the halo
masses. In regards to the halo power spectrum, Figure 7 shows that COLA and
N-body are in agreement within the variance both for GR and for the boost-
factors.

(a)  (b)  (c)  (d)

Figure 6: Halo mass function comparison between COLA (orange dashed lines) and
N-body (blue solid lines) obtained by taking the average over 5 realisations.
The shaded regions represent the standard deviation over the 5 realisations.
The hmf in GR (top panels) show a $\sim 5\%$ accuracy for masses above
$10^{13}\,h^{-1}\,M_{\odot}$. The boost factor in N1 (bottom right) in COLA
agree with that in N-body within the $5\%$ at all masses, while the boost
factor in F5 (bottom left) features an agreement within the $10\%$.

(a)  (b)  (c)  (d)

Figure 7: Halo power spectrum comparison between COLA (orange dashed lines)
and N-body (blue solid lines) obtained by taking the average over 5
realisations. The shaded regions represent the standard deviation over the 5
realisations. The halo power spectrum in GR (top panels) and the boost factors
in F5 and N1 (bottom panels) show agreement within the variance up to $k\sim
0.4\,h\,{\rm Mpc}^{-1}$.

###  4.3 Internal structure of halos

Certain halo properties are needed to create galaxy mock catalogues with the
HOD model discussed in the next section; the phase space distribution of
satellite galaxies is based on the NFW profile, which in turn is completely
determined by the concentration parameter and the mass of the halo.

COLA is purposely designed to accurately determine basic halo quantities
(positions, velocities and masses) but not the internal halo structure, which
would require simulations with a much higher computational cost. Because of
that, in the pipeline developed for this paper, we use several prescriptions
calibrated from high resolution N-body simulations in $\Lambda$CDM to assign
certain halo properties in COLA (see Eq. (4.1) and Eq. (4.2)). Furthermore, MG
may slightly modify some of these relationships (e.g. the concentration
parameter and the velocity dispersion). We use the elephant simulations to
estimate these effects, which we then introduce in the prescription to create
galaxy mock catalogues in the MG models.

In F5 the fifth force is unscreened for halos of mass lower than
$10^{14}\,h^{-1}\,M_{\odot}$, affecting the radial density profile and the
virial equilibrium for such halos. We account for this by implementing
corrections to the concentration parameter and the velocity dispersion.
Detailed studies on the effect of the fifth force on these quantities were
performed in [17] and [18] using higher-resolution simulations, which derived
the fitting formulae for them as a function of halo mass. However, these
studies used a different halo mass definition, thus we use simple hyperbolic
functions to fit concentration parameter and velocity dispersion boost-factors
measured with the rockstar halos in elephant simulations. Although the
concentration parameter measurement for halos with less than $500-1000$
particles, which corresponds to $(4-7)\times 10^{13}\,h^{-1}\,M_{\odot}$ in
elephant, needs to be treated with care, the study using higher-resolution
simulations found a similar increase of the concentration parameter for
unscreened halos [18]. In the case of the concentration parameter we use the
fitting formula:

$\frac{c_{\rm F5}}{c_{\rm GR}}=1+a_{c}\,(1-\tanh{(b_{c}\,\log_{10}(M_{\rm
200c})-M_{c}^{*}})\,.$ (4.3)

The velocity dispersion is obtained by solving the equations for the virial
equilibrium between kinetic energy and gravitational energy of DM in halos:
the solution depends linearly on the square root of the effective
gravitational constant, $G_{\rm eff}$. To fit the transition of the velocity
dispersion between screened and unscreened halos we use the fitting function:

$\frac{V_{\rm rms,F5}}{V_{\rm
rms,GR}}=1+a_{V}\,(1-\tanh{(b_{V}\,\log_{10}(M_{\rm 200c})-M_{V}^{*}})\,.$
(4.4)

The results of the fit compared with the data from rockstar are shown in
Figure 8: the hyperbolic functions correctly represents the transition from
unscreened to screened halos and provides an effective description of the
fifth force effect on these halo properties. The fitted parameters, $a$, $b$
and $M_{*}$, are shown in Figure 8. Note that these fitted parameters are
valid only for $|f_{R0}|=10^{-5}$ and they need to be re-fitted for different
values of $f_{R0}$.

In the case of N1, the Vainshtein mechanism operates efficiently
irrespectively of halo mass and we confirmed that the concentration parameter
and the velocity dispersion in N1 are not modified from GR.

(a)  (b)

Figure 8: Fit of F5 boost factors for concentration parameter (left panel) and
velocity dispersion (right panel). The values measured from rockstar
catalogues (orange dots) are used to fit the functions in Eq. (4.3) and Eq.
(4.4). Only halos with masses above $10^{13}\,h^{-1}\,M_{\odot}$ are used to
perform the fits.

##  5 Galaxy catalogues

The Halo Occupation Distribution (HOD) technique is an empirical model to
connect galaxy and dark matter properties. It allows us to produce galaxy mock
catalogues from low-resolution dark matter simulations that do not resolve the
internal structure of halos.

###  5.1 HOD model

In this paper, due to the limited resolution of elephant simulations, we
create mock galaxies similar to those found in the Baryon Oscillation
Spectroscopic Survey CMASS galaxy sample. These galaxies can be modelled by
the HOD model proposed in [15]999 We used Halotools [55]
(https://halotools.readthedocs.io) to implement the HOD model in our
simulations.. In this model, galaxies are split between centrals and
satellites and the probabilities of a halo of mass $M$ to host either galaxy
type are given by:

$\left\langle
N_{\mathrm{cen}}(M)\right\rangle=\frac{1}{2}\left[1+\operatorname{erf}\left(\frac{\log
M-\log M_{\min}}{\sigma_{\log M}}\right)\right],\quad\left\langle
N_{\mathrm{sat}}(M)\right\rangle=\left\langle
N_{\mathrm{cen}}(M)\right\rangle\left(\frac{M-M_{0}}{M_{1}}\right)^{\alpha}\,,$
(5.1)

where each halo can have only one central galaxy but multiple satellite
galaxies. The parameters $M_{\min}$ and $\sigma_{\log M}$ control respectively
the scale and the width of the transition between probability 0 and 1 to host
a central galaxy. The number of satellites is instead modelled as a power law
with exponent $\alpha$, scale $M_{1}$ and cut-off $M_{0}$, multiplied by the
number of central galaxies.

Once the halo occupation has been determined, we need to assign the position
and velocity for each galaxy. We place central galaxies at the halo centre of
mass position and assign them the velocity of the host halo. On the other
hand, to place satellites we sample a NFW halo density profile [16] centred at
the halo centre of mass and with a concentration parameter determined by Eq.
(4.2). We set the velocities to the halo velocity plus a dispersion term that
corresponds to the halo virial velocity, as calculated in Eq. (24) in [56]. In
the case of F5, we correct the concentration parameter and the velocity
dispersion as discussed in Section 4.3 employing the fitting formulae in Eq.
(4.3) and Eq. (4.4).

###  5.2 HOD parameters fitting

The 5 parameters of this HOD model are to be determined by requiring that the
galaxy clustering reproduces a reference clustering signal, typically coming
from real observations. In this paper, instead, we use as reference the
clustering of the galaxy catalogues measured in the N-body halo catalogues in
GR with a fiducial set of parameters that were constrained in [52] to fit the
clustering of the CMASS galaxies in BOSS. Since we can measure our reference
signal directly from simulations, we can choose the statistics most convenient
for our purpose. The results of [45] show the importance of RSD in breaking
the degeneracy between galaxy bias and MG parameters. Furthermore, it has been
shown that the multipoles of the power spectrum can break the degeneracy
between different HOD parameters better than the projected correlation
function [57]. In light of these works, we employ the multipoles of the galaxy
power spectrum in redshift space to fit the HOD parameters.

| COLA |  N-body
---|---|---
| F5 | GR | N1 | F5 | GR | N1
$\log_{10}\left(M_{\min}\right)$ | 13.205 | 13.155 | 13.170 | 13.174 | 13.090 | 13.116
$\sigma_{\log_{10}M}$ | 0.579 | 0.589 | 0.590 | 0.578 | 0.596 | 0.600
$\log_{10}\left(M_{0}\right)$ | 13.070 | 13.093 | 12.809 | 12.832 | 13.077 | 12.677
$\log_{10}\left(M_{1}\right)$ | 13.969 | 13.955 | 14.056 | 14.030 | 14.000 | 14.107
$\alpha$ | 1.059 | 1.044 | 1.043 | 1.040 | 1.013 | 1.057
$n_{s}[10^{-4}(\,h\,{\rm Mpc}^{-1})^{3}]$ | 3.210 | 3.218 | 3.217 | 3.231 | 3.226 | 3.230
Table 3: HOD parameters and number density of galaxies obtained by minimising
the objective function in Eq. (5.3) with the simplex search algorithm, except
for the N-body GR column where the HOD parameters are the fiducial parameters
taken from [52].

The galaxy power spectrum in redshift space can be decomposed in angular
moments using the projection

$P_{\ell}(k)=(2\ell+1)\int_{0}^{1}d\mu P(k,\mu)\mathcal{L}_{\ell}(\mu)\,,$
(5.2)

where $\mathcal{L}_{\ell}(\mu)$ is the Legendre polynomial of order $\ell$. In
the Kaiser approximation [58], one can see how the quadrupole and the
hexadecapole moments carry information on the redshift space distortion, which
is key to break the degeneracy between cosmological and modified gravity
parameters. The hexadecapole, however, is strongly affected by noise if
estimated in the cosmological volume covered by our simulations. Thus, the
objective function used to fit the HOD parameters includes only the monopole
and the quadrupole. We use a Gaussian model [59] to estimate the covariance
matrix $\text{Cov}_{\ell,\ell^{\prime}}$ in terms of the volume spanned by the
5 realisations ($\sim 5(h^{-1}\,{\rm Gpc})^{3}$), and the galaxy number
density and bias measured from N-body simulations in GR. We add the tuning of
the number density of galaxies to define the objective function as:

$\chi^{2}=\sum_{\ell,\ell^{\prime}\in 0,2}\sum_{i\in k_{\rm
bins}}\left(P_{\ell,i}-P^{\text{ref}}_{\ell,i}\right)\text{Cov}_{\ell,\ell^{\prime},i}^{-1}\left(P_{\ell^{\prime},i}-P^{\text{ref}}_{\ell^{\prime},i}\right)+W_{n_{s}}\left(\frac{n_{s}-n_{s,\text{ref}}}{n_{s,\text{ref}}}\right)^{2}\,.$
(5.3)

Here $W_{n_{s}}$ is the weight to control the importance of the number density
tuning. We will set $W_{n_{s}}=10^{4}$ to enforce roughly a one per cent
agreement in the number density. The multipoles of the power spectrum and the
covariance matrix are evaluated in 25 linearly distributed bins in the range
$k=[0.05,0.3]\,h\,{\rm Mpc}^{-1}$. We use the simplex search algorithm [60] to
find the set of HOD parameters that minimise the objective function (5.3). The
reference power spectra are computed from N-body GR simulations by populating
galaxies using the fiducial HOD parameters of [52], marked in bold in Table 3.
We average over the 5 realisations of halo catalogues in each model to produce
the multipole moments that are used in the objective function. The best fit
HOD parameters resulting from the simplex search are summarised in Table 3
together with the number density of galaxies $n_{s}$.

Using the best-fit parameters, we produce 5 HOD realisations from each halo
catalogue and we average over them to compute galaxy statistics. We add
redshift space distortions along the three different axes, producing 3
redshift space galaxy catalogues for each real space galaxy catalogue, and use
the average over these when computing redshift space clustering statistics.

(a)  (b)  (c)  (d)

Figure 9: Halo mass (top panels) and velocity distributions (bottom panels) of
central (left panels) and satellite (right panels) galaxies in GR. The
velocity distribution shows the number of galaxies in each velocity bin, which
corresponds to peculiar velocities for centrals and velocities with respect to
the host halo rest frame for satellites. The distributions of galaxies in COLA
before (orange dot-dashed lines) and after (blue dashed lines) the fit
compared with the distributions of N-body in GR (black solid lines) show that
the fit is mostly driven by central galaxies.

To interpret the results of the HOD parameters tuning, we measure in the GR
mocks the galaxy number count as a function of the halo mass and the peculiar
velocity (which for satellites we take as the velocity with respect to the
reference frame of the hosting halo). Figure 9 compares the counts
before101010The HOD parameters before the fit are the same as the default
values for the GR case in $N$-body, highlighted in bold in Table 3. and after
fitting the HOD parameters in COLA and $N$-body in GR. We can see how the fit
is driven by central galaxies, due to their higher abundance and the stronger
signal-to-noise ratio in the monopole. The larger value of $M_{\mathrm{min}}$
in COLA places central galaxies in more massive halos, which are less
abundant, compensating for the excess of halos below
$~{}10^{13}\,h^{-1}\,M_{\odot}$. The smaller values of $M_{1}$ compensates for
the increase of $M_{\mathrm{min}}$, while the larger value of $\alpha$ is
responsible for the heavier right tail of the satellite distributions. Figure
10 also shows the host halo mass and velocity distribution but for all the
gravity models (GR, F5 and N1) after the fit of the HOD parameters, which are
summarised in Table 3. We can see that the differences in the velocity
distribution of central galaxies in different models are well captured in
COLA.

(a)  (b)  (c)  (d)

Figure 10:  Same as Figure 9 but for all the gravity models (GR, F5 and N1 in
blue, orange and green respectively) after fitting the HOD parameters, which
are shown in Table 3.

We compare the monopole and the quadrupole of the galaxy power spectrum in
Figures 11 and 12 respectively. The agreement between COLA and N-body for the
monopole is within the variance for all the gravity models up to $k\sim
0.2\,h\,{\rm Mpc}^{-1}$. The GR quadrupole in COLA agrees with that in N-body
within $5\%$ up to $k\sim 0.17\,h\,{\rm Mpc}^{-1}$. The agreement for the
quadrupole boost factors is within $5\%$ up to $k\sim 0.2\,h\,{\rm Mpc}^{-1}$
for both F5 and N1. Additionally, the quadrupole boost factor in F5 shows a
strong scale dependence, with a small departure from GR on linear scales that
increases on non-linear scales. The quadrupole boost factor in N1, instead,
shows a large departure from GR already on linear scales. This different
behaviour of the quadrupole boost factors is due to the difference in the
linear growth rates of these two models, as well as the different screening
mechanisms. This shows that even after the HOD tuning to obtain a similar
monopole, the difference between models shows up in the quadrupole. Thus RSD
provides a powerful way to distinguish between different MG models [45].

We note that even though we fit multipole moments up to $k_{\rm
max}=0.3\,h\,{\rm Mpc}^{-1}$, the agreement between N-body and COLA degrades
substantially at $k>0.2\,h\,{\rm Mpc}^{-1}$ for quadrupole. This is because
the slight difference in satellite distributions gives a large effect on the
higher multipoles at large $k$ due to the Finger-of-God effect [57], which is
highly sensitive to velocity dispersion of galaxies on small scales. This is
not necessarily an issue for the modelling of cosmological and MG effects in
mock galaxies as the velocity dispersion is usually treated as a nuisance
parameter in cosmological analyses.

(a)  (b)  (c)  (d)

Figure 11: Comparison of the monopole of the galaxy power spectrum in redshift
space between the COLA (in orange) and N-body (in blue) catalogues. The power
spectrum in GR (top panels) and the boost factors in F5 and N1 (bottom panels)
show an agreement within the variance between COLA and N-body.

(a)  (b)  (c)  (d)

Figure 12: Comparison of the quadrupole of galaxy power spectrum in redshift
space between COLA (in orange) and N-body (in blue) catalogues. The power
spectrum in GR (top panels) show a $5\%$ agreement up to $k\sim 0.17\,h\,{\rm
Mpc}^{-1}$. The boost factors in F5 and N1 (bottom panels) show show a better
than $5\%$ agreement up to $k\sim 0.2\,h\,{\rm Mpc}^{-1}$, despite $\sim 10\%$
deviations appear above $k\sim 0.17\,h\,{\rm Mpc}^{-1}$ in GR.

To determine the detectability of MG with summary statistics of galaxy
catalogues, we compute $\chi^{2}$ from Eq. (5.3) in F5 and N1 for both the
COLA and N-body catalogues while using the clustering signal in GR as a
reference. To investigate the role of the different scales, we evaluate the
$\chi^{2}$ for different small scale cut-off $k_{\rm max}$ in the range
between $0.1$ and $0.3\,h\,{\rm Mpc}^{-1}$. The results in Figure 13 highlight
the intrinsic difference between F5 and N1, with the former being
significantly different from GR only on non-linear scales while the latter
showing a strong departure from GR already on linear scales. COLA reproduces
the $\chi^{2}(k_{\rm max})$ behaviour of the N-body very well, demonstrating
that it captures the MG effects accurately.

(a)  (b)

Figure 13: Value of $\chi^{2}$ as defined in Eq. (5.3) in COLA (orange dashed
lines) and N-body (blue solid lines) evaluated using the signal in GR for the
corresponding simulation technique as the target function, as a function of
the small scales cut-off $k_{\rm max}$. The $\chi^{2}$ in F5 (left panel)
shows a strong scale dependence, with small values on linear scales and larger
values on mildly non-linear scales, highlighting the importance of mildly non-
linear scales in constraining MG parameters. The $\chi^{2}$ in N1 (right
panel) shows little scale dependence, reflecting the features of Vainshtein
screening. The results in COLA and N-body are in good agreement.

##  6 Conclusions

We have extended, applied and validated the HOD formalism to COLA simulations
in MG, enabling the fast generation of mock galaxy catalogues that are
required to test gravity on cosmological scales in Stage IV surveys.

Using a full N-body simulation suite as a benchmark, we produced an equivalent
simulation suite with mg-picola, focusing on the following gravity models: GR,
braneworld nDGP model (N1) and $n=1$ Hu-Sawicki $f(R)$ gravity (F5). Looking
at the power spectra of the density and velocity divergence fields, we
confirmed the good accuracy of COLA in reproducing the density field
(agreement within the variance up to $k\sim 1\,h\,{\rm Mpc}^{-1}$) and have
given the first measure of its accuracy in the velocity field (within $4\%$ up
to $k\sim 0.4\,h\,{\rm Mpc}^{-1}$).

To apply the HOD formalism to COLA simulations, we first need to employ a halo
finding algorithm to find dark matter halos. We found that the rockstar halo-
finder, which uses the six-dimensional phase-space information to find halos,
did not perform well with our COLA simulations with the default setting,
resulting in a disagreement between COLA and N-body halo statistics ($\sim
25\%$ in the hmf and $10\%$ in the power spectrum at $k\sim 0.4\,h\,{\rm
Mpc}^{-1}$). On the other hand, finding halos with a simple FoF finder and
converting the halo mass with the results of [53] and [54], we have produced
halo catalogues in COLA and N-body which are in good agreement, as shown by
the hmf (better than $10\%$ agreement) and the halo power spectrum (agreement
within the variance).

Since COLA simulations do not resolve the internal halo structure, we have
assumed a NFW profile and the concentration-mass relation of [54] in GR. To
take into account the effect of the fifth force on the internal structure of
halos, we have measured the boost factors of the concentration and velocity
dispersion mass relations in MG with rockstar halo catalogues in N-body. We
have found that no significant difference is present in the N1 boost factors,
which allows us to use the standard NFW profile in N1. In F5, instead, there
is a clear transition from screened to un-screened halos both in the
concentration and the velocity dispersion boost factors. This is in line with
the findings of [17, 18]. Using simple fitting formulae that we calibrated on
these results, we have been able to assign realistic concentration and
velocity dispersion to halos in F5 depending on their mass.

With the HOD model proposed in [15], we have fitted the HOD parameters of the
different models to match the galaxy clustering of N-body catalogues in GR.
The fit performed with a simplex search algorithm is mostly driven by the
central galaxies, whose final distributions in velocity are different in each
model but similar between COLA and N-body. Despite that the objective function
for the fit is designed to minimise the difference of the monopole and the
quadrupole of the galaxy power spectrum with the target signals, we have shown
that MGs leave a characteristic imprint in RSD statistics, especially in the
quadrupole. To estimate the constraining power of galaxy statistics, we have
computed the $\chi^{2}$ of the difference between MG and GR multipoles. We
confirmed that non-linear scales are important in constraining F5, while N1
can be constrained on large scales due to the scale-independent modification
of the growth rate. We also found a good agreement between COLA and N-body
galaxy statistics.

The accuracy of COLA galaxy catalogues may be further improved by tuning the
COLA halo catalogues to reproduce an analytic halo mass function. This can be
achieved by introducing a mass-dependent conversion factor when converting
$M_{\rm FoF}$ to $M_{\rm 200c}$. However, this requires an accurate
measurement of the mass function from N-body simulations in MG, which
undermines the predictability of COLA, one of the major advantages that COLA
has over other approximate methods. In this paper, we focused on the power
spectrum in redshift space. It is an open question whether mg-picola can also
produce higher-order statistics in MG such as the bispectrum (because of the
approximations in the screening, although [40] showed that COLA can reproduce
the bispectrum in GR) as well as other non-linear structures such as voids.
These studies are left for future work.

Using the results of this work, we plan to produce a large number of mock
galaxy catalogues in MG with mg-picola. This will allow us to study the
degeneracy between MG and cosmological parameters by looking at different
probes and formulate the best estimators to use in Stage IV experiments to
constrain gravity on cosmological scales.

## Acknowledgments

KK and AI are supported by the European Research Council under the European
Union’s Horizon 2020 programme (grant agreement No.646702 “CosTesGrav"). KK is
also supported by the UK STFC grant ST/S000550/1. BSW is supported by the
Royal Society grant number RGF\EA\181023. BL is supported by the European
Research Council under the European Union’s Horizon 2020 programme (grant
agreement No. 716532 “PUNCA"), and UK STFC Consolidated Grants ST/L00075X/1,
ST/P000451/1. Numerical computations were done on the Sciama High Performance
Compute (HPC) cluster which is supported by the ICG, SEPNet and the University
of Portsmouth. Some of the simulations described in this work used the DiRAC
Data Centric system at Durham University, operated by the Institute for
Computational Cosmology on behalf of the STFC DiRAC HPC Facility
(www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure
capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1,
STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part
of the National E-Infrastructure. Part of the simulations were run by Wojciech
Hellwing at the Interdisciplinary Center for Mathematical and Computational
Modelling (ICM) at University of Warsaw, under Grants No. GA67-17 and No.
GA65-30.

## References

  * [1] K. Koyama, _Cosmological Tests of Modified Gravity_ , _Rept. Prog. Phys._ 79 (2016) 046902 [1504.04623].
  * [2] S. Alam et al., _Testing the theory of gravity with DESI: estimators, predictions and simulation requirements_ , 2011.05771.
  * [3] H.A. Winther et al., _Modified Gravity N-body Code Comparison Project_ , _Mon. Not. Roy. Astron. Soc._ 454 (2015) 4208 [1506.06384].
  * [4] P. Monaco, _Approximated methods for the generation of dark matter halo catalogs in the age of precision cosmology_ , _ArXiv e-prints_ (2016) [1605.07752].
  * [5] C.-H. Chuang et al., _nIFTy Cosmology: Galaxy/halo mock catalogue comparison project on clustering statistics_ , _Mon. Not. Roy. Astron. Soc._ 452 (2015) 686 [1412.7729].
  * [6] S. Tassev, M. Zaldarriaga and D. Eisenstein, _Solving Large Scale Structure in Ten Easy Steps with COLA_ , _JCAP_ 06 (2013) 036 [1301.0322].
  * [7] J. Koda, C. Blake, F. Beutler, E. Kazin and F. Marin, _Fast and accurate mock catalogue generation for low-mass galaxies_ , _Mon. Not. Roy. Astron. Soc._ 459 (2016) 2118 [1507.05329].
  * [8] A. Izard, M. Crocce and P. Fosalba, _ICE-COLA: Towards fast and accurate synthetic galaxy catalogues optimizing a quasi $N$-body method_, _Mon. Not. Roy. Astron. Soc._ 459 (2016) 2327 [1509.04685].
  * [9] C. Howlett, M. Manera and W.J. Percival, _L-PICOLA: A parallel code for fast dark matter simulation_ , _Astron. Comput._ 12 (2015) 109 [1506.03737].
  * [10] G. Valogiannis and R. Bean, _Efficient simulations of large scale structure in modified gravity cosmologies with comoving Lagrangian acceleration_ , _Phys. Rev. D_ 95 (2017) 103515 [1612.06469].
  * [11] H.A. Winther, K. Koyama, M. Manera, B.S. Wright and G.-B. Zhao, _COLA with scale-dependent growth: applications to screened modified gravity models_ , _JCAP_ 08 (2017) 006 [1703.00879].
  * [12] B.S. Wright, H.A. Winther and K. Koyama, _COLA with massive neutrinos_ , _JCAP_ 10 (2017) 054 [1705.08165].
  * [13] C. Moretti, S. Mozzon, P. Monaco, E. Munari and M. Baldi, _Fast numerical method to generate halo catalogues in modified gravity (part I): second-order Lagrangian perturbation theory_ , _Mon. Not. Roy. Astron. Soc._ 493 (2020) 1153 [1909.06282].
  * [14] A. Knebe et al., _Haloes gone MAD: The Halo-Finder Comparison Project_ , _Mon. Not. Roy. Astron. Soc._ 415 (2011) 2293 [1104.0949].
  * [15] Z. Zheng, A.L. Coil and I. Zehavi, _Galaxy Evolution from Halo Occupation Distribution Modeling of DEEP2 and SDSS Galaxy Clustering_ , _Astrophys. J._ 667 (2007) 760 [astro-ph/0703457].
  * [16] J.F. Navarro, C.S. Frenk and S.D.M. White, _The Structure of cold dark matter halos_ , _Astrophys. J._ 462 (1996) 563 [astro-ph/9508025].
  * [17] M.A. Mitchell, J.-h. He, C. Arnold and B. Li, _A general framework to test gravity using galaxy clusters – I. Modelling the dynamical mass of haloes in $f(R)$ gravity_, _Mon. Not. Roy. Astron. Soc._ 477 (2018) 1133 [1802.02165].
  * [18] M.A. Mitchell, C. Arnold, J.-h. He and B. Li, _A general framework to test gravity using galaxy clusters II: A universal model for the halo concentration in $f(R)$ gravity_, _Mon. Not. Roy. Astron. Soc._ 487 (2019) 1410 [1901.06392].
  * [19] B. Li, G.-B. Zhao, R. Teyssier and K. Koyama, _ECOSMOG: An Efficient Code for Simulating Modified Gravity_ , _JCAP_ 01 (2012) 051 [1110.1379].
  * [20] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [21] A.G. Riess, S. Casertano, W. Yuan, J.B. Bowers, L. Macri, J.C. Zinn et al., _Cosmic Distances Calibrated to 1% Precision with Gaia EDR3 Parallaxes and Hubble Space Telescope Photometry of 75 Milky Way Cepheids Confirm Tension with $\Lambda$CDM_, _Astrophys. J. Lett._ 908 (2021) L6 [2012.08534].
  * [22] T. Tröster et al., _Cosmology from large-scale structure: Constraining $\Lambda$CDM with BOSS_, _Astron. Astrophys._ 633 (2020) L10 [1909.11006].
  * [23] M. Raveri, _Reconstructing Gravity on Cosmological Scales_ , _Phys. Rev. D_ 101 (2020) 083524 [1902.01366].
  * [24] A. Joyce, B. Jain, J. Khoury and M. Trodden, _Beyond the Cosmological Standard Model_ , _Phys. Rept._ 568 (2015) 1 [1407.0059].
  * [25] C.M. Will, _The Confrontation between General Relativity and Experiment_ , _Living Rev. Rel._ 17 (2014) 4 [1403.7377].
  * [26] J. Khoury and A. Weltman, _Chameleon fields: Awaiting surprises for tests of gravity in space_ , _Phys. Rev. Lett._ 93 (2004) 171104 [astro-ph/0309300].
  * [27] J. Khoury and A. Weltman, _Chameleon cosmology_ , _Phys. Rev. D_ 69 (2004) 044026 [astro-ph/0309411].
  * [28] A.I. Vainshtein, _To the problem of nonvanishing gravitation mass_ , _Phys. Lett. B_ 39 (1972) 393.
  * [29] A.A. Starobinsky, _A New Type of Isotropic Cosmological Models Without Singularity_ , _Phys. Lett. B_ 91 (1980) 99.
  * [30] W. Hu and I. Sawicki, _Models of f(R) Cosmic Acceleration that Evade Solar-System Tests_ , _Phys. Rev. D_ 76 (2007) 064004 [0705.1158].
  * [31] L. Pogosian and A. Silvestri, _The pattern of growth in viable f(R) cosmologies_ , _Phys. Rev. D_ 77 (2008) 023503 [0709.0296].
  * [32] H.A. Winther and P.G. Ferreira, _Fast route to nonlinear clustering statistics in modified gravity theories_ , _Phys. Rev. Lett._ 91 (2015) 123507 [1403.6492].
  * [33] G.R. Dvali, G. Gabadadze and M. Porrati, _4-D gravity on a brane in 5-D Minkowski space_ , _Phys. Lett. B_ 485 (2000) 208 [hep-th/0005016].
  * [34] F. Schmidt, _Self-Consistent Cosmological Simulations of DGP Braneworld Gravity_ , _Phys. Rev. D_ 80 (2009) 043001 [0905.0858].
  * [35] S. Bag, S.S. Mishra and V. Sahni, _Emulating a $\Lambda$CDM-like expansion on the phantom brane_, _Phys. Rev. D_ 97 (2018) 123537 [1807.00684].
  * [36] K. Koyama and R. Maartens, _Structure formation in the dgp cosmological model_ , _JCAP_ 01 (2006) 016 [astro-ph/0511634].
  * [37] F. Schmidt, W. Hu and M. Lima, _Spherical Collapse and the Halo Model in Braneworld Gravity_ , _Phys. Rev. D_ 81 (2010) 063005 [0911.5178].
  * [38] M. Lippich, A.G. Sánchez, M. Colavincenzo, E. Sefusatti, P. Monaco, L. Blot et al., _Comparing approximate methods for mock catalogues and covariance matrices – i. correlation function_ , _Monthly Notices of the Royal Astronomical Society_ 482 (2018) 1786–1806.
  * [39] L. Blot, M. Crocce, E. Sefusatti, M. Lippich, A.G. Sánchez, M. Colavincenzo et al., _Comparing approximate methods for mock catalogues and covariance matrices ii: power spectrum multipoles_ , _Monthly Notices of the Royal Astronomical Society_ 485 (2019) 2806–2824.
  * [40] M. Colavincenzo, E. Sefusatti, P. Monaco, L. Blot, M. Crocce, M. Lippich et al., _Comparing approximate methods for mock catalogues and covariance matrices – iii: bispectrum_ , _Monthly Notices of the Royal Astronomical Society_ 482 (2018) 4883–4905.
  * [41] A. Aviles and J.L. Cervantes-Cota, _Lagrangian perturbation theory for modified gravity_ , _Phys. Rev. D_ 96 (2017) 123526 [1705.10719].
  * [42] M. Cautun, E. Paillas, Y.-C. Cai, S. Bose, J. Armijo, B. Li et al., _The Santiago–Harvard–Edinburgh–Durham void comparison – I. SHEDding light on chameleon gravity tests_ , _Mon. Not. Roy. Astron. Soc._ 476 (2018) 3195 [1710.01730].
  * [43] R. Teyssier, _Cosmological hydrodynamics with adaptive mesh refinement: a new high resolution code called ramses_ , _Astron. Astrophys._ 385 (2002) 337 [astro-ph/0111367].
  * [44] C. Hernández-Aguayo, C.M. Baugh and B. Li, _Marked clustering statistics in $f(R)$ gravity cosmologies_, _Mon. Not. Roy. Astron. Soc._ 479 (2018) 4824 [1801.08880].
  * [45] C. Hernández-Aguayo, J. Hou, B. Li, C.M. Baugh and A.G. Sánchez, _Large-scale redshift space distortions in modified gravity theories_ , _Mon. Not. Roy. Astron. Soc._ 485 (2019) 2194 [1811.09197].
  * [46] A. Barreira, A.G. Sánchez and F. Schmidt, _Validating estimates of the growth rate of structure with modified gravity simulations_ , _Phys. Rev. D_ 94 (2016) 084022 [1605.03965].
  * [47] N.C. Devi, A. Rodríguez-Puebla, O. Valenzuela, V. Avila-Reese, C. Hernández-Aguayo and B. Li, _The galaxy–halo connection in modified gravity cosmologies: environment dependence of galaxy luminosity function_ , _Mon. Not. Roy. Astron. Soc._ 488 (2019) 782 [1901.02121].
  * [48] M.C. Cautun and R. van de Weygaert, _The DTFE public software: The Delaunay Tessellation Field Estimator code_ , 1105.0370.
  * [49] P.S. Behroozi, R.H. Wechsler and H.-Y. Wu, _The ROCKSTAR Phase-space Temporal Halo Finder and the Velocity Offsets of Cluster Cores_ , _Astrophys. J._ 762 (2013) 109 [1110.4372].
  * [50] M. Davis, G. Efstathiou, C.S. Frenk and S.D.M. White, _The evolution of large-scale structure in a universe dominated by cold dark matter_ , _Astrophys. J._ 292 (1985) 371.
  * [51] C. Howlett, A. Ross, L. Samushia, W. 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$_, _Mon. Not. Roy. Astron. Soc._ 449 (2015) 848 [1409.3238].
  * [52] M. Manera et al., _The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: a large sample of mock galaxy catalogues_ , _Mon. Not. Roy. Astron. Soc._ 428 (2012) 1036 [1203.6609].
  * [53] Z. Lukic, D. Reed, S. Habib and K. Heitmann, _The Structure of Halos: Implications for Group and Cluster Cosmology_ , _Astrophys. J._ 692 (2009) 217 [0803.3624].
  * [54] A.A. Dutton and A.V. Macciò, _Cold dark matter haloes in the Planck era: evolution of structural parameters for Einasto and NFW profiles_ , _Mon. Not. Roy. Astron. Soc._ 441 (2014) 3359 [1402.7073].
  * [55] A.P. Hearin et al., _Forward Modeling of Large-Scale Structure: An open-source approach with Halotools_ , _Astron. J._ 154 (2017) 190 [1606.04106].
  * [56] S. More, F.C.v.d. Bosch and M. Cacciato, _Satellite Kinematics I: A New Method to Constrain the Halo Mass-Luminosity Relation of Central Galaxies_ , _Mon. Not. Roy. Astron. Soc._ 392 (2009) 917 [0807.4529].
  * [57] C. Hikage, _Constraining Halo Occupation Distribution and Cosmic Growth Rate using Multipole Power Spectrum_ , _Mon. Not. Roy. Astron. Soc._ 441 (2014) L21 [1401.1246].
  * [58] N. Kaiser, _Clustering in real space and in redshift space_ , _Mon. Not. Roy. Astron. Soc._ 227 (1987) 1.
  * [59] A. Taruya, T. Nishimichi and S. Saito, _Baryon Acoustic Oscillations in 2D: Modeling Redshift-space Power Spectrum from Perturbation Theory_ , _Phys. Rev. D_ 82 (2010) 063522 [1006.0699].
  * [60] J.A. Nelder and R. Mead, _A Simplex Method for Function Minimization_ , _Comput. J._ 7 (1965) 308.