# Fast and accurate predictions of the nonlinear matter power spectrum for
general models of Dark Energy and Modified Gravity

B. Bose1,2, M. Tsedrik1, J. Kennedy2,3, L. Lombriser2, A. Pourtsidou1,4, A.
Taylor1,4.
1Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh, EH9 3HJ, U.K.
2Département de Physique Théorique, Université de Genève, 24 quai Ernest
Ansermet, 1211 Genève 4, Switzerland.
3The Global Academy of Agriculture and Food Systems, University of Edinburgh,
Easter Bush Campus, Charnock Bradley Building, EH25 9RG.
4Higgs Centre for Theoretical Physics, School of Physics and Astronomy, The
University of Edinburgh, Edinburgh EH9 3FD, UK<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We embed linear and nonlinear parametrisations of beyond standard cosmological
physics in the halo model reaction framework, providing a model-independent
prescription for the nonlinear matter power spectrum. As an application, we
focus on Horndeski theories, using the Effective Field Theory of Dark Energy
(EFTofDE) to parameterise linear and quasi-nonlinear perturbations. In the
nonlinear regime we investigate both a nonlinear parameterised-post
Friedmannian (nPPF) approach as well as a physically motivated and approximate
phenomenological model based on the error function (Erf). We compare the
parameterised approaches’ predictions of the nonlinear matter power spectrum
to the exact solutions, as well as state-of-the-art emulators, in an evolving
dark energy scenario and two well studied modified gravity models, finding
sub-percent agreement in the reaction using the Erf model at $z\leq 1$ and
$k\leq 5~{}h/{\rm Mpc}$. This suggests only an additional 3 free constants,
above the background and linear theory parameters, are sufficient to model
nonlinear, non-standard cosmology in the matter power spectrum at scales down
to $k\leq 3h~{}/{\rm Mpc}$ within $2\%$ accuracy. We implement the
parametrisations into ver.2.0 of the ReACT code: ACTio et ReACTio.

###### keywords:

cosmology: theory – large-scale structure of the Universe – methods:
analytical – methods: numerical

††pubyear: 2022††pagerange: Fast and accurate predictions of the nonlinear
matter power spectrum for general models of Dark Energy and Modified
Gravity–B.4

##  1 Introduction

Fundamental models of nature generally begin with an action, which when
combined with the principle of least action, gives us the temporal and spatial
dynamics of the system. For the physical system that is our Universe (U), the
action is widely accepted to be the action associated with general relativity
(GR), the Einstein-Hilbert (EH) action, together with a matter contribution
and cosmological constant

$S_{U}=S_{\rm EH}+S_{\rm M}=\int
d^{4}x\sqrt{-g}\left[\frac{R}{2\kappa^{2}}-\frac{\Lambda}{\kappa^{2}}\right]+S_{\rm
M}\,,$ (1)

where $\kappa^{2}=8\pi G_{\rm N}$, $G_{\rm N}$ being Newton’s gravitational
constant and $R$ is the 4-dimensional Ricci scalar that gauges the curvature
of spacetime. $S_{\rm M}$ is the action of the matter content of the Universe,
usually approximated by a perfect, pressureless fluid, but in general will
contain all Standard Model fields. $\Lambda$ is the (cosmological) constant
that can appear naturally in a 4-dimensional action without violating
preferred symmetries (see, for example, Fernandes et al., 2022). This constant
is measured to be non-zero by a suite of cosmological probes such as the
cosmic microwave background (CMB) radiation (Aghanim et al., 2020), type 1a
supernovae (Riess et al., 1998; Perlmutter et al., 1999), and optical galaxy
surveys (see, for example, Alam et al., 2021). This has led to the standard
model of cosmology, $\Lambda$CDM, where CDM stands for cold dark matter111CDM
is the primary matter component in this model, outweighing baryonic matter
five fold according to cosmological and astrophysical measurements such as the
CMB..

Consistently, we would expect a non-zero cosmological constant from quantum
field theory (QFT) predictions, as all vacuum states of standard model
particle fields will contribute an energy density, $\rho_{\rm vac}$, to the
Universe that appears as a constant in the model’s action. Unfortunately, this
results in one of the biggest problems in physics (see Martin, 2012, for a
review). The first aspect of the problem is that our naïve (QFT) predictions
for the energy density of $\Lambda$ is at least 60 orders of magnitude larger
than the (cosmological) measured value222This depends on the energy scales we
are considering in the QFT calculation.. We can still make a (fine) tuning of
the ‘bare’ constant $\Lambda_{\rm bare}$ in the potentials of these fields to
cancel the other vacuum energy contributions to yield the observed value for
$\Lambda$.

One might be fine with this, after all QFTs are used to removing divergences
through renormalisation techniques. The real problem is that we need to
repeatedly fine tune every time a new energy scale or particle field is
considered which changes $\rho_{\rm vac}$ (this can also happen through phase
transitions) (see Padilla, 2015, for an overview). In other words, the value
of $\Lambda$, which is a low energy physics parameter, is incredibly sensitive
to the high energy physics, which is not technically natural and in apparent
opposition to our wide spread employment of effective field theories. These
two aspects of the problem are often referred to collectively as the
‘cosmological constant problem’. We refer the interested reader to the seminal
paper by Weinberg (1989) for a review and the famous no-go theorem which
implicitly delineates possible solutions to the problem.

Prospective solutions to these problems include gravitationally screening the
vacuum energy from our observations by using a scalar field (for example,
Charmousis et al., 2012; Appleby & Linder, 2020; Sobral Blanco & Lombriser,
2020; Khan & Taylor, 2022) or using extra spacetime dimensions (for example,
Burgess, 2004). These solutions would of course also need to produce a small
residual energy that can be used to explain our cosmological measurements, in
particular those associated with an accelerated spatial expansion. This
distinct issue can be called the ‘dark energy problem’, which may be explained
through a variation in the fundamental constants of nature such as Newton’s
gravitational constant, or having the acceleration driven by a scalar field
(see Thomas et al., 2022, for a general parameterisation of such options).

The dark energy and cosmological constant problems motivate a minimal
extension of Equation 1 to include a single extra scalar degree of freedom,
$\phi$, which is both physically and theoretically acceptable, i.e., not
allowing for negative energies for example, and can encapsulate one or more
cosmological constant problem solutions. Such an extension is found in the
well studied Horndeski (H) scalar-tensor theory (Horndeski, 1974). This is the
most general, Lorentz-covariant scalar-tensor theory in 4 spacetime dimensions
that yields second-order equations of motion, a basic condition for the
physical viability of the theory, i.e., it is ghost-free. A universe described
by Horndeski gravity is given as

$\displaystyle S_{U}$ $\displaystyle=S_{\rm H}+S_{\rm M}=\int d^{4}x\sqrt{-g}$
$\displaystyle\big{[}G_{2}(\phi,X)-G_{3}(\phi,X)\Box\phi+G_{4}(\phi,X)R$
$\displaystyle+G_{4,X}(\phi,X)[(\Box\phi)^{2}-(\nabla_{\mu}\nabla_{\nu}\phi)^{2}]$
$\displaystyle+G_{5}(\phi,X)G_{\mu\nu}\nabla^{\mu}\nabla^{\nu}\phi$
$\displaystyle-\frac{1}{6}G_{5,X}(\phi,X)[(\Box\phi)^{3}-3\Box\phi(\nabla_{\mu}\nabla_{\nu}\phi)^{2}+2(\nabla_{\mu}\nabla_{\nu}\phi)^{3}]\Big{]}$
$\displaystyle+S_{\rm M}\,,$ (2)

where each $G_{i}(\phi,X)$, $i=2,3,4,5$ is a free function of the scalar field
$\phi$ and its canonical kinetic term $X=-(\partial\phi)^{2}/2$, and
$G_{i,X}(\phi,X)=\partial G_{i}/\partial X$.

This opens up a very large theory space which needs to be trimmed down with
observational data. We have very strong data constraints at small spatial
scales, i.e., within the Solar System and at astrophysical scales (Will, 1993,
2014), showing gravity is highly consistent with GR in this regime. We also
have high quality observational data from cosmology, primarily from the CMB
which is associated with early cosmological times. This allows new theoretical
models most phenomenological freedom at large temporal and spatial scales as
they must recover CMB and solar system observations. The small spatial scale
constraints can be evaded using so called screening mechanisms (see Koyama,
2018; Burrage & Sakstein, 2018, for reviews) that force predictions of
modified gravity models back to those of GR locally, while early time
measurements like the CMB can easily be recovered through appropriate time
evolution of $\phi$.

An obvious late time cosmological data set directly related to gravity is the
large scale structure of the Universe (LSS). A key summary statistic of this
is the two point correlation function or power spectrum (in Fourier space) of
the cosmological matter field. A prime science goal then becomes the
production of accurate predictions of the matter power spectrum in general
theories beyond-$\Lambda$CDM. For the Horndeski class of models, this is a
nontrivial task as there are an additional four free functions of space and
time to contend with, beyond the matter content and metric freedoms. Of
course, one can always choose particular forms for the $G_{i}(\phi,X)$ and
then produce predictions for the 2-point correlations of matter. This approach
allows one to fully specify how matter should cluster at all physical scales,
and there are many tools and models that do just that to varying degrees of
accuracy (Schmidt et al., 2010; Lombriser, 2014; Arnold et al., 2022; Cataneo
et al., 2019; Bose et al., 2020b; Hernández-Aguayo et al., 2022; Puchwein et
al., 2013; Brax & Valageas, 2013, 2014; Joudaki et al., 2022; Winther et al.,
2017).

If on the other hand we choose not to specify a particular model, we are
required to parameterise both the linear and nonlinear scales i.e., the large
and small physical scales of LSS respectively. At linear scales and for the
Horndeski class of models, we can opt to perform a Taylor expansion of the
$G_{i}$ functions and truncate at some order. Linear theory can then be
characterised by a small number of free functions of time but with no unique
specification in the nonlinear regime. This describes the approach of the
Effective Field Theory of Dark Energy (EFTofDE) by Gubitosi et al. (2013);
Bloomfield et al. (2013) (also see Frusciante & Perenon, 2020, for a review).
Note that if we wish to be even more general than Horndeski we can directly
parametrise the linear relation between matter and the gravitational
potential.

On nonlinear scales, a parameterisation framework one can consider is the
nonlinear parameterised post-Friedmannian (nPPF), which captures modified
gravity or dark energy effects (Lombriser, 2016). Both linear and nonlinear
parameterisations then need to be consistently embedded in some more
comprehensive predictive framework in order to be able to confront theory with
LSS observations.

For past galaxy surveys the precision of the data did not call for high
accuracy in the power spectrum modelling, (as argued in Spurio Mancini et al.,
2019; Traykova et al., 2019). This changes with the next generation of surveys
(Stage-IV) such as Euclid333http://euclid-ec.org (Laureijs et al., 2011) and
the Vera C. Rubin Observatory’s Legacy Survey of Space and Time
(LSST)444https://www.lsst.org/ (Ivezić et al., 2019). These surveys will
provide a significant reduction in statistical errors, errors which will be
lowest in the nonlinear regime. With such precision, we have the opportunity
to greatly constrain deviations to $\Lambda$CDM, including the well defined
model space within Equation 2. This is contingent on whether or not we can
accurately and efficiently map these deviations to the matter power spectrum.
Typically, to remain unbiased in our constraints on Nature, $\mathcal{O}(1)\%$
is quoted as being the target accuracy for theoretical predictions (see
Blanchard et al., 2020a, for example). But this is not sufficient. We also
require this map to be computationally efficient enough to perform data
analyses. Without accuracy, we forfeit trust in our constraints. Without
conciseness and efficiency we face major computational issues.

This paper provides a balance that satisfies these criteria. We mainly focus
on the Horndeski class of models, embedding the EFTofDE and nPPF approaches
into the halo model reaction framework (Cataneo et al., 2019; Giblin et al.,
2019; Cataneo et al., 2020; Bose et al., 2020a, b; Carrilho et al., 2022),
which is able to predict the nonlinear power spectrum for specified theories
beyond-$\Lambda$CDM at $\mathcal{O}(1)\%$-level accuracy. We also present a
completely model independent parametrisation of beyond-$\Lambda$CDM physics at
nonlinear scales, which can be combined with similar parametrisations for the
Universe’s background expansion history and linear structure formation, giving
a parametrisation for general deviations to $\Lambda$CDM. Thus, this work
promotes the halo model reaction framework to being able to perform model
independent predictions, a key step in the search for a more fundamental
description of Nature in the cosmological, low energy regime.

This paper is organised as follows: in section 2 we begin at the observational
end and look how to model the (halo model) reaction. In section 3 we jump to
the theoretical end, and look how we can connect the ingredients of the
reaction to an action of Nature, together with any additional degrees of
freedom characterising nonlinear physics. In section 4 we give an overview of
the mapping between the reaction and the parameterisations of gravity and dark
energy, along with some key simplifying approximations one can consider. We
also perform tests and provide motivations for these approximations. Then, in
section 5 we test the proposed parameterisations by comparing to exact
predictions as well as emulators based on $N$-body simulations in an evolving
dark energy scenario and two representative non-standard models of gravity. In
section 6 we summarise our results and conclude. In the appendices we provide
full expressions for the linear and nonlinear parametrisations as well as
illustrative examples and comparisons within specific non-standard models of
gravity.

##  2 Halo model reaction

The leading order moment of the cosmological matter distribution is the
nonlinear matter power spectrum, $P_{\rm NL}(k,z)$. This Fourier space
quantity captures most of the matter clustering information at all scales (see
Bernardeau et al., 2002, for a review). Following the halo model (see Cooray &
Sheth, 2002, for a review) based approach of Cataneo et al. (2019), in a
target theory of cosmology and gravity this quantity can be modelled as

$P_{\rm NL}(k,z)=\mathcal{R}(k,z)P^{\rm pseudo}_{\rm NL}(k,z)\,,$ (3)

where $P^{\rm pseudo}_{\rm NL}(k,z)$ is called the pseudo power spectrum. This
is defined as the power spectrum of a $\Lambda$CDM universe but whose initial
conditions have been set so as to match the target, beyond-$\Lambda$CDM,
theory’s linear total matter power spectrum $P_{L}(k,z)$ at some target
redshift, $z$. The reason for making such a definition is that it guarantees
the halo mass functions in the target and pseudo universes are similar since
they will have the same linear clustering by definition. This results in a
smoother transition between the clustering statistics in the inter- and intra-
halo regimes. This quantity can be modelled in a number of ways, for example
by using existing halo model based fitting functions such as HMCode (Mead et
al., 2015, 2016, 2021) or for target theories that only predict a redshift
dependent, but scale independent rescaling of the linear spectrum,
$\Lambda$CDM-based emulators such as EuclidEmulator2 (Euclid Collaboration et
al., 2020) or bacco (Angulo et al., 2021) can be used by tuning the spectrum
amplitude parameter to match the modified cosmology’s linear spectrum.

The function $\mathcal{R}(k,z)$ represents all the corrections to the pseudo
spectrum coming from nonlinear beyond-$\Lambda$CDM physics. Following Cataneo
et al. (2020); Bose et al. (2021) we can write this as

$\mathcal{R}(k)=\frac{\left(1-f_{\nu}\right)^{2}P_{\mathrm{hm}}^{(\mathrm{cb})}(k)+2f_{\nu}\left(1-f_{\nu}\right)P_{\mathrm{hm}}^{(\mathrm{cb}\nu)}(k)+f_{\nu}^{2}P_{\mathrm{L}}^{(\nu)}(k)}{P_{\rm
hm}^{\rm pseudo}(k,z)}\,,$ (4)

with the subscript ‘hm’ standing for halo model, $\rm(m)\equiv(cb+\nu)$, cb
for CDM and baryons, $\nu$ for massive neutrinos and
$f_{\nu}=\Omega_{\nu,0}/\Omega_{m,0}$ being the massive neutrino energy
density fraction at $z=0$. The effects of massive neutrinos are included
linearly through the weighted sum of the nonlinear cb halo model and linear
massive neutrino spectra following the findings of Agarwal & Feldman (2011).
We note that we do not consider massive neutrino effects in this work, but
have included them in the expressions to highlight the generality of this
approach (see Bose et al., 2021, for a study with massive neutrinos).

The individual components are given by

$\displaystyle P_{\mathrm{hm}}^{(\mathrm{cb}\nu)}(k)\approx$
$\displaystyle\sqrt{P_{\mathrm{hm}}^{(\mathrm{cb})}(k)P_{\mathrm{L}}^{(\nu)}(k)}\,,$
(5) $\displaystyle P_{\mathrm{hm}}^{(\mathrm{cb})}(k)=$
$\displaystyle\left[(1-\mathcal{E})e^{-k/k_{\star}}+\mathcal{E}\right]P_{\mathrm{L}}^{(\mathrm{cb})}(k)+P_{\mathrm{1h}}^{(\mathrm{cb})}(k)\,,$
(6) $\displaystyle P_{\rm hm}^{\rm pseudo}(k,z)=$ $\displaystyle P_{\rm
L}(k,z)+P_{\rm 1h}^{\rm pseudo}(k,z),$ (7)

where the parameters are given by

$\displaystyle\mathcal{E}(z)=$ $\displaystyle\lim_{k\rightarrow
0}\frac{(1-f_{\nu})^{2}P_{\rm 1h}^{\rm(cb)}(k,z)}{P_{\rm 1h}^{\rm
pseudo}(k,z)},$ (8) $\displaystyle k_{\rm\star}(z)=$
$\displaystyle-\bar{k}\left(\ln\left[\frac{T_{1}(\bar{k},z)\pm
T_{2}(\bar{k},z)}{(1-f_{\nu})^{2}P^{\rm(cb)}_{\rm
L}(\bar{k},z)(1-\mathcal{E}(z))}\right]\right)^{-1}\,.$ (9)

We take the ‘+’ root if $\mathcal{E}>1$, otherwise we take the ‘-’ root. The
$T_{i}$ terms are given by

