Matter Power Spectra in Modified Gravity: A Comparative Study of Approximations and $N$-Body Simulations

Abstract

Testing gravity and the concordance model of cosmology, $\Lambda$CDM, at large scales is a key goal of this decade’s largest galaxy surveys. Here we present a comparative study of dark matter power spectrum predictions from different numerical codes in the context of three popular theories of gravity that induce scale-independent modifications to the linear growth of structure: nDGP, Cubic Galileon and K-mouflage. In particular, we compare the predictions from full $N$-body simulations, two $N$-body codes with approximate time integration schemes, a parametrised modified $N$-body implementation and the analytic halo model reaction approach. We find the modification to the $\Lambda$CDM spectrum is in $2\%$ agreement for $z\leq 1$ and $k\leq 1~{}h/{\rm Mpc}$ over all gravitational models and codes, in accordance with many previous studies, indicating these modelling approaches are robust enough to be used in forthcoming survey analyses under appropriate scale cuts. We further make public the new code implementations presented, specifically the halo model reaction K-mouflage implementation and the relativistic Cubic Galileon implementation.

1 Introduction

Observations of cosmological large-scale structure (LSS) offer a unique laboratory in which to test the concordance cosmological model, $\Lambda$CDM, which assumes General Relativity (GR). Such experiments are highly motivated. Indeed, the nature of the cold dark matter (CDM) and the constant dark energy ($\Lambda$) components, constituting 95% of the Universe’s total energy density (see for example Riess et al., 1998; Perlmutter et al., 1999; Alam et al., 2021; Aghanim et al., 2020), remains elusive. Moreover, $\Lambda$CDM’s inability to reconcile principles of GR with quantum mechanics points to the need for a more unified theory (see Bernardo et al., 2022, for a recent review on gravitational approaches to the cosmological constant problem). These gaps in our understanding motivate the investigation into alternative theories beyond $\Lambda$CDM. By exploring these new frontiers, we hope to uncover a more comprehensive picture of the Universe, potentially leading to groundbreaking insights into its origin, evolution, and ultimate fate.

This decade will provide an immense opportunity for such insights through the efforts of some of the biggest scientific collaborations to date. These include the European Space Agency’s Euclid mission (Barroso et al., 2024), the Vera Rubin Observatory (LSST Dark Energy Science Collaboration, 2012; Ivezić et al., 2019) (LSST), the Dark Energy Survey (Albrecht et al., 2006; Abbott et al., 2016), the Nancy Grace Roman Space Telescope (Akeson et al., 2019) and the Dark Energy Spectroscopic Instrument (Levi et al., 2019). For instance, Euclid and LSST will be measuring up to order 1 billion galaxy shapes (Barroso et al., 2024; Ivezić et al., 2019), 2 orders of magnitude more than previous surveys (see for example Hildebrandt et al., 2017). This means the statistical precision of its resulting weak lensing measurements, such as cosmic shear, will be roughly the same order of magnitude better than previous observations, providing a potentially brilliant probe for new physics.

Consistency tests of $\Lambda$CDM are a primary goal, but these missions are also charged with investigating if there is any statistical preference for new physics. Such beyond-consistency tests require theoretical modelling of any new physics we wish to test. In particular, a key task is to theoretically model key statistical cosmological quantities over a very wide range of physical scales. The 2-point correlation function, or its Fourier analog, the power spectrum, of the cosmological matter distribution is one such summary statistic. At small physical scales, where we have many more galaxy pairs, the measured statistics will be far more precise, potentially providing a heightened signal of any new physics. It should be kept in mind that this work only considers the matter power spectrum, which is a key ingredient for cosmic shear weak lensing analyses.

The small scale precision measurements of forthcoming surveys has forced ambitious accuracy demands on such theoretical predictions (for example $\mathcal{O}(1\%)$ accuracy on the matter power spectrum Hearin et al., 2012; Ivezić et al., 2019; Martinelli et al., 2021). This requires careful consideration of scale cuts. Most Euclid forecasts (Blanchard et al., 2020; Bonici et al., 2023; Frusciante et al., 2023; Casas et al., 2023) consider a ‘pessimistic’ and ‘optimistic’ scale cut in harmonic space, corresponding to a maximum angular multipole of $\ell=1500$ and $\ell=5000$, with the precise value of these cuts in Fourier mode, or $k$-space, varying with redshift. In contrast, LSST applies scale cuts in real space. While percent-level accuracy remains a desirable goal, definitive accuracy up to a specific $k$-cut is only achievable through full parameter inference. This nuanced approach is essential to fully harness the power of these surveys for an exquisite and reliable test of gravity and cosmology.

For these reasons, the community has sought to accurately model these small, nonlinear scales in the matter power spectrum, for beyond-$\Lambda$CDM scenarios. To this end, many methods have been developed to provide such predictions. $N$-body simulations provide our most accurate predictions, and have been extended to many models beyond-$\Lambda$CDM (see for example Li et al., 2012, 2013a; Li et al., 2013b; Puchwein et al., 2013; Llinares et al., 2014; Ruan et al., 2022; Hernández-Aguayo et al., 2022; Hassani et al., 2019; Hassani et al., 2020; Christiansen et al., 2023). This accuracy comes at a large computational cost, making this method inappropriate for expensive data-theory comparisons where we wish to sample a large cosmological and gravitational parameter space. One can alleviate this cost to some extent through approximate methods. For example, Comoving Lagrangian Acceleration (COLA) (Tassev et al., 2013; Howlett et al., 2015) is an $N$-body method that provides a balance between accuracy and speed by reducing the time steps in particle evolution through the perturbative modelling of large scale physics. This method has also been extended to many alternatives to $\Lambda$CDM (Winther et al., 2017; Wright et al., 2023; Brando et al., 2023).

While being faster, COLA methods are still too slow to use directly in data analyses. Despite the computational cost, simulation methods are essential in bench-marking or constructing faster predictive pipelines, such as emulators (Arnold et al., 2021; Harnois-Déraps et al., 2022; Ramachandra et al., 2021; Fiorini et al., 2023; Nouri-Zonoz et al., 2024) or analytic models (Zhao, 2014; Mead et al., 2016). The halo model reaction (Cataneo et al., 2019) is one such analytic method, which can provide a high accuracy at a fraction of the time cost and is theoretically general, allowing its extension to many models of cosmology.

This paper is dedicated to assessing the consistency of these different methods for a few representative beyond-$\Lambda$CDM models of cosmological relevance. The models we consider are the DGP braneworld model (Dvali et al., 2000), the Cubic Galileon model (Nicolis et al., 2009) and the K-mouflage model (Babichev et al., 2009). This work runs in a similar vein to the code comparison project of Ref. Winther et al. (2015), updating the exercise, nearly a decade later, to account for improvements in the codes and methods, as well as approximations and new theoretical models and phenomenology. Such an assessment is vital in modelling the theoretical uncertainty or delimiting the scales of validity of the method under consideration, which will play an important role in forthcoming surveys (Audren et al., 2013; Baldauf et al., 2016). We also present an extension of the halo model reaction code, ReACT, which includes the specific K-mouflage model of gravity considered in this paper.

We outline the paper as follows: In section 2 we briefly introduce the different beyond-$\Lambda$CDM models we consider. In section 3 we outline the different methods we will compare, highlighting the key differences between them and the various approximations they employ. In section 4 we present matter power spectrum boost comparisons of the different methods. We present our conclusions in section 5.

1.1 Notation and conventions

In this work we will use the following definitions and conventions:

1. We use a metric signature of $(-,+,+,+)$.

2. We work in units where $c=\hbar=1$.

3. Jordan frame quantities appear with a hat, e.g. $\hat{q}$.

4. The Planck mass is denoted as $M_{\rm pl}=(8\pi G_{\rm N})^{-1}$, where $G_{\rm N}$ is Newton’s constant.

5. Overdots denote derivatives with respect to cosmic time $t$.

6. Primes denote derivatives with respect to the natural logarithm of the scale factor, $\ln a$, unless otherwise stated.

7. Quantities with a ‘$0$’ subscript denote their value at $z=0$.

8. The canonical scalar field kinetic energy is $X\equiv-(\partial\phi)^{2}/2$.

2 Gravity Beyond General Relativity

The simplest, viable class of alternatives to $\Lambda$CDM can be found by adding a single extra scalar degree of freedom, $\phi$, to GR. Under some basic constraints, such as second-order equations of motion (a generic condition to avoid unbounded negative energies) and four spacetime dimensions, the well-studied Horndeski Lagrangian encompasses all possible scalar-tensor theories with minimally coupled matter (Horndeski, 1974; Deffayet et al., 2011; Kobayashi, 2019). If we accept the speed of light to be the same as gravitational wave propagation, in accordance with the observation of gamma ray burst GRB 170817A (Goldstein et al., 2017) and merger signal of GW170817 (Abbott et al., 2017), the Horndeski Lagrangian reduces to (Lombriser & Taylor, 2016; 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; Quartin et al., 2023)

$$\mathcal{L}_{\rm H}=G_{4}(\phi)\,R+G_{2}(\phi,X)-G_{3}(\phi,X)\Box\phi\,,$$ (1)

where $R$ is the Ricci curvature scalar, $\Box$ is the D’Alembert operator and each $G_{i}(\phi,X)$, $i=2,3,4$ is a free function of the scalar field $\phi$ and its canonical kinetic term $X$. Note that the $G_{4}$ operator is a function of $\phi$ only.

Besides modifying the expansion history of the Universe, modified gravity theories also leave an impact on the growth of structure (see Hou et al., 2023, for a review). This is generally understood by considering linear perturbations on top of a homogeneous and isotropic Friedmann-Lemaître-Robertson-Walker background given by the following line element

$$\mathrm{d}s^{2}=-\left(1+2\Psi\right)\mathrm{d}t^{2}+a^{2}(t)\left(1-2\Phi\right)\delta_{ij}\mathrm{d}x^{i}\mathrm{d}x^{j},$$ (2)

where $\Phi$ is the usual Poisson potential in Newtonian gravity that captures perturbations in the spatial sector of the metric, while $\Psi$ is a gravitational potential corresponding to perturbations in the time-like sector of the line element.

The linear evolution of perturbations of modified gravity theories given by Equation 1 has been thoroughly studied by many different works in the literature (see for example Hu et al., 2014; Zumalacárregui et al., 2017; Frusciante & Perenon, 2020). Within the quasi-static approximation (Sawicki & Bellini, 2015; Winther & Ferreira, 2015b; Pace et al., 2021), the effects of modified gravity on the linear growth of structure in the Universe are encoded in a time- and scale-dependent effective gravitational constant

$$G_{\rm eff,L}(k,a)=G_{\rm N}\left[1+\frac{\Delta G_{\rm eff,L}(k,a)}{G_{N}}\right],$$ (3)

where $k$ is the Fourier mode. In this work, we only consider theories where the linear modification is scale-independent and so we drop the dependency on $k$ for $G_{\rm eff,L}$, where L refers to a linear theory prediction. The Poisson equation at large scales is then written in Fourier space as

$$k^{2}\Phi(k,a)=4\pi G_{\rm eff,L}(a)a^{2}\bar{\rho}_{\mathrm{m}}(a)\delta_{\rm m}(k,a),$$ (4)

where $\bar{\rho}_{\rm m}$ is the background matter density, and $\delta_{\rm m}$ is the corresponding linear matter perturbation.

Another requirement for this class of theories is the inclusion of a theoretical mechanism that prevents large modifications in environments where GR-like physics has been well confirmed by experiment (see Will, 2014; Belgacem et al., 2019, for example). Such mechanisms are known as screening mechanisms (see Brax et al., 2021, for a recent review and experimental tests). The screened environments are thus small scale, dense environments. This means that the modification to Newton’s constant, more generally written as

$$G_{\rm eff}(k,a)=G_{\rm N}\left[1+\frac{\Delta G_{\rm eff}(k,a)}{G_{N}}\right],$$ (5)

requires the condition that $\lim_{k\to\infty}G_{\rm eff}(k,a)\to G_{\rm N}$. In this case $G_{\rm eff}(k,a)$ is the effective gravitational constant valid at all scales - both linear and nonlinear - and it necessarily depends on scale as well as time. In this work we will meet two such screening mechanisms which satisfy this condition: the Vainshtein mechanism and K-mouflage screening.

Returning to Equation 1, we will consider three choices for the Lagrangian functions, each having very particular phenomenological features, including different screening mechanisms and cosmological backgrounds. Where a choice exists, we will give their Lagrangians in the Einstein frame where $G_{4}(\phi)=M_{\rm pl}^{2}/2$, with metric $g_{\mu\nu}$. In this frame the ‘pure gravity’ part of the action resembles the Einstein-Hilbert action for GR, simplifying some computations. However, this frame choice also results in non-minimal coupling of matter to the metric, ensuring the theory behaves very differently to GR.

