v#1

# Fast production of cosmological emulators in modified gravity: the matter
power spectrum

Bartolomeo Fiorini,11footnotetext: Corresponding author. Kazuya Koyama  Tessa
Baker

###### Abstract

We test the convergence of fast simulations based on the COmoving Lagrangian
Acceleration (COLA) method for predictions of the matter power spectrum,
specialising our analysis in the redshift range $1\leq z\leq 1.65$, relevant
to high-redshift spectroscopic galaxy surveys. We then focus on the
enhancement of the matter power spectrum in modified gravity (MG), the boost
factor, using the Dvali-Gabadadze-Porrati (DGP) theory as a test case but
developing a general approach that can be applied to other MG theories. After
identifying the minimal simulation requirements for accurate DGP boost
factors, we design and produce a COLA simulation suite that we use to train a
neural network emulator for the DGP boost factor. Using MG-AREPO simulations
as a reference, we estimate the emulator accuracy to be of $\sim 3\%$ up to
$k=5\,h\,{\rm Mpc}^{-1}$ at $0\leq z\leq 2$. We make the emulator publicly
available at: https://github.com/BartolomeoF/nDGPemu.

##  1 Introduction

With the successful launch of the Euclid222https://www.euclid-ec.org/
satellite and the DESI333https://www.desi.lbl.gov collaboration releasing
their $1\%$ survey data [1], cosmology is living a very exciting time. The
extraordinarily large amount of galaxies that these surveys will observe is
fuelling the efforts of the scientific community aimed at accurately modelling
the non-linear regime of structure formation to obtain precise cosmological
constraints.

One of the key objectives of next-generation surveys is to shed light on the
nature of gravity by constraining alternative gravity models. The simplest
modified gravity (MG) theories can be described by means of a fifth force
mediated by an additional scalar field. Given the stringent constraints coming
from solar system tests of gravity, MG theories of interest to cosmology
incorporate screening mechanisms which suppress the fifth force in specific
environments. Hence, the large-scale structure (LSS) of the Universe provides
a unique laboratory to test the validity of GR and explore alternative gravity
models. By comparing LSS probes with theoretical predictions, we can place
constraints on the deviations from GR and potentially uncover new physics.

Yet, accessing the non-linear regime of structure formation presents
significant challenges. Analytical techniques, such as perturbation theory,
break down in this regime, and numerical simulations become essential.
Cosmological simulations play a vital role in studying the formation and
evolution of structures, enabling us to investigate the complex interplay
between gravity, dark matter (DM), and baryonic physics. To simulate structure
formation on the smallest scales, DM-only N-body simulations are not enough
and it is necessary to run hydrodynamical simulations which incorporate the
effects of baryonic matter. However, despite the advances in simulating
baryonic physics, hydrodynamical simulations show important discrepancies in
their small-scales predictions [2], so weak lensing data below baryonic-
feedback scales is discarded in cosmological analysis (e.g. see [3, 4]).

Furthermore, cosmological simulations are computationally demanding, requiring
substantial computational resources and time. To overcome these computational
challenges, the development of efficient methods to speed up simulations has
been a topic of active research. One such method is the COmoving Lagrangian
Acceleration (COLA) method [5], which combines second-order Lagrangian
Perturbation Theory (2LPT) with the Particle Mesh N-body technique to
accurately reproduce the non-linear dynamics of structure formation while
significantly reducing computational costs. In particular, COLA exploits the
fact that long-range forces perceived by the N-body particles in a frame
comoving with the 2LPT solution are weaker than the ones perceived in the
Eulerian frame. This allows COLA simulations to use a coarser time-resolution
while providing an accurate description of the matter distribution on large
scales [5, 6].

The potential of the COLA method was already exploited in [7] to produce mock
galaxy catalogues for the analysis of the WiggleZ survey [8], where it was
shown that COLA simulations can be used to create synthetic galaxy catalogues
with comparable accuracy to the full N-body code Gadget [9] on scales $k\leq
0.2\,h\,{\rm Mpc}^{-1}$. Since then, the COLA method has seen several
extensions:

  * •

to create more accurate halo catalogues [6],

  * •

to produce matter distribution and mock halo catalogues on the light-cone [10,
11],

  * •

to extend to MG theories the evolution of the matter distribution [12, 13, 14]
and the creation of mock galaxy catalogues [15],

  * •

to include the effects of massive neutrinos [16].

Of particular interest to this work is the extension to MG theories presented
in [13] and as implemented in the publicly available library
FML444https://fml.wintherscoming.no. This relies on the extension to MG of the
2LPT solutions and to an approximation of the screening mechanisms based on
the analytic solution in spherical symmetry [17] which is provided for the
Vainshtein and chameleon screening mechanisms. Recently, an alternative
technique to approximate the Vainshtein screening [18] has been revised and
implemented in a custom version of FML-COLA [19].

Despite the speed-up brought by fast simulation methods like COLA, the
computational cost of simulations is still too high to allow a fine sampling
of the vast theory parameter space that is needed in Monte Carlo searches to
infer precise parameter constraints. One possible solution to this problem is
to approximate the predictions of cosmological simulations for specific
summary statistics making use of machine learning techniques [20, 21] or
fitting functions (see for instance [22]). The resulting approximating
functions are normally referred to as emulators in the field and their key
feature is that they provide theoretical predictions at extremely small
computational cost. Accurate emulators for the matter power spectrum, one of
the most elementary summary statistics that can be predicted thanks to
simulations, are available for standard cosmologies [23, 24]. Matter power
spectrum emulators have been developed also in MG gravity theories such as
[25, 26, 27, 28, 29] in Hu-Sawicki $f(R)$ model [30], and [31] in the Dvali-
Gabadadze-Porrati (DGP) model [32]. An interesting alternative approach based
on the halo model was also developed to predict the matter power spectrum in
MG theories [33].

Due to the many subtleties that characterise N-body simulations and to the
approximate nature of emulators, these tools need to be validated. This can be
done for example by means of comparisons between different approaches [34, 23,
24]. In the case of N-body simulations, it is also important to check that the
agreement with the other codes is not just due to a coincidence and that the
predictions converge to the true answer when increasing resolution.

In this work we push the COLA method to its limits to study the convergence of
COLA predictions for the matter power spectrum and its enhancement in MG, the
boost factor. We employ the results of this analysis to design and produce an
emulator for the matter power spectrum in DGP gravity that we publicly release
together with this paper. While in this work we focus on DGP gravity, the
framework that we develop is quite general and can be applied to other MG
theories.

The remainder of this paper is organised as follows. In section 2 we provide
an in-depth analysis of the convergence and validation of COLA simulations for
prediction of the DM power spectrum in GR. In section 3 we extend the
convergence analysis to the nDGP boost factor for the DM power spectrum and
present an emulator for this signal trained on COLA simulations. Finally, we
come to our conclusions in section 4.

##  2 Convergence tests

Cosmological simulations are subject to errors due to the finite resolution
characterising their numerical implementation. In the case of N-body
simulations, time-resolution, force-resolution and mass-resolution are
important quantities which control the simulation accuracy [35, 36, 23]. To
validate simulation results it is necessary to compare the output of
simulations against well-established results like analytic predictions or
previous simulations already validated (e.g., in [15, 37] we validated COLA
simulations by comparing them with ELEPHANT simulations). Furthermore, the
self-consistency of a simulation technique can be established by performing
relative convergence tests, where simulations with different resolutions are
compared to understand the reliability of specific simulation settings.