$\displaystyle T_{1}(k,z)=$ $\displaystyle
f_{\nu}^{2}P_{\mathrm{L}}^{(\nu)}(k,z)+P_{\rm hm}^{\rm
pseudo}(k,z)\mathcal{R}_{\rm SPT}(k,z)$
$\displaystyle-(1-f_{\nu})^{2}[\mathcal{E}(z)P^{\rm(cb)}_{\rm
L}(k,z)+P_{\mathrm{1h}}^{(\mathrm{cb})}(k,z)]\,,$ (10) $\displaystyle
T_{2}(k,z)=$ $\displaystyle 2\sqrt{f_{\nu}^{2}P_{\rm hm}^{\rm
pseudo}(k,z)P_{\mathrm{L}}^{(\nu)}(k,z)\mathcal{R}_{\rm SPT}(k,z)}\,,$ (11)

where $\mathcal{R}_{\rm SPT}(k,z)$ is the 1-loop standard perturbation theory
(SPT) (Bernardeau et al., 2002) prediction for the reaction given by Equation
4-7 but with the replacements $P_{L}(k,z)\rightarrow P_{\rm 1-loop}(k,z)$ and
$P_{L}(k,z)^{\rm(cb)}\rightarrow P_{\rm 1-loop}(k,z)^{\rm(cb)}$ and
$\mathcal{E}=1$. As in Cataneo et al. (2019) the default scale where we
calculate $k_{\rm\star}$ is set to $\bar{k}=0.06~{}h~{}{\rm Mpc}^{-1}$.

We see that Equation 4 depends on three general predictions for
beyond-$\Lambda$CDM theories: the 2-halo term which we have approximated by
the linear power spectrum $P_{\rm L}$, the quasi-nonlinear power spectrum
given by the 1-loop perturbation theory power spectrum $P_{\rm 1-loop}$, and
the highly nonlinear power spectrum given by the 1-halo term $P_{\rm 1h}$. The
computation of these quantities requires the specification of the matter
density fluctuations at different physical scales. The first two regimes
(linear and quasi-nonlinear) are perturbatively derived up to 3rd order in the
linear density fluctuation $\delta_{L}$, while the fully nonlinear quantity,
$\delta_{\rm NL}$, can be obtained using the assumptions of spherical collapse
(Cooray & Sheth, 2002). Both these routes require us to solve differential
equations representing energy and momentum conservation on a cosmological
background. Our Universe’s spacetime metric is well described by the Friedman-
Lemaître-Robertson-Walker (FLRW) metric, whose background expansion is
described by the Hubble parameter $H(a)\equiv\dot{a}/a$, where $a$ is the
scale factor and an over-dot represents a derivative with respect to the
metric time $t$.

Further, the conservation equations rely on the relation between the
gravitational potential and the matter density fluctuation: the Poisson
equation. In particular, we consider the Poisson equation in the perturbative
limit, only valid up to quasi-nonlinear scales, as well as the fully nonlinear
expression, valid at all scales

$\displaystyle-\left(\frac{k}{aH(a)}\right)^{2}\Phi_{\rm
QNL}(\mbox{\boldmath$k$},a)=$ $\displaystyle\frac{3\Omega_{\rm
m}(a)}{2}\mu(k,a)\,\delta_{\rm
QNL}(\mbox{\boldmath$k$},a)+S(\mbox{\boldmath$k$},a)\,,$ (12)
$\displaystyle-\left(\frac{k}{aH(a)}\right)^{2}\Phi_{\rm
NL}(\mbox{\boldmath$k$},a)=$ $\displaystyle\frac{3\Omega_{\rm
m}(a)}{2}[1+\mathcal{F}(k,a)]\,\delta_{\rm NL}(\mbox{\boldmath$k$},a)\,,$ (13)

where $\Omega_{\rm m}(a)=\Omega_{\rm{m},0}H_{0}^{2}/(H(a)^{2}a^{3})$,
$\Omega_{\rm{m},0}$ being the total matter fraction today. $\Phi$ is the
gravitational potential in the time-time component of the perturbed FLRW
metric. This can be identified with the Newtonian gravitational potential in
the non-relativistic limit, valid for the curvatures and velocities we
consider. The subscripts QNL and NL denote ‘quasi-nonlinear’ and ‘nonlinear’
respectively. One should further note that Equation 12 and Equation 13 also
assume a spherically symmetric density distribution.

The additional functions in Equation 12 and Equation 13 are as follows:
$\mu(k,a)$ characterises the linear modification to GR, $\mathcal{F}(k,a)$ is
the nonlinear modification and $S(\mbox{\boldmath$k$},a)$ is a source term
capturing modifications at 2nd and 3rd order in the linear matter density
perturbations. The source term is given by (Bose & Koyama, 2016)

$\displaystyle S(\mbox{\boldmath$k$},a)$
$\displaystyle=\int\frac{d^{3}\mbox{\boldmath$k$}_{1}d^{3}\mbox{\boldmath$k$}_{2}}{(2\pi)^{3}}\,\delta_{\rm
D}(\mbox{\boldmath$k$}-\mbox{\boldmath$k$}_{12})\gamma_{2}(\mbox{\boldmath$k$}_{1},\mbox{\boldmath$k$}_{2},a)\delta(\mbox{\boldmath$k$}_{1})\,\delta(\mbox{\boldmath$k$}_{2})$
$\displaystyle+\int\frac{d^{3}\mbox{\boldmath$k$}_{1}d^{3}\mbox{\boldmath$k$}_{2}d^{3}\mbox{\boldmath$k$}_{3}}{(2\pi)^{6}}\delta_{\rm
D}(\mbox{\boldmath$k$}-\mbox{\boldmath$k$}_{123})\gamma_{3}(\mbox{\boldmath$k$}_{1},\mbox{\boldmath$k$}_{2},\mbox{\boldmath$k$}_{3},a)$
$\displaystyle\quad\quad\times\delta(\mbox{\boldmath$k$}_{1})\,\delta(\mbox{\boldmath$k$}_{2})\,\delta(\mbox{\boldmath$k$}_{3})\,,$
(14)

introducing two additional functions $\gamma_{2}$ & $\gamma_{3}$
characterising quasi-nonlinear modifications to the Poisson equation (see Bose
& Koyama, 2016, for explicit expressions for these in the Horndeski class of
models). The functions $\gamma_{2}$, $\gamma_{3}$ and $\mathcal{F}$ all encode
details regarding the potential screening mechanism of the theory under
consideration. On this point, it is worth noting that for general theories
beyond-$\Lambda$CDM such mechanisms may not be present, in which case the
spherical density distribution approximation assumed in Equation 12 & Equation
13 may break down (Thomas, 2020). For the modified gravity models considered
in this work, which have some method of screening, this appears to be a
reasonable approximation (Noller et al., 2014). For a study of screened and
unscreened models in the Horndeski class see Noller et al. (2021).

In total, the halo model reaction, and so the nonlinear power spectrum,
requires specification of four functions of space and time - one for the
background $H(a)$, one for the linear regime $\mu(k,a)$, two for the quasi-
nonlinear regime
$\gamma_{2}(\mbox{\boldmath$k$}_{1},\mbox{\boldmath$k$}_{2},a)$ &
$\gamma_{3}(\mbox{\boldmath$k$}_{1},\mbox{\boldmath$k$}_{2},\mbox{\boldmath$k$}_{3},a)$
and finally one for the fully nonlinear regime $\mathcal{F}(k,a)$. In
principle these functions are not completely independent, and one should have
$\mathcal{F}\rightarrow\mu$ in the linear limit. We investigate the importance
of respecting this limit in section 5. Finally, we remind the reader that all
these functions are required to compute the key ingredients of $\mathcal{R}$
(and hence $P_{\rm NL}$): $P_{\rm L}$, $P_{\rm 1-loop}$ and $P_{\rm 1h}$.

The right half of Figure 1 summarises the map from background and Poisson
equations to the halo model reaction as described in this section. The left
half of the figure will be the focus of the next section.

pt$H(a)$pt$\alpha_{M}(a),\alpha_{B}(a),$$\alpha_{K}(a),\alpha_{T}(a)$pt$\xi_{3}^{2}(a),\xi_{4}^{2}(a),\xi_{3}^{3}(a),$$\xi_{4}^{3}(a),\xi_{4}^{4}(a),\xi_{4}^{5}(a)$pt$G_{2}(\phi,X)$$G_{3}(\phi,X)$$G_{4}(\phi,X)$$G_{5}(\phi,X)$pt$q_{i}$
or
$p_{i}(a)$$H(a)$$\mu(k,a)$$\gamma_{2}(\mbox{\boldmath$k$}_{1},\mbox{\boldmath$k$}_{2},a)$$\gamma_{3}(\mbox{\boldmath$k$}_{1},\mbox{\boldmath$k$}_{2},\mbox{\boldmath$k$}_{3},a)$$\mathcal{F}(k,a)$pt$P_{\rm
L}(k,a)$pt$P_{\rm 1-loop}(k,a)$pt$P_{\rm
1h}(k,a)$pt$\mathcal{R}(k,a)$ptKennedy et al (2017)ptLombriser (2016)ptBose et
al (2016)ptCataneo et al (2019)ptAptCptTptIptOptRptEptAptCptTptIptO Figure 1:
A rough schematic of the map from the Horndeski action in the EFTofDE
parametrisation ($H,\alpha_{i},\xi_{i}^{j}$) and nonlinear parameterisations
($p_{i}$ or $q_{i}$) to $\mathcal{R}(k,t)$. The yellow rectangles indicate the
input functions of time (here parametrised by the scale factor $a$) or
constants. The orange rectangles indicate the modifications to the Poisson
equation. $\mu(k,a)$ is bi-coloured indicating we may choose to parametrise it
directly instead of starting at the action level. We provide the main sources
in the literature for each piece of the map along with a dotted box roughly
indicating their associated piece. Note that the solid arrows can only
reconstruct $G_{i}$ to the linear and quasi-nonlinear levels, which can in
turn inform choices for $p_{i}$. The $p_{i}$ provide the nonlinear complement
in the $G_{i}$.

##  3 Parametrisations

We now move away from the observational end and return to the starting point,
the fundamental action of Nature. In particular, here we mostly focus on the
Horndeski action given in Equation 2, but the approach can be trivially
extended to further generality.

As pointed out, given a specific form of the $G_{i}$ functions, the explicit
functional forms of $H$, $\mu$, $\gamma_{2}$, $\gamma_{3}$ and $\mathcal{F}$
can be directly derived. But rather than specifying the full covariant theory,
i.e., 4 free functions of space and time, we ultimately wish to parameterise
the action’s predictions for cosmological matter clustering in terms of a few
free constants.

To do this, we split LSS into three regimes: the background & linear, quasi-
nonlinear and the nonlinear. The background, linear and quasi-nonlinear
regimes will follow the well studied EFTofDE program (Gubitosi et al., 2013;
Bloomfield et al., 2013). For the nonlinear regime we will consider two
different parameterisations. One is the established nonlinear parameterised
post-Friedmannian (nPPF) approach (Lombriser, 2016). The other parametrisation
we propose here is phenomenological and is based on some well known screening
mechanisms. We begin by parameterising the background & linear regime.

###  3.1 Background & Linear: Effective field theory of dark energy

Among the methods to generically parameterise beyond-$\Lambda$CDM physics on
cosmological scales, the methods of Effective Field Theory (EFT) have proven
to be particularly useful. It is simply necessary to determine which
symmetries one wishes the action to have before constructing various operators
out of the fields and derivatives of the fields. One can trust the predictions
made with an EFT as long as it is made at an energy scale below the ‘cutoff’
of the theory, beyond which the validity of the EFT breaks down.

While not being an EFT in this strict sense, the EFTofDE is constructed in a
similar manner and is capable of describing the dynamics of the cosmological
background and perturbations in Horndeski theory in a generic manner. The
EFTofDE approach breaks time diffeomorphism invariance of the cosmological
background by choosing a particular gauge. By doing this one is able to form a
theory out of operators which only respect spatial diffeomorphism invariance.

In constructing the EFTofDE action one begins by foliating spacetime with
constant-time hypersurfaces. Utilising the complete freedom one has in
choosing the coordinates of the theory we can set the scalar field to be only
a function of time such that $\phi(x,t)\rightarrow\phi(t)$. In particular, we
can choose

$\phi=t/\kappa^{2}\,.$ (15)

This choice is called the unitary gauge and in this gauge the scalar field
perturbations vanish, being absorbed into the time-time component of the
metric. The operators in the EFTofDE are the cosmological perturbations
themselves. In the unitary gauge we are free to include operators in the EFT
which are only spatially diffeomorphism invariant, such as $g^{00}$.

Let us denote the normal vector to each spatial hypersurface as

$n_{\mu}=-\frac{\partial_{\mu}\phi}{\sqrt{-(\partial\phi)^{2}}}\,.$ (16)

The induced spatial metric of each hypersurface is then given by
$h_{\mu\nu}=g_{\mu\nu}+n_{\mu}n_{\nu}$. This allows us to include the
extrinsic curvature which is given by the projection of the derivative of the
normal vector along the the hypersurface, onto the hypersurface
$K_{\mu\nu}=h_{\mu\sigma}\nabla^{\sigma}n_{\nu}$. With the induced metric, one
can also compute the intrinsic curvature of each hypersurface given by the
three-dimensional Ricci scalar $R^{(3)}$.

Collecting relevant combinations of the invariants under residual spatial
diffeomorphism symmetry gives the EFTofDE action, which is capable of
describing the dynamics of the background and linear perturbations of
Horndeski theory. The action is given by (see, for example, Kennedy et al.,
2017)

$\displaystyle S_{U,L}=$
$\displaystyle\>S^{(0,1)}+S^{(2)}+S_{M}[g_{\mu\nu},\Psi_{m}]\,,$ (17)
$\displaystyle S^{(0,1)}=$ $\displaystyle\>\int
d^{4}x\sqrt{-g}\left[\frac{\Omega(t)}{2\kappa^{2}}R-\Lambda(t)-c(t)\delta
g^{00}\right]\,,$ (18) $\displaystyle S^{(2)}=$ $\displaystyle\int
d^{4}x\sqrt{-g}\Big{[}\frac{M^{4}_{2}(t)}{2}(\delta
g^{00})^{2}-\frac{\bar{M}^{3}_{1}(t)}{2}\delta K\delta g^{00}$
$\displaystyle-\bar{M}_{2}^{2}(t)\Big{(}\delta K^{2}-\delta K^{\mu\nu}\delta
K_{\mu\nu}-\frac{1}{2}\delta R^{(3)}\delta g^{00}\Big{)}\Big{]}\,.$ (19)

where $S_{U,L}$ represents the action of a Horndeski-universe that describes
field dynamics up to the linear level in the matter and velocity
perturbations. The $(0,1,2)$ represent the order in the perturbed quantities.

In front of each term we include a free function of time called an EFT
coefficient, giving a total of six free functions,
$\\{\Omega(t),\Lambda(t),c(t),M_{2}^{4}(t),\bar{M}_{1}^{3}(t),\bar{M}_{2}^{2}(t)\\}$.
Once we specify a metric, we also introduce any metric degrees of freedom. For
FLRW this is the scale factor $a$, or equivalently the Hubble parameter
$H(a)$. We can then employ the field equation constraints, which in the FLRW
are the Friedmann equations

$\displaystyle 0=$
$\displaystyle\kappa^{2}(2c-\Lambda+\rho_{\mathrm{m}})-3H^{2}(\Omega+a\Omega^{\prime})\,,$
(20) $\displaystyle 0=$
$\displaystyle\kappa^{2}\Lambda+H[aH^{\prime}(2\Omega+a\Omega^{\prime})+H(3\Omega+3a\Omega^{\prime}+a^{2}\Omega^{\prime\prime})]\,,$
(21)

where we have dropped the time dependence in constituent parameters for
compactness, and use the scale factor to parameterise time. A prime denotes a
scale factor derivative and $\rho_{\mathrm{m}}$ is the matter density at $a$.
The Friedmann equations reduce the number of free functions describing the
background and linear perturbations to five. Solving these equations yields

$\displaystyle c(a)$
$\displaystyle=-\frac{\rho_{\mathrm{m}}}{2}-\frac{aH[H^{\prime}(2\Omega+a\Omega^{\prime})+aH\Omega^{\prime\prime}]}{\kappa^{2}}\,,$
(22) $\displaystyle\Lambda(a)=$
$\displaystyle-\frac{H[aH^{\prime}(2\Omega+a\Omega^{\prime})+H(3\Omega+3a\Omega^{\prime}+a^{2}\Omega^{\prime\prime})]}{\kappa^{2}}\,.$
(23)

This means the free functions of the scale factor defining the background and
linear theory would be
$\\{\Omega,H,M_{2}^{4},\bar{M}_{1}^{3},\bar{M}_{2}^{2}\\}$, which we will
refer to as the $M$-basis. We can alternatively write the Hubble function as
the solution to

$H(a):\,\,0=(2\Omega+a\Omega^{\prime})H^{\prime}+aH\Omega^{\prime\prime}+\frac{\kappa\rho_{\mathrm{m}}}{aH}+\frac{2\kappa
c}{aH}\,,$ (24)

if we wish to specify $c$ instead of $H$ for example.

Common in the literature is the $\alpha$-basis
$\\{H,\alpha_{M},\alpha_{B},\alpha_{K},\alpha_{T}\\}$ which has a clearer
physical interpretation of the effects of each function (see, for example,
Bellini & Sawicki, 2014). We provide the map between the $\alpha$\- and
$M$-bases555Note the factor of ‘$-1/2$’ difference in $\alpha_{B}$ between our
expression and that of EFTCAMB (Frusciante et al., 2016) or Kennedy et al.
(2018), for instance.

$\displaystyle\alpha_{M}$ $\displaystyle=\frac{a(M^{2})^{\prime}}{M^{2}}\,,$
(25) $\displaystyle\alpha_{B}$
$\displaystyle=-\frac{aH\Omega^{\prime}+\kappa^{2}\bar{M}^{3}_{1}}{HM^{2}\kappa^{2}}\,,$
(26) $\displaystyle\alpha_{K}$
$\displaystyle=\frac{2c+4\bar{M}_{2}^{4}}{M^{2}H^{2}}\,,$ (27)
$\displaystyle\alpha_{T}$ $\displaystyle=-\frac{\bar{M}_{2}^{2}}{M^{2}}\,,$
(28)