The Einstein frame is obtained by performing a conformal transformation of the Jordan frame. The Jordan frame prioritises use of a metric, $\hat{g}_{\mu\nu}$, which couples minimally to the matter fields but contains the non-trivial $G_{4}$ function. In this frame the gravitational Lagrangian is modified from the the Einstein-Hilbert action. The Jordan-frame metric is related to the Einstein-frame metric, $g_{\mu\nu}$, via a conformal factor $A$ that is a function of the Horndeski scalar:

$$\hat{g}_{\mu\nu}=A^{2}(\phi)g_{\mu\nu}\,.$$ (6)

In what follows, specifically in the case of K-mouflage theories, we will see that some quantities differ between the Jordan and Einstein frame. Though these quantities may be ‘physical’ in nature, they are not directly observable. General coordinate invariance – a key property shared with GR by nearly all modified gravity theories – ensures that observable quantities must be independent of frame choices (see for example Catena et al., 2007; Chiba & Yamaguchi, 2013; Francfort et al., 2019).

We summarise the models considered in this paper, their associated additional parameters and some selected constraints in Table 1.

Table 1: Overview of gravity models considered in this work. Note the K-mouflage kinetic term in Equation 14 does not pass Solar System tests without running into fine-tuning issues (Barreira et al., 2015b). Note the CG has no free parameters with the tracker solution. We have included constraints for the more general GCCG (see subsection 2.2).

model	screening method	free parameters	selected data constraints
nDGP	Vainshtein	$\{\Omega_{\rm rc}\}$	$\Omega_{\rm rc}\leq 0.235\,(2\sigma)$ (LSS) (Barreira et al., 2016)
CG	Vainshtein	$\{s=2,q=0.5\}$	$s=0.05^{+0.08}_{-0.05}$, $q>0.8$ $(2\sigma)$ (Various LSS, GCCG) (Frusciante et al., 2020)
K-mouflage	K-mouflage	$\{n,\beta_{\rm K},K_{0},\lambda\}$	$\beta_{\rm K}\leq 0.1$ (Lunar laser ranging) (Barreira et al., 2015b)

2.1 nDGP

The first model we consider is the Dvali-Gabadadze-Porrati model (Dvali et al., 2000), which does not strictly fall into the Horndeski class, being a five-dimensional braneworld model. It is given by the following action

$$S=\frac{1}{16\pi G_{5}}\int_{\cal M}{\rm d}^{5}x\sqrt{-\gamma}R_{5}+\int_{\partial{\cal M}}{\rm d}^{4}x\sqrt{-g}\left[\frac{M_{\rm pl}^{2}}{2}R+{\cal L}_{\rm m}\right]\,,$$ (7)

where $\gamma$ is the five-dimensional metric and $R_{5}$ its Ricci curvature scalar. $G_{5}$ is the five-dimensional gravitational constant. The matter Lagrangian is restricted to a four-dimensional brane in a five-dimensional Minkowski spacetime. The induced gravity given by the four-dimensional Einstein-Hilbert action is responsible for the recovery of four-dimensional gravity on the brane. The parameter $r_{\rm c}=G_{5}/(2G_{\rm N})$ is called the cross-over scale and is the only free parameter of the model, with its GR limit being $r_{\rm c}\to\infty$.

DGP also exhibits screening coming from higher order derivative terms in the effective 4-dimensional action. Such screening is known as Vainshtein screening (Vainshtein, 1972; Babichev & Deffayet, 2013). The so-called decoupling limit of DGP has the effective action given by (Luty et al., 2003; Gabadadze & Iglesias, 2006)

$$\mathcal{L}_{\rm DGP}=\frac{M_{\rm pl}^{2}}{2}R+\left(3\phi-\frac{1}{\Lambda^{3}_{\rm DGP}}(\partial\phi)^{2}\right)\Box\phi+\frac{1}{2M_{\rm pl}^{2}}\phi T\,,$$ (8)

where $T$ is the trace of the energy momentum tensor and $\Lambda^{3}_{\rm DGP}=M_{\rm pl}^{2}/r_{c}^{2}$. Note that although Equation 7 is not a Horndeski Lagrangian, Equation 8 is (compare to Equation 1). In this case we have

$$G_{3}(\phi,X)=3\phi+\frac{2}{\Lambda_{\rm DGP}^{3}}X\,,$$ (9)

with $G_{2}(\phi,X)=0$ and $G_{4}(\phi)=M_{\rm pl}^{2}/2$.

The literature typically assumes a $\Lambda$CDM background expansion, which is accommodated by introducing an appropriate dark energy contribution (see for example Schmidt, 2009b; Bag et al., 2018) on the stable ‘normal’ branch solution of the Friedmann equations. We follow this here (see Lue, 2006, for more details) and refer to this normal branch as nDGP. We also parametrize the modification to gravity using the energy density fraction $\Omega_{\rm rc}\equiv 1/(4r_{\rm c}^{2}H_{0}^{2})$, where $H_{0}$ is the Hubble constant. The GR-limit is then $\Omega_{\rm rc}\rightarrow 0$.

Although nDGP is now quite strongly constrained by observations (see for example Lombriser et al., 2009; Barreira et al., 2016; Piga et al., 2023), its appeal as a modified gravity model stems from the simplicity of its 4D effective action relative to the new phenomenology it introduces. It is one of the simplest examples of a gravity model that produces Vainshtein screening effects, whilst maintaining scale-independent growth of matter perturbations, and having only one additional parameter relative to $\Lambda$CDM. This has made it a favourite testbed for simulations (Schmidt, 2009a; Khoury & Wyman, 2009; Li et al., 2013a; Winther et al., 2017) and analyses with galaxy surveys (Barreira et al., 2016; Piga et al., 2023; Frusciante et al., 2023). We refer the reader to Refs. Koyama & Maartens (2006); Li et al. (2013a) and Appendix B for details on the modification to the Poisson equation (Equation 4) in linear and nonlinear regimes.

2.2 Cubic Galileon

The Cubic Galileon (CG) model was first derived by Nicolis et al. (2009) as a generalisation of the effective DGP action in 4-dimensions. The Lagrangian is given by (see for example Deffayet et al., 2009; Kobayashi et al., 2010)

$$\mathcal{L}_{\rm CG}=R\frac{M_{\rm pl}^{2}}{2}+c_{2}X+\frac{1}{\Lambda_{3}^{3}}c_{3}X\Box\phi\,,$$ (10)

where $c_{2}$ and $c_{3}$ are dimensionless constants parametrizing the modification to gravity, and the canonical choice for $\Lambda_{3}$ being $\Lambda_{3}^{3}=M_{\rm pl}H_{0}^{2}$, made to give the scalar field non-trivial dynamics on cosmological scales. Comparing with Equation 1 we have

$$G_{2}(\phi,X)=c_{2}X\,,\qquad G_{3}(\phi,X)=-\frac{1}{M_{\rm pl}H_{0}^{2}}c_{3}X\,.$$ (11)

This model also exhibits the Vainshtein mechanism due to the presence of the higher-order derivative terms (see Barreira et al., 2013a, for a derivation in the case of spherical symmetry). In this model, $G_{4}(\phi)=A(\phi)^{-2}/2=1$, and hence there is no difference between Jordan and Einstein frames (see Equation 6). We note that the absence of $G_{4}$ and conformal coupling allows one to interpret this model as a dark energy model with a non-trivial kinetic term.

The Cubic Galileon model is one member of a broader family, the Galileons, which added further derivative terms to Equation 10 (Deffayet et al., 2009). The Galileon family received intense interest from the theoretical physics community due to their shift symmetry properties (the actions are invariant under a shift $\phi(x)\rightarrow\phi(x)+c+b_{\mu}x^{\mu}$, $c$ and $b_{\mu}$ constants); this leads to special properties of the S-matrix. Cosmologically, their impact has been studied on the CMB (for example Barreira et al., 2014; Peirone et al., 2019; Frusciante et al., 2020; Albuquerque et al., 2022), linear matter power spectrum (for example Barreira et al., 2012) and gravitational lensing by voids (for example Baker et al., 2018). See also Refs. Renk et al. (2017); Peirone et al. (2018); Frusciante & Pace (2020) for other observational implications.

The more complex Galileon siblings have been virtually eliminated by their inability to have luminal gravitational waves, leaving behind only the Cubic Galileon (see for example Ezquiaga & Zumalacárregui, 2017; Baker et al., 2017). The Cubic Galileon model can be constrained by considering the integrated Sachs-Wolfe effect cross-correlated with a galaxy sample, as was done in Refs. Renk et al. (2017); Kable et al. (2022). The resulting cross-correlation, however, is shown to be anti-correlated with the expected $\Lambda$CDM signal, which severely constrains this model. It is worth noting, nevertheless, that a broader class of cubic Horndeski theories does not show this anti-correlation (Brando et al., 2019). Similarly to nDGP, it remains a useful testbed displaying Vainshtein screening, with a larger degree of flexbility due to its additional parameters and energy scales. It is a useful starting point from which to investigate the more general model space of the Horndeski Lagrangian (Equation 1).

We also note that the non-zero $G_{2}$ term makes this model phenomenologically distinct from nDGP. Further, in this paper we do not assume a $\Lambda$CDM background as with nDGP, but rather the solution to the Friedmann equations which include the effects of the scalar field (see for example Barreira et al., 2013a). A cosmology with this background evolution but with no further gravitational modification (so the Poisson equation remains as in GR), will be referred to as QCDM as in Ref. Barreira et al. (2013a).

The more general Generalised Covariant Cubic Galileon (GCCG) was recently considered in Ref. Frusciante et al. (2020), which promotes the $G_{i}$ functions to be power law functions of $X$, i.e. $G_{i}\propto X^{p_{i}}$. This model permits a tracker solution at the background level which is given by (De Felice & Tsujikawa, 2012)

$$H^{2q+1}\psi^{2q}=\zeta H_{0}^{2q+1}\,,$$ (12)

where $q\equiv(p_{3}-p_{2})+1/2$ and $\psi=\phi^{\prime}/M_{\rm pl}$. We also have the parameter $s=p_{2}/q$, leaving only 2 additional degrees of freedom for this model over $\Lambda$CDM. The GCCG reverts to the Cubic Galileon model when $q=0.5$ and $s=2$.

The GCCG model has not been ruled out by data, with CMB experiments giving the $2\sigma$ bounds of $q>0$ and $s=0.6^{+1.7}_{-0.6}$, with a slight preference for the model over $\Lambda$CDM (Frusciante et al., 2020). When combined with SN1a and redshift space distortion data sets, the bounds improve to $q>0.8$ and $s=0.05^{+0.08}_{-0.05}$. We note that theoretical stability conditions require both parameters to be positive.

In this paper we will only consider the CG limit of GCCG. We note that we employ the GCCG patch to the ReACT code (Atayde et al., 2024) for those specific predictions. For details on how the Poisson equation is modified in the CG limit, we refer the reader to Refs. Barreira et al. (2013a); Atayde et al. (2024).

2.3 K-mouflage

2.3.1 Lagrangian

The last model we consider is the K-mouflage model (Babichev et al., 2009). This model has the Lagrangian (in the Einstein frame)

$$\mathcal{L}_{\rm K}=R\frac{M_{\rm pl}^{2}}{2}+\mathcal{M}^{4}K(X)\,,$$ (13)

where $K(X)$ is a function of the canonical kinetic term, equivalent to a restricted $G_{2}(\phi,X)$, and $\mathcal{M}^{4}$ is an energy scale of the theory. We will set $\mathcal{M}^{4}=\lambda^{2}H_{0}^{2}M_{\rm pl}^{2}$ as in Ref. Hernández-Aguayo et al. (2022), $\lambda$ being an order $1$ dimensionless constant which can be tuned to give the current accelerated expansion of the Universe today. In this work we will consider a form which has been well studied in the literature (Brax & Valageas, 2014a, b; Barreira et al., 2015a, b; Hernández-Aguayo et al., 2022)

$$K(X)=-1+\frac{1}{H_{0}^{2}\lambda^{2}M_{\rm pl}^{2}}X+K_{0}\frac{1}{H_{0}^{2n}\lambda^{2n}M_{\rm pl}^{2n}}X^{n}\,,$$ (14)

where $K_{0}$ is another dimensionless model parameter and $n\geq 2$ is an integer. For the conformal function, we assume an exponential form

$$A(\phi)=\exp\left(\frac{\beta_{\rm K}\phi}{M_{\rm pl}}\right)\,,$$ (15)

where $\beta_{\rm K}$ is another dimensionless model parameter. In total we then have 4 parameters for this particular model: $\{\lambda,K_{0},n,\beta_{\rm K}\}$.

Unlike the other two models considered, the Jordan and Einstein frames are not set to be identical ($A(\phi)\neq 1$) which distinguishes this model from $k$-essence theories (Armendariz-Picon et al., 1999) where a universal coupling to matter is not present. In this work we will develop predictions for both frames. We provide the transformations of key quantities in the next subsection.

This model exhibits a similar screening mechanism to Vainshtein screening, although quantitatively different due to the absence of the higher-order $G_{3}(\phi,X)\Box\phi$ term, giving it a unique phenomenology. In particular, in dense environments of mass $m$, the K-mouflage radius – the scale below which GR is recovered – goes as $m^{1/2}$, whereas in Vainshtein theories this screening occurs at smaller scales, with a dependence of the Vainshtein radius on the environmental mass being $m^{1/3}$ (Brax & Valageas, 2014a). Vainshtein is also capable of screening large cosmological structures, while K-mouflage is not (Brax et al., 2015).