The COLA method exploits the 2LPT and the Particle-Mesh (PM) method to evolve
the DM distribution. The use of 2LPT allows us to reduce the time resolution
of simulations without losing accuracy on large scales [6]. Conversely,
increasing the time resolution makes the use of 2LPT less important and with
sufficient time resolution, the simulations evolve like PM-only simulations.
In this sense, we will refer to COLA simulations with high time resolution as
PM simulations in the following sections. This equivalence is important since
the accuracy of PM simulations (and therefore of COLA simulations in
$\Lambda$CDM) is only limited by resolution and they converge to full N-body
simulations results with increasing resolution [38]. However, when comparing
COLA with PM-only algorithms, it is not fair to compare their computational
cost using the same number of time-steps as done in [35]. In fact, COLA does
not necessarily need as many time-steps as PM-only simulations to achieve
comparable accuracy. This is evident in the low time-resolution case, where
the 2LPT guarantees that COLA provides accurate results on large scales while
PM-only simulations already show significant deviations [6]. Even in the high
time-resolution case, COLA can achieve the same accuracy as PM-only codes with
fewer time-steps. Indeed, the latter needs a finer time-stepping at high
redshift to deal with the rapidly changing growth factor [35] while the 2LPT
guarantees good accuracy at high redshift in COLA even with larger time-steps
[5, 6]. For instance, the PM code GLAM needs 147 time-steps to have $\Delta
a/a\approx 0.014$ at $z=0$ [35] while COLA needs only $\approx 70$ time-steps
for the same low-redshift time resolution.

The focus of this section is to validate COLA simulations for the predictions
of the matter power spectrum. To do so, we first introduce the cosmological
emulators in subsection 2.1, then we compare COLA results with full N-body
simulations and cosmological emulators in subsection 2.2, and finally, we
study the relative-convergence of COLA simulations for the matter power
spectrum in subsection 2.3.

###  2.1 Cosmological emulators

Cosmological emulators interpolate the results of computationally expensive
cosmological simulations in high-dimensional parameter spaces. Two notable
such emulators are the Bacco emulator [24] and the Euclid Emulator 2 (EE2)
[23]. While these emulators differ on many aspects (e.g. the interpolation
techniques employed, the parameter space spanned and the N-body techniques
used for the training set) they have been shown to be in agreement with each
other within the claimed accuracy [23, 39] for several test cosmologies.

Here in particular we are interested in the agreement between these two
emulators for a specific cosmology, the one that was used to run the N-body
simulation suite that we will discuss in subsection 2.2. In figure 1 we show
the ratio of the matter power spectra estimated with Bacco and EE2 at
redshifts $0$, $0.5$ and $1$ for the above-mentioned cosmology (see eq 2.1).
The agreement is within $2\%$ at all the scales considered. At $z=1$, the
agreement is well within $1\%$. This is consistent with the accuracy claimed
by the two emulators and gives additional reliability to the emulators’
results, in particular at redshift $z=1$, that we use in the next subsection
to prove the accuracy of COLA simulations for predictions of the matter power
spectrum.

Figure 1:  Ratio of Bacco and EE2 power spectra for the cosmology defined in
eq (2.1).

###  2.2 Comparison of COLA with emulators and N-body simulations

As reference results to compare COLA with, we make use of the power spectra
obtained from a suite of DM-only simulations in GR and nDGP gravity555The nDGP
results (for $H_{0}r_{c}$ values of 0.5, 1, 2, 5) are used in section 3. ran
with MG-AREPO [40, 41] which uses the tree-particle-mesh technique for
gravitational interactions and a moving mesh to solve the hydrodynamic
equations (see references [42, 9] for more on the N-body techniques at the
base of this code). The simulations were started at redshift $z=127$ using the
cosmological parameters

$\begin{array}[]{ccc}\Omega_{m,0}=0.3089\,,&\Omega_{\Lambda,0}=0.6911\,,&\Omega_{b,0}=0.0486\,,\\\
n_{s}=0.9667\,,&\sigma_{8}=0.8159\,,&h=0.6774\,,\end{array}$ (2.1)

and evolved up to redshift $z=0$ using $1024^{3}$ particles in a $(1{\rm
Gpc}/h)^{3}$ box [43].

To assess the accuracy of COLA simulations we produce a suite of simulations
in a $L=512$Mpc$/h$ box, with the same cosmology of the Arepo simulations and
with different mass resolutions ($N_{\rm part}^{1/3}=512,1024,2048$), time
resolutions ($N_{\rm step}=25\times[1,2,3]=25,50,75$) and with force
resolutions depending on the mass resolutions ($N_{\rm
mesh}^{1/3}=[1,2,3,4]\times N_{\rm part}^{1/3}$ for $N_{\rm
part}^{1/3}=512,1024$ and $N_{\rm mesh}^{1/3}=[1,1.5,2]\times N_{\rm
part}^{1/3}$ for $N_{\rm part}^{1/3}=2048$). We start the simulations at
redshift $z=19$ and stop them at redshift $z=1$. The choice of 25 time-steps
as the low time-resolution case is motivated by the time-step size of typical
COLA configurations ($\Delta a\approx 0.02$). In the high time-resolution
case, where the late-time evolution is PM-like, 75 time-steps produce a time-
step size of $\Delta a\approx 0.006$. This is the same step size of the
reference simulations used in the convergence tests of GLAM at $z<3$ [35].

In the next subsection, we will use the highest resolution PM run ($N_{\rm
part}^{1/3}=2048$, $N_{\rm mesh}^{1/3}=4096$, $N_{\rm step}=75$) as the
reference to assess the accuracy of all the other runs. Before that, we
validate the highest resolution run by comparing its matter power spectrum
with the ones from the Bacco emulator and from the Arepo simulations. In
Figure 2

Figure 2:  Comparison of PM simulations matter power spectra with Bacco and
Arepo power spectra. The high mass-resolution PM run agrees with Bacco up to
$k\sim 2h/$Mpc within $1\%$. Arepo does not agree with Bacco but is in better
agreement with the low mass-resolution PM run which has approximately the same
mass resolution as the Arepo simulation.

we plot the ratio of the matter power spectra of PM and Arepo simulations with
the Bacco power spectrum for the cosmology of eq. (2.1), in GR at redshift 1.
The first thing that comes to the eye, is that the Arepo power spectrum (green
line) departs from Bacco already at $k\sim 0.6h$/Mpc, and at $k\sim 1h/$Mpc
the discrepancy is already at $4\%$. We attribute this to the low mass
resolutions of these Arepo simulations ($M_{\rm part}\approx 8\cdot
10^{10}\,h^{-1}\,M_{\odot}$) [23], and to motivate this claim we add to the
comparison the PM simulation for $N_{\rm part}^{1/3}=512$, $N_{\rm
mesh}^{1/3}=2048$, $N_{\rm step}=75$ (blue line) which has a similar mass
resolution of Arepo simulations. This is in good agreement with Arepo up to
$k\sim 2h/$Mpc. The blue line instead shows that the higher resolution PM run
is in sub-percent agreement with Bacco up to $k\sim 2\,h\,{\rm Mpc}^{-1}$.

###  2.3 Relative convergence of COLA simulations

With all the simulations of the convergence suite, we compare the accuracy of
all the different settings using the high-resolution run as the reference. To
summarise the results, we show in different figures different mass
resolutions, figure 3 for $N_{\rm part}^{1/3}=512$, figure 4 for $N_{\rm
part}^{1/3}=1024$, figure 5 for $N_{\rm part}^{1/3}=2048$, and split each
figure into 2 panels for two redshift values $z=1.65,1$666 This choice of
redshift values is motivated by the expected redshift localisation of
H$\alpha$ emitters in the Euclid spectroscopic survey [44]. respectively on
the left and on the right. In each panel, we show the ratio of the power
spectra in the COLA run with that of the high-resolution run for the various
force and time resolutions described in the legends.

Figure 3:  Ratio of the matter power spectrum of simulations with 512
particles for the force and time resolutions shown in the legend with that of
the high-resolution run. Figure 4:  Ratio of the matter power spectrum of
simulations with 1024 particles for the force and time resolutions shown in
the legend with that of the high-resolution run. Figure 5:  Ratio of the
matter power spectrum of simulations with 2048 particles for the force and
time resolutions shown in the legend with that of the high-resolution run.