where $M^{2}=\Omega~{}\kappa^{-2}+\bar{M}_{2}^{2}$. Note that one can
alternatively specify $M^{2}$ and solve for $H$.

To end this section, another basis worth considering is the basis introduced
in Kennedy et al. (2018): $\\{H,M^{2},c_{s}^{2},\alpha,\alpha_{B0}\\}$ (also
see Lombriser et al., 2019), which implicitly assumes $\alpha_{T}=0$ (see
subsection 4.3 for motivation). This basis allows for some simple priors on
the functions that ensure the theory has no ghost or gradient instabilities,
i.e., negative energies or imaginary sound speeds. We will refer to this basis
as the $s$-basis. The priors to ensure stability on these functions are then
simply $M^{2},c_{s}^{2},\alpha>0$, and $\alpha_{B0}$ is constant666We note
that this basis does not ensure the absence of a tachyonic instability
(Gsponer & Noller, 2022).. The map between the $s$\- and $\alpha$-bases is
given by

$\displaystyle c_{s}^{2}=$
$\displaystyle\frac{2}{\alpha}\Big{[}\frac{a\alpha_{B}^{\prime}}{2}-(1+\alpha_{T})\left(1-\frac{\alpha_{B}}{2}\right)^{2}$
$\displaystyle+\left(1+\alpha_{M}-\frac{aH^{\prime}}{H}\right)\left(1-\frac{\alpha_{B}}{2}\right)-\frac{\rho_{\mathrm{m}}}{2H^{2}M^{2}}\Big{]}\,,$
(29) $\displaystyle\alpha=$
$\displaystyle\alpha_{K}+\frac{3}{2}\alpha_{B}^{2}\,,$ (30)

where $c_{s}$ is the speed of sound, while $\alpha_{B0}=\alpha_{B}(a=1)$ is
the boundary condition ($\alpha_{B}$’s value today) specified to solve the
differential equation given by Equation 29.

In what follows we will stick with the $\alpha$-basis and implement this as
the default basis in the accompanying code ACTio et ReACTio. We provide the
explicit form of the linear modification to the Poisson equation in this basis
in subsection A.1. We leave it to the user to perform the transformation from
their preferred basis to the $\alpha$-basis, and provide an accompanying
notebook GtoPT.nb that performs some of these transformations.

###  3.2 Quasi-nonlinear: Covariant theory map

To fully specify the halo model reaction, we need to go beyond the linear
matter perturbations. In particular, we also require the 2nd and 3rd order
density perturbations to solve for the 1-loop power spectrum entering
$\mathcal{R}_{\rm SPT}$ in section 2-11. This requires us to expand to fourth
order in the metric perturbation $\delta g^{00}$ and extrinsic curvature
$\delta K^{\mu\nu}$ in Equation 17. This has been done in Cusin et al. (2018b)
and has been used to calculate the 1-loop spectrum in Cusin et al. (2018a).
Further, in Kennedy et al. (2017) the authors relate the EFTofDE functions up
to a given order to the corresponding covariant theory’s Lagrangian $G_{i}$
functions as

$G_{i}(\phi,X)=g_{i}(\phi,X)+\Delta G_{i}(\phi,X)\,,$ (31)

where $g_{i}$, $i\in\\{2,3,4,5\\}$, are well-defined functions of $\phi$, $X$
and the lower order EFTofDE parameters, e.g.,
$\\{H,\alpha_{M},\alpha_{B},\alpha_{K},\alpha_{T}\\}$. The other terms are
given as

$\displaystyle\Delta G_{2,3}=$
$\displaystyle\sum_{n>2}\xi_{n}^{(2,3)}(\phi)\Big{(}1+X\kappa^{4}\Big{)}^{n}\,,$
(32) $\displaystyle\Delta G_{4,5}=$
$\displaystyle\sum_{n>3}\xi_{n}^{(4,5)}(\phi)\Big{(}1+X\kappa^{4}\Big{)}^{n}\,,$
(33)

where $\Delta G_{i}$ are higher order corrections to the covariant action and
$\xi_{n}^{i}(\phi)$ are higher order EFTofDE functions, $X$ again being the
scalar field canonical kinetic energy term.

A particular covariant theory is specified once $\xi_{n}^{i}$ are given for
all $n\in\mathbb{N}$, but if we truncate at some order $n_{t}$, we specify the
subset of Horndeski theories which are identical on scales described by the
EFTofDE up to $\xi_{n_{t}}^{i}$. Up to 3rd order in the matter density
perturbation, we introduce 6 new functions with $n_{t}=4$. Together with the
background and linear order functions, this gives a total of 11 free functions
of time for the quasi-nonlinear scales. The $G_{i}$ given in Equation 31 can
then be related to $\mu$, $\gamma_{2}$ and $\gamma_{3}$ by the map provided in
the Appendices of Bose & Koyama (2016); Takushima et al. (2015).

In Appendix A we provide the map between the 5 linear EFTofDE functions in the
$\alpha$-basis and the linear modification to the Poisson equation, $\mu$,
used in Equation 12. The 2nd and 3rd order functions $\gamma_{2}$ and
$\gamma_{3}$ (see Equation 14) are significantly more complicated but can be
derived by using the map from the EFTofDE to $G_{i}(\phi,X)$ provided in
Kennedy et al. (2017) and then the $G_{i}(\phi,X)$ to $\gamma_{2}$ &
$\gamma_{3}$ given in Bose & Koyama (2016). The map, although not reproduced
here in full, is given in detail in a Mathematica notebook provided in the
ACTio et ReACTio repository, GtoPT.nb. This being said, in section 4 we give
support for the omission of $\gamma_{2}$ and $\gamma_{3}$ in the calculation
of $\mathcal{R}$ for moderate to low modifications to gravity, and given the
additional degrees of freedom we will introduce in the nonlinear regime.

Having specified a route between the Horndeski action of nature and the linear
and 1-loop power spectra, $P_{\rm L}(k,a)$ & $P_{\rm 1-loop}(k,a)$, we now
look at two methods of parameterising clustering in the highly nonlinear
regime, characterised by the 1-halo term, $P_{\rm 1h}(k,a)$. This will then
specify a full parameterisation of the halo model reaction $\mathcal{R}$, and
consequently the nonlinear power spectrum, $P_{\rm NL}(k,a)$.

###  3.3 Nonlinear

The effects of modified gravity on the nonlinear cosmic structure formation
are captured by the effective deviation $\mathcal{F}$ from the gravitational
constant in the nonlinear Poisson equation given in Equation 13 and the
cosmological background evolution. Specifically, the modified Poisson equation
alters the evolution equation for the halo top-hat radius $R_{\rm TH}$ (see,
for example, Schmidt et al., 2009). This quantity gives an estimate for
$\delta_{\rm NL}$, needed to compute the 1-halo power spectrum. Here we
discuss two parameterisations of $\mathcal{F}$.

####  3.3.1 Nonlinear parametrised post-Friedmannian framework

Following the nPPF approach of Lombriser (2016), the effective gravitational
coupling for generic screening mechanisms and other suppression effects can be
decomposed as a function of scale $r$

$1+\mathcal{F}(a,r)=A+\sum_{i}^{N_{0}}B_{i}\prod_{j}^{N_{i}}\mathcal{F}_{ij}\,,$
(34)

where $\mathcal{F}_{ij}$ are some transition functions encapsulating screening
or other suppression effects such as a Yukawa suppression. $N_{0}$ and $N_{i}$
characterise their respective number. In the fully screened limit, the
effective coupling reduces to $A$, typically unity, whereas it becomes $B_{i}$
in the fully unscreened limit, matching linear theory. To parameterise these
transitions, Lombriser (2016) adopted a generalised form of the Vainshtein
screening effect in the DGP braneworld model (Dvali et al., 2000)

$\mathcal{F}\sim b\left(\frac{r}{r_{\rm
scr}}\right)^{a_{f}}\left\\{\left[1+\left(\frac{r_{\rm
scr}}{r}\right)^{a_{f}}\right]^{1/b}-1\right\\}\,,$ (35)

where $r_{\rm scr}$ denotes the screening scale, which in general can be time,
mass, and environment dependent. The parameter $a_{f}$ (not to be confused
with the scale factor) determines the radial dependence of the coupling in the
screening limit along with $b$ that characterises an interpolation rate
between the screened and unscreened limits.

Screening effects such as the chameleon (Li & Efstathiou, 2012; Khoury &
Weltman, 2004; Lombriser et al., 2014) symmmetron (Hinterbichler & Khoury,
2010; Taddei et al., 2014), k-mouflage (Babichev et al., 2009; Brax &
Valageas, 2014), and Vainshtein (Vainshtein, 1972; Schmidt et al., 2010; Dvali
et al., 2000) mechanisms as well as other suppression effects such as the
linear shielding mechanism (Lombriser & Taylor, 2015b) or Yukawa suppression,
can be analytically mapped onto Equation 35 by matching expressions in the
limits of large and small $r$ and $r\rightarrow r_{\rm scr}$. The relevant
expressions may be found in Lombriser (2016). It is worth highlighting that
the parameters of Equation 35 for a given screening model may in principle be
directly read off from Equation 2 by employing the scaling method of McManus
et al. (2016) (also see Renevey et al., 2020) and counting the powers of
second and first spatial derivatives and the scalar field potential. Note that
the parameter $b$ may be understood as the choice of transition template used
to approximately cast the screening effect into. Alternatively to Equation 35,
one could also adopt other transition functions such as a hyperbolic tangent,
a sigmoid or an error function as we will propose in subsubsection 3.3.2. For
DGP, the choice of Equation 35 with $b=2$ becomes exact.

To implement Equation 35 in the spherical collapse model, one replaces
$r/r_{\rm scr}\rightarrow y/y_{\rm scr}$, where $y$ is the normalised top-hat
radius (Equation 103). A single general element $N_{0}=N_{1}=1$ can then be
described by seven parameters (or functions) $p_{1-7}$ in addition to
$p_{0}=A$ (typically $=1$). The first three, $p_{1-3}$, determine $a_{f}$,
$b$, and $B$. The other four are used to generally capture possible time,
mass, and environmental dependencies of the dimensionless screening scale,
which can be modelled as (Lombriser, 2016)

$y_{\rm scr}=p_{4}a^{p_{5}}\left(2G_{N}\,H_{0}M_{\rm
vir}\right)^{p_{6}}\left(\frac{y_{\rm env}}{y_{\rm h}}\right)^{p_{7}}\,,$ (36)

where $y_{\rm h}$ and $y_{\rm env}$ refer to the normalised radii of the the
halo and the environment respectively, $H_{0}$ is the Hubble constant and
$M_{\rm vir}$ is the virial mass of the halo777Note that in ReACT we use the
initial comoving top-hat radius, $R_{\rm th}$ (see subsection B.4), as an
input parameter instead of mass, related as
$M_{\mathrm{vir}}=4\pi\bar{\rho}_{\mathrm{m},i}(1+\delta_{i})(a_{i}R_{\rm
th})^{3}/3\approx 4\pi\Omega_{\mathrm{m,0}}\rho_{\mathrm{crit}}R_{\rm
th}^{3}/3$ with the critical density $\rho_{\rm crit}$ and
$1+\delta_{i}\approx 1$.. In this way, we can simplify Equation 35 to
(Lombriser, 2016)

$\mathcal{F}_{\rm
nPPF}=p_{1}p_{2}\frac{(1+s^{a_{f}})^{\frac{1}{p_{1}}}-1}{s^{a_{f}}}\,,$ (37)

where

$a_{f}=\frac{p_{1}}{p_{1}-1}p_{3}\,$ (38)

and $s=y_{\rm scr}/y_{\rm h}$. Note we have set $p_{0}=1$. The parameters
$p_{1-7}$ can be computed from theory and in many cases take on trivial values
(see subsection B.2). It is worth highlighting here that the nPPF formalism
has also been implemented in $N$-body simulations and cast into Fourier space
(Hassani & Lombriser, 2020), where it was shown to accurately match
simulations of exact model implementations.

Finally, we consider the large, linear scale limit of $\mathcal{F}$. Equation
39 provides a parametrised function for the screening regime, where we have a
transition to GR from some large scale modification. In this form, it does not
capture any additional effects coming from say Yukawa suppression, typical of
chameleon theories. Such phenomena may become relevant for the spherical
collapse calculation at early times or for very large halo masses. In order to
correctly capture this, we could either model the Yukawa suppression as
another transition cast into Equation 37 or simply augment Equation 37 with
the linear modification $\mu(k,a)$ as

$\mathcal{F}_{\rm
nPPF}=p_{1}p_{2}\frac{\left(1+s^{a_{f}}\right)^{\frac{1}{p_{1}}}-1}{s^{a_{f}}}\times(1-\mu(\hat{k},a))\,.$
(39)

In this case we also need to perform the Fourier transform of
$\mu(\hat{k},a)$, which is non-trivial. As a first order approximation, we
parametrise this with a simple scaling of the inverse of the comoving initial
top-hat radius $R_{\rm th}$ as

$\hat{k}=\frac{10^{p_{8}}}{a^{2}y_{\rm h}R_{\rm th}}\,,$ (40)

where the dimensionless constant $p_{8}$ calibrates the Yukawa suppression.
The Fourier transform can be made more sophisticated (see, for example,
Hassani & Lombriser, 2020) but in section 5 we find the impact of Yukawa
suppression is negligible for the $f(R)$ models we consider, and so only
include this augmentation for completeness. Further, Equation 39 would only be
meaningful for a non-trivial scale dependent $\mu(k,a)$. For scale-independent
theories one can absorb the scaling provided by $\mu(a)$ in the $p_{2}(a)$
parameter of Equation 37.

####  3.3.2 Phenomenological parameterisation

With its full freedom, the nPPF parameterisation is a very flexible way of
modelling the nonlinear scales. It is able to capture various specific
covariant theories exactly or to high accuracy (see subsection B.2 and section
5), and given a covariant theory, say from the Horndeski class, we can map its
nonlinear Poisson modification to the $p_{i}$ parameters. On the other hand,
if we remain agnostic about the covariant theory, 8 additional parameters,
some of which may also be time dependent, poses computational issues as well
as degrades the amount of cosmological and gravitational information we can
extract due to degeneracies between these nuisance and the physical parameters
of interest.

With this in mind, we propose the following general and reduced
parameterisation of $\mathcal{F}$ based on the error function (Erf). We have
found this mimics the general profile of the effective gravitational constant
in various modified gravity theories. Essentially we wish to capture a basic
transition from unscreened to screened regimes. The simple form we adopt is
given by

$\mathcal{F}_{\rm Erf}={\rm Erf}[ay_{\rm
h}10^{\bar{J}}]\times(1-\mu(\hat{k},a))\,,$ (41)

where as in the nPPF case, we use

$\hat{k}=\frac{10^{q_{4}}}{a^{2}y_{\rm h}R_{\rm th}}\,,$ (42)

and

$\bar{J}=q_{1}-q_{2}\log(R_{\rm th})+q_{3}\log(ay_{\rm env})\,.$ (43)

$\mu$ is the linear modification to gravity. In the EFTofDE parameterisation
$\mu$ is given in Equation 60, but this can also be parametrised more
generally (see, for example, Silvestri et al., 2013; Kennedy et al., 2018;
Srinivasan et al., 2021).

The Erf model introduces 4 free constants:

  * $\mathbf{q_{1}}$:

This parametrises the screening scale and goes as its inverse.

  * $\mathbf{q_{2}}$:

This gives the halo mass dependency of the screening scale.

  * $\mathbf{q_{3}}$:

This gives the environment dependency of the screening scale.

  * $\mathbf{q_{4}}$:

This calibrates any existing Yukawa suppression scale.

The time dependence of $\mathcal{F}_{\rm Erf}$ is fixed and so for a specified
cosmology and set of EFTofDE parameters, we only need to adjust the constants
$\\{q_{1},q_{2},q_{3},q_{4}\\}$. To provide some insight, we note the
following limits

$\displaystyle\lim_{q_{1}\rightarrow\infty}1+\mathcal{F}_{\rm Erf}$
$\displaystyle=\mu\rightarrow{\rm\bf Unscreened\,limit}\,,$ (44)
$\displaystyle\lim_{q_{1}\rightarrow-\infty}1+\mathcal{F}_{\rm Erf}$
$\displaystyle=1\rightarrow{\rm\bf GR\,limit}\,,$ (45)
$\displaystyle\lim_{q_{2},q_{3}\rightarrow 0}1+\mathcal{F}_{\rm Erf}$
$\displaystyle\rightarrow{\rm\bf Vainshtein\,type\,models}\,,$ (46)
$\displaystyle\lim_{q_{3}\rightarrow 0}1+\mathcal{F}_{\rm Erf}$
$\displaystyle\rightarrow{\rm\bf k-mouflage\,type\,models}\,,$ (47)
$\displaystyle q_{3}>0:\,1+\mathcal{F}_{\rm Erf}$
$\displaystyle\rightarrow{\rm\bf chameleon\,type\,models}\,,$ (48)

where we refer to the main types of screening mechanisms typical of scalar-
tensor theories (see subsubsection 3.3.1). Note that all parameters lose their
meaning as $\mu(k,a)\rightarrow 1$, which in the EFTofDE case is when the
relevant parameters assume their GR values.

Given this, we can take $q_{2}$ and $q_{3}$ to be positive. Being exponents of
the top-hat radius and environment parameter, they are also not expected to be
very large, and as we will see in section 5, they turn out to be
$\mathcal{O}(1)$. Further, since in the GR limit $\mu\rightarrow 1$, and so
$\mathcal{F}_{\rm Erf}\rightarrow 0$ irrespective of the value of $q_{1}$, we
can also take $q_{1}$ to be positive. We also find $q_{1}$ to be an
$\mathcal{O}(1)$ parameter.