K-mouflage has been confronted with a number of cosmological data sets in Refs. Barreira et al. (2015b); Benevento et al. (2019), with a review of current constraints given in Ref. Brax et al. (2021) and forecasts using spectroscopic and photometric primary probes by Euclid given in Ref. Frusciante et al. (2023). In particular, in Ref. Barreira et al. (2015b), the authors place a Solar System constraint on the coupling parameter $\beta_{\rm K}\leq 0.1$, and argue that the power law form for $K(X)$ as chosen here will necessarily require a degree of fine tuning to avoid constraints. Despite this, this model is a good test case for implementation as it has been well studied in the literature and there are available $N$-body simulations with which to compare to (Hernández-Aguayo et al., 2022). More viable non-canonical kinetic terms can easily be implemented following the current implementations.

We alert the reader that we have made public a Mathematica notebook with some key Einstein frame quantities and derivations for the model along with this work. This contains useful expressions such as the exact solutions for the Einstein frame background $H(a)$ in the $n=2$ and $n=3$ cases.

2.3.2 Transformation to Jordan Frame

In this section we provide some basic translations between Einstein and Jordan frames which will be useful for our comparisons of the K-mouflage model. We follow Ref. Francfort et al. (2019) for these expressions. We use subscripts ‘J’ and ‘E’ to denote Jordan and Einstein frame quantities respectively.

The scale factor transforms as

$$a_{\rm J}=\bar{A}\,a_{\rm E}\,,$$ (16)

where $\bar{A}$ is the conformal factor evaluated at the background level (see Equation 15). The Hubble rate transforms as

$$H_{\rm J}(a)=\frac{H_{\rm E}}{\bar{A}}\left[1+\frac{\beta_{\rm K}}{M_{\rm pl}}\frac{{\rm d}\phi}{{\rm d}\ln a_{\rm E}}\right]\,.$$ (17)

The matter power spectrum transforms as (Francfort et al., 2019)

$$(2\pi)^{3}\delta_{\rm D}(\mbox{\boldmath$k$}_{1}+\mbox{\boldmath$k$}_{2})P_{\rm J}(k_{1}) =\langle\delta_{\rm J}(\mbox{\boldmath$k$}_{1})\delta_{\rm J}(\mbox{\boldmath$k$}_{2})\rangle$$

$$=\langle\delta_{\rm E}(\mbox{\boldmath$k$}_{1})\delta_{\rm E}(\mbox{\boldmath$k$}_{2})\rangle$$

$$-4\frac{\bar{A}_{\phi}}{\bar{A}}\langle\delta_{\rm E}(\mbox{\boldmath$k$}_{1})\delta\phi(\mbox{\boldmath$k$}_{2})\rangle$$

$$-4\frac{\bar{A}_{\phi}}{\bar{A}}\langle\delta\phi(\mbox{\boldmath$k$}_{1})\delta_{\rm E}(\mbox{\boldmath$k$}_{2})\rangle$$

$$+16\left(\frac{\bar{A}_{\phi}}{\bar{A}}\right)^{2}\langle\delta\phi(\mbox{\boldmath$k$}_{1})\delta\phi(\mbox{\boldmath$k$}_{2})\rangle\,,$$ (18)

where and we used $\delta_{\rm J}=\delta_{\rm E}-4\delta\phi\bar{A}_{\phi}/\bar{A}$, with $\bar{A}_{\phi}={\rm d}\bar{A}(\phi)/{\rm d}\phi$, $\delta$ is shorthand for the matter density field perturbation $\delta_{\rm m}$, $\delta\phi$ is the scalar field perturbation, $\phi=\bar{\phi}+\delta\phi$, and $k$ is the comoving Fourier mode in $h/{\rm Mpc}$. Angular brackets denote an ensemble average. The linear order Klein-Gordan equation for the scalar field perturbation in Fourier space under the quasi-static approximation is (Brax & Valageas, 2014b)

$$\delta\phi(\mbox{\boldmath$k$})=-\frac{\bar{A}\beta_{\rm K}a^{2}}{M_{\rm pl}K_{X}k^{2}}\bar{\rho}_{\rm m}\delta_{\rm E}(\mbox{\boldmath$k$})\,,$$ (19)

where $K_{X}={\rm d}\,K(X)/{\rm d}X$. Substituting $\delta\phi$ into Equation 18 gives us the following relationship between linear matter power spectra predictions

$$P_{\rm L,J}(k)=P_{\rm L,E}(k)\left[1+2\mathcal{J}(k,a)+\mathcal{J}(k,a)^{2}\right]\,,$$ (20)

where

$$\mathcal{J}(k,a)=\frac{12\bar{A}_{\phi}\beta_{\rm K}H_{0}^{2}M_{\rm pl}\Omega_{\rm m,0}}{ak^{2}K_{X}}\,,$$ (21)

where we have used the relation $\bar{\rho}_{\rm m}=3H_{0}^{2}M_{\rm pl}^{2}\Omega_{\rm m,0}a^{-3}$. We see that the linear matter power spectra in both frames are identical up to corrections that are suppressed by powers of $\sim H_{0}^{2}/k^{2}$.

It was argued in Ref. Francfort et al. (2019) that this correction to the matter power spectrum at nonlinear scales continues to go as $\sim H^{4}/k^{4}$, and so becomes negligible on all sub-horizon scales. This argument hinged on a number of assumptions, including $A_{\phi}\sim-A(\phi)/(2\phi)$. We will show in section 4 that the corrections are indeed small at nonlinear scales for the K-mouflage model, using a conformal factor given in Equation 15.

3 Tools & Methods

In this section we give an overview of the methods developed to give predictions for the large-scale structure in all modified gravity scenarios considered. After explaining details of how we compute matter power spectra, we describe the methods we will compare in this work. Most of these are $N$-body simulation-based approaches with various degrees of approximation. The halo model reaction (Cataneo et al., 2019) is also considered, which is an analytic method based on the halo model and perturbation theory. Table 2 gives an overview of these methods.

Table 2: Overview of the numerical codes employed in this comparison. For more information on screening approximations see main text. PPF: parametrised post-Friedmannian; PM: particle mesh; AMR: adaptive mesh refinement; 2LPT: 2nd order Lagrangian perturbation theory; K-G: Klein-Gordon.

code	type	screening approximation	reference(s)
ECOSMOG	$N$-body (AMR)	full K-G solution	Li et al. (2012)
MG-GLAM	$N$-body (uniform PM)	full K-G solution	Hernández-Aguayo et al. (2022)
MG-evolution	$N$-body (uniform PM)	PPF with free parameter $k_{*}$	Hassani & Lombriser (2020); Adamek et al. (2016b, a)
Hi-COLA	$N$-body (PM in the 2LPT frame)	screening factor	Wright et al. (2023), Sen Gupta et al. (in prep.)
COLA-FML	$N$-body (PM in the 2LPT frame)	linear K-G equation in Fourier space	Winther et al. (2017); Brando et al. (2023); Scoccimarro (2009)
ReACT	Halo model and perturbation theory	spherical collapse	Bose et al. (2022); Atayde et al. (2024)

3.1 $P(k)$ estimation

$N$-body simulations track the time-evolution of the matter distribution in the simulation box (of side $L_{\rm box}$) by means of a number of $N$-body particles ($N_{\rm P}$). To estimate the matter power spectrum from these sort of discrete distributions it is necessary to deal with some subtleties. The number of particles used in $N$-body simulations is often large (i.e. $10^{8}-10^{12}$) so that it would be computationally impractical to estimate the matter power spectrum by computing the distances between each pair of particles. Hence, the particles are normally interpolated on a regular grid using mass assignment schemes (MAS). Then the matter power spectrum is estimated exploiting the Fast Fourier Transform (FFT) algorithm. However, the modes close to the Nyquist frequency of the FFT grid can be significantly affected by aliasing (Jing, 2005; Sefusatti et al., 2016). To avoid this problem we use the interlacing technique with the triangular-shaped-cloud MAS (Sefusatti et al., 2016) to compute the matter power spectra from the simulations. Aiming to compare our matter power spectra deep in the nonlinear regime but mindful of the limited mass-resolution of our simulations, we use a FFT grid of size $N_{\rm mesh,1D}=L_{\rm box}/({\rm d}x)$ where ${\rm d}x$ is the domain grid resolution of the simulation, and use a simple linear binning with $k_{\rm min}=k_{\rm f}/2$ and $\Delta k=k_{\rm f}$, where $k_{\rm f}\equiv\frac{2\pi}{L_{\rm box}}$ is the fundamental frequency of the box.

3.2 Full $N$-body

Our reference predictions will come from numerical simulations that solve the nonlinear Klein-Gordon equation with multi-grid relaxation to get the precise modified force law. They also employ a large number of time steps over which the particles are evolved, ensuring the accuracy of the resulting predictions. We consider two variants of these codes.

3.2.1 ECOSMOG

The ECOSMOG simulation code (Li et al., 2012, 2013a) is a modified gravity extension of the adaptive mesh refinement (AMR) code Ramses (Teyssier, 2002). ECOSMOG relies on multigrid relaxation techniques to solve the nonlinear Klein-Gordon equations for the additional scalar fields that appear in some modified gravity theories (such as those considered here). This code was used to simulate several gravity models in the literature:

• f(R) (Li et al., 2012);

• nDGP (Li et al., 2013a);

• symmetron (Davis et al., 2012; Brax et al., 2013);

• dilaton (Brax et al., 2011);

• galileon (cubic, quartic, cubic vector) (Barreira et al., 2013b; Barreira et al., 2013a; Becker et al., 2020).

The accuracy of this code for predictions of $f(R)$ (Hu & Sawicki, 2007) effects on the matter power spectrum has been estimated to be of $\sim 1\%$ up to $k\sim 7h/{\rm Mpc}$ in the code comparison paper Ref. Winther et al. (2015). This code constitutes the highest precision predictive tool to be considered in this work.

3.2.2 MG-GLAM

MG-GLAM (Hernández-Aguayo et al., 2022; Ruan et al., 2022) extends the Particle Mesh (PM) code GLAM (Klypin & Prada, 2018) to a general class of modified gravity theories (including the K-mouflage model) by adding extra modules for solving the Klein-Gordon equations, using the multigrid relaxation algorithm. It uses a regularly spaced 3D mesh covering the cubic simulation box, solves the Poisson equation for the Newtonian potential using the Fast Fourier Transform (FFT) algorithm, and adopts the Cloud-In-Cell (CIC) scheme to implement the matter density assignment and force interpolation.

MG-GLAM has been tested with the results from other high-precision modified gravity codes, such as ECOSMOG (Li et al., 2012, 2013a), MG-GADGET (Puchwein et al., 2013), and MG-AREPO (Arnold et al., 2019; Hernández-Aguayo et al., 2021). For example, using $1024^{3}$ particles in a box of size $512\,\mathrm{Mpc}/h$, MG-GLAM simulations can accurately predict the matter power boost, $P_{\mathrm{MG}}/P_{\mathrm{\Lambda CDM}}$ at $k\lesssim 3\,h/\mathrm{Mpc}$, with about $1\%$ of the computational costs of the high-fidelity code ECOSMOG. Being the only code that has been used in the literature to simulate K-mouflage cosmologies, an estimate of its accuracy for the K-mouflage boost factor is not available. However it has been compared to the tree-PM code MG-Arepo for another derivative coupling model (nDGP) where it showed an agreement of $\sim 2\%$ up to $k=3\,h/\mathrm{Mpc}$, with deviations of $\sim 1\%$ from MG-Arepo (and theory predictions) already present on linear scales.

3.3 MG-evolution

We further consider the relativistic $N$-body code, MG-evolution (github Hassani & Lombriser, 2020). This code is based on gevolution (Adamek et al., 2016b), integrating parametrised modifications of gravity for various dark energy scenarios. MG-evolution is implemented based on a parametrisation framework that includes both linear and deeply nonlinear scales, where the nonlinear parametrisation is based on modified spherical collapse computations and a parametrised post-Friedmannian expression.

MG-evolution has been tested for a number of well-studied modified gravity models encompassing $f(R)$ and nDGP gravity that include large-field value and derivative screening effects (Hassani & Lombriser, 2020). Unlike most modified gravity $N$-body implementations, MG-evolution is as fast as the $\Lambda$CDM simulations as it does not need to deal with solving computationally expensive scalar field equations.

In Appendix B we discuss the nDGP and CG implementations in MG-evolution through a parametrisation with one screening transition, $k_{*}$, which is treated as a free parameter (see Appendix B). The effective gravitational constant is expressed as

$$\frac{\Delta{G}_{\text{eff }}(k,a)}{G_{\rm N}}=\frac{\Delta G_{\text{eff, L}}(a)}{G_{\rm N}}\times\frac{\Delta G_{\text{eff, NL}}(k,a)}{G_{\rm N}}\,,$$ (22)