For the $N_{\rm part}^{1/3}=512$ case, we notice that the results are almost
insensitive to the time resolution and do not converge towards the high-
resolution result. This is attributed to the fact that the low mass resolution
has a strong effect on the scales considered in the comparison, especially for
higher redshifts. In addition to this, we notice that the convergence is non-
monotonic: the fact that $N_{\rm mesh}^{1/3}=1024$ results have significantly
larger power than the ones of $N_{\rm mesh}^{1/3}=1536$ and $N_{\rm
mesh}^{1/3}=2048$ is somehow surprising. From additional tests that we have
run, we could identify this behaviour as being due to discreteness effects
happening at high redshifts when the particles’ displacements from the Fourier
grid are small in combination with the continuous Fourier space kernels used
for the force computation in the PM scheme777In the FML library it is now
possible to choose the technique adopted by the PM solver for the force
computation. In particular, it is possible to choose the same finite
differencing scheme used in codes like FastPM [45] and GLAM [35] which seems
to avoid these problems of non-monotonic convergence in low mass-resolution PM
simulations.. While this effect is moderate when starting the simulations at
$z=19$, it becomes more important when starting the simulations at higher
redshift.

For the $N_{\rm part}^{1/3}=1024$ and $2048$ cases, the convergence towards
the high-resolution PM result happens only if the force and time resolutions
are increased accordingly. This condition is particularly evident in the
simulations with $N_{\rm part}^{1/3}=1024$, where increasing the number of
time steps from 25 to 75 keeping the fixed force resolution of $N_{\rm
mesh}^{1/3}=1024$ does not improve the convergence toward the high-resolution
run. Similarly, increasing the force resolution above $N_{\rm
mesh}^{1/3}=2048$ but using only 25 time-steps does not improve the agreement
with the highest-resolution run. We notice that, in the $N_{\rm
part}^{1/3}=1024$ case, the results for high force and time resolution appear
to be more converged at low redshift than at high redshift: this is due to the
fact that discreteness effects are more important at early times while force
resolution becomes more important at late times.

To assess the convergence of the highest-resolution run used as a reference in
the previous figures, we compare in figure 6 the power spectra of lower-
resolution simulations with the highest-resolution one, as done previously,
but now increasing mass, force and time resolution simultaneously. The
agreement between the $N_{\rm part}^{1/3}=1024$, $N_{\rm mesh}^{1/3}=3072$,
$N_{\rm step}=50$ run and the highest-resolution run indicates that the latter
is converged at sub-percent level up to $k\sim 1\,h\,{\rm Mpc}^{-1}$. This is
in line with the result of the comparison with the Bacco emulator.

Figure 6:  Convergence of matter power spectra towards the high-resolution
result as mass-, force- and time-resolution are simultaneously increased (as
indicated by the colour legend). Different line styles are used to distinguish
results at different redshift values as indicated in the line-style legend.

##  3 Modified gravity boost factor

Another interesting application of COLA simulations is to predict the matter
power spectrum in MG [25]. As non-linear matter power spectrum emulators are
already available, it is possible to obtain the matter power spectrum in MG by
multiplying the matter power spectrum in GR with the MG boost factor

$B_{\rm MG}(k,z)\equiv\frac{P_{\rm MG}(k,z)}{P_{\rm GR}(k,z)}\,.$ (3.1)

In general, the MG boost factor depends on MG parameters as well as on
cosmological parameters and redshift. In the following, we focus on the DGP
model and use COLA simulations to train an emulator for the matter power
spectrum in this gravity theory. Nonetheless, the method that we develop is
not specific to the DGP model and can be used to produce emulators for other
MG theories.

The DGP gravity is a braneworld model in 5-dimensions, which has been
thoroughly studied by the modified gravity community due to its interesting
cosmological features [46, 47, 48]. This theory has 2 branches of solutions
for the Hubble expansion on the brane [46]: a self-accelerating branch with
ghost instabilities on a de Sitter background [47, 49], and a normal branch
which requires fine-tuning in the form of an effective dark energy to recover
the observed $\Lambda$CDM-like expansion history [50, 51]. Nonetheless, the
normal branch of the DGP theory (nDGP) also presents a non-trivial fifth force
dynamics which can be described by the Vainshtein mechanism [52], so this
theory was often used as a test-case representative for the class of
Vainshtein-screened theories. In particular, several N-body codes were
developed to simulate non-linear structure formation in nDGP theory (e.g. DGPM
[50], ECOSMOG [53], MG-AREPO [54], MG-GLAM [38]) and two of them, DGPM and
ECOSMOG, were showed to give consistent prediction for the matter power
spectrum boost-factor in [55].

To develop the nDGP emulator, we start by studying the simulation requirements
to accurately predict the nDGP boost factor with COLA and investigate the
sensitivity of the boost factor to the theory parameters in subsection 3.1.
Then we present the simulations suite and the data processing used to create
the emulator’s data sets in subsection 3.2 and finally, we discuss the
emulator’s training and performance in subsection 3.3.

###  3.1 Accuracy and sensitivity

The computational cost of cosmological simulations depends on several
parameters. In particular, it is very sensitive to volume and resolution. With
this in mind, we want to find the optimal simulation settings that allow us to
have the target accuracy while probing large enough scales but without
requiring excessive computational resources.

To define the simulation requirements for the nDGP emulator, we start by
assessing the accuracy of COLA in predicting the boost factor of nDGP theories
against Arepo simulations. To this extent, we run high-resolution PM
simulations888We emphasize that PM here refers to the high time-resolution
used by these COLA simulations at late time, which makes the 2LPT contribution
important only at early times as discussed in section 2. It is important to
note that the screening mechanism is implemented only approximately in COLA
simulations for DGP gravity [17, 13, 15, 56]. For an alternative screening
approximation for general Vainshtein mechanism see [19]. On the other hand,
Arepo simulations solve the non-linear scalar field exactly. with the same
cosmology of Arepo simulations for several values of the nDGP parameter,
$H_{0}r_{c}$. We simulate the evolution of $2048^{3}$ particles in a box of
side $L=1024\,h^{-1}\,{\rm Mpc}$ using $6048^{3}$ mesh-grids and 75 time-steps
from redshift $z_{\rm ini}=127$ to $z_{\rm fin}=1$. The mass resolutions that
we use in the PM simulations is $M_{\rm part}\sim 1\cdot
10^{10}\,h^{-1}\,M_{\odot}$, roughly 8 times higher than the one of the Arepo
simulations. The force resolution is $\ell_{\rm Force}\equiv L/N_{\rm
mesh}^{1/3}\approx 0.17\,h^{-1}\,{\rm Mpc}$ and the time resolution is $\Delta
a\approx 0.006$. We run the simulations in GR and nDGP from the same IC. For
nDGP simulations we use the 4 values of the $H_{0}r_{c}$ parameter
($H_{0}r_{c}=$0.5, 1, 2, 5) used in the Arepo simulations suite.

Figure 7 shows a comparison between the nDGP boost factors of PM and Arepo
simulations for the different values of the $H_{0}r_{c}$ parameter. We also
include in the comparison the nDGP boost factor estimated from linear theory.
In the top panel, we show the nDGP boost factors, while in the bottom we show
the relative difference between the PM and linear boost factors from Arepo
boost factors. On large scales, PM results seem to be more consistent with
linear theory than the ones from Arepo simulations. Overall, we find $\lesssim
0.5\%$ agreement in all gravity models up to $k\sim 5\,h\,{\rm Mpc}^{-1}$
between the high-resolution simulations and Arepo simulations. This means that
the screening approximation used in COLA simulations in nDGP theory [13, 15],
once tuned on full N-body simulations for a single value of $H_{0}r_{c}$, is
sub-percent accurate in reproducing the power spectrum boost factor of nDGP
theories for a wide range of the parameter $H_{0}r_{c}$ and for different
redshift values 999We are using here the same screening settings adopted in
[15] and tuned against Ecosmog [53] results at $z\sim 0.5$..

Figure 7:  Ratio of the nDGP boost factors at redshift $z=1$ between PM
simulations and Arepo simulations in 4 nDGP gravity models as indicated in the
legend.

The COLA simulations that we run in nDGP use IC produced with the back-scaling
approach. The IC are set in such a way that the IC power spectrum is the same
for nDGP and GR simulations. This is done by integrating the growth factors
from matter domination at redshift $z=499$ and by normalising them at the
redshift of IC $z_{\rm ini}$

$D^{z_{\rm ini}}(z_{\rm ini})=1\,.$ (3.2)