Parameter $q_{4}$, which calibrates the Yukawa suppression scale, is generally
only relevant for theories where the linear growth factor, or Poisson
modification $\mu$, is scale-dependent. As we will show in subsection 5.3,
$q_{4}$ does not appear to be relevant for the scales associated with
spherical collapse. We note $q_{4}$ can in principle take on negative values,
pushing the Yukawa suppression to smaller scales. As $q_{4}\rightarrow\infty$
the Yukawa suppression scale also goes to infinity. We leave its relevance for
more general theories for a future work.

We provide a Mathematica notebook, Nonlinear.nb, with all the forms of
$\mathcal{F}$ considered in this paper along with comparisons.

Finally, the left half of Figure 1 summarises the map from the parametrised
action, together with additional parameters, to the Poisson equation
modifications as described in this section, completing the map from action to
reaction.

##  4 Approximations and Overview

We have outlined a map that goes from the parameterised action of nature and
structure formation $\\{H\\}_{b}$,
$\\{\alpha_{M},\alpha_{B},\alpha_{K},\alpha_{T}\\}_{\rm L}$,
$\\{\xi_{3}^{2},\xi_{4}^{2},\xi_{3}^{3},\xi_{4}^{3},\xi_{4}^{4},\xi_{4}^{5}\\}_{\rm
QNL}$ & $\\{p_{1},p_{2},p_{3},p_{4},p_{5},p_{6},p_{7},p_{8}\\}_{\rm NL}$ or
$\\{q_{1},q_{2},q_{3},q_{4}\\}_{\rm NL}$ to the nonlinear effects on the power
spectrum $\mathcal{R}(k,a)$, where ‘b’ stands for background, ‘L’ for linear,
‘QNL’ for quasi-nonlinear and ‘NL’ for nonlinear. A schematic of this map is
given in Figure 1. An important point worth stressing is that our nonlinear
parametrisations are completely general, and not specific to the Horndeski
class of theories. They do however rely on $\mu(k,a)$, which one can always
choose to parametrise in a model independent way.

Considering the Horndeski class for concreteness, the EFTofDE and nonlinear
parametrisations constitute a very large set of arbitrary functions of time
and constants. Despite it being significantly less than the infinite number of
theories contained within the Horndeski class, it is still arguably too many
for statistical data analyses, both on computational and scientific grounds.
Thankfully, as we will shortly motivate, these sets can be yet reduced
significantly.

To reduce or optimise the parameter space, we consider the following:

  1. (1)

We assume the quasi-static approximation (QS) for all perturbative
calculations (see Sawicki & Bellini, 2015; Pace et al., 2021, for example).

  2. (2)

We assume $\gamma_{2}=\gamma_{3}=0$.

  3. (3)

Observational and theoretical constraints.

  4. (4)

Time parameterisations of EFTofDE functions, $\alpha_{i}(a)$.

  5. (5)

The parameterised nPPF (see Equation 39) or phenomenological (see Equation 41)
form of $\mathcal{F}$ is flexible enough to capture general modifications to
gravity.

In this section we will motivate approximations (1) \- (4) with direct
reference to the accompanying code ACTio et ReACTio. Assumption (5) will be
addressed separately in section 5.

###  4.1 Quasi-static approximation

We begin by noting that the QS in linear theory can be easily avoided by using
a Boltzmann code such as EFTCAMB (Hu et al., 2014; Raveri et al., 2014) to
calculate the linear input spectrum or transfer function888The QS can be
partly circumvented in the nonlinear regime, Equation 39 and Equation 41, by
also using the prediction of $\mu(k,a)$ taken from say EFTCAMB. . This option
is available in our code, but the default setting assumes a $\Lambda$CDM
linear spectrum or transfer function at $z=0$ and rescales it using the
internally calculated growth functions of the desired theory. This is done
using the linear form of Equation 12 (see Equation 60) which assumes the
quasi-static approximation. Being able to use a $\Lambda$CDM linear spectrum
enhances the computational efficiency of our code as it avoids a call to
EFTCAMB. EFTCAMB is significantly slower than CAMB (Lewis et al., 2000), which
already takes $\mathcal{O}(1)$ seconds to produce a linear spectrum. In this
case one can also use a linear spectrum emulator like CosmoPower (Spurio
Mancini et al., 2022) or bacco (Aricò et al., 2021), which takes
$\mathcal{O}(0.1)$ seconds to produce the linear spectrum. Note that one can
also employ CosmoPower to construct an emulator for the linear power spectrum
in the EFTofDE based on EFTCAMB output, overcoming the QS and computational
inefficiency issues.

Given the utility in using the QS, we want to get an idea of its validity. In
Figure 2 we show the effects of the QS at $z=0$ and $z=1$ for models with non-
zero $\alpha_{K}$ and $\alpha_{B}$ (KGB Deffayet et al., 2010), on the
nonlinear spectrum as given by Equation 3. We use the halofit (Takahashi et
al., 2012) formula for $P_{\rm NL}^{\rm pseudo}$ and assume a $\Lambda$CDM
background expansion, $H(a)=H_{\rm\Lambda CDM}(a)$ as well as no screening
effects, i.e., $\mathcal{F}=\mu-1$ and $\gamma_{2}=\gamma_{3}=0$.

We find that the QS is valid for these mild to moderate parameter choices on
scales of $k\geq 0.1h/{\rm Mpc}$. Upcoming surveys will probe scales larger
than this which may be an issue. Taking into account cosmic variance assuming
a galaxy survey volume similar to the effective volume of forthcoming surveys,
$V_{\rm eff}=20~{}{\rm Gpc}^{3}/h^{3}$ (Laureijs et al., 2011; Aghamousa et
al., 2016; Blanchard et al., 2020b), the QS is still a sub-dominant source of
error for even extreme choices of $\alpha_{B}$ and $\alpha_{K}$ (see
subsection 4.3). Note that time derivatives of the fields drop out from the
calculation of $\mu$ for $k\rightarrow\infty$ in Horndeski theories (Lombriser
& Taylor, 2015a; Pace et al., 2021). We note at small scales, modelling
inaccuracies and shot noise errors will arguably dominate any inaccuracies
incurred from using the QS.

We do however warn that the QS begins to break significantly for beyond
Horndeski theories (Lombriser & Taylor, 2015a). For large modifications to GR
within Horndeski, we advise comparing the resulting nonlinear spectrum with
and without the QS against the predicted errors on the specific data that is
being analysed. Further, we have implemented the following necessary condition
for the QS to hold in our code (Peirone et al., 2018)

$\frac{k}{aH(a)}>c_{s}^{2}(a)\,,$ (49)

where $c_{s}^{2}$ is given by Equation 29, with its violation producing a
warning prompt.

Figure 2: Ratio of the quasi-static approximated (QS) nonlinear spectrum to
the exact calculation. We show the results for four EFTofDE models with
$\\{\alpha_{K},\alpha_{B}\\}$ non-zero and all other $\alpha$ parameters set
to 0 and a $\Lambda$CDM background expansion, at $z=0$ (left) and $z=1$
(right). The exact calculation uses Equation 4 with an EFTCAMB linear spectrum
while the QS uses Equation 4 with a rescaled $\Lambda$CDM linear spectrum
using the modified growth equations, making use of Equation 60. The orange
band represents the error coming from cosmic variance assuming an effective
survey volume of $V_{\rm eff}=20~{}{\rm Gpc}^{3}/h^{3}$. We assume
$\mathcal{F}=\mu-1$ and $\mathcal{E}=1$ in all nonlinear computations. The
dotted lines mark 1% deviations which is an optimistic estimate on the
modelling errors of the halo model reaction framework.

###  4.2 $\gamma_{2}=\gamma_{3}=0$ approximation

We begin by noting that setting $\gamma_{2}=\gamma_{3}=0$ implies we have
$R_{\rm SPT}\approx 1$ in section 2 as the 1-halo terms are subdominant. This
forces the argument of the logarithm in Equation 9 to be very close to unity,
giving a very large $k_{\rm star}$. Effectively, this is the same as setting
$\mathcal{E}=1$ in Equation 6. This is the choice we take when adopting this
approximation. We should remark that simply setting $\gamma_{2}=\gamma_{3}=0$
leaves one slightly sensitive to the $k_{\rm star}$ correction through the
1-halo terms and consequently on the particular choice of halo mass function.

Using the exact forms of $\gamma_{2}$ & $\gamma_{3}$ as described in
subsection 3.2 is a big challenge. This is primarily for computational reasons
as it involves numerical time derivatives. Smoothness of such derivatives is
difficult to ensure and can affect results. In particular, the exponential
dependence of $\mathcal{R}$ on $k_{\star}$ (see Equation 4) makes it very
sensitive to inaccuracies in the 1-loop calculation. Further, the full map to
$\gamma_{2}$ & $\gamma_{3}$ from the EFTofDE would increase computational time
significantly, degrading our code’s ability to perform statistical analyses on
data.

To test the impact of setting $\gamma_{2}=\gamma_{3}=0$ we compare Equation 4
with and without these terms switched on for two different theories of
gravity, DGP and the Hu-Sawicki $f(R)$ model (Hu & Sawicki, 2007). The former
is an instance of derivative or Vainshtein screening and the latter of
potential or chameleon screening, covering two main types of screening
mechanism.

This comparison is shown in Figure 3. We find that in the case of DGP, the
correction coming from the 1-loop computation is negligible for small and
moderate modifications to GR at all scales. On the other hand, the corrections
to the $f(R)$ theory can be up to $1.5\%$ at $z=0$ for moderate modifications
to GR. This may be acceptable if these inaccuracies can be partially absorbed
into the nonlinear degrees of freedom. We explore this in section 5.

Figure 3: Ratio of the approximated reaction to the full calculation. The
approximation assumes $\mathcal{E}=1$ in Equation 6 which is approximately
equivalent to no higher order perturbative, screening terms i.e,
$\gamma_{2}=\gamma_{3}=0$. We show Hu-Sawicki $f(R)$ gravity on the left and
the normal branch of DGP on the right for varying modifications to GR. For
$f(R)$ we show predictions for when the model parameter takes the value of
$|f_{R0}|=10^{-5}$ (moderate modification, blue), $|f_{R0}|=10^{-6}$ (low
modification, red) and $|f_{R0}|=10^{-7}$ (very low modification, green). For
DGP we show predictions for when the model parameter is $\Omega_{\rm rc}=0.25$
(moderate modification, blue) and $\Omega_{\rm rc}=0.01$ (low modification,
red). We also show the comparison over two redshifts, $z=0$ (solid lines) and
$z=1$ (dotted lines).

###  4.3 Observational and theoretical constraints

Table 1: Theoretical and observational constraints on $\alpha$-parameters from references in the right column. Theoretical constraints are coming from low-energy (EFT) and high-energy (positivity bounds) physics. Note that $\alpha_{K}$ is not constrained by data, hence the subluminality condition does not impose any constraining power on the scalar mode perturbations. Also note the contradiction in the condition for GW propagation: subluminal versus superluminal speed. The positivity bounds do not hold in general, they are derived for a quadratic subclass of Horndeski theories with $G_{3}=G_{5}=G_{5,X}=0$ in Equation 2. Data driven constraints strongly depend on the imposed theoretical priors and time-dependent parametrisation of $\alpha$-functions. Here we mention only two prior-independent observational constraints. |  | scalar | tensor |
---|---|---|---|---
| no ghost | $\alpha_{K}+\frac{3}{2}\alpha_{B}^{2}>0$ | $M^{2}>0$ | Bellini & Sawicki (2014)
Low | gradient stability | $c_{s}^{2}\geq 0$ | $\alpha_{T}\geq-1$
Energy | (sub)luminality |  large $\alpha_{K}$ | $\alpha_{T}\leq 0$ | de Rham & Melville (2018)
| no GW-induced instability | $|\alpha_{M}+\alpha_{B}|\lesssim 10^{-2}$ |  | Creminelli et al. (2020)
High | scalar-scalar scattering | $\alpha_{B}\leq\frac{2\alpha_{T}}{1+\alpha_{T}}$ | Melville & Noller (2020)
Energy | scalar-matter scattering | $\alpha_{T}\geq 0$ | de Rham et al. (2021)
Data | GW propagation speed | $|\alpha_{T}|\leq 10^{-15}$ | Abbott et al. (2017)
| CMB and LSS | $|\alpha_{M}|,|\alpha_{B}|\leq\mathcal{O}(0.1)$ | Spurio Mancini et al. (2019)

Firstly, we want to eliminate a range of $\alpha$-parameter values that leads
to two pathological instabilities: ghost (i.e., negative kinetic energy) and
gradient (i.e., imaginary speed of sound). These constraints for the Horndeski
theories were first derived in De Felice & Tsujikawa (2012). In terms of the
$\alpha$-functions, Bellini & Sawicki (2014) found that the stability of the
background requires

$\alpha>0\,,\,\,\,\,\,c_{s}^{2}\geq 0\,,$ (50)

from Equation 29 and Equation 30 for scalar modes, and

$M^{2}>0\,,\,\,\,\,\,c_{T}^{2}=1+\alpha_{T}\geq 0\,,$ (51)

for tensor modes of perturbations. An additional theoretical constraint is the
stability of scalar modes in the presence of gravitational waves of large
amplitude, for instance, sourced by massive binary systems (Creminelli et al.,
2020). Mapped to the parameterisation used in this work this requires the
following bound (Noller, 2020):

$|\alpha_{M}+\alpha_{B}|\lesssim 10^{-2}\,.$ (52)

Previously, it was argued that the constraining power of upcoming cosmological
surveys will allow us to pin down the $\alpha$-parameters at the
$\mathcal{O}(0.1)$-level (e.g., Frusciante et al., 2019). For the condition
above this implies that $\alpha_{M}\approx-\alpha_{B}$. However, in such
forecasts nonlinear scales were ignored with a typical highest mode around
$k_{\mathrm{max}}\approx 0.15$ $h$ Mpc-1. We speculate that this constraint
may be improved upon by inclusion of the nonlinear scales. Therefore, in our
code we treat $\alpha_{B}$ and $\alpha_{M}$ independently.

Secondly, one may consider that the new physics should not modify the speed of
gravitational wave propagation (Lombriser & Taylor, 2016; Abbott et al., 2017;
Lombriser & Lima, 2017; Creminelli & Vernizzi, 2017; Ezquiaga &
Zumalacárregui, 2017; Baker et al., 2017; Sakstein & Jain, 2017; Battye et
al., 2018; de Rham & Melville, 2018; Creminelli et al., 2018), and so
$\alpha_{T}=\bar{M}_{2}^{2}=0$. This luminality condition has been argued to
not be as clear cut a constraint through EFT considerations (de Rham &
Melville, 2018; Baker et al., 2022) as well as through the positivity bounds
from high energy physics (de Rham et al., 2021), so in our code we keep the
$\alpha_{T}$ dependence in $\mu$. Subluminality, stated in the former
references, follows from the existence of a Wilsonian UV completion (Adams et
al., 2006) and dependence on the theory’s ‘cutoff’ scale. From Equation 29 it
can be seen that subliminality of scalar modes is guaranteed for large values
of $\alpha_{K}$, while for tensor modes subluminality requires $\alpha_{T}<0$.
Superluminality, stated in de Rham et al. (2021), is a consequence of the
positivity bounds for scattering between scalar and matter fields. Such
positivity bounds require a unitary, causal, local UV completion of our low-
energy EFT theory. However, superluminality does not necessary result in
casual paradoxes (Babichev et al., 2008; Burrage et al., 2012). In general,
the notion of causality in terms of the low-energy EFT is a rather subtle
topic (for instance, see de Rham & Tolley, 2020; Reall, 2021).

Thirdly, in the QS $\alpha_{K}$ does not enter the equations of motion
(Bellini & Sawicki, 2014). Therefore, it is completely unconstrained in our
approach, or for any model with $c^{2}_{s}\approx 1$. However, in the exact
computation $\alpha_{K}$ affects only the largest scales (see Figure 2), which
are dominated by cosmic variance. This can be a motivation to not consider
$\alpha_{K}$ in data analyses, leaving only $\alpha_{M}$ and $\alpha_{B}$ in a
‘bare-bones’ case. We do not impose any of these reductions in our code and
leave it to the user to specify well motivated priors on the full set of
EFTofDE parameters in their analyses.

Lastly, we note that there are a host of data driven constraints that one can
put on the EFTofDE parameters (Huang, 2016; Bellini et al., 2016; Noller &
Nicola, 2019a, b; Spurio Mancini et al., 2019; Melville & Noller, 2020;
Noller, 2020; de Rham et al., 2021). Such constraints strongly depend on the
imposed theoretical priors and time-dependent parameterisation of the
$\alpha$-functions (see subsection 4.4). However, they all agree that the
uncertainties and values of the $\alpha$-parameters are of order
$\mathcal{O}(0.1)$. The future CMB and LSS surveys promise to improve the
constraints up to at least one order of magnitude
$\sigma(|\alpha_{i}|)\sim\mathcal{O}(0.01)$ (see, for example, Abazajian et
al., 2016). One may also assume a $\Lambda$CDM background, well motivated by
CMB data (e.g., Aghanim et al., 2020), and so set $H(a)=H_{\rm\Lambda
CDM}$(a)999Our code defaults to this assumption, but there is the option to
parameterise the background too.. We summarize the constraints discussed above
in Table 1.

###  4.4 Parameterising time dependence

Here, we look at how one can parameterise the time dependence of the EFTofDE
functions. To first order this can be approximated by a Taylor expansion,
$\theta_{i}(a)\approx\theta_{i,0}+\theta_{i,p}(1-a)$, leaving at least 6 free
constants characterising deviations from $\Lambda$CDM. In typical data
analyses, only a 1-parameter time dependence is considered. For example, in
Noller & Nicola (2019a) the authors consider the following three
parameterisations for the $\alpha_{i}$, $i\in\\{M,B,K,T\\}$

$\displaystyle\mathbf{(1):}\,\alpha_{i}(a)$
$\displaystyle=c_{i}\Omega_{\rm\Lambda}(a)\,,$ (53)
$\displaystyle\mathbf{(2):}\,\alpha_{i}(a)$ $\displaystyle=c_{i}a\,,$ (54)
$\displaystyle\mathbf{(3):}\,\alpha_{i}(a)$ $\displaystyle=c_{i}a^{n_{i}}\,,$
(55)