where we recall $G_{\rm eff}(k,a)=G_{\rm N}[1+\Delta G_{\rm eff}(k,a)/G_{\rm N}]$, $G_{\text{eff, L}}$ denoting the linear regime parametrisation and $G_{\text{eff, NL}}$ refers to the parametrisation of the nonlinear regime that includes the screening or other suppression effects. The expressions for $G_{\text{eff, NL}}$ are given in Appendix B. MG-evolution then solves the modified Poisson equation (Equation 4) based on $G_{\rm eff}$ obtained from Equation 22. It is worth noting that this parametrization of gravitational modification is done in Fourier space. As detailed in Ref. Hassani & Lombriser (2020), this transformation yields an effective screening wavenumber $k_{*}$, which can be modeled (Lombriser, 2016) for different screening types. Currently, as mentioned, we treat $k_{*}$ as a free parameter to be set by the user. In this work we tune the values of $k_{*}$ in order to optimize the agreement with the reference predictions in each model and at each redshift considered. The resulting values of $k_{*}$ are presented in section 4.

3.4 COmoving Lagrangian Acceleration

The COmoving Lagrangian Acceleration (COLA) method (Tassev et al., 2013) is a hybrid $N$-body approach to performing dark matter simulations to study the effects of gravity on the formation of large-scale structure. It leverages the fact that the growth of structure on large scales can be computed analytically through Lagrangian Perturbation Theory (LPT). This informs the small-scale $N$-body part of COLA codes, thereby allowing for a significant speedup in the production of results at the cost of a modest loss of accuracy at small scales. In short, the COLA approach is a method well-suited for producing large-scale structure results on mildly nonlinear scales much faster than traditional $N$-body codes.

The COLA method introduced by Tassev et al. (2013) was originally developed for use with GR $N$-body codes. However, being able to probe nonlinear scales is particularly useful in the study of modified gravity theories, as key phenomenology, such as the effects of screening mechanisms, become apparent on these scales. Since Tassev et al., implementations of COLA codes for modified gravity theories have followed for specific theories, such as $f($R) and nDGP (Winther et al., 2017; Valogiannis & Bean, 2017). Below we describe two branches of work that extend the COLA method to more general families of gravity models.

3.4.1 Hi-COLA

Horndeski-in-COLA (Hi-COLA) (github Wright et al., 2023) is an implementation of the COLA methodology for a broad Horndeski class of scalar-tensor theories. Hi-COLA aims not to carry hard-coded theory-specific implementations, but instead receives as input the Lagrangian functions for a given theory of interest, making it generic. The action of the new scalar degree of freedom, $\phi$, is included as a fifth force in the COLA simulation. The class of theories Hi-COLA is designed to work with are reduced-Horndeski theories, where gravitational waves travel at the speed of light. Such theories are described by Equation 1. This restriction follows the constraints on the theory space suggested by the analysis of GW170817.

After receiving inputs for the forms of the Horndeski functions, $G_{2}$, $G_{3}$ and $G_{4}$, the symbolic manipulation modules of Hi-COLA construct the appropriate background equations of motion and background-dependent fifth force expressions and solves them. These are used to handle the expansion of the simulation box, compute 2nd-order LPT factors (2LPT) and construct the total force experienced by dark matter particles. This force can be schematically written as

$$F_{\textrm{total}}=G_{\rm eff}F_{\rm N}\,,$$ (23)

where

$$G_{\textrm{eff}}=\frac{G_{\rm G4}}{G_{\rm N}}\big{\{}1+\beta_{\rm HC}(z)S_{\textrm{HC}}(z,\delta_{m})\big{\}}.$$ (24)

$F_{\rm N}$ is the regular Newtonian force which is present in GR, and the multiplicative factor in braces represents the extra force contributions from $\phi$. $G_{G4}$ is the effective gravitational constant, which can differ from $G_{\rm N}$ in a time-dependent manner if $G_{4}$ in Equation 1 is non-trivial. This term will play a role in the results of subsubsection 4.2.3.

$\beta_{\textrm{HC}}$ is a background-dependent function known as the coupling factor; it controls the total possible strength of the fifth force at a given point in time. $S_{\textrm{HC}}$ is a background and density-dependent function called the screening factor. On linear scales $S_{\textrm{HC}}\rightarrow 1$, whilst in screened regimes $S_{\textrm{HC}}\rightarrow 0$. Hence this factor is responsible for the suppression of the fifth force in on small scales, returning the theory’s behaviour to GR.

$S_{\textrm{HC}}$ is derived under a quasi-linear perturbative treatment, where the metric perturbations are considered to first order, whilst the scalar field derivative perturbations are kept up to third order, following Ref. Kimura et al. (2012). Combined with the assumptions that the quasi-static approximation holds and that the matter over-density is distributed spherically in space leads to the analytic form of $S_{\textrm{HC}}$ (see Equation 3.15 in Wright et al., 2023). These assumptions in the derivation of $S_{\textrm{HC}}$ lead to a caveat: that in its current public state, Hi-COLA is designed to work with theories that exhibit Vainshtein screening. However, recent development of Hi-COLA has focused on extending the formalism to other screening mechanisms like K-mouflage, and these results are presented in subsubsection 4.2.3. The full details of K-mouflage in Hi-COLA will be the subject of an upcoming publication, Sen Gupta et al. (in prep.).

3.4.2 COLA-FML

In this subsection we describe another approximate simulation method to modified gravity theories endowed with the Vainshtein mechanism, such as nDGP and the Galileon theory family. This method was initially proposed in Ref. Scoccimarro (2009), and later revisited in Ref. Brando et al. (2023). It consists of linearizing the Klein-Gordon equation in Fourier space, and implementing a resummation scheme to find a function, $G_{\rm eff}(k,a)$, defined in the same way as Equation 22, that approximately captures the nonlinear corrections introduced by the Vainshtein mechanism on small scales. Specifically, this function transitions between an unscreened regime at large scales, where $G_{\rm eff}(k,a)\to G_{\rm eff,L}(a)$, to the small scale regime where GR is recovered, $G_{\rm eff}(k,a)\to G_{\rm N}$.

Consequently, in order to do so, in Refs. Scoccimarro (2009); Brando et al. (2023) the authors require the nonlinear function, $\Delta G_{\text{eff, NL}}(k,a)/G_{\rm N}$, has the screening property, i.e.

$$\frac{\Delta G_{\text{eff, NL}}(k,a)}{G_{\rm N}}\left(k/k_{*}\ll 1\right) \to 1,$$

$$\frac{\Delta G_{\text{eff, NL}}(k,a)}{G_{\rm N}}\left(k/k_{*}\gg 1\right) \to 0\,,$$ (25)

where $k_{*}$ is the wavenumber associated with the Vainshtein radius, defined in Eq. 53. The specifics behind the computation of the function $\Delta G_{\text{eff, NL}}(k,a)/G_{\rm N}$ is explicitly shown in Ref. Brando et al. (2023). This screening approximation scheme has the advantage of not introducing additional screening parameters used to tune the approximate results with results from N-body simulations that consistently solve the full Klein-Gordon equation at each timestep of the simulation. The whole dependence of the gravity theory is encoded in the $\Delta G_{\text{eff, NL}}(k,a)/G_{\rm N}$ function.

The methodology of this approximate method for Vainshtein screening is computed using an external Python notebook, where one can follow the steps outlined in Ref. Brando et al. (2023) to compute $G_{\rm eff}(k,a)$ externally. With the tabulated function computed, the results are then implemented in the COLA-FML (github) N-body solver, that implements the COLA method in a parallelised manner, ideal for fast and approximate simulations. The COLA-FML library also has different screening approximations for theories other than the ones considered here, and are presented in Ref. Winther & Ferreira (2015a). Importantly for this paper, our results for the $G_{\rm eff}(k,a)$ screening case will be different than the ones of Hi-COLA at nonlinear scales, however, at linear scales the two codes are identical.

3.5 Halo Model Reaction

The halo model reaction (Cataneo et al., 2019) is a flexible, accurate and fast means to model the nonlinear matter power spectrum. This model has been demonstrated to align with $N$-body simulations at the 2% level down to $k=3~{}h/{\rm Mpc}$, with minor variations depending on redshift, the extent of modification to GR, and the mass of neutrinos (Cataneo et al., 2020; Bose et al., 2021). The method aims to model nonlinear corrections to the matter power spectrum resulting from modified gravity through the reaction $\mathcal{R}(k,z)$, which incorporates both 1-loop perturbation theory and the halo model (see Bernardeau et al., 2002; Cooray & Sheth, 2002, for reviews). In this framework, the nonlinear matter power spectrum is expressed as the product

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

where the pseudo power spectrum is defined such that all nonlinear physics are modeled using GR but the initial conditions are adjusted to mimic the modified linear clustering at the target redshift.

The halo model reaction without massive neutrinos, $\mathcal{R}(k,z)$, is given as a corrected ratio of target-to-pseudo halo model spectra

$$\mathcal{R}(k,z)=\frac{\{[1-\mathcal{E}(z)]\,e^{-k/k_{\star}(z)}+\mathcal{E}(z)\}\,P_{\rm 2H}(k,z)+P_{\rm 1H}(k,z)}{P_{\rm hm}^{\rm pseudo}(k,z)}\,.$$ (27)

The components are explicitly given as

$$P_{\rm hm}^{\rm pseudo}(k,z)= P_{\rm 2H}(k,z)+P_{\rm 1H}^{\rm pseudo}(k,z),$$ (28)

$$\mathcal{E}(z)= \lim_{k\rightarrow 0}\frac{P_{\rm 1H}(k,z)}{P_{\rm 1H}^{\rm pseudo}(k,z)},$$ (29)

$$k_{\rm\star}(z)= -\bar{k}\left\{\ln\left[\frac{A(\bar{k},z)}{P_{\rm 2H}(\bar{k},z)}-\mathcal{E}(z)\right]-\ln\left[1-\mathcal{E}(z)\right]\right\}^{-1}\,,$$ (30)

with

$$A(k,z)=\frac{P_{\rm 1-loop}(k,z)+P_{\rm 1H}(k,z)}{P^{\rm pseudo}_{\rm 1-loop}(k,z)+P_{\rm 1H}^{\rm pseudo}(k,z)}P_{\rm hm}^{\rm pseudo}(k,z)-P_{\rm 1H}(k,z)\,.$$ (31)

$P_{\rm 2H}(k,z)$ is the 2-halo term which we approximate with the linear matter power spectrum, $P_{\rm L}(k,z)$. $P_{\rm 1H}(k,z)$ and $P_{\rm 1H}^{\rm pseudo}(k,z)$ are the 1-halo terms as predicted by the halo model, with and without modifications to the standard spherical collapse equations, respectively. Recall that by definition, the pseudo cosmology has no nonlinear beyond-$\Lambda$CDM modifications. Similarly, $P_{\rm 1-loop}(k,z)$ and $P_{\rm 1-loop}^{\rm pseudo}(k,z)$ are the standard perturbation theory 1-loop matter power spectra with and without nonlinear modifications to $\Lambda$CDM, respectively. As in the literature, Equation 29’s limit is taken to be at $k=0.01\,h/{\rm Mpc}$ and $k_{\star}$ is computed using $\bar{k}=0.06\,h/{\rm Mpc}$.

The nDGP model was part of the initial release of the publicly available halo model reaction code, ReACT (Bose et al., 2020). This code has been updated to include massive neutrinos in Ref. Bose et al. (2021) and model independent parametrisations in Ref. Bose et al. (2022), which constituted version 2 of the code (github). The GCCG model was recently implemented in this version of ReACT (Atayde et al., 2024), which is employed in section 4 in the CG limit. The K-mouflage patch is being made public with this work and we give all the relevant expressions in Appendix A.

For the pseudo spectrum appearing in Equation 26 we use HMCode2020 (Mead et al., 2020). This is currently the most accurate and flexible prescription for the pseudo spectrum and has been tested in a number of works (see for example Cataneo et al., 2019; Bose et al., 2021; Bose et al., 2022). It is more accurate than the halofit prescription of Ref. Takahashi et al. (2012), quoting a $2.5-5\%$ accuracy for $k\leq 1~{}h/{\rm Mpc}$. It can also accommodate modifications that induce an additional scale dependence in the linear matter power spectrum. For modifications that only introduce a scale-independent shift in the linear spectrum amplitude, more accurate emulators can be used, such as the EuclidEmulator2 (Knabenhans et al., 2019), which are quoted to be $1\%$ accurate when compared to high fidelity $N$-body simulations down to $k\leq 10~{}h/{\rm Mpc}$. Despite this, the reaction function $R(k,z)$ is only expected to be $1\%$ accurate for $k\leq 1~{}h/{\rm Mpc}$ (Cataneo et al., 2019).

It is also worth noting that EuclidEmulator2’s internal accuracy is restricted to a hyperspheroidal region of their parameter space. Points outside this region might have considerable degradation in accuracy. This is considerably important in the context of beyond-$\Lambda$CDM scenarios as we need tools that work in extreme regions of the parameter space. For these reasons, work is currently being undertaken to build a pseudo spectrum emulator based on appropriate numerical simulations (Giblin et al., 2019).