Due to this choice of normalisation, COLA simulations ignore the fact that
nDGP can be different from GR also at higher redshift. In particular, COLA
simulations are usually started at redshift $z_{\rm ini}=19$, while nDGP
effects can be significant up to redshift $z\sim 100$. We can correct this by
multiplying the final power spectrum by

$\mathcal{A}_{\rm lin}(z)=\left(\frac{D_{\rm nDGP}^{z_{\rm ini}=127}(z)\cdot
D_{\rm GR}^{z_{\rm ini}=19}(z)}{D_{\rm GR}^{z_{\rm ini}=127}(z)\cdot D_{\rm
nDGP}^{z_{\rm ini}=19}(z)}\right)^{2}\,,$ (3.3)

which has the effect of re-scaling the amplitude of the power spectrum as if
the simulations were started at redshift $z=127$. This approach neglects the
effects of non-linearity before $z_{\rm ini}=19$, but these are sub-dominant
due to the weakness of the nDGP fifth force at early times. The importance of
this correction is evident from figure 8, where we show that the agreement
with Arepo of the nDGP boost factors from PM simulations started at
$z_{ini}=19$ is improved when applying the linear correction. In fact, the raw
results of $z_{ini}=19$ PM simulations show $1\%$ deviations from Arepo (left
panel), while they are in better than $0.5\%$ agreement with Arepo once
corrected using eq (3.3) (right panel).

Figure 8:  Ratio of the nDGP boost factors at redshift $z=1$ between PM
simulations started at redshift 19 and Arepo simulations in 4 nDGP gravity
models as indicated in the legend. The left panel shows the raw results from
simulations, and the right panel shows the results after applying the linear
correction in eq (3.3).

Having tested the accuracy of COLA simulations in the PM-limit in reproducing
the nDGP boost factor, we want to test the impact of force and mass resolution
on the nDGP boost factor in COLA simulations to find the optimal setup to span
the theory parameter space running a large number of simulations. To do so, we
run more COLA simulations in a box of side $L=512\,h\,{\rm Mpc}^{-1}$ with
lower mass-resolution ($M_{\rm part}\sim 8\cdot 10^{10}\,h^{-1}\,M_{\odot}$,
similarly to Arepo simulations), with time-resolution $\Delta a\approx 0.02$
and with 4 force resolutions: $N_{\rm mesh}^{1/3}=512$ (low force resolution),
$N_{\rm mesh}^{1/3}=1024$, $N_{\rm mesh}^{1/3}=1536$, $N_{\rm
mesh}^{1/3}=2048$ (high force resolution). We compare the results for the nDGP
boost factor of low mass-resolutions COLA simulations with high-resolutions PM
simulations in figure 9 where we can see that the results are converged at
sub-percent level up to $k\sim 1\,h\,{\rm Mpc}^{-1}$ and the force resolution
has a significant impact on the boost factor only on very small scales,
$k>2\,h\,{\rm Mpc}^{-1}$. The results of the intermediate cases, $N_{\rm
mesh}^{1/3}=1024$ and $N_{\rm mesh}^{1/3}=1536$, are consistent with the high-
resolution PM simulations within $2\%$ up to $k\sim 5\,h\,{\rm Mpc}^{-1}$.

Figure 9:  Ratio of the nDGP boost factors in low mass-resolution COLA
simulations with that in high-resolution PM simulations for different force
resolutions, $N_{\rm mesh}^{1/3}=512$, $N_{\rm mesh}^{1/3}=1024$, $N_{\rm
mesh}^{1/3}=1536$ and $N_{\rm mesh}^{1/3}=2048$. In different panels, we show
results for different nDGP gravity models as indicated in the top-left corner
of each panel. The black dashed lines show the typical amplitude of the boost
factor in the different gravity models for visual reference.

Before choosing a strategy to sample the parameter space for producing the
emulator’s training set, we explore the impact that each parameter has on the
boost factor. To do so we first investigate the range of values that the
parameter $H_{0}r_{c}$ should span: for stronger deviations from GR we choose
$\sim 50\%$ deviations as the upper limit, while for the other limit, we want
to "smoothly" connect to GR so we look for models where the deviation is
$<1\%$ (i.e., below our target accuracy). Based on previous results in the
literature [57], we select $H_{0}r_{c}=0.2$ as the strong deviation case and
$H_{0}r_{c}=20$ as the weak deviation case.

The results of [39] show that the response function of the DM power spectrum
for changes in cosmological parameters estimated with COLA agrees with
cosmological emulators in a wide range of cosmologies when using $\ell_{\rm
Force}\equiv L/N_{\rm part}^{1/3}=0.5\,h^{-1}\,{\rm Mpc}$. We shall use at
least the same force resolution here for the boost factor since we want to
predict it for a similar range of cosmologies. As an additional test to check
the accuracy of the boost factor for the wider range of gravity parameters, we
produce high-resolution PM simulations and use them to make a comparison with
the boost factors of $\ell_{\rm Force}=0.5\,h^{-1}\,{\rm Mpc}$ simulations.
Using a default box side of $512\,h^{-1}\,{\rm Mpc}$, we in fact compare the
boost factors obtained using $512^{3}$ particles and $N_{\rm mesh}^{1/3}=1024$
with the ones obtained using $1024^{3}$ particles and $N_{\rm
mesh}^{1/3}=3072$. We show the ratio of the two boost factors for 4 redshift
values (legend) and 3 values of the parameter $H_{0}r_{c}$ (plot titles) in
figure 10. In the most extreme case, the difference in the two boost factors
is $\lesssim 1\%$ up to $k=1\,h\,{\rm Mpc}^{-1}$ ($\lesssim 3\%$ up to
$k=5\,h\,{\rm Mpc}^{-1}$). This is in line with our target accuracy so we rely
on these low-resolution settings for the simulations that we use to conduct
the rest of our analysis.

Figure 10:  Ratio of the nDGP boost factors obtained with low-resolution COLA
simulations ($N_{\rm part}^{1/3}=512$, $N_{\rm mesh}^{1/3}=1024$) with respect
to the same boost factors but obtained with high-resolution PM simulations
($N_{\rm part}^{1/3}=1024$, $N_{\rm mesh}^{1/3}=3072$). The legend describes
the redshift value of each line.

The boost factor of nDGP is strongly dependent on the theory parameter
$H_{0}r_{c}$ and on the redshift as shown in figure 11, where we plot the
boost factors for three values of $H_{0}r_{c}$ (titles) and 4 redshift values
(legend). However, the cosmological dependence of the nDGP boost factor is
quite weak. In support of this statement, we show the effect of changing each
cosmological parameter on the nDGP boost factor for $H_{0}r_{c}=0.2$ in figure
12, $H_{0}r_{c}=2$ in figure 13 and $H_{0}r_{c}=20$ in figure 14. These
figures not only show that the response is limited to $\sim 7\%$ in the most
extreme case (variations of $\Omega_{m}$ for $H_{0}r_{c}=0.2$ at redshift
$z=0$) but also that its scale dependence is similar across the different
cosmological parameters and values of $H_{0}r_{c}$. The amplitude of the
effect is smaller for larger $H_{0}r_{c}$ values.

Figure 11:  nDGP boost factor for three values of the parameter $H_{0}r_{c}$.
The legend in figure 10 describes the redshift value of each line.  Figure 12:
Change on the nDGP boost factor for $H_{0}r_{c}=0.2$ due to the variation of
each cosmological parameter (specified in each panel) with respect to the
reference cosmology. The redshift values of each line are as described in the
legend. Figure 13:  Same as figure 12 but for $H_{0}r_{c}=2$ Figure 14:  Same
as figure 12 but for $H_{0}r_{c}=20$

###  3.2 Simulations and data sets