where $c_{i}$ and $n_{i}$ are free constants and $\Omega_{\rm\Lambda}(a)$ is
the $\Lambda$CDM cosmological constant energy density fraction as a function
of time. For a comprehensive list of various other time parameterisations see
Appendix B of Frusciante & Perenon (2020). These all draw on the motivation
that modifications should only become relevant at late times. In our code, the
default is set to (2) for all $\alpha_{i}$. We note that such parametrisations
may exclude well-known theories as shown in Kennedy et al. (2018), which
motivated the $s$-basis introduced in subsection 3.1.

We can also adopt similar parametrisations for the background $H(a)$, but a
more general choice would be for example the Chevalier-Polarski-Linder (CPL)
parametrisation (Chevallier & Polarski, 2001; Linder, 2003), which
parametrises the dark energy equation of state $w(a)$ in terms of two free
constants, $\\{w_{0},w_{a}\\}$ as

$w(a)=w_{0}+w_{a}(1-a)\,,$ (56)

which gives the following form for $H(a)$

$H^{2}(a)=H^{2}_{0}\left(\Omega_{{\rm
m},0}a^{-3}+\Omega_{\Lambda}e^{-3\int(1+w(a)){\rm d}\log a}\right)\,.$ (57)

Table 2: A maximal, reduced and minimal set of parameters needed for a comprehensive nonlinear power spectrum analysis of the unrestricted theory space of Equation 2, together with a proposed minimal set for more general theories. The Horndeski minimal set assumes $\alpha_{T}=\alpha_{K}=0$ and Equation 43, while the maximal and the reduced assume full freedom of Equation 39. Both reduced and minimal assume $\gamma_{2}=\gamma_{3}=0$. The totals show the number of free functions of time plus any free constants. We note $q_{4}$ and $p_{8}$ are found in section 5 to be likely irrelevant for the 1-halo computation, and so we do not consider them in the reduced or minimal cases. We also note very tight constraints on $\alpha_{B}-\alpha_{M}$ (see subsection 4.3) relevant to the minimal case.  | Maximal | Reduced | Minimal (Horndeski) | Minimal (general)
---|---|---|---|---
Background | $H(a)$ | $H(a)$ | $H(a)$ | $w_{0},w_{a}$
Linear | $\alpha_{M}(a),\alpha_{B}(a),\alpha_{K}(a),\alpha_{T}(a)$ | $\alpha_{M}(a),\alpha_{B}(a),\alpha_{K}(a),\alpha_{T}(a)$ | $\alpha_{M}(a),\alpha_{B}(a)$ | $\gamma$
Quasi-nonlinear | $\xi_{3}^{2}(a),\xi_{4}^{2}(a),\xi_{3}^{3}(a),\xi_{4}^{3}(a),\xi_{4}^{4}(a),\xi_{4}^{5}(a)$ | - | - | -
Nonlinear | $p_{1-7}(a)+p_{8}$ | $p_{1-7}(a)$ | $q_{1},q_{2},q_{3}$ | $q_{1},q_{2},q_{3}$
Total | 18+1 | 12 | 3 + 3 constants | 6 constants

###  4.5 Parametrisation of $\mathcal{F}$

The nPPF form for $\mathcal{F}$ given in Equation 39 captures dependencies of
the nonlinear modification to the Poisson equation on the relevant variables,
namely $\\{y_{\rm h},a,M_{vir},y_{\rm env}\\}$. Being motivated by the form of
$\mathcal{F}$ in DGP (Equation 100), it can recover the DGP form given
appropriate choices for $p_{i}$ albeit with a non-trivial dependency of
$p_{2}$ on $a$ (see Equation 108). Equation 39 becomes approximate when moving
beyond DGP. On the other hand, the Erf form, Equation 41, is completely
phenomenological and is an approximation even in DGP.

Note that the nPPF is also more directly relatable to specific actions and
gravity models, in which case its degrees of freedom can be significantly
restricted. It is thus far more suitable when particular models are being
targeted for analysis. The Erf model on the other hand is completely general
and has no direct relation to specific actions of gravity. It is thus more
suitable when no specific model is being targeted and we want to place
constraints on general models of gravity. In section 5 we test these two
approximations in both DGP and $f(R)$ gravity.

###  4.6 Overview

With all these approximations and constraints, the arguable minimal parameter
space characterising deviations to $\Lambda$CDM is 3 free functions of time
and 4 constants. Without approximations or constraints, the maximal is 18 free
functions of time and a constant. Of course we can also find intermediate
reduced sets, such as using the nPPF but with $\gamma_{2}=\gamma_{3}=0$. Given
we need to parameterise these functions of time, the maximal set is currently
an unfeasible parameter space to probe comprehensively, both in terms of data
processing as well as parameter degeneracies which limits the amount of useful
physical information one can extract from the data.

Finally, we have focused on the Horndeski class of models, but one can extend
this to larger generality by considering for example the growth index $\gamma$
parametrisation for $\mu(a)$ (Peebles, 1980; Linder & Cahn, 2007) (explicitly,
see Eq. 47 of Kennedy et al., 2018) and Equation 57 for $H(a)$. Combined with
the Erf model, this would constitute a minimal set of 6 free constants for
general modifications to $\Lambda$CDM. This minimal model-independent
parametrisation has also been implemented into the code.

We summarise these parameterisations in Table 2.

##  5 Testing the nonlinear parameterisations

In this section we compare the predictions for the halo model reaction
$\mathcal{R}$, using the various nonlinear parameterisations of modifications
to the Poisson equation outlined in subsection 3.3, to exact solutions as well
as state-of-the-art emulators within an evolving dark energy scenario (wCDM),
DGP and Hu-Sawicki $f(R)$ gravity. Note that the exact solutions for the
reaction have in turn been themselves compared to full $N$-body simulations in
other works (see Cataneo et al., 2019, for example), exhibiting
$\mathcal{O}(1)\%$ agreement. These models cover a fair range of theoretical
and phenomenological features typical of modified gravity and dark energy
models, making them good representatives and test cases.

We look to test predictions for $\mathcal{R}$ using Equation 39 (nPPF) and
Equation 41 (Erf) with $\gamma_{2}=\gamma_{3}=0$ against the full calculation
which computes $\mathcal{R}$ using exact forms for $\gamma_{2}$, $\gamma_{3}$
and $\mathcal{F}$ (see Appendices of Bose et al., 2020b, for all relevant
expressions). We further employ the EuclidEmulator2 emulator (Knabenhans et
al., 2021) and the fofr emulator (Winther et al., 2019) for the wCDM and
$f(R)$ cases respectively. These emulators have been trained on high quality
$N$-body simulations and are 1-2% accurate within the scales we examine,
providing a good benchmark for our predictions. One should keep in mind that
the halo model reaction approach’s accuracy is limited by the pseudo power
spectrum employed. For example, if we use HMCode2020 (Mead et al., 2021) for
the pseudo, which is claimed to be $2.5\%$ accurate down to $k\leq 10~{}h/{\rm
Mpc}$, we then expect any power spectrum comparisons to then be $\sim 4.5\%$
consistent with $N$-body at $k\leq 3~{}h/{\rm Mpc}$, which assumes the result
of Cataneo et al. (2016), i.e., that the exact solution for $\mathcal{R}$ is
$\sim 2\%$ accurate at these scales. In the wCDM and DGP cases, both
$\Lambda$CDM and $P_{\rm NL}^{\rm pseudo}$ are computed using the halofit
fitting function (Takahashi et al., 2012), but the $f(R)$ case uses
HMCode2020.

The computation of $\mathcal{R}$ requires us to solve the evolution equations
for the spherical top-hat radius parametrised by $y_{\rm h}$ (Equation 103).
This necessitates the specification of $\mathcal{F}$ at all redshifts up to
the target redshift. We then should test approximations for $\mathcal{F}$ even
at high redshifts, which is done in subsection B.3, where we compare
$\mathcal{F}$ at $z=0,1,4$. For comparisons of the halo model reaction, we
only consider $z=0,1$ which are more observationally relevant.

We fit $\\{q_{1},q_{2},q_{3},q_{4}\\}$ for the Erf model, $\mathcal{F}_{\rm
Erf}$. In the nPPF case, we do not fit all the 8 free parameters of
$\mathcal{F}_{\rm nPPF}$, and only consider $p_{1}$ and $p_{8}$, treating both
as constants. In principle, and indeed for unspecified theories of gravity,
all 8 parameters will be fit to the data. For the comparisons made here,
$p_{2-7}$ are fixed to the theoretically predicted values quoted in Appendix
B. Fitting such a high dimensional parameter space is beyond the scope of this
paper.

In what follows we fit the free parameters by performing a least square fit to
the exact $\mathcal{R}$ prediction. We choose to fit our parametrised models
to the exact predictions for $\mathcal{R}$, rather than the emulator
predictions for $P_{\rm NL}(k)$ for two reasons. First so as to test the
ansatz for the phenomenological screening and the consistency of the
predictions (see Appendix B). Second, we do not want to assume anything about
the pseudo spectrum in these fits. To fit we minimise the following merit
function

$s^{2}=\sum^{\rm max}_{j={\rm min}}\sum^{\rm max}_{i={\rm
min}}\frac{[\mathcal{R}_{\rm exact}(k_{i},z_{j})-\mathcal{R}_{\rm
approx}(k_{i},z_{j})]^{2}}{\sigma_{i,j}^{2}}\,,$ (58)

where we assume error bars on $\mathcal{R}_{\rm exact}$ coming from cosmic
variance (Zhao, 2014; Blanchard et al., 2020a; Mancarella et al., 2022) and a
constant systematic error added in quadrature

$\sigma_{i,j}^{2}(k,z)=\frac{4\pi^{2}}{k_{i}^{2}\Delta
k_{i}V_{s,j}}+\sigma_{\rm sys}^{2}\,,$ (59)

where $V_{s,j}\in\\{0.3,8\\}~{}{\rm Gpc}^{3}/h^{3}$ is taken to be a stage IV
survey-like volume for each bin $z_{j}\in\\{0,1\\}$ respectively (Laureijs et
al., 2011; Aghamousa et al., 2016; Mancarella et al., 2022; Blanchard et al.,
2020b). We fit in the range $k_{\rm min}=0.1\leq k_{i}\leq 3=k_{\rm max}$
which is the range over which the exact computation of $\mathcal{R}$ is $2\%$
accurate (Cataneo et al., 2019), sampling logarithmically, with $\Delta k_{i}$
being the bin width. We take $\sigma_{\rm sys}=0.02$ to reflect the systematic
error in the parametrised reaction when compared to simulations by proxy of
the exact solution. The best fit parameter values are shown in Table 3.

Table 3: Best fit parameter values for the DGP and $f(R)$ models. The fit is performed to the exact solution for $\mathcal{R}$ in the range $0.1~{}h/{\rm Mpc}\leq k\leq 3~{}h/{\rm Mpc}$ and at $z=0,1$ as described in the main text. For the Erf model, we do not fit $q_{2}$, $q_{3}$ and $q_{4}$ for DGP and for $f(R)$ we find the quality of fit with and without $q_{4}$ is similar. For all fits we thus set $q_{4}=0$. The nPPF is exact for DGP and so we only consider $f(R)$, fixing all $p_{2-7}$ to the values given in Equation 109.  | nPPF | Erf
---|---|---
$\Omega_{\rm rc}$ | $p_{1}$ | $p_{8}$ | $q_{1}$ | $q_{2}$ | $q_{3}$
0.25 | - | - | 0.76 | 0 | 0
0.01 | - | - | 0.71 | 0 | 0
$|f_{\rm R0}|$ | $p_{1}$ | $p_{8}$ | $q_{1}$ | $q_{2}$ | $q_{3}$
$10^{-5}$ | 3 | -0.8 | 0.9 | 0.35 | 0.65
$10^{-6}$ | 8.5 | -0.5 | 1.65 | 0.7 | 2.45
$10^{-7}$ | 5.65 | -0.45 | 0.6 | 0.8 | 2.15

###  5.1 Evolving dark energy example: wCDM

Here we perform a sanity check that the general minimal model outlined in
Table 2 produces consistent results for a wCDM cosmology, and is at least as
accurate as the exact solution. To do this we compare a minimal model with CPL
parameters $w_{0}=-1.2$ and $w_{a}=0.4$, and a growth index of $\gamma=0.55$
to the exact solution as well as predictions from EuclidEmulator2 using the
same CPL parameters. We further set the nonlinear parameters of the Erf model
($q_{i}$) to unity, but check that they have no impact on the results as
expected from Equation 41 ($\mu\approx 1$ for $\gamma=0.55$).

We show our results in Figure 4. We see that the minimal model is both
completely consistent with the exact solution which has no nonlinear or linear
modification to the Poisson equation, as well as $1\%$ consistent with the
emulator down to $k\leq 2~{}h/{\rm Mpc}$ and $2\%$ down to $k\leq 3~{}h/{\rm
Mpc}$. The minimal general model could feasibly outperform the exact solution
given its degrees of freedom. In a future work we plan to check forecasted
constraints and possible biases on cosmological parameters for the minimal
general model, in full posterior estimation analyses employing $N$-body
simulation measurements.

Figure 4: Top panel: The ratio of the $w$CDM nonlinear power spectrum to the
$\Lambda$CDM nonlinear power spectrum computed using the EuclidEmulator2
(black) and halofit together with the halo model reaction (see Equation 3) for
the minimal general model (blue) as outlined in the right most column of Table
2. Middle panel: The ratio between theoretical and emulator predictions for
the ratio between wCDM to $\Lambda$CDM spectra, i.e., the ratio of blue to
black top panels curves. Bottom panel: The ratio of the exact halo model
reaction to the minimal general model. We plot the ratio for two
observationally relevant redshifts, $z=0$ (solid) and $z=1$ (dotted). We show
these results for $w_{0}=-1.2$ and $w_{a}=0.4$. The minimal general model also
has $\gamma=0.55$ and $q_{1}=q_{2}=q_{3}=q_{4}=1$.

###  5.2 Vainshtein example: DGP

For DGP the nPPF parameterisation reproduces the exact form of $\mathcal{F}$
(Equation 100) for specific choices of the $p_{i}$ parameters (Equation 108).
On the other hand, the Erf parametrisation (Equation 41) is approximate and we
fit the associated parameters. We note that DGP has no Yukawa suppression at
large scales and produces a constant enhancement of the $\Lambda$CDM linear
growth factor. This enhancement is controlled by the DGP degree of freedom
$\Omega_{\rm rc}\equiv 1/(4H_{0}^{2}r_{c}^{2})$ where $r_{c}$ is the cross-
over scale dictating where gravity goes from behaving 4-dimensionally to
5-dimensionally. We consider two levels of deviation to $\Lambda$CDM: a
moderate modification given by $\Omega_{\rm rc}=0.25$ and a small modification
given by $\Omega_{\rm rc}=0.01$.

We only fit $q_{1}$ as we do not have any mass, environment or Yukawa-
suppression scale dependence, and so we set $q_{2}=q_{3}=q_{4}=0$ in this
case. The best fit values of $q_{1}$ are given in Table 3. Further, we employ
the exact form of $\mu(a)$ in Equation 41 (see appendices of Bose et al.,
2020b, for the explicit expression).

In the top panels of Figure 5 we show the ratio of a DGP power spectrum to a
$\Lambda$CDM spectrum with the same background expansion history, normalised
to unity at linear scales for easier comparisons of nonlinear effects. The DGP
spectrum is given by Equation 3. We see the moderate modification gives up to
a $6\%$ deviation from $\Lambda$CDM (above the linear growth enhancement) for
$k\leq 3h/{\rm Mpc}$ while the small modification can reach $2\%$ over the
same range of scales. Reassuringly, in the bottom panels we find sub-percent
agreement between the Erf and exact predictions down to $k=5~{}h/{\rm Mpc}$,
with a smaller disagreement for the smaller deviation from $\Lambda$CDM.

One can further parameterise the time dependence of $q_{1}$ which would
alleviate some of these deviations, but we find these differences to be more
than acceptable given the relative size compared to the modification to
$\Lambda$CDM shown in the top panels. Moreover, a large number of additional
degrees of freedom will be introduced in real data analyses such as intrinsic
alignments and parameterisations of baryonic physics. These will be degenerate
to some level with modified gravity effects (see Schneider et al., 2020, for
example), allowing lower accuracy demands in the modelling of $\mathcal{R}$.

This additional time dependence is highlighted in Figure 8, where we find that
the Erf model can match the exact form of $\mathcal{F}$ extremely well at a
fixed redshift. Upon investigation, we found this dependence to be highly
degenerate with $q_{2}$ which prompted us to not introduce new freedom to the
model, especially because we can achieve very good fits already, even without
$q_{2}$.

Note we have not compared the parametrised model to an emulator nor
simulations in this case. Given the excellent agreement with the exact
solution we can infer its accuracy is at least as good as the exact solution,
given it employs 3 additional degrees of freedom. We remind the reader that
the exact reaction was found to be $2\%$ accurate when compared to $N$-body
simulations in Cataneo et al. (2019).