The choice of HMCode2020 keeps in line with the halo model reaction’s claim of generality, while maintaining competitive accuracy within the reaction’s percent-level accuracy range, especially when taking the ratio of modified to unmodified spectra, i.e. the matter power spectrum boost (see Equation 32).

4 Results

Our main results are the comparisons of the nonlinear matter power spectrum between the different codes. Specifically, we consider the models described in Table 3 for which we have full $N$-body simulations available to use as benchmarks. We list the specifications of each simulation ran for these comparisons below. These include: box size ($L_{\rm box}$), number of particles ($N_{\rm P}$), particle mass ($m_{\rm P}$), grid cells ($N_{\rm g}$), initial redshift ($z_{\rm ini}$) and force resolution.

• ECOSMOG runs: $L_{\rm box}=1024~{}\textrm{Mpc}/h$,   $N_{\rm P}=1024^{3}$,   $m_{\rm P}\simeq 7.8\times 10^{10}~{}M_{\odot}/h$. The initial conditions are generated at $z_{\rm ini}=49$ by MPGrafic (Prunet et al., 2008) using the Zel’dovich approximation. It uses a force resolution of $\sim 15.6\textrm{kpc}/h$.

• MG-GLAM runs: $L_{\rm box}=512\,\mathrm{Mpc}/h,\,N_{\rm P}=1024^{3},\,N_{\rm g}=2048^{3},\,m_{\rm P}=1.07\times 10^{10}\,M_{\odot},/h$, where $N_{\rm g}$ is the number of grid cells. Initial conditions are generated at $z_{\rm ini}=100$ using GLAM’s own initial condition generator. It uses a fixed force resolution of $250\,\mathrm{kpc}/h$ with an adaptive time-stepping described in the original GLAM paper (Klypin & Prada, 2018).

• MG-evolution runs: The nDGP simulation runs use $L_{\rm box}=1000~{}{\rm Mpc}/h$ with $N_{\rm g}=N_{\rm p}=1024^{3}$. The initial conditions are generated at $z_{\rm ini}=49$. For the CG case, the initial conditions are the same as in the COLA runs. These runs use $L_{\rm box}=400~{}{\rm Mpc}/h$, $N_{\rm P}=N_{\rm g}=512^{3}$.

• COLA runs: $L_{\rm box}=400~{}\textrm{Mpc}/h$ and $N_{\rm P}=512^{3}$ with initial conditions generated using 2LPT for all simulations. For K-mouflage: $m_{\rm P}\simeq 4.1\times 10^{10}~{}{\rm M_{\odot}}/h$. Initial conditions are generated at $z_{\rm ini}=19$. For nDGP: $m_{\rm P}\simeq 3.7\times 10^{10}~{}M_{\odot}/h$. Initial conditions are generated at $z_{\rm ini}=49$. For CG and QCDM: $m_{\rm P}\simeq 4.1\times 10^{10}~{}{\rm M_{\odot}}/h$. Initial conditions are generated at $z_{\rm ini}=49$.

Before presenting the spectra comparisons, we take a look at how each model presented in section 2 modifies the standard $\Lambda$CDM background evolution. This background evolution is adopted for each of the different codes and so differences seen in the following section only arise from how the perturbations are treated.

Table 3: Models considered in this work. The $\Lambda$CDM $\sigma_{8}(z=0)=0.851,0.805,0.815$ for nDGP, CG and K-mouflage cosmologies respectively. We remind the reader that the values for $\{s,q\}$ are fixed in CG, while adopting the tracker solution for the scalar field imposes the values for $c_{2}$ and $c_{3}$ quoted in the table.

Parameter	nDGP-N1	nDGP-N5	CG	QCDM	K-mouflage - A	K-mouflage - B	K-mouflage - C
$\Omega_{\rm m,0}$	0.281		0.313		0.3089
$\Omega_{\rm b,0}$	0.046		0.049		0.0486
$H_{0}$	69.7		67.32		67.74
$n_{s}$	0.971		0.9655		0.9667
$A_{s}$	$2.297\times 10^{-9}$		$2.010\times 10^{-9}$		$2.064\times 10^{-9}$
$\sigma_{8}(z=0)$	0.912	0.865	0.884	0.865	0.881	0.852	0.837
$\Omega_{\rm rc}$	0.25	0.01	-	-	-	-	-
$c_{2}/c_{3}^{2/3}$	-	-	-5.378	-	-	-	-
$c_{3}$	-	-	10	-	-	-	-
$s$	-	-	2.0	-	-	-	-
$q$	-	-	0.5	-	-	-	-
$n$	-	-	-	-	2	2	2
$\lambda$	-	-	-	-	1.475	1.460	1.452
$K_{0}$	-	-	-	-	1	10	1
$\beta_{\rm K}$	-	-	-	-	0.2	0.2	0.1

Figure 1: The ratio of the normalised Hubble expansion rate ($E(a)=H(a)/H_{0}$) between the modified gravity and GR models. The left panel shows the K-mouflage models shown in Table 3 in the Einstein frame, while the middle panel shows the same models in the Jordan frame, with $a$ now being the Jordan frame scale factor. The right panel shows the QCDM model, which has the same background expansion as the CG model. The model parameters for K-mouflage are defined in subsection 2.3.

4.1 Background Evolution

In Figure 1 we show the modification to the standard $\Lambda$CDM background expansion for the models described in Table 3. We remind the reader that we assume a $\Lambda$CDM expansion for the nDGP models and so this is not shown. We see that the QCDM and CG cases give a much larger modification at late times than any of the K-mouflage models in the Einstein frame. In all models, we see a slower expansion rate at roughly $a>0.2$ which acts to enhance structure formation. Indeed, the $\sigma_{8}$ is larger for the QCDM model than for both K-mouflage models B and C (see Table 3), despite having a lower $A_{s}$ (although the QCDM cosmology has a slightly larger $\Omega_{\rm m,0}$). In all cases, the maximum modification is $\sim 8\%$ (QCDM), with the K-mouflage models giving a maximum modification of $3\%$ at $a=0.5$.

In the same figure we also show the modification in the Jordan frame for the K-mouflage models (middle panel). We see here that relative to $\Lambda$CDM, we have a significantly slower expansion at $a>0.03$, with a maxmimum modification of $11\%$ at $a\sim 0.4$. Further, the current day expansion is larger than the one expected from $\Lambda$CDM by $5\%$ under the strongest modification considered here. We do remind the reader that the free parameter $\lambda$ has been tuned to match the current day expansion rate in $\Lambda$CDM in the Einstein frame. These panels show that relatively large differences can be observed at the background level when switching frames, which we will see in the next subsection are not evident at the level of the perturbations (also see subsubsection 2.3.2).

4.2 Matter power spectrum boost

Next we take a look at the perturbations, specifically how the matter power spectrum is modified over $\Lambda$CDM. For this, we consider the modified gravity boost, defined as

$$B(k,z)\equiv\frac{P_{\rm NL}(k,z)}{P^{\Lambda\rm CDM}_{\rm NL}(k,z)}\,.$$ (32)

Figure 2: Comparing boost factors for the various codes listed in Table 2 for nDGP at $z=0,0.5,1$ (from left to right) with ECOSMOG as benchmark. The upper panels show the results for the nDGP-N5 (low modification) model and the lower panels for the nDGP-N1 (high modification) model (see Table 3).

4.2.1 nDGP

In Figure 2 we show how the various predictions for the boost compare, using the ECOSMOG measurements as a reference, for the nDGP cosmologies found in Table 3. Boost comparisons for nDGP amongst different codes have already been performed extensively in the literature, and so this case is shown mainly as a consistency check, but also to compare the Hi-COLA implementation which has not yet been tested before.

We find that for the low modification case, N5, all predictions remain within $1\%$ of each other for $k\leq 3~{}h/{\rm Mpc}$, including the linear prediction, which for $z=0$ gives a modification of $B(k\rightarrow 0,z=0)=1.033$. For N1, $B(k\rightarrow 0,z=0)=1.149\%$. In this case, all predictions except the linear remain within $2\%$ of the ECOSMOG reference for $k\leq 1~{}h/{\rm Mpc}$. MG-evolution performs the best as expected, having an additional free parameter giving the screening transition, $k_{*}$. We have found $k_{*}=2~{}h/{\rm Mpc}$ and $k_{*}=1~{}h/{\rm Mpc}$ give a good overall agreement with the ECOSMOG simulations for the N1 and N5 models respectively. Using these values the MG-evolution boost remains within $1\%$ up to $k\leq 3~{}h/{\rm Mpc}$ except for the largest modes at $z=1$ where it worsens to $2\%$, consistent with what was found in Ref. Hassani & Lombriser (2020). Similarly, the halo model reaction remains within $2\%$ for $k\leq 3~{}h/{\rm Mpc}$ except for the largest modes at $z=1$, where it worsens to $3\%$ agreement, in accordance with Ref. Cataneo et al. (2019).

The two COLA methods show similar agreement, but deviate the most on average from the reference boost measurements. COLA-FML performs slightly better at $z=0$ while Hi-COLA does better at higher $z$, with deviations up to $4\%$ at $k=3~{}h/{\rm Mpc}$. This is very consistent with the results of Ref. Winther et al. (2017).

Figure 3: Comparing boost factors for the various codes listed in Table 2 for the CG (upper) and QCDM (middle) and QCDM-based boost (bottom) at $z=0,0.5,1$ (from left to right) with ECOSMOG as the benchmark. Note COLA-FML and Hi-COLA’s results for QCDM are identical and so we only show the Hi-COLA QCDM ratio in the middle panels.

Figure 4: Comparing boost factors for the various codes listed in Table 2 for the K-mouflage models listed in Table 3 with $\{K_{0},\beta_{\rm K}\}=\{1,0.2\},\{10,0.2\},\{1,0.1\}$ (from top to bottom) at $z=0,0.5,1$ (from left to right) with MG-GLAM as the benchmark. All models assume $n=2$.

4.2.2 Cubic Galileon (CG)

Figure 3 shows the boost comparisons between the various codes for the CG and QCDM cases, again using ECOSMOG as a reference. These ECOSMOG simulations were ran using the same code as presented in Ref. Barreira et al. (2013a). We have changed the baseline cosmology for these new runs, particularly lowering the value of $A_{s}$ and $H_{0}$. We also run a $\Lambda$CDM counterpart with which to calculate Equation 32. Previous works have compared the boost ratio of the CG spectrum to that in QCDM (Wright et al., 2023), or have performed direct spectra comparisons (Atayde et al., 2024). Further, in Ref. Atayde et al. (2024) the authors found significant disagreement when using an HMCode2020 prescription, which was outperformed by the halofit pseudo spectrum prescription. This was being caused by the $\sigma_{8}$-dependent damping introduced into HMCode2020 (Mead et al., 2020), which was not calibrated for particularly high values of $\sigma_{8}$ as that of the simulations found in Ref. Barreira et al. (2013a). The lower value of $\sigma_{8}$ in our simulations was found to greatly improve the performance of HMCode2020 over halofit. For comparison with previous work, we also show the comparisons for the ratio of CG to QCDM power spectra, or QCDM-based boost, in the bottom panels of Figure 3.

The MG-evolution predictions again give the best agreement, remaining within $1\%$ in the CG case for $z\geq 0.5$ down to $k=3~{}h/{\rm Mpc}$. For $z=1$, the linear implementation, or equivalently $k_{*}\to\infty$, provides the best match. However, in the figure, we have plotted the case $k_{*}(z\geq 0.5)=6~{}h/{\rm Mpc}$ as it appears to work well given the resolution of the simulation. Adopting the value $k_{*}=6~{}h/{\rm Mpc}$ at $z=0$ causes a quick divergence of the predictions as expected from the nature of $k_{*}$, amounting to a $8\%$ disagreement at $k=1~{}h/{\rm Mpc}$. Adopting $k_{*}=0.4~{}h/{\rm Mpc}$ at $z=0$ brings the predictions to within $2\%$ agreement in the same range of scales. Interestingly, in the QCDM-based boost case we can adopt the same value of $k_{*}=0.3~{}h/{\rm Mpc}$ for all redshifts considered while keeping a good fit to the ECOSMOG measurements. This seems to indicate that $k_{*}$ is also degenerate to some extent with the background modification.

In the QCDM case, the predictions are consistent within $1\%$ down to $k=1~{}h/{\rm Mpc}$. The disagreement for $k>1~{}h/{\rm Mpc}$ arises from resolution effects, as evidenced by the agreement between MG-evolution and Hi-COLA. This suggests that the tuning of $k_{*}$ performed to match the reference boost factor in CG is compensating for the resolution-induced loss of boost due to the background by suppressing the small-scale screening.