We produce two independent simulation suites, one for the training set and the
other for the test set. Both training and test set simulations use $512^{3}$
particles to track the evolution of the matter density field in a
$(512\,h^{-1}\,{\rm Mpc})^{3}$ box. The forces are calculated using $1024^{3}$
mesh grids. The simulations are started at $z_{\rm ini}=19$ and evolved up to
$z_{\rm fin}=0$ using 50 time-steps linearly spaced in the scale factor. We
sample the cosmological parameter space using the Latin Hypercube Sampling
(LHS) technique [58], a space-filling sampling method widely used in computer
experiments to produce near-random sets of points which represent the real
variability of multi-dimensional sample spaces. The range of cosmological
parameters spanned by the training and test sets is described in table 1. We
use 20 cosmologies101010In this work we assume $k_{\rm pivot}=0.05\,{\rm
Mpc}^{-1}$ for the primordial power spectrum. for the training set and 10 for
the test set, which are listed in table 2 and table 3 respectively. For the
modified gravity parameter $H_{0}r_{c}$ we opt for a logarithmic sampling with
21 points in the range $[0.2,20]$ to produce the training set, and with 10
log-random points in the same range for the test set. The redshift sampling is
determined by the time-stepping of COLA simulations. For each cosmology in the
training set (test set), we run a vanilla GR simulation and nDGP simulations
in all the gravity models of the training set (test set) using the same IC.

Parameter | Min. | Max.
---|---|---
$h$ | $0.61$ | $0.73$
$\Omega_{\rm b}$ | $0.04$ | $0.06$
$\Omega_{\rm m}$ | $0.28$ | $0.36$
$n_{\rm s}$ | $0.92$ | $1.0$
$A_{\rm s}$ | $1.7\times 10^{-9}$ | $2.5\times 10^{-9}$
Table 1: Range of cosmological parameters spanned by the nDGP emulator.

Each simulation takes roughly 28 minutes on a single node of the Sciama High-
Performance Compute (HPC) cluster using 32 tasks111111The specific partition
used to run the simulations mounts Intel Xeon E5-2683v4 processors.. This
corresponds to $\approx 15$ CPU hours. In total, we ran $22\times 20=440$
simulations for the training set and $11\times 10=110$ simulations for the
test set. This means that the total computational cost for the simulations is
$\approx 8\cdot 10^{3}$ CPU hours. In particular, we have used a partition
counting 12 nodes with 32 cores each for a total of 384 cores (a typical value
for an HPC cluster) obtaining all the power spectra in less than one wall-
clock day.

$\Omega_{\rm m}$ | $\Omega_{\rm b}$ | $h$ | $n_{s}$ |  $A_{s}$ [$\cdot 10^{-9}$]
---|---|---|---|---
0.282 | 0.0515 | 0.667 | 0.950 | 2.12
0.286 | 0.0445 | 0.655 | 0.990 | 2.44
0.294 | 0.0585 | 0.697 | 0.938 | 2.00
0.298 | 0.0415 | 0.619 | 0.978 | 2.04
0.290 | 0.0545 | 0.625 | 0.986 | 1.80
0.302 | 0.0485 | 0.709 | 0.922 | 1.72
0.306 | 0.0435 | 0.649 | 0.942 | 2.36
0.314 | 0.0595 | 0.721 | 0.958 | 2.28
0.318 | 0.0575 | 0.661 | 0.962 | 2.32
0.310 | 0.0405 | 0.673 | 0.966 | 1.88

$\Omega_{\rm m}$ | $\Omega_{\rm b}$ | $h$ | $n_{s}$ |  $A_{s}$ [$\cdot 10^{-9}$]
---|---|---|---|---
0.322 | 0.0465 | 0.727 | 0.982 | 1.96
0.326 | 0.0455 | 0.637 | 0.994 | 2.08
0.334 | 0.0535 | 0.631 | 0.974 | 2.16
0.338 | 0.0555 | 0.703 | 0.954 | 1.92
0.330 | 0.0525 | 0.613 | 0.934 | 2.24
0.342 | 0.0495 | 0.685 | 0.998 | 2.20
0.346 | 0.0565 | 0.691 | 0.970 | 2.48
0.354 | 0.0505 | 0.643 | 0.930 | 1.76
0.358 | 0.0425 | 0.679 | 0.926 | 2.40
0.350 | 0.0475 | 0.715 | 0.946 | 1.84

Table 2: List of cosmologies used to create the training set for the nDGP
emulator.

$\Omega_{\rm m}$ | $\Omega_{\rm b}$ | $h$ | $n_{s}$ |  $A_{s}$ [$\cdot 10^{-9}$]
---|---|---|---|---
0.284 | 0.0410 | 0.700 | 0.956 | 1.82
0.292 | 0.0490 | 0.628 | 0.940 | 1.74
0.308 | 0.0510 | 0.652 | 0.964 | 2.22
0.316 | 0.0430 | 0.688 | 0.996 | 2.30
0.324 | 0.0470 | 0.712 | 0.932 | 2.06

$\Omega_{\rm m}$ | $\Omega_{\rm b}$ | $h$ | $n_{s}$ |  $A_{s}$ [$\cdot 10^{-9}$]
---|---|---|---|---
0.332 | 0.0530 | 0.724 | 0.980 | 2.38
0.348 | 0.0450 | 0.616 | 0.972 | 1.98
0.340 | 0.0550 | 0.664 | 0.948 | 1.90
0.356 | 0.0570 | 0.640 | 0.988 | 2.46
0.300 | 0.0590 | 0.676 | 0.924 | 2.14

Table 3: List of cosmologies used to create the test set for the nDGP
emulator.

The power spectra are computed by interpolating the DM distribution on a
$1024^{3}$ mesh grid with the cloud-in-cell mass-assignment scheme. The power
spectra are corrected for the window function used in the interpolation and
for shot noise. We use bins of width $\Delta k=k_{f}$ between $k_{\rm
min}=\frac{1}{2}k_{f}$ and $k_{\rm Nyquist}=5\,h\,{\rm Mpc}^{-1}$, where
$k_{f}$ is the fundamental frequency of the box. Our target accuracy is $1\%$
up to $k=1\,h\,{\rm Mpc}^{-1}$, scales up to which we have already proven COLA
simulations with these specifications to have $\sim 1\%$ reliable
determination of the cosmological response function [39]. However for
practical applications (e.g. weak lensing predictions) we extend the emulation
range to $k=5\,h\,{\rm Mpc}^{-1}$, where the analysis discussed in the
previous section suggests an accuracy of our simulations of $\sim 3\%$ for the
nDGP boost factor. Beyond $k=5\,h\,{\rm Mpc}^{-1}$, baryonic effects are
dominant over MG effects and it is still unclear whether we can extract
additional information from these small scales [3]. As we have seen that mass
resolution is a limiting factor at high redshift, we make the conservative
choice of restricting our data sets to the power spectra at $z_{\rm max}=2$.
We take the ratio of the power spectra in nDGP theories with the respective
power spectrum in GR for each redshift and cosmology to obtain the nDGP boost
factors. As the nDGP and GR simulations are run from the same IC, the sample
variance is largely cancelled out when computing the nDGP boost factors. Since
our simulations are started at redshift $z_{\rm ini}=19$ we multiply the nDGP
boost factors with the linear theory correction in eq. (3.3).

Due to the numerical methods used to estimate the power spectrum of the N-body
particles, the nDGP boost factors that we obtain are sampled with a linear
binning in $k$-space and are affected by finite resolution noise at small
scales. This noise is non-physical and may worsen the interpolation accuracy,
therefore we decide to smooth the boost factors. We compare two smoothing
techniques: one based on a logarithmic re-sampling of the data, the other
based on the Savitzky-Golay filter [59]. In the left plot of Fig 15 is shown a
comparison of the two smoothing techniques for the central case of
$H_{0}r_{c}=2$ and $z=0$. The two techniques give consistent results so we
decided to perform the rest of the analysis with the Savitzky-Golay filter.

Multi-layer perceptrons (MLP), simple neural networks that we want to employ
for the parameter space interpolation, work better when they deal with
normalised input values [60], hence we re-scale our parameters accordingly. We
map the cosmological parameters in such a way that the minimum and maximum
values they assume in parameter space are 0 and 1. The redshift is naturally
mapped in the range $[\frac{1}{3},1]$ by computing the scale factor. We
finally convert the MG parameter $H_{0}r_{c}$ into
$W_{r_{c}}\equiv\frac{0.2}{H_{0}r_{c}}$ which spans $[0.01,1]$ for our
$H_{0}r_{c}$ range (or more precisely it spans $[0,1]$ as we include the GR
trivial results). In the right plot of Fig 15, the distribution of the
normalised parameters is summarised in 10 linear bins between 0 and 1.