Figure 5: Top panels: The ratio of the DGP nonlinear power spectrum to the
$\Lambda$CDM nonlinear power spectrum computed using halofit and the halo
model reaction (see Equation 3) for the exact (black) and Erf (blue) cases. We
do not show the nPPF case as it reduces to the exact solution for specific
choices of its parameters. The Erf model assumes $\gamma_{2}=\gamma_{3}=0$. We
have normalised the ratio to unity at large scales for easier comparisons.
Bottom panels: The ratio of halo model reactions; the Erf model
$\mathcal{R}_{\rm Erf}$ to the exact solution. This is equivalent to the ratio
of the top panel blue to black curves. We show these results for a moderate
modification, $\Omega_{\rm rc}=0.25$ (left) and a low modification,
$\Omega_{\rm rc}=0.01$ (right). We plot the ratio for two observationally
relevant redshifts, $z=0$ (solid) and $z=1$ (dotted).  Figure 6: Top panel:
The ratio of the $f(R)$ nonlinear power spectrum to the $\Lambda$CDM nonlinear
power spectrum computed using HMCode2020 and the halo model reaction (see
Equation 3) for the exact (grey), Erf (blue) and nPPF (red) cases. The fofr
emulator is also shown in black. Middle panel: The ratio between theoretical
and emulator predictions for the ratio between $f(R)$ to $\Lambda$CDM spectra,
i.e., the ratio of grey, blue and red to black top panels curves. Bottom
panel: The ratio of halo model reactions; the parameterised models to the
exact solution. This is equivalent to the ratio of the top panel coloured
curves to the gray curves. Note that both parameterised models have
$\gamma_{2}=\gamma_{3}=0$. We plot the ratio for two observationally relevant
redshifts, $z=0$ (solid) and $z=1$ (dotted). We show these results for a
moderate modification, $|f_{\rm R0}|=10^{-5}$. The orange bands indicate the
$2\%$ region which is the current absolute accuracy of the exact
$\mathcal{R}$.

Figure 7: Same as Figure 6 for a low modification, $|f_{\rm R0}|=10^{-6}$
(top) and a very low modification, $|f_{\rm R0}|=10^{-7}$ (bottom). Note that
the fofr emulator for $|f_{\rm R0}|=10^{-7}$ gives the $\Lambda$CDM prediction
and so we omit the middle panel.

###  5.3 Chameleon example: Hu-Sawicki $f(R)$

For this theory we consider both the nPPF and Erf models for $\mathcal{R}$,
and compare them to the exact solution (Equation 105) as well as at the power
spectrum level to the fofr emulator of Winther et al. (2019). This model makes
use of the chameleon screening mechanism which exhibits an environmental and
mass dependence. It also has a Yukawa suppression which returns it to GR at
large scales. The additional degree of freedom is the value of the background
scalar field at $z=0$, $f_{\rm R0}$, which controls the level of deviation
from GR. We consider three levels of deviation from $\Lambda$CDM, $|f_{\rm
R0}|=10^{-5}$ (moderate modification), $|f_{\rm R0}|=10^{-6}$ (low
modification) and $|f_{\rm R0}|=10^{-7}$ (very low modification). We note that
the moderate $f(R)$ modification is already ruled out by data (see Cataneo et
al., 2015; Desmond & Ferreira, 2020; Lombriser, 2014; Brax et al., 2021, for
example), but provides a good flexibility test of the parameterisation.

In the nPPF case, we choose the theoretically motivated parameters given in
Equation 109. These emerge from a parameterised form of $f(R)$ gravity
(Lombriser et al., 2014) and so are approximate. $p_{1}$ and our new parameter
$p_{8}$ remain free. Treating them both as constants, we fit them in the same
way that we fit the Erf model’s parameters, by minimising Equation 58. We note
that the other nPPF parameters, $p_{2}-p_{7}$, take on different forms for the
chameleon screening and Yukawa suppression regimes. We only consider the
screening regime which is more relevant for the spherical collapse
calculation, and rely on $\mu(\hat{k},a)$ appearing in Equation 39 to take
care of the Yukawa suppression.

Yukawa suppression is relevant for large masses, large $y_{\rm env}$ or small
values of $f_{\rm R0}$. Given this, we do not expect $p_{8}$ or $q_{4}$ to be
relevant for spherical collapse where $y_{\rm h}\leq y_{\rm env}\leq 1$, and
even less for the 1-halo spectrum where the Sheth-Torman mass function down-
weights large masses (see, for example, Schmidt et al., 2009). We verify this
by performing two separate fits: the first only including the parameter sets
$\\{p_{1}\\}$ and $\\{q_{1},q_{2},q_{3}\\}$ for the nPPF and Erf model
respectively, while the second extending these sets to include $p_{8}$ and
$q_{4}$ respectively.

We find that values of $q_{4},p_{8}\geq 0$ negligibly change the goodness of
fit for the low and very low modification strengths, while sufficiently
negative values degrade the fit, which is expected as the Yukawa scale begins
to overlap with the screening scale. Further, we observe only a marginal
improvement at $z=0$ for $|f_{\rm R0}|=10^{-5}$ in the Erf case. Given this,
all fits shown and quoted here set $q_{4}=0$ for the Erf case. In the nPPF
case, we observe a moderate improvement for $|f_{\rm R0}|=10^{-5}$ and so keep
$p_{8}$. We report the best-fit parameters in Table 3.

The $f(R)$ results are shown in Figure 6 and Figure 7. We see the moderate
modification can reach a $20\%$ deviation from $\Lambda$CDM for $k\leq 3h/{\rm
Mpc}$ while the low and very low modifications reach $10\%$ and $3\%$
respectively. Both parameterisations do well in modelling the moderate
modification case $|f_{\rm R0}|=10^{-5}$, shown in Figure 6. The Erf model
prediction for $\mathcal{R}$ stays within $1\%$ of the exact solution for
$k\leq 3h/{\rm Mpc}$. Similarly, the nPPF remains within $2\%$ for $k\leq
2h/{\rm Mpc}$. The situation improves for the lower modification cases, shown
in Figure 7. These comparisons exhibit sub-$1\%$ agreement between the Erf
(nPPF) model and exact solution for $k\leq 5(3)~{}h/{\rm Mpc}$ at $z=0$ and
$z=1$.

All power spectra predictions are $\sim 3\%$ consistent with the fofr emulator
which mainly demonstrates the accuracy of HMCode2020. Interestingly, we find
that the additional degrees of freedom within the nPPF and Erf models are
degenerate with possible inaccuracies in the pseudo, even down to
$k=5~{}h/{\rm Mpc}$. Again, we leave it to a future work to see if these
additional degrees of freedom can improve constraining power on cosmological
and gravitational parameters while remaining unbiased.

Our comparisons indicate that for the Erf model, degeneracies between
$q_{1-3}$ and $q_{4}$ make the latter parameter unnecessary. We note that the
fit of $q_{1-3}$ becomes insensitive to the value of $q_{4}$ if it is
sufficiently large, here found to be $q_{4}=0$. For the nPPF model, the
additional freedom provided by $p_{8}$ is necessary to improve the fit, but it
does not help substantially for observationally viable values of $f_{\rm R0}$.
Further, we remind the reader that we do not know $p_{2-7}$ a priori for
unspecified theories of gravity, and so the importance of $p_{8}$ is likely to
be minimal when considering these additional degrees of freedom.

Lastly, we remark that the Erf model gives a good fit for a range of values
for $q_{1-3}$101010Similar fits were found for $\mathcal{O}(0.1)$ values for
these parameters.. The values quoted in Table 3 are only the best fit values,
which are also very dependent on Equation 59. This makes it hard to extract
any further dependence on $f_{\rm R0}$ in Equation 41 (note this already
depends on $\mu(k,a)$), which is also beyond the scope of this parametrisation
which aims to be general in terms of gravitational degrees of freedom.

##  6 Summary

In this paper we have presented a significant extension of the code described
in Bose et al. (2020b) which produces nonlinear corrections to the matter
power spectrum coming from beyond-$\Lambda$CDM physics in the form of the halo
model reaction $\mathcal{R}$. In particular, we have focused on implementing
parameterisations of key equations, in particular the background expansion
history and the linear and nonlinear Poisson equations.

For the linear scales and background we have considered the effective field
theory of dark energy (EFTofDE) while for the nonlinear scales we have
considered two distinct parameterisations, a nonlinear parameterised post-
Friedmannian (nPPF) based model and a more phenomenological model based on the
error function (Erf). Together, these give a general parameterisation of the
nonlinear matter power spectrum in Horndeski models. We neglect loop
corrections in these parameterisations but leave these as viable additions and
we provide theoretical and numerical means of deriving these for the Horndeski
class of theories. This being said, we remark that the nonlinear
parametrisations are completely general, and so to move beyond the Horndeski
class it is sufficient to parametrise only the background expansion history
and the linear modification to the Poisson equation. Further, the nonlinear
parametrisations also have unscreened limits, and so we are not restricted to
theories exhibiting screening. In summary, this work presents a fast, accurate
and highly general nonlinear power spectrum predictor for non-standard models
of gravity and cosmology including massive neutrinos, parameterised with a
minimal set of free, physically meaningful constants.

We have tested these parameterisations against the full solutions for
$\mathcal{R}$ in three beyond-$\Lambda$CDM models, wCDM, Hu-Sawicki $f(R)$ and
DGP gravity. This has identified a minimal set of 3 free functions of time and
3 dimensionless, positive, $\mathcal{O}(1)$ dimensionless constants, which can
replicate the exact solutions to within $1\%$ at $k\leq 5h/{\rm Mpc}$ and at
$z\leq 1$ for modifications to GR within current data constraints and within
the Horndeski class. This level of imprecision is sub-dominant to the $2\%$
accuracy currently achieved by the reaction method at these scales (Cataneo et
al., 2019, 2020), and further to the inaccuracies in current pseudo spectrum
prescriptions (Bose et al., 2021; Carrilho et al., 2022). We have seen that
the additional parameters have some degree of degeneracy with pseudo spectrum
inaccuracies, which may improve the scales of validity for the nonlinear power
spectrum $P_{\rm NL}$ as predicted within the halo model reaction framework.
We thus suspect that this minimal parametrisation is acceptable for upcoming
Stage IV cosmic shear analyses given the flexibility of the nonlinear
parameterisation and the many other nuisance degrees of freedom entering a
real data analyses, such as those characterising baryonic physics or intrinsic
galaxy alignments (see, for example, Tröster et al., 2021).

The Erf model is also highly model independent, capturing the basic
phenomenology of screening mechanisms. It can thus be suitable for analyses
targeting general deviations from $\Lambda$CDM. For example, one may perform a
model independent analysis combining the Erf parametrisation with the linear
theory growth index $\gamma$-parametrisation (Peebles, 1980; Linder & Cahn,
2007) (also see Eq. 47 of Kennedy et al., 2018) and say the background
parametrisation of Chevallier & Polarski (2001); Linder (2003), giving 6 free
constants characteristing general deviations from $\Lambda$CDM in the matter
power spectrum at a wide range of scales. On the other hand, the nPPF approach
is complementary as it can be directly related to specific actions of Nature,
making it very suitable when we look to constrain more specific classes of
theories.

In future work we will test the robustness of the minimal parameterisation,
and forecast constraints on deviations to $\Lambda$CDM by performing full
Markov chain Monte Carlo (MCMC) analyses on mock data of the cosmic shear
spectrum. Consistency and accuracy checks can also be performed using recently
developed parametrised modified gravity simulations (Hassani & Lombriser,
2020; Srinivasan et al., 2021; Fiorini et al., 2021; Wright et al., 2022;
Brando et al., 2022). On this note, our code is as fast as the original ReACT
and so is capable of running MCMC analyses. Despite its appreciable baseline
speed, we aim to make this even faster by creating emulators based off halo
model reaction predictions using the recently released CosmoPower code (Spurio
Mancini et al., 2022) which will highly optimise such analyses. It is a future
plan to also perform real data analyses on currently available cosmic shear
data to constrain deviations to $\Lambda$CDM using the general minimal
parametrisation given in Table 2.

It is currently an ongoing project to also extend the halo model reaction to
redshift space and biased tracers in a vein similar to Bose et al. (2020a). We
also plan to include interacting dark energy parametrisations (Gleyzes et al.,
2015; Skordis et al., 2015), a scenario where essentially one decouples the
baryons from $\Lambda$CDM modifications, contrary to the scenario considered
in this paper where all matter is coupled to the scalar field.

## Acknowledgments

The authors would like to thank Matteo Cataneo, Daniel B Thomas, Tessa Baker
and Filippo Vernizzi for useful comments and suggestions. They further thank
the referee for their useful comments. Lastly, they thank all fautors of
fautor.org/papers/0005. BB, JK and LL acknowledge support from the Swiss
National Science Foundation (SNSF) Professorship grant Nos. 170547 & 202671\.
BB was supported by a UK Research and Innovation Stephen Hawking Fellowship
(EP/W005654/1). ANT acknowledges support from a STFC Consolidated Grant. MT’s
research is supported by a doctoral studentship in the School of Physics and
Astronomy, University of Edinburgh. AP is a UK Research and Innovation Future
Leaders Fellow [grant MR/S016066/1]. For the purpose of open access, the
author has applied a Creative Commons Attribution (CC BY) licence to any
Author Accepted Manuscript version arising from this submission.

## Data Availability

The software used in this article is publicly available in the ACTio-ReACTio
repository at https://github.com/nebblu/ACTio-ReACTio. In the same repository
we also provide two Mathematica notebooks: GtoPT.nb that explicitly calculates
the modified gravity 1st, 2nd and 3rd order Poisson equation modifications for
particular covariant theories of gravity as well as provides maps between the
EFTofDE $M$\- and $\alpha$-bases, and Nonlinear.nb which provides expressions,
tests and comparisons of the nonlinear Poisson equation modifications.