The halo model reaction remains within $1\%$ for $k\leq 1~{}h/{\rm Mpc}$ for both QCDM and CG cases, with the exception of the CG at $z=0$. Here we find up to $4\%$ disagreement with ECOSMOG. This is an atypically large disagreement given the similarity of CG to nDGP, for which the halo model reaction performs significantly better. To investigate this, we have tested different pseudo spectra prescriptions, specifically halofit and EuclidEmulator2, neither offering significant improvement for the matter power spectrum boost. We have also tried omitting the 1-loop correction (see Equation 31) with little change to the predictions as found in Ref. Bose et al. (2022). The excellent agreement in the QCDM case at $z=0$, with $1\%$ agreement beyond $k=3~{}h/{\rm Mpc}$, indicates no issue in the background implementation. Further, the QCDM-based boost comparisons show the same disagreement at $z=0$, but the same or better agreement at higher redshifts, which is just a partial cancellation of inaccuracies in the QCDM and CG $\Lambda$CDM-based boost cases.

Lastly, we also checked the behaviour of the reaction function $\mathcal{R}$ for varying GCCG modification strengths by changing the value of $s$. We compared these to corresponding nDGP predictions for $\mathcal{R}$ such that the nDGP models gave the same linear enhancement of structure as the GCCG cases, making their pseudo spectra identical. We observed significantly more suppression coming from $\mathcal{R}$ in the GCCG than nDGP, especially for large modifications (large $s$ or large $\Omega_{\rm rc}$), which may be due to the $G_{2}$ term present in the GCCG. We do note that the CG has a very large linear enhancement of clustering at $z=0$, equivalent to a nDGP model with $\Omega_{\rm rc}=0.6$. This may indicate a break down of the halo model reaction’s assumptions, specifically the $\Lambda$CDM fits it assumes for the halo mass function and virial concentration. The latter has been shown to significantly impact its accuracy (Cataneo et al., 2019; Srinivasan et al., 2021; Srinivasan et al., 2024), especially when the modification to gravity is large. To further pin the $z=0$ CG disagreement down, we would need to run a CG pseudo cosmology simulation which would make it clear whether or not the reaction modelling or $\Lambda$CDM-fits in the halo model components are failing. GCCG simulations with a smaller modification will also indicate this validity range. This will be the focus of future work.

Finally, both COLA implementations remain $2\%$ consistent with ECOSMOG in the CG case at scales $k\leq 1~{}h/{\rm Mpc}$. COLA-FML performs slightly better at low $z$ while Hi-COLA shows better agreement at higher $z$. The implementations differ only in their approach to screening and so we only show the Hi-COLA results for QCDM, where it is similarly consistent as in the CG case. We note that all codes tend to under-predict the boost at small scales. Part of this difference surely comes from the fact that while the ECOSMOG code consistently solves the full Klein-Gordon equation, the other codes implement the screening mechanism in an approximate way, making use of the spherical approximation in one way or another. Therefore, at smaller scales these approximate methods are not guaranteed to be valid. A better test of the accuracy of screening is provided by the QCDM-based boost in the bottom panels, where we see far better agreement between the COLA methods and ECOSMOG.

On this note, we remark that both COLA and MG-evolution’s disagreement with the benchmarks in both nDGP, CG and QCDM cases is also partially due to a low force resolution which can lead to a loss of power on small scales (see Brando et al., 2022, for example). By increasing the force resolution, and time steps in the COLA cases, we expect to find much better agreement above $k=1~{}h/{\rm Mpc}$, particularly in the QCDM case which does not have screening. We note that the limited force accuracy will affect all particle mesh codes, including MG-GLAM, and the most efficient and sure way to go to smaller scales would be to use Tree-particle mesh or AMR codes like ECOSMOG.

Table 4: Maximal percent deviation of the nonlinear matter power spectrum boost under various modelling approaches against benchmarks, at different redshifts for $k\leq 1(3)~{}h/{\rm Mpc}$.

	MG-evolution			COLA-FML			Hi-COLA			ReACT
Model	z=0	z=0.5	z=1	z=0	z=0.5	z=1	z=0	z=0.5	z=1	z=0	z=0.5	z=1
N1	1(1)%	1(1)%	1(2)%	1(2)%	1(4)%	1(4)%	2(4)%	1(4)%	1(3)%	2(2)%	2(2) %	2(3)%
N5	1(1)%	1(1)%	1(1)%	1(1)%	1(1)%	1(1)%	1(1)%	1(1)%	1(1)%	1(1)%	1(1)%	1(1)%
CG	1(1)%	1(1)%	1(2)%	1(6)%	1(5)%	1(5)%	1(10)%	2(6)%	1(3)%	4(8)%	2(5)%	1(3)%
QCDM	1(5)%	1(3)%	1(2)%	-	-	-	2(6)%	1(3)%	1(3)%	1(1)%	1(4)%	1(5)%
K-mouflage A	-	-	-	-	-	-	1(4)%	1(3)%	1(3)%	2(17)%	1(8)%	1(4)%
K-mouflage B	-	-	-	-	-	-	1(2)%	1(1)%	1(1)%	2(10)%	1(5)%	1(2)%
K-mouflage C	-	-	-	-	-	-	1(1)%	1(1)%	1(1)%	1(4)%	1(3)%	1(1)%

4.2.3 K-mouflage

For K-mouflage, we restrict our comparisons to MG-GLAM, Hi-COLA and ReACT with the MG-evolution implementation to be the focus of an upcoming work. We expect the same level of accuracy as exhibited in the CG and nDGP cases, especially given the freedom imparted by $k_{*}$.

In Figure 4 we show the results for the K-mouflage model. As a reference we use the MG-GLAM simulations, ran for the purpose of this comparison. We compare the K-mouflage boost for the three models listed in Table 3, all of which assume $n=2$ in Equation 14. We begin by noting that the coupling of matter to the scalar field is proportional to $\beta_{\rm K}/K_{0}$ (see Equation. 81 of Brax & Valageas, 2014b, for example), and so large positive $K_{0}$ decreases the fifth force while large $\beta_{\rm K}$ increases it. We find the larger the modification, the worse the agreement between ReACT, Hi-COLA and MG-GLAM. We can see this by moving from top to bottom panels in Figure 4. Further, we note for the largest modification (top panels), there is a $~{}1\%$ offset between MG-GLAM and linear theory (as well as the other codes). This was also seen in Figure. 10 of Ref. Hernández-Aguayo et al. (2022) but not seen in the linearised simulations presented in that reference, suggesting this is a consequence of the nonlinear treatment of MG-GLAM. We also note much smaller linear theory offsets at large scales for the weaker modifications.

For the strongest modification, K-mouflage A in Table 3, at low $z$, all codes are consistent within $2\%$ for $k\leq 1~{}h/{\rm Mpc}$. This agreement improves for the halo model reaction to $1\%$ agreement for $k\leq 3~{}h/{\rm Mpc}$ for $z=1$ and the weakest modification, K-mouflage C in Table 3. Overall, Hi-COLA does not show significant improvement or degradation with redshift or modification strength, consistently remaining within $2\%$ for $k\leq 3~{}h/{\rm Mpc}$. The exception is K-mouflage A for $z=0$ (upper left panel), where it degrades to 4% discrepancy at $k=3~{}h/{\rm Mpc}$. The Hi-COLA predictions are all made in the Jordan frame while MG-GLAM and ReACT produce predictions in the Einstein frame. It is here we note the consistency of the nonlinear matter power spectrum in both frames, confirming the claim of Ref. Francfort et al. (2019).

Before concluding we make some technical notes on the comparisons. In the case of the Jordan frame predictions from Hi-COLA, the boost is taken with the K-mouflage spectrum measured at $a_{\rm J}$, calculated using Equation 16. Finally, we note that ReACT has the option to use the PPF screening formalism for K-mouflage as derived in Ref. Lombriser (2016), and which we present in Appendix B for completeness. This framework comes with an additional degree of freedom and so we have chosen not to use this in our comparisons. The inclusion of K-mouflage theories in Hi-COLA will be the subject of an upcoming publication, Sen Gupta et al. (in prep.).

5 Conclusions

High quality $N$-body codes for modified gravity are essential in order to place reliable constraints on gravity using large-scale structure (LSS) observations. Ongoing galaxy surveys such as Euclid or the Dark Energy Survey will heighten their necessity by beating down the statistical uncertainty on our measurements, making theoretical accuracy essential. Benchmarking the accuracy of approximate but computationally efficient numerical methods against these high quality simulations is an important step towards reliable constraints from the forthcoming data.

In this paper we have performed comparisons of the matter power spectrum modification induced by three distinct theories of modified gravity, each of which induces a scale-independent enhancement of the linear growth of structure: the normal branch of the DGP braneworld model, the Cubic Galileon and K-mouflage. The former two employ the Vainshtein screening mechanism, while the latter employs the K-mouflage screening mechanism. For similar comparisons with scale-dependent modifications to the linear growth and the chameleon (Khoury & Weltman, 2004) or symmmetron (Hinterbichler & Khoury, 2010) screening mechanisms, we refer the reader to Refs. Winther et al. (2015); Winther et al. (2017); Cataneo et al. (2019); Hassani & Lombriser (2020).

We compare the matter power spectrum boost predicted by 6 different numerical codes, each of which has a varying approach to the nonlinear gravitational coupling: full $N$-body (ECOSMOG and MG-GLAM), Comoving Lagrangian Acceleration (COLA) of which we compare two distinct codes, Hi-COLA and COLA-FML, the relativistic parametrised $N$-body code, MG-evolution, and the semi-analytic halo model reaction approach expressed by the ReACT code. We summarise the distinctions of each code below:

• Full $N$-body: Solves the Klein-Gordan equation exactly to get the force applied to particles in a box. Serves as an accuracy benchmark.

• Hi-COLA: Includes a fifth force in the COLA formalism via a screening factor, as well as consistently solving the modified cosmological expansion history. Screening factors are derived using a quasi-linear treatment of the metric and scalar field perturbations, along with assuming the quasi-static approximation and spherically distributed over-densities.

• COLA-FML: Introduces the Vainshtein mechanism by evaluating a function, $G_{\rm eff}(k,a)$, that captures on average nonlinear corrections from the screening mechanism. This method is performed by linearizing the Klein-Gordon equation in Fourier space, and the full function is found by an iterative process.

• MG-evolution: Employs a parametrised ansatz for the nonlinear force law which comes with a screening parameter.

• ReACT: Uses spherical collapse, the halo model and 1-loop perturbation theory to predict the matter power spectrum.

We summarize the overall accuracy exhibited by each approach in Table 4 with respect to the full $N$-body benchmark. We remark that $N$-body codes solving the full Klein-Gordon equation in modified gravity are $1\%$ consistent (Winther et al., 2015) for $k\lesssim 7~{}h/{\rm Mpc}$ in their prediction for the boost.

We find that all approaches considered here are overall $2\%$ consistent with the benchmark $N$-body boost at scales $k\leq 1~{}h/{\rm Mpc}$ and for $z\leq 1$. The only exceptions are ReACT for the strongest modifications to $\Lambda$CDM and at $z=0$. MG-evolution performs the best, with a general accuracy of $1\%$ at all scales considered ($k\leq 3~{}h/{\rm Mpc}$), but this accuracy comes at the cost of tuning the screening parameter depending on the output redshift or modification strength, which might undermine the predictivity of the code.

We thus can advocate the safe use of these codes, and any emulators based upon them (see Tsedrik et al., 2024; Carrion et al., 2024; Gordon et al., 2024, for example), at fairly nonlinear scales for scale-independent models. We note the caveat that emulation error should be quantified and appropriately accounted for.

For a more concrete estimate on the validity of these methods, we can consider a Euclid-like survey whose weak lensing analysis will have a signal to noise peaking at (conservatively) $z\approx 0.7$ (see Lepori et al., 2022, for example). Imposing a $2\%$ accuracy demand on the matter power spectrum model, and assuming a $\Lambda$CDM fiducial background cosmology, we can arguably trust all method predictions for $\ell_{\rm max}\lesssim 1800$. This roughly corresponds to the pessimistic scenario described in Ref. Blanchard et al. (2020).

At scales $k>1~{}h/{\rm Mpc}$ we find all codes begin to diverge by more than $2\%$ for the strongest modifications considered. They should thus not be used to model the highly nonlinear scales of structure formation in the context of forthcoming LSS analyses without considering an appropriate theoretical error contribution to the error budget (see Audren et al., 2013, for example).

The goal of this work was to validate different methods to compute the nonlinear matter power spectrum boost (see Equation 32). This function inherently depends on the nonlinear matter power spectrum of $\Lambda$CDM. But further, the boost must be applied to an accurate $\Lambda$CDM spectrum prediction in order to get a nonlinear modified matter power spectrum prediction. Therefore, the final modified gravity prescription inherits a dependence on predictions of the standard model. While we now have state-of-the-art high resolution tools to evaluate $P_{\rm NL}^{\rm\Lambda CDM}(k,z)$, the region in which these tools have internal accuracy within $1\%-2\%$ may not be as broad as we need for extracting unbiased constraints on cosmological parameters for Stage-IV LSS surveys (see Gordon et al., 2024, for a more in depth discussion). Furhtermore it is expected that in beyond-$\Lambda$CDM analysis, extreme regions of the parameter space need to be sampled, which heightens the need for the development of more comprehensive emulators in $\Lambda$CDM, as these are also imperative for the study of beyond-$\Lambda$CDM models.