Figure 15:  On the left. Comparison of the two smoothing techniques for the
boost factor of nDGP with $H_{0}r_{c}=2$. On the right. Distribution of the
normalised parameters across the training set in 10 bins between 0 and 1: all
but $W_{r_{c}}$ and $a$ are uniformly distributed, simplifying the emulation
problem.

We convert the boost factor with the formula

$B_{\rm Norm}(k)=-\frac{\log_{10}(B(k)-0.999)}{3}\,,$ (3.4)

to give more importance to small deviations from GR than to larger
deviations121212Notice that the boost for $W_{r_{c}}=0$ (limit in which we
recover GR) is trivial, i.e. $B(k)=1$. In this case, we get $B_{\rm Norm}=1$
which explains the choice of the factor $1/3$ in eq (3.4). The subtraction of
the value $0.999$ inside the logarithm is instead performed to place the GR
result at a reasonable distance from the nDGP weakest deviation result in the
$B_{\rm Norm}$-space.. The simplicity and self-similarity of the nDGP boost
factor curve suggest that the 407 $k$-modes that we are using to represent the
boost factor are probably unnecessary. After subtracting the mean value of the
boost factor from each element of the dataset, we perform a Principal
Component Analysis (PCA) [61] to study if the variance of the training set is
significantly larger in some $k$-modes-space directions than in others. We
find that the component with the greatest variance accounts for more than
$99\%$ of the variance, while by using the first three components it is
possible to describe $\sim 99.98\%$ of the variance. This means that we can
reduce the dimensionality of the problem from 407 output values (one for each
$k$ bin) to 3 output values losing only $\sim 0.02\%$ of the information.

###  3.3 Emulator

We train an MLP with one hidden layer of 100 nodes activated by hyperbolic
tangent functions with the limited-memory Broyden–Fletcher–Goldfarb–Shanno
[62] optimiser131313This iterative optimisation technique exploits a numerical
estimate of the Hessian based only on gradient evaluations to determine the
descent direction, which makes the technique computationally more efficient
than Newton’s method. on the training set. This converges in less than 5000
iterations, producing the nDGP emulator that we benchmark on the test set. In
the top panel of figure 16 we show some boost factors examples taken from the
test set (solid lines) and compared with the emulator’s predictions (dashed
lines). It is worth noticing that the emulated boost factors are free of high-
frequency noise thanks to the smoothing discussed in sub-section 3.2. In the
bottom panel of figure 16 we show the mean (black dashed line) and variance
($1\sigma$ and $2\sigma$ shaded regions) of the emulation error on the test
set together with the residuals for the 30 examples shown in the top panel.
The $2\sigma$ contours are well within the $1\%$ threshold at all scales.

Figure 16:  Comparison between boost factors of the test set (solid lines) and
emulated boost factors (dashed lines) for 30 examples randomly chosen in the
test set. The residuals between emulated and true boost factors (coloured
solid lines) are shown in the bottom panel together with the mean (dashed
black line) and variance ($1\sigma$ and $2\sigma$ shaded regions) of the
emulation error on the test-set.

As an additional test of the emulation accuracy, we compare the emulator’s
predictions with the boost computed from Arepo simulations in nDGP gravity
using 4 values of $H_{0}r_{c}$: 0.5, 1, 2, 5. Thanks to the flexibility of the
emulator, unlike in the comparison between PM and Arepo simulations of figure
7, we are able in this case to compare the boosts at each of the 9 output
redshifts of the Arepo simulation suite within the interpolation range. The
results are shown in figure 17. The emulator predicts the nDGP boost factor of
Arepo with $\sim 2\%$ accuracy in all cases and at all the scales considered.
Additionally, we have checked that the emulator performance reaches $\sim 1\%$
accuracy if restricting the redshift range to $z\in[0.5,2]$. This is made
possible by the excellent agreement between the nDGP boost factors of COLA and
Arepo simulations but also by the linear-theory correction of eq (3.3) applied
to the nDGP boost factors of COLA simulations started at redshift
$z_{ini}=19$.

Figure 17:  Comparison between emulator predictions (dashed lines) and Arepo
results (solid lines) of the nDGP boost factor at 9 redshifts $z\in[0,2]$ for
4 values of the parameter $H_{0}r_{c}$: 0.5 (in red), 1 (in orange), 2 (in
blue) and 5 (in green). The bottom panel shows the residuals.

It is important to notice that even the results of these Arepo simulations
might be not perfectly converged, but it is beyond the scope of this paper to
assess the accuracy of these full N-body simulations. Furthermore, we have
Arepo results for specific values of the modified gravity parameter
$H_{0}r_{c}$, so the accuracy of the emulator at the edges of the parameter
space cannot be fully assessed. With these caveats, and considering the level
of convergence of the COLA predictions for the nDGP boost factor discussed in
subsection 3.1, we estimate the accuracy of the emulator to be of $\sim 2\%$
up to $k=1\,h\,{\rm Mpc}^{-1}$ and of $\sim 3\%$ up to $k=5\,h\,{\rm
Mpc}^{-1}$.

The biggest computational cost behind the nDGP emulator comes from running the
COLA simulations necessary to produce the power spectra in GR and nDGP
theories ($\sim 10^{4}$ CPU hours). The data processing applied to the power
spectra in order to produce the training set has a minor impact on the
computational cost ($\sim 10$ CPU minutes). Also, the training of the emulator
is relatively fast ($\sim 1$ CPU hour). Once the model is trained, predictions
are extremely fast. On an average single-core CPU, the model can produce $\sim
5000$ predictions in a second. The prediction speed of emulators is key to
enabling fast cosmological inference with Monte Carlo techniques which require
$10^{5}-10^{6}$ evaluations [63, 64, 65].

As we have shown here, DM power spectrum emulators can be extended to MG
models, and in particular to nDGP gravity, with modest computational cost
thanks to COLA simulations, enabling accurate and extremely fast predictions
of the matter power spectrum up to $k=5\,h^{-1}\,{\rm Mpc}$.

##  4 Conclusions

This article focused on the application of the COLA simulations for
predictions of the matter power spectrum in GR and Modified Gravity (MG)
theories. We addressed the issue of errors arising from finite resolution in
numerical implementations of simulations, including time-resolution, force-
resolution, and mass-resolution. Based on our analysis we designed and
exploited a COLA simulation suite to extend cosmological emulators to nDGP
gravity. Together with this paper, we release the nDGP emulator that will be
publicly available at the address: https://github.com/BartolomeoF/nDGPemu.

To study the convergence of COLA simulations in $\Lambda$CDM, we argued that
COLA converges to PM-only simulations increasing the time resolution, and
referred to results obtained in this limit as PM results. By comparing high
force-resolution PM simulations with Bacco emulator (trained with high mass-
resolution simulations) and low mass-resolution N-bodysimulations (Arepo) for
predictions of the matter power spectrum at redshift $z=1$, we highlighted the
importance of mass resolution, responsible for $\sim 4\%$ differences at
$k=1\,h\,{\rm Mpc}^{-1}$, while showing that PM simulations reproduce the
reference results in both the low mass-resolution and high-mass-resolution
cases with $\sim 1\%$ accuracy up to $k\sim 2\,h\,{\rm Mpc}^{-1}$. After that,
we produced a suite of simulations in $\Lambda$CDM with the same cosmology but
with different numbers of particles, force grids and time-steps, and using
high-resolution PM simulations as a reference, we tested the relative
convergence of COLA simulations. We concluded that the power spectra converge
towards the reference only when mass-resolution, force-resolution and time-
resolution are increased accordingly.