## References

  * Abazajian et al. (2016) Abazajian K. N., et al., 2016, CMB-S4 Science Book, First Edition (arXiv:1610.02743)
  * Abbott et al. (2017) Abbott B. P., et al., 2017, Astrophys. J., 848, L13
  * Adams et al. (2006) Adams A., Arkani-Hamed N., Dubovsky S., Nicolis A., Rattazzi R., 2006, JHEP, 10, 014
  * Agarwal & Feldman (2011) Agarwal S., Feldman H. A., 2011, Mon. Not. Roy. Astron. Soc., 410, 1647
  * Aghamousa et al. (2016) Aghamousa A., et al., 2016, The DESI Experiment Part I: Science,Targeting, and Survey Design (arXiv:1611.00036)
  * Aghanim et al. (2020) Aghanim N., et al., 2020, Astron. Astrophys., 641, A6
  * Alam et al. (2021) Alam S., et al., 2021, Phys. Rev. D, 103, 083533
  * Angulo et al. (2021) Angulo R. E., Zennaro M., Contreras S., Aricò G., Pellejero-Ibañez M., Stücker J., 2021, MNRAS, 507, 5869
  * Appleby & Linder (2020) Appleby S., Linder E. V., 2020, JCAP, 12, 036
  * Aricò et al. (2021) Aricò G., Angulo R. E., Zennaro M., 2021, doi:10.12688/openreseurope.14310.2
  * Arnold et al. (2022) Arnold C., Li B., Giblin B., Harnois-Déraps J., Cai Y.-C., 2022, Mon. Not. Roy. Astron. Soc., 515, 4161
  * Babichev et al. (2008) Babichev E., Mukhanov V., Vikman A., 2008, JHEP, 02, 101
  * Babichev et al. (2009) Babichev E., Deffayet C., Ziour R., 2009, Int. J. Mod. Phys. D, 18, 2147
  * Baker et al. (2017) Baker T., Bellini E., Ferreira P. G., Lagos M., Noller J., Sawicki I., 2017, Phys. Rev. Lett., 119, 251301
  * Baker et al. (2022) Baker T., et al., 2022, JCAP, 08, 031
  * Battye et al. (2018) Battye R. A., Pace F., Trinh D., 2018, Phys. Rev., D98, 023504
  * Bellini & Sawicki (2014) Bellini E., Sawicki I., 2014, JCAP, 07, 050
  * Bellini et al. (2016) Bellini E., Cuesta A. J., Jimenez R., Verde L., 2016, J. Cosmology Astropart. Phys., 2016, 053
  * Bernardeau et al. (2002) Bernardeau F., Colombi S., Gaztanaga E., Scoccimarro R., 2002, Phys. Rept., 367, 1
  * Blanchard et al. (2020a) Blanchard A., et al., 2020a, Astron. Astrophys., 642, A191
  * Blanchard et al. (2020b) Blanchard A., et al., 2020b, Astron. Astrophys., 642, A191
  * Bloomfield et al. (2013) Bloomfield J. K., Flanagan E. E., Park M., Watson S., 2013, JCAP, 1308, 010
  * Bose & Koyama (2016) Bose B., Koyama K., 2016, JCAP, 1608, 032
  * Bose et al. (2020a) Bose B., Winther H. A., Pourtsidou A., Casas S., Lombriser L., Xia Q., Cataneo M., 2020a, JCAP, 09, 001
  * Bose et al. (2020b) Bose B., Cataneo M., Tröster T., Xia Q., Heymans C., Lombriser L., 2020b, Mon. Not. Roy. Astron. Soc., 498, 4650
  * Bose et al. (2021) Bose B., et al., 2021, Mon. Not. Roy. Astron. Soc., 508, 2479
  * Brando et al. (2022) Brando G., Fiorini B., Koyama K., Winther H. A., 2022, JCAP, 09, 051
  * Brax & Valageas (2013) Brax P., Valageas P., 2013, Phys. Rev. D, 88, 023527
  * Brax & Valageas (2014) Brax P., Valageas P., 2014, Phys. Rev. D, 90, 023508
  * Brax et al. (2021) Brax P., Casas S., Desmond H., Elder B., 2021, Universe, 8, 11
  * Burgess (2004) Burgess C. P., 2004, Annals Phys., 313, 283
  * Burrage & Sakstein (2018) Burrage C., Sakstein J., 2018, Living Rev. Rel., 21, 1
  * Burrage et al. (2012) Burrage C., de Rham C., Heisenberg L., Tolley A. J., 2012, JCAP, 07, 004
  * Carrilho et al. (2022) Carrilho P., Carrion K., Bose B., Pourtsidou A., Hidalgo J. C., Lombriser L., Baldi M., 2022, Mon. Not. Roy. Astron. Soc., 512, 3691
  * Cataneo et al. (2015) Cataneo M., et al., 2015, Phys. Rev. D, 92, 044009
  * Cataneo et al. (2016) Cataneo M., Rapetti D., Lombriser L., Li B., 2016, JCAP, 12, 024
  * Cataneo et al. (2019) Cataneo M., Lombriser L., Heymans C., Mead A., Barreira A., Bose S., Li B., 2019, Mon. Not. Roy. Astron. Soc., 488, 2121
  * Cataneo et al. (2020) Cataneo M., Emberson J., Inman D., Harnois-Deraps J., Heymans C., 2020, Mon. Not. Roy. Astron. Soc., 491, 3101
  * Charmousis et al. (2012) Charmousis C., Copeland E. J., Padilla A., Saffin P. M., 2012, Phys.Rev.Lett., 108, 051101
  * Chevallier & Polarski (2001) Chevallier M., Polarski D., 2001, Int. J. Mod. Phys., D10, 213
  * Cooray & Sheth (2002) Cooray A., Sheth R. K., 2002, Phys. Rept., 372, 1
  * Creminelli & Vernizzi (2017) Creminelli P., Vernizzi F., 2017, Phys. Rev. Lett., 119, 251302
  * Creminelli et al. (2018) Creminelli P., Lewandowski M., Tambalo G., Vernizzi F., 2018, JCAP, 1812, 025
  * Creminelli et al. (2020) Creminelli P., Tambalo G., Vernizzi F., Yingcharoenrat V., 2020, J. Cosmology Astropart. Phys., 2020, 002
  * Cusin et al. (2018a) Cusin G., Lewandowski M., Vernizzi F., 2018a, JCAP, 04, 005
  * Cusin et al. (2018b) Cusin G., Lewandowski M., Vernizzi F., 2018b, JCAP, 04, 061
  * De Felice & Tsujikawa (2012) De Felice A., Tsujikawa S., 2012, JCAP, 1202, 007
  * Deffayet et al. (2010) Deffayet C., Pujolas O., Sawicki I., Vikman A., 2010, JCAP, 10, 026
  * Desmond & Ferreira (2020) Desmond H., Ferreira P. G., 2020, Phys. Rev. D, 102, 104060
  * Dvali et al. (2000) Dvali G., Gabadadze G., Porrati M., 2000, Phys.Lett., B485, 208
  * Euclid Collaboration et al. (2020) Euclid Collaboration et al., 2020, arXiv e-prints, p. arXiv:2010.11288
  * Ezquiaga & Zumalacárregui (2017) Ezquiaga J. M., Zumalacárregui M., 2017, Phys. Rev. Lett., 119, 251304
  * Fernandes et al. (2022) Fernandes P. G. S., Carrilho P., Clifton T., Mulryne D. J., 2022, Class. Quant. Grav., 39, 063001
  * Fiorini et al. (2021) Fiorini B., Koyama K., Izard A., Winther H. A., Wright B. S., Li B., 2021, JCAP, 09, 021
  * Frusciante & Perenon (2020) Frusciante N., Perenon L., 2020, Phys. Rept., 857, 1
  * Frusciante et al. (2016) Frusciante N., Papadomanolakis G., Silvestri A., 2016, JCAP, 07, 018
  * Frusciante et al. (2019) Frusciante N., Peirone S., Casas S., Lima N. A., 2019, Phys. Rev. D, 99, 063538
  * Giblin et al. (2019) Giblin B., Cataneo M., Moews B., Heymans C., 2019, Mon. Not. Roy. Astron. Soc., 490, 4826
  * Gleyzes et al. (2015) Gleyzes J., Langlois D., Mancarella M., Vernizzi F., 2015, JCAP, 08, 054
  * Gsponer & Noller (2022) Gsponer R., Noller J., 2022, Phys. Rev. D, 105, 064002
  * Gubitosi et al. (2013) Gubitosi G., Piazza F., Vernizzi F., 2013, JCAP, 1302, 032
  * Hassani & Lombriser (2020) Hassani F., Lombriser L., 2020, Mon. Not. Roy. Astron. Soc., 497, 1885
  * Hernández-Aguayo et al. (2022) Hernández-Aguayo C., Ruan C.-Z., Li B., Arnold C., Baugh C. M., Klypin A., Prada F., 2022, JCAP, 01, 048
  * Hinterbichler & Khoury (2010) Hinterbichler K., Khoury J., 2010, Phys.Rev.Lett., 104, 231301
  * Horndeski (1974) Horndeski G. W., 1974, Int.J.Theor.Phys., 10, 363
  * Hu & Sawicki (2007) Hu W., Sawicki I., 2007, Phys.Rev., D76, 064004
  * Hu et al. (2014) Hu B., Raveri M., Frusciante N., Silvestri A., 2014, Phys.Rev., D89, 103530
  * Huang (2016) Huang Z., 2016, Phys. Rev. D, 93, 043538
  * Ivezić et al. (2019) Ivezić v., et al., 2019, Astrophys. J., 873, 111
  * Joudaki et al. (2022) Joudaki S., Ferreira P. G., Lima N. A., Winther H. A., 2022, Phys. Rev. D, 105, 043522
  * Kennedy et al. (2017) Kennedy J., Lombriser L., Taylor A., 2017, Phys. Rev., D96, 084051
  * Kennedy et al. (2018) Kennedy J., Lombriser L., Taylor A., 2018, Phys. Rev., D98, 044051
  * Khan & Taylor (2022) Khan A., Taylor A., 2022, JCAP, 10, 075
  * Khoury & Weltman (2004) Khoury J., Weltman A., 2004, Phys.Rev.Lett., 93, 171104
  * Knabenhans et al. (2021) Knabenhans M., et al., 2021, Mon. Not. Roy. Astron. Soc., 505, 2840
  * Koyama (2018) Koyama K., 2018, Int. J. Mod. Phys., D27, 1848001
  * Laureijs et al. (2011) Laureijs R., et al., 2011,
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, Astrophys. J., 538, 473
  * Li & Efstathiou (2012) Li B., Efstathiou G., 2012, Mon. Not. Roy. Astron. Soc., 421, 1431
  * Linder (2003) Linder E. V., 2003, Phys. Rev. Lett., 90, 091301
  * Linder & Cahn (2007) Linder E. V., Cahn R. N., 2007, Astropart. Phys., 28, 481
  * Lombriser (2014) Lombriser L., 2014, Annalen Phys., 526, 259
  * Lombriser (2016) Lombriser L., 2016, JCAP, 11, 039
  * Lombriser & Lima (2017) Lombriser L., Lima N. A., 2017, Phys. Lett., B765, 382
  * Lombriser & Taylor (2015a) Lombriser L., Taylor A., 2015a, JCAP, 11, 040
  * Lombriser & Taylor (2015b) Lombriser L., Taylor A., 2015b, Phys.Rev.Lett., 114, 031101
  * Lombriser & Taylor (2016) Lombriser L., Taylor A., 2016, JCAP, 1603, 031
  * Lombriser et al. (2014) Lombriser L., Koyama K., Li B., 2014, JCAP, 03, 021
  * Lombriser et al. (2019) Lombriser L., Dalang C., Kennedy J., Taylor A., 2019, JCAP, 01, 041
  * Mancarella et al. (2022) Mancarella M., Kennedy J., Bose B., Lombriser L., 2022, Phys. Rev. D, 105, 023531
  * Martin (2012) Martin J., 2012, Comptes Rendus Physique, 13, 566
  * McManus et al. (2016) McManus R., Lombriser L., Peñarrubia J., 2016, JCAP, 11, 006
  * Mead et al. (2015) Mead A., Peacock J., Heymans C., Joudaki S., Heavens A., 2015, Mon. Not. Roy. Astron. Soc., 454, 1958
  * Mead et al. (2016) Mead A., Heymans C., Lombriser L., Peacock J., Steele O., Winther H., 2016, Mon. Not. Roy. Astron. Soc., 459, 1468
  * Mead et al. (2021) Mead A. J., Brieden S., Tröster T., Heymans C., 2021, MNRAS, 502, 1401
  * Melville & Noller (2020) Melville S., Noller J., 2020, Phys. Rev. D, 101, 021502
  * Noller (2020) Noller J., 2020, Phys. Rev. D, 101, 063524
  * Noller & Nicola (2019a) Noller J., Nicola A., 2019a, Phys. Rev. D, 99, 103502
  * Noller & Nicola (2019b) Noller J., Nicola A., 2019b, Phys. Rev. D, 99, 103502
  * Noller et al. (2014) Noller J., von Braun-Bates F., Ferreira P. G., 2014, Phys. Rev. D, 89, 023521
  * Noller et al. (2021) Noller J., Santoni L., Trincherini E., Trombetta L. G., 2021, JCAP, 01, 045
  * Pace et al. (2021) Pace F., Battye R., Bellini E., Lombriser L., Vernizzi F., Bolliet B., 2021, JCAP, 06, 017
  * Padilla (2015) Padilla A., 2015, Lectures on the Cosmological Constant Problem
  * Peebles (1980) Peebles P. J. E., 1980, The large-scale structure of the universe. Research supported by the National Science Foundation. Princeton, N.J., Princeton University Press, 1980. 435 p.
  * Peirone et al. (2018) Peirone S., Koyama K., Pogosian L., Raveri M., Silvestri A., 2018, Phys. Rev. D, 97, 043519
  * Perlmutter et al. (1999) Perlmutter S., et al., 1999, Astrophys. J., 517, 565
  * Pogosian & Silvestri (2016) Pogosian L., Silvestri A., 2016, Phys. Rev. D, 94, 104014
  * Puchwein et al. (2013) Puchwein E., Baldi M., Springel V., 2013, Mon.Not.Roy.Astron.Soc., 436, 348
  * Raveri et al. (2014) Raveri M., Hu B., Frusciante N., Silvestri A., 2014, Phys.Rev., D90, 043513
  * Reall (2021) Reall H. S., 2021, Phys. Rev. D, 103, 084027
  * Renevey et al. (2020) Renevey C., Kennedy J., Lombriser L., 2020, JCAP, 12, 032
  * Riess et al. (1998) Riess A. G., et al., 1998, Astron. J., 116, 1009
  * Sakstein & Jain (2017) Sakstein J., Jain B., 2017, Phys. Rev. Lett., 119, 251303
  * Sawicki & Bellini (2015) Sawicki I., Bellini E., 2015, Phys. Rev. D, 92, 084061
  * Schmidt (2010) Schmidt F., 2010, Phys.Rev., D81, 103002
  * Schmidt et al. (2009) Schmidt F., Lima M. V., Oyaizu H., Hu W., 2009, Phys.Rev., D79, 083518
  * Schmidt et al. (2010) Schmidt F., Hu W., Lima M., 2010, Phys.Rev., D81, 063005
  * Schneider et al. (2020) Schneider A., et al., 2020, JCAP, 04, 020
  * Silvestri et al. (2013) Silvestri A., Pogosian L., Buniy R. V., 2013, Phys. Rev. D, 87, 104015
  * Skordis et al. (2015) Skordis C., Pourtsidou A., Copeland E. J., 2015, Phys. Rev., D91, 083537
  * Sobral Blanco & Lombriser (2020) Sobral Blanco D., Lombriser L., 2020, Phys. Rev. D, 102, 043506
  * Spurio Mancini et al. (2019) Spurio Mancini A., et al., 2019, Mon. Not. Roy. Astron. Soc., 490, 2155
  * Spurio Mancini et al. (2022) Spurio Mancini A., Piras D., Alsing J., Joachimi B., Hobson M. P., 2022, Mon. Not. Roy. Astron. Soc., 511, 1771
  * Srinivasan et al. (2021) Srinivasan S., Thomas D. B., Pace F., Battye R., 2021, JCAP, 06, 016
  * Taddei et al. (2014) Taddei L., Catena R., Pietroni M., 2014, Phys. Rev. D, 89, 023523
  * Takahashi et al. (2012) Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, Astrophys. J., 761, 152
  * Takushima et al. (2015) Takushima Y., Terukina A., Yamamoto K., 2015, Phys. Rev., D92, 104033
  * Thomas (2020) Thomas D. B., 2020, Phys. Rev. D, 101, 123517
  * Thomas et al. (2022) Thomas D. B., Clifton T., Anton T., 2022, Scale-Dependent Gravitational Couplings in Parameterised Post-Newtonian Cosmology (arXiv:2207.14713)
  * Traykova et al. (2019) Traykova D., Bellini E., Ferreira P. G., 2019, JCAP, 08, 035
  * Tröster et al. (2021) Tröster T., et al., 2021, Astron. Astrophys., 649, A88
  * Vainshtein (1972) Vainshtein A., 1972, Phys.Lett., B39, 393
  * Weinberg (1989) Weinberg S., 1989, Rev.Mod.Phys., 61, 1
  * Will (1993) Will C. M., 1993, Theory and Experiment in Gravitational Physics
  * Will (2014) Will C. M., 2014, Living Rev. Rel., 17, 4
  * Winther et al. (2017) Winther H. A., Koyama K., Manera M., Wright B. S., Zhao G.-B., 2017, JCAP, 08, 006
  * Winther et al. (2019) Winther H., Casas S., Baldi M., Koyama K., Li B., Lombriser L., Zhao G.-B., 2019, Phys. Rev. D, 100, 123540
  * Wright et al. (2022) Wright B. S., Gupta A. S., Baker T., Valogiannis G., 2022, Hi-COLA: Fast, approximate simulations of structure formation in Horndeski gravity (arXiv:2209.01666)
  * Zhao (2014) Zhao G.-B., 2014, Astrophys. J. Suppl., 211, 23
  * de Felice et al. (2011) de Felice A., Kobayashi T., Tsujikawa S., 2011, Physics Letters B, 706, 123
  * de Rham & Melville (2018) de Rham C., Melville S., 2018, Phys. Rev. Lett., 121, 221101
  * de Rham & Tolley (2020) de Rham C., Tolley A. J., 2020, Phys. Rev. D, 102, 084048
  * de Rham et al. (2021) de Rham C., Melville S., Noller J., 2021, JCAP, 08, 018

##  Appendix A Linear Poisson modifications

###  A.1 Horndeski in the $\alpha$-basis

The linear modification to the Poisson equation in the Horndeski class of
EFTofDE, under the quasi-static approximation, and in the $\alpha$-basis,
$\\{H,\alpha_{M},\alpha_{B},\alpha_{K},\alpha_{T}\\}$, is given by

$\mu(k,a)=\frac{2}{\kappa^{2}}\frac{f_{1}(a)+f_{2}(a)a^{2}/k^{2}}{f_{3}(a)+f_{4}(a)a^{2}/k^{2}}\,,$
(60)

where the constituent functions are given by

$\displaystyle f_{1}=$ $\displaystyle B_{2}C_{3}-C_{1}B_{3}\,,$ (61)
$\displaystyle f_{2}=$ $\displaystyle B_{2}C_{\pi}\,,$ (62) $\displaystyle
f_{3}=$ $\displaystyle A_{1}[B_{3}C_{2}-B_{1}C_{3}(k,a)]$
$\displaystyle+A_{2}[B_{1}C_{1}-B_{2}C_{2}]$
$\displaystyle+A_{3}[B_{2}C_{3}-B_{3}C_{1}]\,,$ (63) $\displaystyle f_{4}=$
$\displaystyle[A_{3}B_{2}-A_{1}B_{1}]C_{\pi}\,.$ (64)

Finally, we give the $A$, $B$ and $C$ functions in terms of the $\alpha$-basis
as (Pogosian & Silvestri, 2016)

$\displaystyle A_{1}=$ $\displaystyle 2M^{2}\,,$ (65) $\displaystyle A_{2}=$
$\displaystyle\alpha_{B}HM^{2}\,,$ (66) $\displaystyle A_{3}=$ $\displaystyle
0\,,$ (67) $\displaystyle B_{1}=$ $\displaystyle-\frac{1}{c_{T}^{2}}\,,$ (68)
$\displaystyle B_{2}=$ $\displaystyle 1\,,$ (69) $\displaystyle B_{3}=$
$\displaystyle\frac{(-\alpha_{M}+\alpha_{T})H}{c_{T}^{2}}\,,$ (70)
$\displaystyle C_{1}=$ $\displaystyle-B_{3}M^{2}c_{T}^{2}\,,$ (71)
$\displaystyle C_{2}=$ $\displaystyle\frac{A_{2}}{2}\,,$ (72) $\displaystyle
C_{3}=$ $\displaystyle
c+\frac{HM^{2}}{2}\Big{[}-2\alpha_{T}H+\alpha_{M}^{2}c_{T}^{2}H+aH\alpha_{B}^{\prime}+aH\alpha_{M}^{\prime}$
$\displaystyle+a\alpha_{T}H\alpha_{M}^{\prime}+2a\alpha_{T}H^{\prime}+a^{2}\alpha_{T}^{\prime}H^{\prime}+\alpha_{B}[(1+\alpha_{M})H+aH^{\prime}]$
$\displaystyle+\alpha_{M}[H(1-\alpha_{T}+2a\alpha_{T}^{\prime})+ac_{T}^{2}H^{\prime}]+a^{2}H\alpha_{T}^{\prime\prime}\Big{]}\,,$
(73) $\displaystyle C_{\pi}=$
$\displaystyle-\frac{1}{4}aH\Big{[}12cH^{\prime}+HM^{2}(6\alpha_{M}^{2}c_{T}^{2}HH^{\prime}$
$\displaystyle+6\alpha_{B}(2a(H^{\prime})^{2}+H[(4+\alpha_{M})H^{\prime}+aH^{\prime\prime}])$
$\displaystyle+\alpha_{M}(c_{T}^{2}(12a(H^{\prime})^{2}-\bar{R}^{\prime})+6H(2(2c_{T}^{2}+a\alpha_{T}^{\prime})H^{\prime}+ac_{T}^{2}H^{\prime\prime}))$
$\displaystyle+a[12(\alpha_{T}+a\alpha_{T}^{\prime})(H^{\prime})^{2}$
$\displaystyle-\alpha_{T}^{\prime}\bar{R}^{\prime}+6H(H^{\prime}(\alpha_{B}^{\prime}+c_{T}^{2}\alpha_{M}^{\prime}+5\alpha_{T}^{\prime}+a\alpha_{T}^{\prime\prime})+a\alpha_{T}^{\prime}H^{\prime\prime})])\Big{]}\,,$
(74)

where $c_{T}^{2}=(1+\alpha_{T})$ and $\bar{R}$ is the background Ricci scalar.
The $c(a)$ parameter (Equation 22) in the $\alpha$-basis is given by