In a similar vein, a further investigation of the impact of baryons in a full parameter inference scenario remains imperative to perform using the codes validated in this work. It has been shown that the interplay between baryonic physics and cosmology exhibit some dependence at small scales (Elbers et al., 2024). However, it is unknown to what extent in the nonlinear regime we can still extract relevant cosmological information to improve our constraints, i.e., if we need to model baryonic physics deep inside the nonlinear regime, $k\sim 10~{}h/$Mpc or not. Alternatively to modelling baryonic physics at the level of the power spectrum, it would be interesting to investigate the performance of procedures that mitigate the impact of baryons in the parameter constrains, such as the methods described in Refs. Eifler et al. (2015); Huang et al. (2019); Huang et al. (2021).

To conclude, let us highlight that the methods compared in this work have been designed with an element of theoretical flexibility in mind. There is a general shift to move beyond hard-coded codes designed to run with only one gravity model, and instead build more general tools that can be calibrated to a range of different models. This is an essential step forward to streamline the testing of new theoretical ideas with data from Stage IV surveys.

Data Availability

The halo model reaction 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 a Mathematica notebook, kmouflage.nb, which contains relevant calculations for the K-mouflage model. $N$-body matter power spectra measurements are available upon request.

Appendix A K-mouflage ReACT patch

Here we present the expressions needed to calculate the halo model reaction (see Equation 27) in the K-mouflage model. The halo model reaction relies on both the halo model (see Cooray & Sheth, 2002, for a review) and 1-loop perturbation theory (see Bernardeau et al., 2002, for a review). In particular, besides the background expansion $H(a)$, we require the modifications to the 1st, 2nd, 3rd order perturbative and nonlinear Poisson equations, as well as contributions to the potential energy of halos in order to solve the virial theorem (see Cataneo et al., 2019; Bose et al., 2020, for more details). K-mouflage also comes with a friction term correction to the Euler equation (Brax & Valageas, 2014b).

A.1 Background

For the background expansion we must solve the Klein-Gordon and Friedmann equations simultaneously. We do this numerically in ReACT as done in Ref. Hernández-Aguayo et al. (2022). The Friedmann equations are

$$\frac{\mathcal{H}^{2}}{H_{0}^{2}}\left[1-\frac{\varphi^{\prime}{}^{2}}{6}\right]= \frac{A(\varphi)\Omega_{\rm m,0}}{a}+\frac{1}{3}\lambda^{2}a^{2}$$

$$+\frac{(2n-1)}{3}\lambda^{2}a^{2}K_{0}\left(\frac{\varphi^{\prime}{}^{2}}{2\lambda^{2}a^{2}}\right)^{n}\frac{\mathcal{H}^{2n}}{H_{0}^{2n}}\,,$$ (33)

$$\frac{{\rm d}\mathcal{H}}{{\rm d\tau}}\frac{1}{H_{0}^{2}}= -\frac{A(\varphi)\Omega_{\rm m,0}}{2a}+\frac{1}{3}\lambda^{2}a^{2}-\frac{1}{3}\varphi^{\prime}{}^{2}\frac{\mathcal{H}^{2}}{H_{0}^{2}}$$

$$-\frac{(n+1)}{3}\lambda^{2}a^{2}K_{0}\left(\frac{\varphi^{\prime}{}^{2}}{2\lambda^{2}a^{2}}\right)^{n}\frac{\mathcal{H}^{2n}}{H_{0}^{2n}}\,,$$ (34)

while the Klein-Gordan equation is given as

$$(K_{X}+2\bar{X}K_{XX})\left[\frac{\mathcal{H}^{2}}{H_{0}^{2}}\varphi^{\prime\prime}+\frac{{\rm d}\mathcal{H}}{{\rm d\tau}}\frac{1}{H_{0}^{2}}\varphi^{\prime}\right]$$

$$+2(K_{X}-\bar{X}K_{XX})\frac{\mathcal{H}^{2}}{H_{0}^{2}}\varphi^{\prime}+3\frac{{\rm d}A(\varphi)}{\rm d\varphi}\frac{\Omega_{\rm m,0}}{a}=0\,,$$ (35)

where $K_{XX}={\rm d}^{2}K/{\rm d}X^{2}$ and we recall primes denote derivatives with respect to $\ln a$. In these equations we have defined the normalised scalar field $\varphi\equiv\phi/M_{\rm pl}$ and used the conformal Hubble rate $\mathcal{H}(a)=H(a)a$. $\tau$ is conformal time. We indicate that a few typographical errors did occur in Ref. Hernández-Aguayo et al. (2022) which have been corrected in the above equations.

To solve these equations we first find the analytic solutions to Equation 33 for a given value of $n$ . For $n=2$ this is a quadratic equation in $\mathcal{H}^{2}/H_{0}^{2}$. Then, for a given value of $a$ (or $\ln a$) we can substitute $\mathcal{H}$ and Equation 34 in Equation 35, enabling us to solve for the entire evolution of $\varphi$ (and $\varphi^{\prime}$), and consequently $H(a)$.

A.2 Perturbations

The linear modification to the Poisson equation is given by (Brax & Valageas, 2014b)

$$\frac{G_{\rm eff,L}}{G_{\rm N}}=A(\varphi)\left(1+\frac{2\beta_{\rm K}^{2}}{K_{X}}\right)\,.$$ (36)

Here we have included the conformal factor $A(\varphi)$, that comes along with $\rho_{\rm m}$ in the Poisson equation, Equation 4. Note that $G_{\rm eff,L}/G_{\rm N}=\mu(k,a)$ in the ReACT standard notation of Refs. Bose & Koyama (2016); Cataneo et al. (2019); Bose et al. (2020, 2022) for example.

The 2nd and 3rd order symmetrised modifications to the Poisson equation, in the same notation of Ref. Bose & Koyama (2016); Bose et al. (2020), are (Brax & Valageas, 2014b)

$$\gamma_{2}(\mbox{\boldmath$k$}_{1},\mbox{\boldmath$k$}_{2},a)= 0\,,$$

$$\gamma_{3}(\mbox{\boldmath$k$}_{1},\mbox{\boldmath$k$}_{2},\mbox{\boldmath$k$}_{3},a)= -\frac{9}{2}K_{XX}\left(\frac{A(\varphi)\Omega_{\rm m,0}}{a}\frac{H_{0}^{2}}{\mathcal{H}^{2}}\right)^{3}\left(\frac{\beta_{\rm K}}{K_{X}}\right)^{4}\frac{\mathcal{H}^{4}}{H_{0}^{2}}\frac{1}{a^{2}\lambda^{2}}$$

$$\times\left[\frac{(\mu_{12}+2\mu_{13}\mu_{23})}{k_{1}k_{2}}+\frac{(\mu_{13}+2\mu_{23}\mu_{12})}{k_{1}k_{3}}\right.$$

$$\left.+\frac{(\mu_{23}+2\mu_{13}\mu_{12})}{k_{2}k_{3}}\right]\,,$$ (37)

where we have defined $\mu_{ij}\equiv\hat{k}_{i}\cdot\hat{k}_{j}$ and $k_{i}=|\mbox{\boldmath$k$}_{i}|$.

Lastly, we also have a modification to the Euler equation in the form of a friction term (Brax & Valageas, 2014b). Similar terms have been included in ReACT in the context of interacting dark energy models (Simpson, 2010; Baldi & Simpson, 2015; Bose et al., 2018; Carrilho et al., 2022). In the K-mouflage model considered here, this term is given as

$$A_{\rm friction}=\beta_{\rm K}\varphi^{\prime}\,.$$ (38)

This term enters the Euler equation as expressed in Equation. 2.10 of Ref. Bose et al. (2018) for example.

A.3 Spherical collapse

The halo model reaction also requires us to solve for the spherical top-hat overdensity. This involves solving the evolution equation for the top-hat radius which requires specification of the nonlinear Poisson equation. The modification to this equation is to a good approximation equal to the linear modification at extra galactic scales given the smallness of the K-mouflage radius (Brax & Valageas, 2014b)

$$\frac{G_{\text{eff}}(k,a)}{G_{\rm N}}=\frac{G_{\text{eff, L}}(a)}{G_{\rm N}}\,.$$ (39)

We note in the notation of Ref. Cataneo et al. (2019), $\mathcal{F}=G_{\rm eff}/G_{\rm N}-1=\Delta G_{\rm eff}/G_{\rm N}$. In ReACT $\mathcal{F}$ appears as $1+\mathcal{F}$ in the Poisson equation. This yields the correct conformal factor accounting for the Einstein-frame transformation of the background density in the nonlinear Poisson equation, as it is already explicit in Equation 36.

Lastly, we note that the top-hat radius evolution also must include the friction term Equation 38.

A.4 Virial theorem

Here we present the potential energy contributions to the virial theorem in the K-mouflage model considered. This is needed to calculate the virial mass in the halo model reaction calculations. The specific components we require are (Schmidt et al., 2010; Cataneo et al., 2019)

$$\frac{W_{\rm N}}{E_{0}} =-\Omega_{\rm m,0}\frac{a^{-1}}{a_{\rm i}^{2}}y^{2}(1+\delta)\,;$$ (40)

$$\frac{W_{\rm\phi}}{E_{0}} =-\Omega_{\rm m,0}\mathcal{F}\,\frac{a^{-1}}{a_{\rm i}^{2}}y^{2}\delta\,;$$ (41)

$$\frac{W_{\rm eff}}{E_{0}} =-\frac{1}{3M_{\rm pl}^{2}H_{0}^{2}}(1+3w_{\rm eff})\,\bar{\rho}_{\rm eff}\frac{a^{2}}{a_{\rm i}^{2}}y^{2}\,;$$ (42)

$$\frac{W_{\rm fric}}{E_{0}} =-2A_{\rm friction}\frac{H^{2}}{H_{0}^{2}}\frac{a^{2}}{a_{\rm i}^{2}}\,y\,y^{\prime}\,,$$ (43)

where $y\equiv\frac{R_{\rm TH}}{R_{\rm i}}\frac{a_{\rm i}}{a}$, $R_{\rm TH}$ being the comoving top-hat radius, $R_{\rm i}$ the initial top-hat radius and $E_{0}$ is a normalisation. These represent the Newtonian contribution, the scalar field contribution, the effective dark energy contribution and a frictional force contribution as derived in Ref. Carrilho et al. (2022). In the K-mouflage model the scalar field affects both force enhancement and acts as an effective dark energy component.

We recall that $\mathcal{F}=G_{\rm eff}/G_{\rm N}-1$ which does not account for the correct conformal factor to appear in Equation 41 in the K-mouflage model. In this case we should have

$$\frac{W_{\rm\phi}}{E_{0}} =-\Omega_{\rm m,0}\left[G_{\rm eff,L}/G_{\rm N}-A(\varphi)\right]\,\frac{a^{-1}}{a_{\rm i}^{2}}y^{2}\delta\,,$$

$$=-\Omega_{\rm m,0}\left[A(\varphi)\frac{2\beta_{\rm K}^{2}}{K_{X}}\right]\,\frac{a^{-1}}{a_{\rm i}^{2}}y^{2}\delta\,,$$ (44)

where we used Equation 39 and Equation 36. $A_{\rm friction}$ is given by Equation 38. $w_{\rm eff}=\bar{p}_{\rm eff}/\bar{\rho}_{\rm eff}$ and $\bar{\rho}_{\rm eff}$ are the equation of state and energy density of the effective dark energy fluid component, with $\bar{p}_{\rm eff}$ being the fluid’s pressure. These are given in the Einstein frame by (Brax & Valageas, 2014a, b):

$$\bar{\rho}_{\rm eff} =-M_{\rm pl}^{2}H_{0}^{2}\lambda^{2}(K-M_{\rm pl}^{2}H^{2}\bar{\phi}^{\prime}{}^{2}K_{X})\,;$$ (45)

$$\bar{p}_{\rm eff} =M_{\rm pl}^{2}H_{0}^{2}\lambda^{2}K\,.$$ (46)

We then get

$$w_{\rm eff}=-\frac{K}{K-M_{\rm pl}^{2}H^{2}\bar{\phi}^{\prime}{}^{2}K_{X}}\,.$$ (47)

We can simplify Equation 42 further by noting that when adopting the model in Equation 14 we have

$$K_{X}=\frac{1}{H_{0}^{2}\lambda^{2}M_{\rm pl}^{2}}+K_{0}\frac{1}{H_{0}^{2n}\lambda^{2n}M_{\rm pl}^{2n}}n\,X^{n-1}\,.$$ (48)

Substituting this into Equation 47, we find the effective dark energy contribution to the potential energy is given by

$$\frac{W_{\rm eff}}{E_{0}}=-\frac{\lambda^{2}}{3}\left[2K+\frac{H^{2}}{H_{0}^{2}}\frac{\varphi^{\prime}{}^{2}}{\lambda^{2}}\left(1+K_{0}nX^{n-1}\right)\right]\frac{a^{2}}{a_{\rm i}^{2}}y^{2}\,.$$ (49)