To showcase the potentiality of COLA for extending cosmological emulators to
MG theories, we developed an nDGP boost-factor emulator of the matter power
spectrum in section 3. By using Arepo predictions as a reference, we
investigated the accuracy of COLA simulations for the nDGP boost factor in
relation to the resolution settings and identified minimal requirements for
our target accuracy. To include MG effects at $z>z_{\rm ini}$ in the back-
scaling approach, we proposed and validated the linear theory correction of eq
(3.3). We also studied the response of the nDGP boost factors to the change of
cosmological parameters finding weak cosmological dependence of the nDGP boost
factor and confirming that the cosmological dependence becomes weaker for
smaller deviations from $\Lambda$CDM (larger values of $H_{0}r_{c}$). In sub-
section 3.2, we detailed our simulation strategy to compute the nDGP boost
factors and examined the computational cost consisting of a total of $\approx
8\cdot 10^{3}$ CPU hours. After smoothing the signals and remapping the
parameters in a unitary interval, we performed a PCA and selected only 3
components while preserving $99.98\%$ of the information. We then used the
resulting data set to train a neural network and bench-marked its accuracy on
an independent test set ($<1\%$ at all scales considered). As an additional
test, we compared the emulator predictions with the Arepo results and found
better than $2\%$ agreement at all scales considered. We concluded that COLA
simulations can be used to accurately extend matter power spectrum emulators
to MG theories up to $k=1\,h\,{\rm Mpc}^{-1}$ ($k=5\,h\,{\rm Mpc}^{-1}$) with
$2\%$ accuracy ($3\%$ accuracy).

The approach developed in this work can be extended to other gravity models.
In particular it would be interesting to train an emulator for the matter
power spectrum in Horndeski gravity [66, 67] which is made possible thanks to
the recent extension of the COLA method [14, 39, 19] to the sub-class of
Horndeski theories incorporating the Vainshtein screening mechanism.
Furthermore, we plan to extend the convergence of COLA simulations for
predictions of halo statistics for which we have already collected promising
results. Finally, we aim to develop a suite of COLA-based emulators for
several LSS probes, such as the ones studied in [15, 37].

## Acknowledgments

We thank Albert Izard, Daniela Saadeh and Jamie Donald-McCann for useful
insight and stimulating discussions on COLA simulations and cosmological
emulators. We thank Hans Winther for creating the FML library and continuously
improving it with new features. We thank César Hernández-Aguayo for providing
access to the MG-Arepo simulations. For the purpose of open access, the
authors have applied a Creative Commons Attribution (CC BY) licence to any
Author Accepted Manuscript version arising from this work. Numerical
computations were done on the Sciama High Performance Compute (HPC) cluster
which is supported by the ICG, SEPNet and the University of Portsmouth. BF is
supported by a Royal Society Enhancement Award (grant no.
RF$\backslash$ERE$\backslash$210304). KK is supported by the STFC grant
ST/W001225/1. TB is supported by ERC Starting Grant SHADE (grant no. StG
949572) and a Royal Society University Research Fellowship (grant no.
URF$\backslash$R1$\backslash$180009).

#### Data availability

Supporting research data are available on reasonable request from the
corresponding author.