$\displaystyle c(a)$
$\displaystyle=-M^{2}\Big{[}\frac{3H_{0}^{2}\Omega_{m,0}}{2a^{3}\kappa^{2}M^{2}}$
$\displaystyle+\frac{1}{2}H\Big{(}aH^{\prime}[(2+\alpha_{M})c_{T}^{2}+a\alpha_{T}^{\prime}]$
$\displaystyle+H(c_{T}^{2}[(\alpha_{M}-1)\alpha_{M}+a\alpha_{M}^{\prime}]+2a\alpha_{M}\alpha_{T}^{\prime}+a^{2}\alpha_{T}^{\prime\prime})\Big{)}\Big{]}\,,$
(75)

where we have used $\rho_{m}=3H_{0}^{2}\Omega_{m,0}/(\kappa^{2}a^{3})$, with
$\Omega_{m,0}$ being the matter density fraction today and $H_{0}=H(a=1)$ is
the Hubble constant.

We note that in our code we make the redefinition
$M^{2}\kappa^{2}=M^{2}/m_{0}^{2}\rightarrow M^{2}$ where $m_{0}^{2}($ is the
Planck mass. Further, we comment on the flexibility offered here. One can
choose to specify any two of $\\{\alpha_{M},M^{2},H\\}$. If $H$ is specified
then either $\alpha_{M}$ or $M^{2}$ must also be specified, with the third
function given by the relation in Equation 25. If $H$ is not specified, then
we must solve the Friedmann equations to obtain $H$. As a default in our code,
$H$ is specified and it is assumed that the specified expressions for
$\alpha_{M}$ and $M^{2}$ are consistent with Equation 25.

We can also take the small scale ($k\rightarrow\infty$) limit of Equation 60
to get a simpler expression valid at scales where the QS is a safer
approximation and for models exhibiting negligible scale dependence in the
linear growth. This is given by

$\mu_{\infty}=\frac{1}{M^{2}\kappa^{2}}\Big{(}1+\alpha_{T}+\beta_{\xi}\Big{)}\,,$
(76)

where

$\beta_{\xi}=\frac{2}{c_{s}^{2}\alpha}\Big{(}c_{T}^{2}\frac{\alpha_{B}}{2}+\alpha_{M}-\alpha_{T}\Big{)}^{2}\,$
(77)

with $c_{s}^{2}$ and $\alpha$ and $\alpha$ given by Equation 29 and Equation
30.

###  A.2 Example: Hu-Sawicki $f(R)$

In this section we derive the relevant EFTofDE parameters and linear Poisson
modification for the Hu-Sawicki form of $f(R)$ gravity (Hu & Sawicki, 2007).
This exercise is also performed in the GtoPT notebook provided in the ACTio-
ReACTio github repository.

The action in $f(R)$ gravity is given by

$\displaystyle S=$ $\displaystyle\int
d^{4}x\sqrt{-g}\frac{1}{2\kappa^{2}}\Big{(}R+f(R)\Big{)}$
$\displaystyle\approx$ $\displaystyle\int
d^{4}x\sqrt{-g}\frac{1}{2\kappa^{2}}\Big{(}R+f(\bar{R})+f_{R}(\bar{R})(R-\bar{R})\Big{)}$
$\displaystyle=$ $\displaystyle\int
d^{4}x\sqrt{-g}\frac{1}{2\kappa^{2}}\Big{(}(1+f_{R})R+f-f_{R}\bar{R}\Big{)}\,,$
(78)

where $f_{R}=df(R)/dR$ and we have performed a Taylor expansion in the second
line. We can then map this action onto the functions given in Equation 2
together with an identification of the scalar degree of freedom
$\phi\equiv(1+f_{R})/\kappa^{2}$ (de Felice et al., 2011):

$\displaystyle
G_{2}=-\frac{1}{2\kappa^{2}}(\bar{R}f_{R}-f),\,\,\,\,G_{3}=0,\,\,\,\,\mathrm{and}\,\,\,\,G_{4}=\frac{1}{2\kappa^{2}}(1+f_{R})\,.$
(79)

If we now write down the action in the ADM formalism and use the Gauss-Codazzi
relation, we can compare to Equation 18 and Equation 19 to get

$\displaystyle\Omega=(1+f_{R}),\,\Lambda=\frac{1}{2\kappa^{2}}(f-\bar{R}f_{R}),\,\
c=\bar{M}^{2}_{2}=\bar{M}_{1}^{3}=M_{2}^{4}=0\,.$ (80)

Using Equation 25 \- Equation 27 we have

$\displaystyle\alpha_{M}$ $\displaystyle=\frac{af_{R}^{\prime}}{1+f_{R}}\,,$
(81) $\displaystyle\alpha_{T}$ $\displaystyle=0\,,$ (82)
$\displaystyle\alpha_{B}$ $\displaystyle=-\frac{af_{R}^{\prime}}{1+f_{R}}\,,$
(83) $\displaystyle\alpha_{K}$ $\displaystyle=0\,,$ (84)

with

$M^{2}=\frac{(1+f_{R})}{\kappa^{2}}\,,\\\ $ (85)

where a prime denotes a scale factor derivative. When substituting into the
expressions in Appendix A.1 we get the following solution for $\mu$ (Equation
60)

$\mu=\frac{1}{1+f_{R}}\Big{[}1+\Big{(}\frac{k}{a}\Big{)}^{2}\frac{1}{3\tilde{\Pi}(k,a)}\Big{]}\,,$
(86)

where

$\tilde{\Pi}(k,a)=\left(\frac{k}{a}\right)^{2}+(1+f_{R})\frac{\bar{R}_{f}}{3}\,,$
(87)

and

$\bar{R}_{f}\equiv\frac{d\bar{R}}{df_{R}}=f_{RR}^{-1}\Big{(}=\frac{\bar{R}^{\prime}}{f_{R}^{\prime}}\Big{)}\,.$
(88)

In the Hu-Sawicki model we have the following choice for $f(R)$

$f(R)=-m^{2}\frac{c_{1}(R/m^{2})^{n}}{c_{2}(R/m^{2})^{n}+1}\,,$ (89)

where in this work we set the index $n=1$ and the mass scale $m^{2}$, $c_{1}$
and $c_{2}$ are free parameters to be constrained by data. Taking the
derivative of Equation 89 with respect to $R$ and the high curvature limit
($R\gg m^{2}$) gives

$f_{R}=-\frac{c_{1}}{c_{2}^{2}}\Big{(}\frac{m^{2}}{R}\Big{)}^{2}\,.$ (90)

By rearranging this equation and evaluating at the background level at $a=1$
(today), we can apply the following standard reparameterisation

$\frac{c_{1}}{c_{2}^{2}}=-\bar{f}_{R0}\Big{(}\frac{\bar{R}_{0}}{m^{2}}\Big{)}^{2}\,,$
(91)

where $\bar{f}_{R0}$ is the background value of $f_{R}$ evaluated today and is
a free parameter governing the level of modification to $\Lambda$CDM at the
level of structure formation. Substituting Equation 91 into Equation 90 gives

$f_{R}=\bar{f}_{R0}\Big{(}\frac{\bar{R}_{0}}{R}\Big{)}^{2}\,.$ (92)

Further, we have

$f_{RR}=-2\bar{f}_{R0}\Big{(}\frac{\bar{R}_{0}}{R}\Big{)}^{2}\frac{1}{R}\,.$
(93)

Using the background expression for $f_{RR}$ in Equation 88 and substituting
into Equation 87 gives

$\tilde{\Pi}(k,a)=\left(\frac{k}{a}\right)^{2}+(1+f_{R})\frac{1}{6|f_{R0}|}\frac{\bar{R}^{3}}{\bar{R}_{0}^{2}}\,.$
(94)

Now if we approximate the background to be close to $\Lambda$CDM, as supported
by observations and by construction for $|f_{\rm R0}|\ll 1$, we have

$\bar{R}\approx
3\frac{H_{0}^{2}}{a^{3}}\Big{(}\Omega_{m,0}+4a^{3}\Omega_{\Lambda,0}\Big{)}\,,$
(95)

where $\Omega_{\Lambda,0}=1-\Omega_{m,0}$ for a flat $\Lambda$CDM universe.
Taking $a=1$ we have the curvature today

$\bar{R}_{0}\approx 3H_{0}^{2}\Big{(}4-3\Omega_{m,0}\Big{)}\,.$ (96)

Finally substituting $\bar{R}$ and $\bar{R}_{0}$ in Equation 94 we get the
expression for $\mu$ as it appears in ACTio-ReACTio (Bose et al., 2020b)

$\mu=1+\Big{(}\frac{k}{a}\Big{)}^{2}\frac{1}{3\Pi(k,a)}\,,$ (97)

with

$\displaystyle\Pi(k,a)=$
$\displaystyle\left(\frac{k}{a}\right)^{2}+\frac{\Xi(a)^{3}}{2f_{0}(4-3\Omega_{\rm
m,0})^{2}},$ (98) $\displaystyle\Xi(a)=$ $\displaystyle\frac{\Omega_{\rm
m,0}+4a^{3}\Omega_{\Lambda,0}}{a^{3}},$ (99)

where $f_{0}=|{f}_{\rm R0}|/H_{0}^{2}$.

We make the crucial note here that in the derivation above we have over-
constrained our system. Namely we have specified all of $\Omega,c$ as well as
set $H=H_{\rm\Lambda CDM}$. If we use Equation 81-Equation 85 together with
$H=H_{\rm\Lambda CDM}$ we find that $c\neq 0$ and we don’t recover Equation
97. This follows directly from the fact that $H_{\rm\Lambda CDM}$ is not an
exact solution for the Friedmann equations in $f(R)$ which implicitly assumes
$c=0$.

In our code we give the option to over constrain by specifying $c(a)$.
Alternatively, one can code in the Friedmann equations and solve for $H(a)$.
This of course increases computational inefficiency. To partially alleviate
this issue, we can also place the additional constraint on $c_{s}^{2}(a)$
instead of $c$ which has well motivated physical priors (see text around
Equation 29). The relationship between $c_{s}^{2}(a)$ and $c(a)$ is derived
from Equation 29 and Equation 75. It is given explicitly in the GtoPT notebook
as well as the Actio-Reactio source code.

The derivation of the 2nd and 3rd order modifications to the Poisson equation
(see Equation 14) from the ADM decomposed action requires us to go to higher
order in the metric perturbations. We do not do this here as we omit these
corrections from our code due to computational difficulty and the low level of
impact they have on the final nonlinear power spectrum (see Section 5).
However, in the provided Mathematica notebook, GtoPT, one can go from a
specified $G_{i}$ of the Horndeski Lagrangian to $\mu$, $\gamma_{2}$ and
$\gamma_{3}$ following the map given in Bose & Koyama (2016). We provide a
number of examples in that notebook and refer the reader to Bose et al.
(2020b); Bose & Koyama (2016) for the forms of $\mu$, $\gamma_{2}$ and
$\gamma_{3}$ in DGP and Hu-Sawicki $f(R)$ gravity.

##  Appendix B Nonlinear Poisson modifications

###  B.1 Exact forms

We provide the exact forms for the nonlinear modification to the Poisson
equation (see Equation 13) in DGP and $f(R)$ gravity, which are reproduced
from Bose et al. (2020b).

The modification in DGP is given by (Schmidt et al., 2010)

$\mathcal{F}_{\rm DGP}=\frac{2}{3\beta(a)}\frac{\sqrt{1+s^{3}}-1}{s^{3}},$
(100)

where111111We note a typo appearing in Eq. C7 of Bose et al. (2020b) where
$\delta$ should have been $(\delta+1)$.

$s=\left[\frac{2\Omega_{m,0}(\delta+1)}{9a^{3}\beta(a)^{2}\Omega_{rc}}\right]^{\frac{1}{3}},$
(101)

$\delta$ being the nonlinear over-density given by

$\delta=y^{-3}(1+\delta_{i})-1\,,$ (102)

with $\delta_{i}$ being the initial over-density and

$y\equiv\frac{R_{\rm TH}/a}{R_{i}/a_{i}}\,,$ (103)

$R_{\rm TH}$ and $R_{i}$ being the physical halo top-hat radius at the target
scale factor $a$ and the initial scale factor $a_{i}$ respectively.
$\Omega_{\rm rc}\equiv 1/(4H_{0}^{2}r_{c}^{2})$ where $r_{c}$ is the cross-
over scale and is the free parameter of the theory governing the level of
modification. Finally, $\beta(a)$ is given by

$\beta(a)\equiv 1+\frac{H}{H_{0}}\frac{1}{\sqrt{\Omega_{\rm
rc}}}\left(1+\frac{aH^{\prime}}{3H}\right)\,.$ (104)

The fully nonlinear modification in Hu-Sawicki $f(R)$ is given by (Lombriser
et al., 2014)

$\mathcal{F}_{\rm fR}={\rm
min}\left[O-O^{2}+\frac{O^{3}}{3},\frac{1}{3}\right],$ (105)

where

$O=\frac{f_{0}y_{h}a(3\Omega_{m,0}-4)^{2}}{\Omega_{m,0}(R_{i}/a_{i})^{2}}\times\left[\tilde{G}(y_{\rm
env})-\tilde{G}(y_{\rm h})\right]\,,$ (106)

and

$\tilde{G}(y)=\left[\frac{\Omega_{m,0}}{(ya)^{3}}+4-4\Omega_{m,0}\right]^{-2}\,,$
(107)

where $y_{h}$ is the quantity solved for using $f(R)$ halos whereas $y_{\rm
env}$ is that solved for in the environment, which is approximated by
performing the same calculation but with $f_{0}=0$.

###  B.2 nPPF forms

We also reproduce the nPPF expressions for both of these theories from
Lombriser (2016). In DGP we have the following values for the $p_{i}$
parameters in Equation 39

$\displaystyle p_{1}$ $\displaystyle=2,\qquad p_{2}=1,\qquad
p_{3}=\frac{3}{2},$ $\displaystyle p_{4}(a)$
$\displaystyle=2\Big{(}\frac{\Omega_{\rm m,0}}{4\Omega_{\rm
rc}}\frac{1}{9\beta(a)^{2}}\Big{)}^{1/3},\qquad p_{5}=-1,\qquad p_{6}=0,$
$\displaystyle p_{7}$ $\displaystyle=0\,,$ (108)

which reproduce Equation 100 exactly. Note if using Equation 37 we simply set
$p_{2}=\frac{1}{3\beta(a)}$.

On the other hand, the Hu-Sawicki $f(R)$ parameterisation is not exact but is
closely matched by the following parameters in the screening regime
(Lombriser, 2016) for a given choice of $p_{1}$ (using Equation 39)

$\displaystyle p_{2}$ $\displaystyle=1,\qquad p_{3}=7,$ $\displaystyle p_{4}$
$\displaystyle=2\Omega_{\rm m,0}^{1/3}\Big{[}(\Omega_{\rm m,0}+4(1-\Omega_{\rm
m,0})^{-2}\frac{p_{1}}{3|f_{\rm R0}|}\Big{]}^{1/p_{3}},$ $\displaystyle p_{5}$
$\displaystyle=-1,\qquad p_{6}=\frac{2}{3p_{3}},\qquad p_{7}=-\frac{6}{7}\,,$
(109)

where we used $\alpha=1/(n+1)=0.5$ (Lombriser et al., 2014) in Equation 5.6 of
Lombriser (2016). Again, if using Equation 37 we set $p_{2}=\frac{1}{3}$.

###  B.3 Comparisons

Here we provide some comparisons of the approximate expressions for
$\mathcal{F}$ given by the nPPF model (Equation 39) and the Erf model
(Equation 41) against the exact expressions in DGP (Equation 100) and $f(R)$
(Equation 105). Since the nPPF form is exact for DGP, we only compare it in
the $f(R)$ case. Unless otherwise stated, the fits are performed as described
in section 5 and shown in Table 3.

####  B.3.1 DGP

In Figure 8 we show the nonlinear modification to the Poisson equation,
$1+\mathcal{F}$, in DGP for $\Omega_{\rm rc}=0.25$ and $\Omega_{\rm rc}=0.01$
as a function of top-hat radius parameter $y_{\rm h}$. We plot the exact
solution given by Equation 100 given as solid curves to the best fit Erf model
given as dashed curves. We see an additional redshift dependence of the
screening scale becoming important for high $z$. As modifications to GR are
expected to be small at high redshift, this deviation may not be so important,
which is supported by Figure 5. We have performed a fit of this redshift
dependence and find it behaves very well as a power law, with an
$\mathcal{O}(0.1)$ exponent (specifically $\leq 0.15$), which is first of all
small and second of all degenerate with $q_{2}$ and $q_{3}$, meaning the model
likely has sufficient freedom to very well capture a DGP type of modification
to gravity without biasing cosmological or gravitational constraints.

####  B.3.2 Hu-Sawicki $f(R)$

Here we check that the nPPF and Erf models can qualitatively reproduce the
exact form of $\mathcal{F}$ (Equation 105) across all scales, masses and
environments for Hu-Sawicki $f(R)$ gravity.

Before showing the results, we make a note on the best fitting parameters for
the Erf model. We find that the fits in Table 3, performed by fitting the
exact prediction for the reaction $\mathcal{R}$, do not give a very good
agreement when comparing to the exact form of $\mathcal{F}$. In particular, we
find that the mass dependence parameter, $q_{2}$, seems to be underestimated
when fitting $\mathcal{R}$. This parameter dictates the left hand slope in the
contour plots in this section. Such a discrepancy may be due to a number of
factors including a missing redshift dependency, the details of the fit,
degeneracies with $q_{3}$ and failings of the power law description. We find a
better by-eye fit for $\mathcal{F}_{\rm Erf}$ across redshifts $z=0,1,4$ and
all values of $f_{\rm R0}$ is $q_{2}=0.85$. All other parameters are as in
Table 3 unless otherwise stated.

In Figure 9 we show $1+\mathcal{F}$ for the exact (top panels), the Erf
(middle panels) and the nPPF (lower panels) cases with $|f_{\rm R0}|=10^{-5}$,