Finally, we should note that the Newtonian contribution also should have a conformal factor along with $\Omega_{\rm m,0}$

$$\frac{W_{\rm N}}{E_{0}}=-A(\varphi)\Omega_{\rm m,0}\frac{a^{-1}}{a_{\rm i}^{2}}y^{2}(1+\delta)\,.$$ (50)

Appendix B Parametrised post-Friedmannian expressions

Here we review expressions for the general parametrisation of the effective gravitational constant appearing in the nonlinear Poisson equation as described in Ref. Lombriser (2016). This is based on the parametrised post-Friedmannian framework and is the means of modelling modifications to nonlinear structure formation in the MG-evolution code. This parametrisation has also been implemented in the ReACT code (Bose et al., 2022).

MG-evolution adopts a generalised form of the Vainshtein screening effect given by (Lombriser, 2016)

$$\frac{\Delta G_{\text{eff, NL}}}{G}=b\left(\frac{k_{*}}{k}\right)^{a_{f}}\left\{\left[1+\left(\frac{k}{k_{*}}\right)^{a_{f}}\right]^{1/b}-1\right\},$$ (51)

where NL stands for nonlinear, and $k_{*}$ and $b$, respectively, characterise the effective screening wavenumber and the interpolation rate between the screened and unscreened regimes. This expression augments the linear theory prediction as given in Equation 22 to give the full solution for $G_{\rm eff}$. We shall briefly provide the particular form of this expression for the three models considered in this work and refer the reader to Ref. Lombriser (2016) for more details.

B.1 nDGP

To parametrise nDGP gravity we consider (see Lombriser, 2016)

$$\frac{\Delta G_{\rm nDGP,NL}}{G_{\rm N}}=\frac{1}{3\beta(a)}\left(\frac{k_{*}}{k}\right)^{3}\left\{\left[1+\left(\frac{k}{k_{*}}\right)^{3}\right]^{\frac{1}{2}}-1\right\},$$ (52)

where $k_{*}$ corresponds approximately to the Vainshtein radius:

$$r_{*}=\frac{16\pi G_{\rm N}\delta Mr_{c}^{2}}{9\beta^{2}},$$ (53)

where $\delta M$ is the mass enclosed by a spherical region, $r_{c}$ is the crossover scale in nDGP theories, and $\beta(a)$ is given below. The Vainshtein radius effectively defines a region where the fifth force introduced by the scalar field gets shielded. The effective screening wavenumber $k_{*}$ can in principle be modelled. However, it is treated as a free parameter in MG-evolution. The function $\beta(a)$ reads as

$$\beta(a)=1+2Hr_{c}\Big{(}1+\frac{\dot{H}}{3H^{2}}\Big{)}\,,$$ (54)

and the linear effective gravitational constant in nDGP is given by

$$\frac{G_{\rm eff,\ L}}{G_{\rm N}}=1+\frac{1}{3\beta(a)}.$$ (55)

We remind the reader that a cosmological background that matches that of $\Lambda$CDM is assumed.

B.2 Cubic Galileon

To parametrise the Cubic Galileon we adopt Equation 22, with the nonlinear parametrisation

$$\frac{\Delta G_{\text{CG, NL}}}{G_{\rm N}}=\left(\frac{k_{*}}{k}\right)^{3}\left\{\left[1+\left(\frac{k}{k_{*}}\right)^{3}\right]^{\frac{1}{2}}-1\right\}\,.$$ (56)

To obtain the linear regime parametrisation we use the effective gravitational potential in Cubic Galileon theory in the linear regime, which reads as (Barreira et al., 2013a)

$$\frac{\Delta G_{\text{CG, L}}}{G_{\rm N}}=-\frac{2}{3}\frac{c_{3}\dot{\phi}^{2}}{M_{\mathrm{pl}}\mathcal{M}^{3}\beta_{2}}\,,$$ (57)

where $\phi$ is the Galileon scalar field. $\beta_{2}$ and $\mathcal{M}^{3}$ read,

$$\beta_{2} \equiv 2\frac{\mathcal{M}^{3}M_{\mathrm{pl}}}{\dot{\phi}^{2}}\beta_{1},$$ (58)

$$\mathcal{M}^{3} \equiv M_{\mathrm{pl}}H_{0}^{2}\,,$$ (59)

where

$$\beta_{1}\equiv\frac{1}{6c_{3}}\left[-c_{2}-\frac{4c_{3}}{\mathcal{M}^{3}}(\ddot{\phi}+2H\dot{\phi})+2\frac{c_{3}^{2}}{M_{\mathrm{pl}}^{2}\mathcal{M}^{6}}\dot{\phi}^{4}\right]\,.$$ (60)

We consider $c_{2}=-1$ and we use the tracker solution (Bellini et al., 2018),

$$\xi\equiv\frac{\dot{\phi}H}{M_{\mathrm{pl}}H_{0}^{2}}\,,$$ (61)

where $\xi$ is a constant and can be written in terms of $c_{2},c_{3}$,

$$\xi=-\frac{c_{2}}{6c_{3}}=\frac{1}{6c_{3}}\,.$$ (62)

As a result, we have the following solutions for $\dot{\phi}$ and $\ddot{\phi}$

$$\dot{\phi}=\frac{\xi M_{\mathrm{pl}}H_{0}^{2}}{H},\;\ddot{\phi}=-\frac{\xi M_{\mathrm{pl}}H_{0}^{2}\dot{H}}{H^{2}}\,.$$ (63)

Following the discussion presented in Ref. Barreira et al. (2013b) for the background tracker solution, we can derive the Hubble expansion rate as a function of scale factor

$$\frac{H^{2}}{H_{0}^{2}} =\frac{1}{2}\Big{[}\left(\Omega_{\rm m,0}a^{-3}+\Omega_{\rm r,0}a^{-4}\right)$$ (64)

$$+\sqrt{\left(\Omega_{\rm m,0}a^{-3}+\Omega_{\rm r,0}a^{-4}\right)^{2}+4\left(1-\Omega_{\rm m,0}-\Omega_{\rm r,0}\right)}\Big{]}\,,$$

where $H_{0}^{2}=\frac{8\pi G}{3}$ in MG-evolution units and $\Omega_{\rm r,0}$ is the radiation energy density fraction today. Computing the cosmic time derivative results in,

$$\frac{\dot{H}+H^{2}}{H_{0}^{2}}= -\frac{(a\Omega_{\rm m,0}+2\Omega_{\rm r,0})}{4a^{2}}$$

$$-\frac{(a\Omega_{\rm m,0}+\Omega_{\rm r,0})(3a\Omega_{\rm m,0}+4\Omega_{\rm r,0})}{4a^{6}\sqrt{4(1-\Omega_{\rm m,0}-\Omega_{\rm r,0})+\frac{(a\Omega_{\rm m,0}+\Omega_{\rm r,0})^{2}}{a^{8}}}}$$

$$+\frac{a^{2}}{2}\sqrt{4(1-\Omega_{\rm m,0}-\Omega_{\rm r,0})+\frac{(a\Omega_{\rm m,0}+\Omega_{\rm r,0})^{2}}{a^{8}}}\,.$$

B.3 K-mouflage

Here we derive the effective gravitational constant appearing in the nonlinear Poisson equation for the K-mouflage model described in subsection 2.3. We follow Ref. Lombriser (2016). We also note that a conformal factor, $A(\varphi)$, still needs to be applied to transform the density appearing in the nonlinear Poisson equation which is not included in the $G_{\rm KM,eff}$ expressions below.

The effective modification assuming a spherically symmetric matter distribution is given as (Winther & Ferreira, 2015a; Brax & Valageas, 2014b)

$$\frac{\Delta G_{\rm KM,eff}}{G_{\rm N}}=\frac{2\,\beta_{\rm K}^{2}}{K_{X}H_{0}^{2}\lambda^{2}M_{\rm pl}^{2}}\,.$$ (65)

Using the Klein-Gordon equation for a spherically symmetric matter distribution we can write $K_{X}$ as

$$K_{X}^{2}X=-\frac{2\,\beta_{\rm K}^{2}}{H_{0}^{4}\lambda^{4}M_{\rm pl}^{2}}F_{\rm N}^{2}\,,$$ (66)

where the Newtonian force is just $F_{\rm N}=G_{\rm N}M(<r)/r^{2}$, $r$ being the physical radial coordinate and $M(<r)$ being the mass enclosed in radius $r$. Substituting for $K_{X}$ in Equation 65 gives

$$\frac{\Delta G_{\rm KM,eff}}{G_{\rm N}}=\frac{\beta_{\rm K}}{M_{\rm pl}}\frac{1}{F_{\rm N}}\sqrt{-2X}\,.$$ (67)

Now to solve for $X$ we can adopt the model in Equation 14. By using Equation 48 to solve Equation 66 we get

$$X=\frac{H_{0}^{2}\lambda^{2}M_{\rm pl}^{2}}{6K_{0}}\frac{\left[1-f(x)\right]^{2}}{f(x)}\,,$$ (68)

where

$$f(x)=\left(1+x+\sqrt{x(x+2)}\right)^{\frac{1}{3}}\,,$$ (69)

and we have defined $x\equiv-C_{A}/r^{4}$, $C_{A}$ being a parameter proportional to $K_{0}$, defined as

$$C_{A}\equiv\frac{54\beta_{\rm K}^{2}G_{\rm N}^{2}M^{2}}{H_{0}^{2}\lambda^{2}}K_{0}\,.$$ (70)

We have written $M=M(<r)$ for compactness. It should be pointed out that $x\in(-2,0)$ yields no solution for $X$ which can be a problem for for very large $r$ and $K_{0}>0$. This will not generally be an issue as we look for solutions in the nonlinear regime.

Substituting Equation 68 into Equation 67 gives

$$\frac{\Delta G_{\rm KM,eff}}{G_{\rm N}}=C_{B}\frac{1-f(x)}{\sqrt{xf(x)}}\,,$$ (71)

where

$$C_{B}\equiv 3\sqrt{2}\beta_{\rm K}^{2}\,.$$ (72)

We note that in Ref. Lombriser (2016) there seems to be a missing factor of $1/(2\sqrt{2})$ in Equation. 3.26 in order to have the identification $C_{A}=C_{2}^{2}$. We allow here the case when $x\leq 0$ which can occur for $K_{0}>0$. Further, we note $C_{1}=-C_{B}$, $C_{1}$ and $C_{2}$ being the equivalent quantities for $C_{A}$ and $C_{B}$ in Ref. Lombriser (2016). We include a Mathematica notebook with our derivations.

We now derive the PPF expression from the limits of Equation 71:

$$\frac{\Delta G_{\rm KM,eff}}{G_{\rm N}}\rightarrow 2\beta_{\rm K}^{2}\,\qquad\quad\qquad {\rm for}\qquad|x|\ll 1\,\,({\rm i.e.}\,r^{4}\gg|C_{A}|)\,,$$ (73)

$$\frac{\Delta G_{\rm KM,eff}}{G_{\rm N}}\rightarrow\frac{C_{B}r^{4/3}}{(-C_{A})^{1/3}}\,\qquad {\rm for}\qquad|x|\gg 1\,\,({\rm i.e.}\,r^{4}\ll|C_{A}|)\,,$$ (74)

which are the same limits obtained by Ref. Lombriser (2016).

We now map these onto the (real space) PPF expression, Equation. 5.3 of Ref. Lombriser (2016)

$$\frac{\Delta G_{\rm eff}}{G_{\rm N}}=p_{1}p_{2}\frac{(1+s^{a_{f}})^{\frac{1}{p_{1}}}-1}{s^{a_{f}}}\,,$$ (75)

where

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

and $s=y_{\rm scr}/y_{\rm h}$. $y$ is the normalised top-hat radius

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

$R_{\rm TH}$ and $R_{i}$ are the comoving halo top-hat radius and $a_{i}$ the initial scale factor. The dimensionless screening scale is given by

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

$y_{\rm env}$ refers to the normalised radius of the environment and $M_{\rm vir}$ is the virial mass of the halo.

Comparing Equation 75 and Equation 78 with Equation 73 and Equation 74 we find, for a choice of $p_{1}$,

$$p_{2} =\frac{2\beta_{\rm K}^{2}}{p_{1}},\qquad p_{3}=\frac{4}{3}\frac{p_{1}-1}{p_{1}},$$

$$p_{4} =\left[\frac{-\sqrt{2}K_{0}p_{1}^{3}\beta_{\rm K}^{2}}{\lambda^{2}}\right]^{\frac{1}{4}}\Omega_{\rm m,0}^{\frac{1}{3}},$$

$$p_{5} =-1,\qquad p_{6}=1/6,\qquad p_{7}=0\,.$$ (79)

We note that whether there’s a $p_{1}$ in $p_{3}$ depends on whether $p_{1}$ is positive or negative (see Equation. 2.14 of Ref. Lombriser (2016)). We have also used $M_{\rm vir}\approx 4\pi\Omega_{\rm m,0}\rho_{\rm crit}R_{\rm th}^{3}/3$.