## References

  * [1] DESI collaboration, _Validation of the Scientific Program for the Dark Energy Spectroscopic Instrument_ , 2306.06307.
  * [2] N.E. Chisari et al., _Modelling baryonic feedback for survey cosmology_ , _Open J. Astrophys._ 2 (2019) 4 [1905.06082].
  * [3] DES collaboration, _Dark Energy Survey Year 3 results: Cosmological constraints from galaxy clustering and weak lensing_ , _Phys. Rev. D_ 105 (2022) 023520 [2105.13549].
  * [4] J.U. Lange, A.P. Hearin, A. Leauthaud, F.C. van den Bosch, E. Xhakaj, H. Guo et al., _Constraints on $S_{8}$ from a full-scale and full-shape analysis of redshift-space clustering and galaxy-galaxy lensing in BOSS_, 2301.08692.
  * [5] S. Tassev, M. Zaldarriaga and D. Eisenstein, _Solving Large Scale Structure in Ten Easy Steps with COLA_ , _JCAP_ 06 (2013) 036 [1301.0322].
  * [6] 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].
  * [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] M.J. Drinkwater et al., _The WiggleZ Dark Energy Survey: Survey Design and First Data Release_ , _Mon. Not. Roy. Astron. Soc._ 401 (2010) 1429 [0911.4246].
  * [9] V. Springel, _The Cosmological simulation code GADGET-2_ , _Mon. Not. Roy. Astron. Soc._ 364 (2005) 1105 [astro-ph/0505010].
  * [10] 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].
  * [11] A. Izard, P. Fosalba and M. Crocce, _ICE-COLA: fast simulations for weak lensing observables_ , _Mon. Not. Roy. Astron. Soc._ 473 (2018) 3051 [1707.06312].
  * [12] 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].
  * [13] 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].
  * [14] LSST Dark Energy Science collaboration, _Hi-COLA: fast, approximate simulations of structure formation in Horndeski gravity_ , _JCAP_ 03 (2023) 040 [2209.01666].
  * [15] B. Fiorini, K. Koyama, A. Izard, H.A. Winther, B.S. Wright and B. Li, _Fast generation of mock galaxy catalogues in modified gravity models with COLA_ , _JCAP_ 09 (2021) 021 [2106.05197].
  * [16] B.S. Wright, H.A. Winther and K. Koyama, _COLA with massive neutrinos_ , _JCAP_ 10 (2017) 054 [1705.08165].
  * [17] H.A. Winther and P.G. Ferreira, _Fast route to nonlinear clustering statistics in modified gravity theories_ , _Phys. Rev. D_ 91 (2015) 123507 [1403.6492].
  * [18] R. Scoccimarro, _Large-Scale Structure in Brane-Induced Gravity I. Perturbation Theory_ , _Phys. Rev. D_ 80 (2009) 104006 [0906.4545].
  * [19] G. Brando, K. Koyama and H.A. Winther, _Revisiting Vainshtein Screening for fast N-body simulations_ , 2303.09549.
  * [20] S. Agarwal, F.B. Abdalla, H.A. Feldman, O. Lahav and S.A. Thomas, _PkANN - I. Non-linear matter power spectrum interpolation through artificial neural networks_ , _"Mon. Not. Roy. Astron. Soc."_ 424 (2012) 1409 [1203.1695].
  * [21] S. Habib, K. Heitmann, D. Higdon, C. Nakhleh and B. Williams, _Cosmic Calibration: Constraints from the Matter Power Spectrum and the Cosmic Microwave Background_ , _Phys. Rev. D_ 76 (2007) 083503 [astro-ph/0702348].
  * [22] VIRGO Consortium collaboration, _Stable clustering, the halo model and nonlinear cosmological power spectra_ , _Mon. Not. Roy. Astron. Soc._ 341 (2003) 1311 [astro-ph/0207664].
  * [23] Euclid collaboration, _Euclid preparation: IX. EuclidEmulator2 – power spectrum emulation with massive neutrinos and self-consistent dark energy perturbations_ , _Mon. Not. Roy. Astron. Soc._ 505 (2021) 2840 [2010.11288].
  * [24] R.E. Angulo, M. Zennaro, S. Contreras, G. Aricò, M. Pellejero-Ibañez and J. Stücker, _The BACCO simulation project: exploiting the full power of large-scale structure for cosmology_ , _Mon. Not. Roy. Astron. Soc._ 507 (2021) 5869 [2004.06245].
  * [25] H. Winther, S. Casas, M. Baldi, K. Koyama, B. Li, L. Lombriser et al., _Emulators for the nonlinear matter power spectrum beyond $\Lambda$CDM_, _Phys. Rev. D_ 100 (2019) 123540 [1903.08798].
  * [26] LSST Dark Energy Science collaboration, _Matter Power Spectrum Emulator for f(R) Modified Gravity Cosmologies_ , _Phys. Rev. D_ 103 (2021) 123525 [2010.00596].
  * [27] C. Arnold, B. Li, B. Giblin, J. Harnois-Déraps and Y.-C. Cai, _FORGE – the f(R) gravity cosmic emulator project I: Introduction and matter power spectrum emulator_ , 2109.04984.
  * [28] I.n. Sáez-Casares, Y. Rasera and B. Li, _The e-MANTIS emulator: fast predictions of the non-linear matter power spectrum in $f(R)$CDM cosmology_, 2303.08899.
  * [29] R. Mauland, H.A. Winther and C.-Z. Ruan, _Sesame: A power spectrum emulator pipeline for beyond- $\Lambda$CDM models_, 2309.13295.
  * [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] J. Harnois-Déraps, C. Hernandez-Aguayo, C. Cuesta-Lazaro, C. Arnold, B. Li, C.T. Davies et al., _mglens: Modified gravity weak lensing simulations for emulation-based cosmological inference_ , _Mon. Not. Roy. Astron. Soc._ 525 (2023) 6336 [2211.05779].
  * [32] 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].
  * [33] B. Bose, M. Tsedrik, J. Kennedy, L. Lombriser, A. Pourtsidou and A. Taylor, _Fast and accurate predictions of the nonlinear matter power spectrum for general models of Dark Energy and Modified Gravity_ , 2210.01094.
  * [34] H.A. Winther et al., _Modified Gravity N-body Code Comparison Project_ , _Mon. Not. Roy. Astron. Soc._ 454 (2015) 4208 [1506.06384].
  * [35] A. Klypin and F. Prada, _Dark matter statistics for large galaxy catalogues: power spectra and covariance matrices_ , _Mon. Not. Roy. Astron. Soc._ 478 (2018) 4602 [1701.05690].
  * [36] J. DeRose, R.H. Wechsler, J.L. Tinker, M.R. Becker, Y.-Y. Mao, T. McClintock et al., _The Aemulus Project I: Numerical Simulations for Precision Cosmology_ , _Astrophys. J._ 875 (2019) 69 [1804.05865].
  * [37] B. Fiorini, K. Koyama and A. Izard, _Studying large-scale structure probes of modified gravity with COLA_ , 2208.01345.
  * [38] C. Hernández-Aguayo, C.-Z. Ruan, B. Li, C. Arnold, C.M. Baugh, A. Klypin et al., _Fast full N-body simulations of generic modified gravity: derivative coupling models_ , _JCAP_ 01 (2022) 048 [2110.00566].
  * [39] G. Brando, B. Fiorini, K. Koyama and H.A. Winther, _Enabling matter power spectrum emulation in beyond- $\Lambda$CDM cosmologies with COLA_, _JCAP_ 09 (2022) 051 [2203.11120].
  * [40] C. Arnold, M. Leo and B. Li, _Realistic simulations of galaxy formation in $f(R)$ modified gravity_, _Nature Astron._ 3 (2019) 945 [1907.02977].
  * [41] C. Hernández-Aguayo, C. Arnold, B. Li and C.M. Baugh, _Galaxy formation in the brane world I: overview and first results_ , _Mon. Not. Roy. Astron. Soc._ 503 (2021) 3867 [2006.15467].
  * [42] V. Springel, _E pur si muove: Galiliean-invariant cosmological hydrodynamical simulations on a moving mesh_ , _Mon. Not. Roy. Astron. Soc._ 401 (2010) 791 [0901.4107].
  * [43] M.A. Mitchell, C. Hernández-Aguayo, C. Arnold and B. Li, _A general framework to test gravity using galaxy clusters IV: cluster and halo properties in DGP gravity_ , _Mon. Not. Roy. Astron. Soc._ 508 (2021) 4140 [2106.13815].
  * [44] Euclid collaboration, _Euclid preparation: VII. Forecast validation for Euclid cosmological probes_ , _Astron. Astrophys._ 642 (2020) A191 [1910.09273].
  * [45] Y. Feng, M.-Y. Chu, U. Seljak and P. McDonald, _FastPM: a new scheme for fast simulations of dark matter and haloes_ , _Mon. Not. Roy. Astron. Soc._ 463 (2016) 2273 [1603.00476].
  * [46] C. Deffayet, _Cosmology on a brane in Minkowski bulk_ , _Phys. Lett. B_ 502 (2001) 199 [hep-th/0010186].
  * [47] M.A. Luty, M. Porrati and R. Rattazzi, _Strong interactions and stability in the DGP model_ , _JHEP_ 09 (2003) 029 [hep-th/0303116].
  * [48] K. Koyama and F.P. Silva, _Non-linear interactions in a cosmological background in the DGP braneworld_ , _Phys. Rev. D_ 75 (2007) 084040 [hep-th/0702169].
  * [49] K. Koyama, _Ghosts in the self-accelerating universe_ , _Class. Quant. Grav._ 24 (2007) R231 [0709.2399].
  * [50] F. Schmidt, _Self-Consistent Cosmological Simulations of DGP Braneworld Gravity_ , _Phys. Rev. D_ 80 (2009) 043001 [0905.0858].
  * [51] L. Lombriser, W. Hu, W. Fang and U. Seljak, _Cosmological Constraints on DGP Braneworld Gravity with Brane Tension_ , _Phys. Rev. D_ 80 (2009) 063536 [0905.1112].
  * [52] C. Deffayet, G.R. Dvali, G. Gabadadze and A.I. Vainshtein, _Nonperturbative continuity in graviton mass versus perturbative discontinuity_ , _Phys. Rev. D_ 65 (2002) 044026 [hep-th/0106001].
  * [53] B. Li, G.-B. Zhao, R. Teyssier and K. Koyama, _ECOSMOG: An Efficient Code for Simulating Modified Gravity_ , _JCAP_ 01 (2012) 051 [1110.1379].
  * [54] C. Hernández-Aguayo, C. Arnold, B. Li and C.M. Baugh, _Galaxy formation in the brane world I: overview and first results_ , _Mon. Not. Roy. Astron. Soc._ 503 (2021) 3867 [2006.15467].
  * [55] 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].
  * [56] B. Fiorini, _Testing gravity on cosmological scales : theoretical predictions with the COLA method_ , Ph.D. thesis, University of Portsmouth, 2023\.  2303.07121.
  * [57] C.-Z. Ruan, C. Hernández-Aguayo, B. Li, C. Arnold, C.M. Baugh, A. Klypin et al., _Fast full N-body simulations of generic modified gravity: conformal coupling models_ , _JCAP_ 05 (2022) 018 [2110.00328].
  * [58] M.D. McKay, R.J. Beckman and W.J. Conover, _Comparison of three methods for selecting values of input variables in the analysis of output from a computer code_ , _Technometrics_ 21 (1979) 239 [https://doi.org/10.1080/00401706.1979.10489755].
  * [59] A. Savitzky and M.J.E. Golay, _Smoothing and differentiation of data by simplified least squares procedures._ , _Analytical Chemistry_ 36 (1964) 1627 [https://doi.org/10.1021/ac60214a047].
  * [60] C.M. Bishop, _Neural Networks for Pattern Recognition_ , Oxford University Press, Inc., USA (1995).
  * [61] I.T. Jolliffe, _Principal Component Analysis_ , Springer Series in Statistics, Springer-Verlag (2002), 10.1007/b98835.
  * [62] R. Fletcher, _Practical Methods of Optimization_ , Wiley (2000).
  * [63] J. Donald-McCann, F. Beutler, K. Koyama and M. Karamanis, _matryoshka: halo model emulator for the galaxy power spectrum_ , _Mon. Not. Roy. Astron. Soc._ 511 (2022) 3768 [2109.15236].
  * [64] A. Spurio Mancini, D. Piras, J. Alsing, B. Joachimi and M.P. Hobson, _CosmoPower: emulating cosmological power spectra for accelerated Bayesian inference from next-generation surveys_ , _Mon. Not. Roy. Astron. Soc._ 511 (2022) 1771 [2106.03846].
  * [65] J. Donald-McCann, K. Koyama and F. Beutler, _matryoshka II: Accelerating Effective Field Theory Analyses of the Galaxy Power Spectrum_, 2202.07557.
  * [66] G.W. Horndeski, _Second-order scalar-tensor field equations in a four-dimensional space_ , _Int. J. Theor. Phys._ 10 (1974) 363.
  * [67] C. Deffayet, X. Gao, D.A. Steer and G. Zahariade, _From k-essence to generalised Galileons_ , _Phys. Rev. D_ 84 (2011) 064039 [1103.3260].