On the road to percent accuracy: nonlinear reaction of the matter power spectrum to dark energy and modified gravity

TL;DR

This paper introduces a method combining the halo model and perturbation theory to accurately predict the nonlinear matter power spectrum for dark energy and modified gravity models, achieving better than 1% accuracy across various scenarios.

Contribution

The authors develop a parameter-free approach that accurately predicts the nonlinear matter power spectrum for diverse dark energy and modified gravity models using existing simulations.

Findings

Achieves better than 1% accuracy for dark energy models up to k ≈ 1 h/Mpc.

Predicts the nonlinear matter power spectrum for DGP and f(R) gravity within 3% accuracy.

Including halo mass function data improves accuracy to 1%.

Abstract

We present a general method to compute the nonlinear matter power spectrum for dark energy and modified gravity scenarios with percent-level accuracy. By adopting the halo model and nonlinear perturbation theory, we predict the reaction of a $\Lambda$CDM matter power spectrum to the physics of an extended cosmological parameter space. By comparing our predictions to $N$-body simulations we demonstrate that with no-free parameters we can recover the nonlinear matter power spectrum for a wide range of different $w_0$-$w_a$ dark energy models to better than 1% accuracy out to $k \approx 1 \, h \, {\rm Mpc}^{-1}$. We obtain a similar performance for both DGP and $f(R)$ gravity, with the nonlinear matter power spectrum predicted to better than 3% accuracy over the same range of scales. When including direct measurements of the halo mass function from the simulations, this accuracy improves to 1%. With a single suite of standard $\Lambda$CDM $N$-body simulations, our methodology provides a direct route to constrain a wide range of non-standard extensions to the concordance cosmology in the high signal-to-noise nonlinear regime.

Figures

Tables

Table 1. Table 1: Equation of state parameters defining the dark energy models used in this work.

Model	$w_0$	$w_a$
DE1	$-0.7$	$0$
DE2	$-1.3$	$0$
DE3	$-1$	$0.5$
DE4	$-1$	$-0.5$
DE5	$-0.7$	$-1.5$
DE6	$-1.3$	$0.5$

Table 2. Table 2: Main technical details of the simulations employed in this work. The Nyquist frequency k Ny = π ​ N p / L box subscript 𝑘 Ny 𝜋 subscript 𝑁 𝑝 subscript 𝐿 box k_{\rm Ny}=\pi N_{p}/L_{\rm box} , ϵ italic-ϵ \epsilon is the force resolution, and m p subscript 𝑚 𝑝 m_{p} the particle mass. The standard cosmological parameters for the real f ​ ( R ) 𝑓 𝑅 f(R) and nDGP simulations, as well as for their standard Λ Λ \Lambda CDM counterparts, are ω b ≡ Ω b ​ h 2 = 0.02225 subscript 𝜔 b subscript Ω b superscript ℎ 2 0.02225 \omega_{\rm b}\equiv\Omega_{\rm b}h^{2}=0.02225 , ω c ≡ Ω c ​ h 2 = 0.1198 subscript 𝜔 𝑐 subscript Ω c superscript ℎ 2 0.1198 \omega_{c}\equiv\Omega_{\rm c}h^{2}=0.1198 , H 0 = 100 ​ h = 68 subscript 𝐻 0 100 ℎ 68 H_{0}=100h=68 km s -1 Mpc -1 , A s = 2.085 × 10 − 9 subscript 𝐴 s 2.085 superscript 10 9 A_{\rm s}=2.085\times 10^{-9} , n s = 0.9645 subscript 𝑛 s 0.9645 n_{\rm s}=0.9645 . For the evolving dark energy models we have instead ω b = 0.0245 subscript 𝜔 b 0.0245 \omega_{\rm b}=0.0245 , ω c = 0.1225 subscript 𝜔 c 0.1225 \omega_{\rm c}=0.1225 , H 0 = 70 subscript 𝐻 0 70 H_{0}=70 km s -1 Mpc -1 , σ 8 = 0.8 subscript 𝜎 8 0.8 \sigma_{8}=0.8 , n s = 0.96 subscript 𝑛 s 0.96 n_{\rm s}=0.96 .

model	$L_{\mathrm{box}}$	$N_p^3$	$k_{\mathrm{Ny}}$	$ϵ$	$m_p$	realisations	code
$f\left(R\right)$	512 $\mathrm{Mpc}{h}^{-1}$	${1024}^3$	6.3 $h{\mathrm{Mpc}}^{-1}$	15.6 kpc$h^{-1}$	$1.1\times {10}^{10}{M}_{\odot }{h}^{-1}$	1	ecosmog
nDGP	512 $\mathrm{Mpc}{h}^{-1}$	${1024}^3$	6.3 $h{\mathrm{Mpc}}^{-1}$	15.6 kpc$h^{-1}$	$1.1\times {10}^{10}{M}_{\odot }{h}^{-1}$	1	ecosmog
DE	200 $\mathrm{Mpc}{h}^{-1}$	${512}^3$	8 $h{\mathrm{Mpc}}^{-1}$	7.8 kpc$h^{-1}$	$5\times {10}^9{M}_{\odot }{h}^{-1}$	3	gadget-2

Table 3. Table 3: Refitted Sheth-Tormen halo mass function parameters obtained from the simulation data shown in Figs.  13 and 14 . Fits are consistent with p = 0 𝑝 0 p=0 in all cases. The standard values are { A , q , p } = { 0.3222 , 0.75 , 0.3 } 𝐴 𝑞 𝑝 0.3222 0.75 0.3 \{A,q,p\}=\{0.3222,0.75,0.3\} .

Model	$A$	$q$
nDGPm ($z=0$)	0.3427	0.819
nDGPm ($z=1$)	0.3067	0.757
nDGPw ($z=0$)	0.3347	0.819
nDGPw ($z=1$)	0.3023	0.754
pseudo-nDGPm ($z=0$)	0.3438	0.829
pseudo-nDGPm ($z=1$)	0.3057	0.761
pseudo-nDGPw ($z=0$)	0.3332	0.823
pseudo-nDGPw ($z=1$)	0.3013	0.753

Equations

$$S=\int {\mathrm{d}}^{4}x\sqrt{-g}\left\{\sum _{i=2}^5{\mathcal{L}}_i[\varphi ,g_{\mu \nu }]+{\mathcal{L}}_{\mathrm{m}}[\psi ,g_{\mu \nu }]\right\},$$

$${\mathcal{L}}_2$$

$${\mathcal{L}}_3$$

$${\mathcal{L}}_4$$

$${\mathcal{L}}_5$$

$${\mathcal{L}}_{\mathrm{H}}=K(\varphi ,X)+G_4(\varphi )R-G_3(\varphi ,X)\square \varphi .$$

$$\mathrm{d}{s}^2=g_{\mu \nu }\mathrm{d}{x}^{\mu }\mathrm{d}{x}^{\nu }=-(1+2\Psi )\mathrm{d}{t}^2+a^2(1+2\Phi )\mathrm{d}{\mathbf{x}}^2,$$

$$S_{\mathrm{EH}}=\int {\mathrm{d}}^{4}x\sqrt{-g}\frac{1}{16\pi G}\left[R+f(R)\right].$$

$$K$$

$$G_3$$

$$G_4$$

$${\nabla }^2\Psi$$

$${\nabla }^2\delta f_R$$

$$\nabla \Psi =\nabla {\Psi }_{\mathrm{N}}-\frac{1}{2}\nabla \delta f_R,$$

$$|\delta f_R|\ll \frac{2}{3}|{\Psi }_{\mathrm{N}}|,$$

$$f(R)=-2\Lambda -{\bar{f}}_{R0}\frac{{\bar{R}}_0^2}{R},$$

$$H^2=\frac{8\pi G}{3}({\bar{\rho }}_{\mathrm{m}}+{\bar{\rho }}_{\Lambda }),$$

$$K$$

$$G_3$$

$$G_4$$

$$\nabla \Psi =\nabla {\Psi }_{\mathrm{N}}+\frac{1}{2}\nabla \phi ,$$

$${\nabla }^2\phi +\frac{r_{\mathrm{c}}^2}{3\beta }\left[({\nabla }^2\phi {)}^2-({\nabla }_i{\nabla }_j\phi )({\nabla }^i{\nabla }^j\phi )\right]=\frac{8\pi G}{3\beta }\delta {\rho }_{\mathrm{m}},$$

$$\beta (a)\equiv 1+2H{r}_{\mathrm{c}}\left[1+\frac{\dot{H}}{3H^2}\right],$$

$$\delta M(r)=4\pi {\int }_0^r\mathrm{d}{r}^{\prime }{r}^{\prime 2}\delta {\rho }_{\mathrm{m}}(r^{\prime }).$$

$$\frac{\mathrm{d}\phi }{\mathrm{d}r}=\frac{4}{3\beta }{\left(\frac{r}{r_{\mathrm{V}}}\right)}^3\left[\sqrt{1+{\left(\frac{r_{\mathrm{V}}}{r}\right)}^3}-1\right]\frac{G\delta M(r)}{r^2},$$

$$r_{\mathrm{V}}(r)\equiv {\left[\frac{16r_{\mathrm{c}}^{2}G\delta M(r)}{9{\beta }^2}\right]}^{1/3}.$$

$$\left\{\begin{array}[]{ll}\frac{{\rm d}\varphi}{{\rm d}r}=\frac{2}{3\beta}\frac{{\rm d}\Psi_{\rm N}}{{\rm d}r}&\quad r\gg r_{\rm V}>R_{\rm\scriptscriptstyle TH}\,,\\ \frac{{\rm d}\varphi}{{\rm d}r}\approx 0&\quad R_{\rm\scriptscriptstyle TH}<r\ll r_{\rm V}\,.\end{array}\right.$$

$$K$$

$$G_3$$

$$G_4$$

$${\rho }_{\mathrm{DE}}$$

$$P_{\mathrm{DE}}$$

$$\frac{\mathrm{d}{\bar{\rho }}_{\mathrm{DE}}}{\mathrm{d}t}+3H{\bar{\rho }}_{\mathrm{DE}}(1+w)=0$$

$${\bar{\rho }}_{\mathrm{DE}}(a)={\bar{\rho }}_{\mathrm{DE},0}\exp \left[3{\int }_a^1\frac{\mathrm{d}{a}^{\prime }}{a^{\prime }}(1+w)\right],$$

$$w(a)=w_0+(1-a)w_a,$$

$$\frac{{\ddot{R}}_{\mathrm{TH}}}{R_{\mathrm{TH}}}=-\frac{4\pi G}{3}[{\bar{\rho }}_{\mathrm{m}}+(1+3w){\bar{\rho }}_{\mathrm{eff}}]-\frac{1}{3}{\nabla }^2\Psi .$$

Full text

††pubyear: 2018††pagerange: On the road to percent accuracy: nonlinear reaction of the matter power spectrum to dark energy and modified gravity–C

1 Introduction

General Relativity (GR) has been put under intense scrutiny in the Solar System, where it has successfully passed all tests (Will, 2014). Its application to cosmology, however, involves vastly different length scales and is comparable in orders of magnitude to an extrapolation from an atomic nucleus to the scale of human experience. It is therefore important to perform independent tests of our Theory of Gravity in the cosmological regime. Further motivation for a thorough inspection of cosmological gravity can be drawn from the necessity of a large dark sector in the energy budget of our Universe to explain large-scale observations with GR (Riess et al., 1998; Perlmutter et al., 1999; Aghanim et al., 2018; Hildebrandt et al., 2017; Abbott et al., 2018). In particular the late-time accelerated expansion of the cosmos has traditionally been an important driver for the development of alternative theories of gravity, a concept that has however become strongly challenged with the confirmation of the luminal speed of gravity (Lombriser & Taylor, 2016; Abbott et al., 2017b; Lombriser & Lima, 2017; Creminelli & Vernizzi, 2017; María Ezquiaga & Zumalacárregui, 2017; Baker et al., 2017; Sakstein & Jain, 2017; Battye et al., 2018; de Rham & Melville, 2018; Creminelli et al., 2018). Nevertheless, cosmic acceleration could be the result of a dark energy field permeating the Universe that may well be coupled to matter with an observable impact on cosmological scales. Importantly, should that coupling be universal, i.e. affecting baryons and dark matter equally, the corresponding models must then rely on the employment of screening mechanisms to comply with the stringent Solar-System bounds (Vainshtein, 1972; Khoury & Weltman, 2004; Babichev et al., 2009; Hinterbichler & Khoury, 2010). Signatures of screening are naturally to be expected in the nonlinear cosmological small-scale structure, where modified gravity transitions to GR, and for some models are even exclusively confined to these scales (Wang et al., 2012; Heymans & Zhao, 2018). The increasing wealth of high-quality data at these scales (Laureijs et al., 2011; LSST Dark Energy Science Collaboration, 2012; Hildebrandt et al., 2017; Abbott et al., 2018) renders cosmological tests of gravity a very timely enterprise. At the same time, cosmological structure formation proves notoriously difficult to model to sufficient accuracy in this regime, where high signal-to-noise measurements have the potential to distinguish a few percent deviation from GR (Heymans & Zhao, 2018).

For any given theory of gravity or dark energy model, our current best predictions for the statistical properties of the resulting matter distribution come from large-volume high-resolution $N$-body simulations (Oyaizu, 2008; Schmidt, 2009b; Zhao et al., 2011; Li et al., 2012; Brax et al., 2012; Baldi, 2012; Puchwein et al., 2013; Wyman et al., 2013; Barreira et al., 2013; Li et al., 2013b; Llinares et al., 2014; Winther et al., 2015). Running these, however, can take up to thousands of node-hours on dedicated cluster facilities, and although methods to partially alleviate this drawback exist (see, e.g., Barreira et al., 2015; Mead et al., 2015a; Valogiannis & Bean, 2017; Winther et al., 2017; Bose et al., 2017; Llinares, 2017) exploring vast swathes of the theory space remains currently unfeasible. Alternatively, analytical and semi-analytical methods can be used to swiftly predict specific large-scale structure observables, such as the matter power spectrum (Koyama et al., 2009; Schmidt et al., 2009; Schmidt et al., 2010; Li & Hu, 2011; Fedeli et al., 2012; Brax & Valageas, 2013; Lombriser et al., 2014; Barreira et al., 2014a, b; Zhao, 2014; Achitouv et al., 2016; Mead et al., 2016; Aviles & Cervantes-Cota, 2017; Cusin et al., 2018; Bose et al., 2018; Hu et al., 2018), with the important caveat that they have limited accuracy in the nonlinear regime of structure formation, and often involve some level of fitting to the same quantity measured in simulations. These approaches are therefore inadequate for future applications to high-quality data from Stage IV surveys (Laureijs et al., 2011; LSST Dark Energy Science Collaboration, 2012; Levi et al., 2013; Koopmans et al., 2015), where percent level accuracy over a wide range of scales will be necessary to obtain tight and unbiased constraints on departures from GR (Alonso et al., 2017; Casas et al., 2017; Reischke et al., 2018; Spurio Mancini et al., 2018; Taylor et al., 2018b) and the nature of dark energy (Albrecht et al., 2006). Matter power spectrum emulators can provide a solution to this problem for particular modified gravity or dark energy models (Heitmann et al., 2014; Lawrence et al., 2017; Euclid Collaboration et al., 2018; Winther et al., 2019), but still rely on the availability of large quantities of computational resources to the determine the properties of the matter power spectrum at the location of the emulator nodes. The absence of a clear attractive alternative to the $\Lambda$CDM paradigm calls for a more general framework, one easily adaptable to non-standard cosmologies beyond the handful of well studied cases.

Here we take an important step in this direction by extending the method proposed in Mead (2017), where the halo model is used to compute matter power spectrum ratios with respect to a convenient baseline cosmology. Mead (2017) showed that by determining these ratios, rather than the absolute value of the matter power spectrum, the shortcomings of the halo model are mitigated. The initial conditions of the baseline cosmology are designed so that under GR+$\Lambda$CDM evolution the linear clustering of matter at some given redshift exactly reproduces that of the target cosmology of interest, whose evolution is instead governed by non-standard laws of gravity and/or background expansion. Assuming one can generate an accurate nonlinear matter power spectrum for the reference cosmology (e.g with a suitable emulator), recovery of the target power spectrum then hinges on the computation of a ‘correction’ factor that incorporates the nonlinear effects of fifth forces, screening mechanisms and deviations from the cosmological constant. We use the halo model and nonlinear perturbation theory to obtain such corrections, and refer to this quantity as the reaction.

The paper is organised as follows. In Sec. 2 we briefly describe popular modified gravity and dark energy models used here as testbeds for our methodology. Sec. 3 reviews the halo model formalism and introduces the matter power spectrum reactions. The cosmological simulations used to validate our approach are described in Sec. 4, and Sec. 5 presents the capability of the halo model reactions to predict the nonlinear matter power spectrum. We summarise our conclusions in Sec. 6. We provide details of the spherical collapse and perturbation theory calculations employed in this work in App. A and App. B, respectively. Additional tests to gauge the importance of the halo mass function and halo concentration in our predictions are presented in App. C.

2 Dark energy and modified gravity theory

The most general four dimensional scalar-tensor theory with second-order equations of motion is described by the action (Horndeski, 1974; Deffayet et al., 2011; Kobayashi et al., 2011)

[TABLE]

where $g$ is the determinant of the metric $g_{\mu\nu}$ minimally coupled to a generic matter field $\psi$ (Jordan frame), ${\cal L}_{\rm m}$ is the matter Lagrangian, $\phi$ is the scalar degree of freedom, and the terms entering the Einstein-Hilbert Lagrangian are

[TABLE]

Here, $K$ and $G_{i}$ are arbitrary functions of $\phi$ and $X\equiv-\nabla^{\nu}\phi\nabla_{\nu}\phi/2$, and the subscripts $X$ and $\phi$ denote derivatives.

The nearly simultaneous detection of gravitational waves and electromagnetic signals emitted from two colliding neutron stars (Abbott et al., 2017a) imposes tight constraints on the present-day speed of gravitational waves $c_{T}$, i.e. $|c_{T}/c-1|\lesssim 10^{-15}$, where $c$ is the speed of light (Abbott et al., 2017b). Restricting ourselves to theories of gravity with non-evolving $c_{T}$, and requiring this not to be achieved by extreme fine tuning of the $G_{4}$ and $G_{5}$ functions (Lombriser & Taylor, 2016; Lombriser & Lima, 2017; Creminelli & Vernizzi, 2017; María Ezquiaga & Zumalacárregui, 2017; Baker et al., 2017; Sakstein & Jain, 2017; Battye et al., 2018; de Rham & Melville, 2018), implies that the remaining Horndeski Lagrangian takes the form (McManus et al., 2016)

[TABLE]

In this paper we focus on well-studied models of modified gravity and dark energy, each exploring the effects introduced by the individual terms in Eq. (3). Quintessence and k-essence dark energy models are described by a contribution of $K$ only (Sec. 2.3). $G_{4}$ introduces a coupling of this field to the metric that modifies gravity. The class of models described by this term encompasses the chameleon (Khoury & Weltman, 2004), symmetron (Hinterbichler & Khoury, 2010), and k-mouflage (Babichev et al., 2009) screening mechanisms, and we will study a particular example of this action with chameleon screening in Sec. 2.1 when considering a realisation in $f(R)$ gravity. The $G_{3}$ term appears, for instance, in the four-dimensional effective scalar-tensor theory of DGP braneworld gravity (Sec. 2.2) and gives rise to the Vainshtein screening mechanism (Vainshtein, 1972). Either $G_{3}$ or noncanonical kinetic contributions in $K$ produce a nonluminal sound speed of the scalar field fluctuations that can yield observable scale-dependent effects beyond the sound horizon. The mass scale associated with a scalar field potential in $K$ can furthermore introduce a scale-dependent growth of structure below the sound horizon. With $c_{T}=1$, a genuine self-acceleration of the cosmological background that is directly attributed to modified gravity must arise from $G_{4}$ (Lombriser & Taylor, 2016), which is however in tension with observations (Lombriser & Lima, 2017). A self-acceleration from $K$ or $G_{3}$ that dispenses with the need of a cosmological constant is, in contrast, still observationally feasible.

Throughout, we assume a flat Friedmann-Robertson-Walker (FRW) background, and the perturbed metric in the Newtonian gauge reads111Here and throughout we work in natural units, and set $c=1$.

[TABLE]

where $\Psi$ and $\Phi$ denote the two gravitational potentials, and $a$ is the scale factor. The evolution of non-relativistic matter perturbations is determined by $\Psi$, whereas photons follow the null geodesics defined by the lensing potential $\Phi_{-}=(\Psi-\Phi)/2$ (see, e.g., Carroll, 2004). For all models considered here $\Phi_{-}=\Psi_{N}$, where $\Psi_{N}$ is the standard Newtonian potential.

2.1 $f(R)$ gravity

In $f(R)$ gravity the Einstein-Hilbert action is modified to contain an additional non-linear function of the Ricci scalar $R$, that is

[TABLE]

The $f(R)$ Lagrangian is a particular case of Eq. (3) with the Horndeski functions (see, e.g., de Felice et al., 2011)

[TABLE]

where we defined the scalaron field $f_{R}\equiv{\rm d}f/{\rm d}R$, and used $\phi=(1+f_{R})/\sqrt{8\pi G}$. In the quasi-static regime222See, e.g., Noller et al. (2014), Bose et al. (2015) and Lagos et al. (2018) for a detailed discussion on the validity of the quasi-static approximation in modified gravity., structure formation is governed by the following coupled equations (e.g., Oyaizu, 2008)

[TABLE]

where $\delta\rho_{\rm m}=\rho_{\rm m}-\bar{\rho}_{\rm m}$, $\delta R=R-\bar{R}$ and $\delta f_{R}=f_{R}-\bar{f}_{R}$ are the matter density, curvature and scalaron perturbations with respect to their background averaged values. Eqs. (9) and (10) can be combined to give

[TABLE]

which explicitly shows that the scalar field fluctuations source an additional fifth force.

Since GR accurately describes gravity in our Solar System, viable modifications must also be compatible with local constraints. In $f(R)$ gravity this is achieved by means of the chameleon screening mechanism (Khoury & Weltman, 2004), which suppresses departures from standard gravity for large enough potential wells $\Psi_{\rm N}$. In practice, structures are screened if the thin shell condition,

[TABLE]

is satisfied. Stable theories require $f_{R}<0$ (Hu & Sawicki, 2007), thus the chameleon screening activates throughout an isolated object if $|\delta f_{R}|\leq|\bar{f}_{R}|\ll|\Psi_{\rm N}|$. Assuming that the Milky Way is placed in the cosmological background, and knowing that $|\Psi_{\rm MW}|\sim 10^{-6}$, this in turn imposes $|\bar{f}_{R0}|<10^{-6}$ for the present value of the background scalaron field.

Hereafter, we adopt the following $f(R)$ functional form (Hu & Sawicki, 2007)

[TABLE]

where $\Lambda$ is an effective cosmological constant driving the background cosmic acceleration, and $\bar{R}_{0}$ corresponds to the background Ricci scalar today. We will work with values $|\bar{f}_{R0}|=10^{-5}$ (F5) and $|\bar{f}_{R0}|=10^{-6}$ (F6), for which cosmological structures are, respectively, partially unscreened or largely screened throughout the cosmic history. Note that deviations from the $\Lambda$CDM expansion history are of order $\bar{f}_{R0}$ (Hu & Sawicki, 2007). Hence, for the $f(R)$ models considered here the background evolution is in effect equivalent to that of the concordance cosmology, with the Hubble parameter given by

[TABLE]

where $\bar{\rho}_{\Lambda}$ is the energy density of the cosmological constant. The large-scale structure data currently available allows amplitudes $|\bar{f}_{R0}|\lesssim 10^{-5}$ (Terukina et al., 2014; Lombriser, 2014; Cataneo et al., 2015; Liu et al., 2016; Alam et al., 2016), thus placing the F5 model on the edge of the region of parameter space still relevant for cosmological applications333See, however, the recent work by He et al. (2018) where it was showed that deviations as small as $|\bar{f}_{R0}|=10^{-6}$ could already be in strong tension with redshift space distortions data. At the present time, the tightest constraints on $f(R)$ gravity come from the analysis of kinematic data for the gaseous and stellar components in nearby galaxies, which only allows $|\bar{f}_{R0}|\lesssim 10^{-8}$ (Desmond et al., 2018). if other effects degenerate with the enhanced growth of structure are ignored. Accounting for massive neutrinos (Baldi et al., 2014) and baryonic feedback (Puchwein et al., 2013; Hammami et al., 2015; Arnold et al., 2018) will loosen the existing constraints (see, e.g., Hagstotz et al., 2018; Giocoli et al., 2018). In addition, alternative functional forms to Eq. (13) can lead to different upper bounds on $|\bar{f}_{R0}|$ (see, e.g., Cataneo et al., 2015).

2.2 DGP

In the DGP braneworld model the matter fields live on a four-dimensional brane embedded in a five-dimensional Minkowski space (Dvali et al., 2000). In this model the dimensionality of the gravitational interaction is controlled by the crossover scale parameter $r_{\rm c}$, such that on scales smaller than $r_{\rm c}$ DGP becomes a four-dimensional scalar-tensor theory described by an effective Lagrangian with terms (Nicolis & Rattazzi, 2004; Park et al., 2010)

[TABLE]

where the brane-bending mode $\varphi$ represents the scalar field. Hereafter we will be working with the normal branch DGP model (nDGP), which despite being a stable solution of the theory is also incompatible with the observed late-time cosmic acceleration. To obviate this problem the Lagrangian given by Eqs. (15)-(17) is extended to include a smooth, quintessence-type dark energy with a potential conveniently designed to match the expansion history of a flat $\Lambda$CDM cosmology (Schmidt, 2009b)444This is an assumption made to ease comparisons to $\Lambda$CDM simulations, and is not a strict observational requirement (cf. Lombriser et al., 2009).. Therefore, the Friedmann equation Eq. (14) applies here as well.

The scalar field $\varphi$ couples to non-relativistic matter by sourcing the dynamical potential $\Psi$, which in turn produces a gravitational force given by

[TABLE]

where the second term on the right hand side is the attractive fifth force contribution. On length scales $\lambda\ll H^{-1},r_{\rm c}$, and in the quasi-static regime (Schmidt, 2009a; Brito et al., 2014; Winther & Ferreira, 2015), the evolution of the brane-bending mode is described by (Koyama & Silva, 2007)

[TABLE]

with the function $\beta(a)$ defined as

[TABLE]

where overdots denote derivatives with respect to cosmic time. The derivative self-interactions in Eq. (19) suppress the field in high-density regions, where the matter density field is nonlinear, effectively restoring GR. This is the so-called Vainshtein screening. To explicitly illustrate how this mechanism works we shall consider a spherically symmetric overdensity with mass

[TABLE]

Then, for this system the gradient of $\varphi$ reads (Koyama & Silva, 2007; Schmidt et al., 2010)

[TABLE]

where we defined the Vainshtein radius

[TABLE]

The scale introduced in Eq. (23) sets the distance from the centre of the spherical mass distribution above which fifth force effects are observable. For instance, for a top-hat overdensity of radius $R_{\rm\scriptscriptstyle TH}$ one has the two limiting cases

[TABLE]

In the following we will consider the medium and weak nDGP variants used in Barreira et al. (2016), with crossover scales $r_{\rm c}H_{0}=0.5$ (nDGPm) and $r_{\rm c}H_{0}=2$ (nDGPw) in units of the present-day Hubble horizon $H_{0}^{-1}$. Note that, at present, of these two cases nDGPw is the only one compatible with growth rate data (Barreira et al., 2016). Hence, similarly to F5, the use of nDGPm will serve as a testbed for our methodology in conditions of relatively strong departures from GR.

2.3 Dark energy

The simplest models described by the Lagrangian in Eq. (3) are those in which the scalar field is minimally coupled to gravity, that is

[TABLE]

In this scenario, the field $\phi$ is associated with a fluid called dark energy (DE) with energy density and pressure (see, e.g., Amendola & Tsujikawa, 2010)

[TABLE]

respectively. Its background evolution is controlled by the equation of state parameter $w\equiv\bar{P}_{\rm\scriptscriptstyle DE}/\bar{\rho}_{\rm\scriptscriptstyle DE}$, and the solution to the continuity equation

[TABLE]

is given by

[TABLE]

where $\bar{\rho}_{{\rm\scriptscriptstyle DE},0}$ is the present-day dark energy density.

Popular models of dark energy, such as quintessence (Wetterich, 1988; Ratra & Peebles, 1988), k-essence (Armendariz-Picon et al., 2000) and clustering quintessence (Creminelli et al., 2009), belong to this subclass of theories. In this paper we restrict our discussion to a quintessence-like dark fluid with rest-frame sound speed $c_{\rm s}^{2}=1$ and equation of state (Chevallier & Polarski, 2001; Linder, 2003)

[TABLE]

where $\{w_{0},w_{a}\}$ are free phenomenological parameters. The relativistic sound speed washes out the dark energy perturbations on sub-horizon scales, resulting in modifications to the growth of structure tied solely to the different expansion history compared to $\Lambda$CDM. Our methodology can, in principle, also be applied to forms of dark energy clustering on small scales.

Table 1 summarises the dark energy models selected for this work, which have been chosen to roughly enclose the 2$\sigma$ region of parameter space allowed by the Planck 2015 temperature and polarization data in combination with baryon acoustic oscillations, supernova Ia and $H_{0}$ measurements (Planck Collaboration et al., 2016).

3 Matter power spectrum reaction

In Sec. 3.1 and Sec. 3.2 we briefly review the spherical collapse model and the halo model formalism, which we use to predict the nonlinear matter power spectrum for the range of cosmological models listed in Sec. 2. The halo model assumes that all matter in the Universe is localized in virialized structures, called halos. In this approach, the spatial distribution of these objects and their density profiles determine the statistics of the matter density field on all scales. It is typically assumed that each mass element belongs to one halo only, i.e. halos are spatially exclusive. Below we introduce the ingredients entering the halo model prescription, and refer the interested reader to the Cooray & Sheth (2002) review on the topic for more details. In Sec. 3.3 we then detail our new approach to reach percent level accuracy on these power spectra over a range of scales where the halo model alone is known to fail.

3.1 Spherical collapse model

The Press-Schechter formalism (Press & Schechter, 1974) approximates halo formation following the evolution of a spherical top-hat overdensity of radius $R_{\rm\scriptscriptstyle TH}$ and mass $M=4\pi R_{\rm\scriptscriptstyle TH}^{3}\bar{\rho}_{\rm m}(1+\delta)/3$ in an otherwise homogenous background. Mass conservation and the Euler equations imply (see, e.g., Schmidt et al., 2009)

[TABLE]

Here, $\bar{\rho}_{\rm eff}$ and $w$ are, respectively, the background energy density and equation of state of an effective dark energy component causing the late-time cosmic acceleration. Hence, in $f(R)$ gravity and nDGP, $\bar{\rho}_{\rm eff}=\bar{\rho}_{\Lambda}$ and $w=-1$. For the smooth dark energy models in Sec. 2.3 we have $\bar{\rho}_{\rm eff}=\bar{\rho}_{\rm\scriptscriptstyle DE}$ and $w$ given by Eq. (31). Modifications of gravity enter through the potential term in Eq. (32), which we parametrize as

[TABLE]

where $\mathcal{F}$ can depend on time, mass and environment. Eq. (33) reduces to the standard Poisson equation for $\mathcal{F}=0$, and expressions for $\mathcal{F}$ in $f(R)$ and nDGP cosmologies are given in App. A.

The mass fluctuation $\delta_{\rm i}$ at the initial time $a_{\rm i}$ within $R_{\rm\scriptscriptstyle TH}$ evolves as

[TABLE]

where $R_{\rm i}$ is the initial top-hat radius. Using Eqs. (32)-(34) we then find $\delta_{\rm i}$ such that collapse (i.e. $R_{\rm\scriptscriptstyle TH}=0$) occurs at a chosen time $a=a_{\rm coll}$. The Press-Schechter approach assumes that all regions in the initial density field with overdensities larger than $\delta_{i}$ have collapsed into halos by $a_{\rm coll}$. Equivalently, one can compare the linearly evolved initial fluctuations to the linearly extrapolated collapse overdensity $\delta_{\rm c}(a_{\rm coll})\equiv D_{\Lambda}(a_{\rm coll})\delta_{i}/a_{i}$, with (see, e.g., Dodelson, 2003)

[TABLE]

being the linear growth factor in $\Lambda$CDM555For cosmologies with a scale-independent linear growth, such as nDGP and $w$CDM, using the $\Lambda$CDM growth is simply a matter of convenience. In $f(R)$ gravity this approach has the advantage of preserving the statistics of the initial mass fluctuations. For more details see Cataneo et al. (2016)., and $\Omega_{\rm m}\equiv 8\pi G\bar{\rho}_{\rm m,0}/3H_{0}^{2}$.

In the idealised top-hat scenario, the spherical mass collapses to a point of infinite density. However, processes in the real Universe act so that, after turnaround, the overdensity eventually reaches virial equilibrium (see, e.g., Mo et al., 2010). Following Schmidt et al. (2010) (for an earlier work see also Maor & Lahav, 2005), we do not assume energy conservation during collapse, and compute the time of virialization, $a_{\rm vir}$, from the virial theorem alone (see App. A for details). This approach differs from previous works where changes induced by dark energy (Mead, 2017) or modified gravity (Lombriser et al., 2014) were neglected. The virial comoving radius $R_{\rm vir}$ of the formed halo can be derived from its virial mass

[TABLE]

knowing that the virial overdensity is given by

[TABLE]

with the mass fluctuation $\delta$ obtained from Eq. (34).

3.2 Halo model

The simplest statistics describing the clustering properties of the matter density field $\rho_{\rm m}({\bf x})$ is the 2-point correlation function or, its Fourier transform, the power spectrum $P(k)$ defined as

[TABLE]

where $\delta_{\rm D}$ denotes the Dirac delta function, and $\tilde{\delta}({\bf k})$ represents the Fourier transform of the matter density fluctuations relative to the background mean density, $\delta({\bf x})=\rho_{\rm m}({\bf x})/\bar{\rho}_{\rm m}-1$. Note that Eq. (38) assumes statistical homogeneity and isotropy.

In the halo model the matter power spectrum results from the contribution of correlations between halos ($P_{2\rm h}$) and those within halos ($P_{1\rm h}$), and can be written as666In this instance, and whenever the context is clear, we omit the time-dependence from our notation. However, we reintroduce it any time this can become a source of ambiguity.

[TABLE]

To properly account for these correlations we need to know the abundance of such halos. For any redshift $z$, the halo mass function provides the comoving number density of halos of mass $M_{\rm vir}$, and it is defined as

[TABLE]

where the peak height $\nu\equiv\delta_{\rm c}/\sigma$, and we adopt the Sheth-Tormen (ST) multiplicity function (Sheth & Tormen, 1999, 2002)

[TABLE]

Here, the normalization constant $A$ is found imposing that all mass in the Universe is confined into halos, i.e. $\int{\rm d}\nu f(\nu)=1$, and the remaining parameters take the $\Lambda$CDM standard values $q=0.75$ and $p=0.3$, unless stated otherwise. The variance of the linear density field smoothed with a top-hat filter of comoving radius $R$ enclosing a mass $M=4\pi R^{3}\bar{\rho}_{\rm m,0}/3$ is given by

[TABLE]

where $\tilde{W}$ is the Fourier transform of the top-hat filter, and $P_{\rm\scriptscriptstyle L}(k,z)$ is the $\Lambda$CDM linear power spectrum. At this point it is worth emphasising that in some GR extensions, besides its usual dependence on background cosmology and redshift, the spherical collapse threshold $\delta_{\rm c}$ can also vary with halo mass and environment (Li & Efstathiou, 2012; Li & Lam, 2012; Lam & Li, 2012; Lombriser et al., 2013, 2014). When appropriate we include both these dependencies in our modelling by following the approach of Cataneo et al. (2016), where the initial value of the environmental overdensity is derived from the peak of the environment probability distribution.

Halos are biased tracers of the underlying dark matter density field, and at the linear level the halo and matter density fields are connected by the relation $\delta_{\rm h}=b_{\rm\scriptscriptstyle L}\delta$. Adopting the ST mass function, the peak-background split formalism predicts the linear halo bias777Valogiannis & Bean (2019) recently found that in $f(R)$ gravity the linear halo bias contains an additional term accounting for the environmental dependence, which we omit in Eq. (43). Given the relative unimportance of the bias for our halo model reactions (see Secs. 3.3 and 5), this choice is, in effect, inconsequential for the accuracy of our predictions. (Sheth & Tormen, 1999)

[TABLE]

The last piece of information required by the halo model is a description of the matter distribution within halos. We adopt Navarro-Frenk-White (NFW) halo profiles (Navarro et al., 1996)

[TABLE]

where the scale radius $r_{\rm s}$ is parametrized through the virial concentration $c_{\rm vir}\equiv R_{\rm vir}/r_{\rm s}$, and the normalization $\rho_{\rm s}$ follows from the virial mass as

[TABLE]

Inside the virial radius, and for all cosmological models studied here, the NFW profiles are a good representation of the averaged halo profiles measured in simulations (Schmidt et al., 2009; Schmidt, 2009b; Zhao et al., 2011; Lombriser et al., 2012; Kwan et al., 2013; Shi et al., 2015; Achitouv et al., 2016).

In $\Lambda$CDM, $f(R)$ gravity and nDGP we model the $c$-$M$ relation as the power law

[TABLE]

fixing $c_{0}=9$ and $\alpha=0.13$ (Bullock et al., 2001), and $M_{\ast}$ is defined by $\nu(M_{\ast})=1$. In particular, for $f(R)$ gravity $M_{\ast}$ depends itself on the halo mass (Lombriser et al., 2014), which means the $c$-$M$ relation for these models is no longer described by a simple power law (Shi et al., 2015). For the smooth dark energy models in Sec. 2.3 we correct for the different expansion histories following Dolag et al. (2004), that is

[TABLE]

where $g_{\rm\scriptscriptstyle X}$ is the linear growth factor normalized to $z=0$ (see App. B). This correction reflects that halos collapse at different times in cosmological models with different growth histories. In cosmological models where halos collapse earlier these halos will be more concentrated compared to the same mass halos if they form later. In App. C we demonstrate that our results are insensitive to the correct shape of the $c$-$M$ relation on scales $k\lesssim 0.5$ $h\,{\rm Mpc}^{-1}$.

We can now predict the nonlinear matter power spectrum, and rewrite Eq. (39) as

[TABLE]

where, more explicitly,

[TABLE]

In the equations above, $u(k,M)$ corresponds to the Fourier transform of an NFW profile truncated at $R_{\rm vir}$, normalized such that $u(k\rightarrow 0,M)\rightarrow 1$. Note that from Eqs. (41) and (43) it follows that $\lim_{k\rightarrow 0}I(k)=1$.

3.3 Halo model reactions

The apparent simplicity and versatility of the halo model has contributed to its widespread use as a method to predict the nonlinear matter power spectrum in diverse scenarios. Examples include the $\Lambda$CDM cosmology (Seljak, 2000; Peacock & Smith, 2000; Giocoli et al., 2010; Valageas & Nishimichi, 2011; Valageas et al., 2013; Mohammed & Seljak, 2014; Seljak & Vlah, 2015; van Daalen & Schaye, 2015; Mead et al., 2015b; Schmidt, 2016), dark energy and modified gravity models (Schmidt et al., 2009; Schmidt et al., 2010; Li & Hu, 2011; Fedeli et al., 2012; Brax & Valageas, 2013; Lombriser et al., 2014; Barreira et al., 2014a, b; Achitouv et al., 2016; Mead et al., 2016; Hu et al., 2018), massive neutrinos (Abazajian et al., 2005; Massara et al., 2014; Mead et al., 2016), baryonic physics (Mohammed & Seljak, 2014; Fedeli, 2014; Fedeli et al., 2014; Mead et al., 2015b), alternatives to cold dark matter (Dunstan et al., 2011; Schneider et al., 2012; Marsh, 2016), and primordial non-Gaussianity (Smith et al., 2011). Its imperfect underlying assumptions are however responsible for inaccuracies that limit its applicability to future high-quality data (see, e.g., Figure 1 in Massara et al., 2014), where percent level accuracy is required in order to obtain unbiased cosmological constraints (Huterer & Takada, 2005; Eifler, 2011; Hearin et al., 2012; Taylor et al., 2018a).

To mitigate these downsides one can add complexity to the model at the expense of introducing new free parameters (see, e.g., Seljak & Vlah, 2015), fitting the existing ones to the matter power spectrum measured in simulations (see, e.g., Mead et al., 2015b), or sensibly increasing the computational costs by going beyond linear order in perturbation theory (see, e.g., Valageas & Nishimichi, 2011). Here, instead, we follow and extend the approach presented in Mead (2017), which we shall refer to as halo model reactions888Note that in Mead (2017) this is referred to as response. We use the term reaction to distinguish it from the quantities studied in Neyrinck & Yang (2013), Nishimichi et al. (2016) or Barreira & Schmidt (2017). Our and these other definitions are all conceptually analogous, in the sense that they describe how the nonlinear power spectrum responds to changes in some feature, which in our case is physics beyond the vanilla $\Lambda$CDM cosmology, e.g. fifth forces, evolving dark energy, massive neutrinos, baryons etc..

Our goal is to model the nonlinear power spectrum of fairly general extensions to the standard cosmology, a flat $\Lambda$CDM Universe with massless neutrinos. These cosmologies equipped with beyond-$\Lambda$CDM physics are what we will call real cosmologies. We use the halo model to determine the change (i.e. the reaction) that this new physics induces in a reference $\Lambda$CDM cosmology, for which simulations are considerably cheaper. Key to the success of our method is how this reference cosmology is defined, which is what we will call the pseudo cosmology. Essentially, this is a $\Lambda$CDM cosmology evolved with standard gravity up to a final redshift $z_{\rm f}$, with the additional property that its linear clustering of matter exactly matches that of the target real cosmology of interest at $z_{\rm f}$. In other words, the cold dark matter and the cosmological constant determine the expansion history and growth of structure of the pseudo cosmology, but the initial conditions (see Sec. 4) are adjusted so that

[TABLE]

The reaction function is then defined as the ratio of the nonlinear matter power spectrum in the real cosmology to that in the pseudo cosmology,

[TABLE]

and our corresponding halo model prediction takes the heuristic form999Note that we neglect the integral factor Eq. (50) in our two-halo terms. We checked that setting $I^{2}(k)=1$ for all scales has no measurable impact on our halo model reactions.

[TABLE]

Here, $\mathcal{E}\sim 1$ and $k_{\star}>0$ are parameters introduced to improve the accuracy of the halo model reactions in modified gravity theories. These are not free parameters, and we shall derive them using the halo model and standard perturbation theory below. However, let us first examine the general behaviour of Eq. (53):

On large linear scales $\mathcal{R}\rightarrow 1$ by definition101010This is not strictly true in the traditional halo model implementation we adopt in this work. In fact, the one-halo terms have a constant tail in the low-$k$ limit (see Eq. 54) that dominates the two-halo contributions on very large scales. In a consistent formulation of the halo model, however, where mass and momentum conservation are enforced, this tail disappears leaving only the theoretically motivated two-halo term (Schmidt, 2016). For our purposes we can simply ignore this inconsistency, and restrict the use of the reaction function to scales $k>0.01$ $h\,{\rm Mpc}^{-1}$.; 2. 2.

On small nonlinear scales $\mathcal{R}\approx P_{1\rm h}^{\rm real}/P_{1\rm h}^{\rm pseudo}$; 3. 3.

Quasi-linear scales $0.01\lesssim k\,{\rm Mpc}\,h^{-1}\lesssim 0.1$ are well described by perturbation theory, while intermediate scales $0.1\lesssim k\,{\rm Mpc}\,h^{-1}\lesssim 1$ are primarily controlled by the halo mass function ratio $n_{\rm vir}^{\rm real}/n_{\rm vir}^{\rm pseudo}$.

Fixing the real and pseudo linear power spectra to be identical (as in Eq. 51) forces the corresponding mass functions to be somewhat similar. Therefore, owing to (iii), reaction functions can overcome the typical inaccuracies that plague the halo model in the transition region between large and small scales.

To assign a value to the boost/suppression term $\mathcal{E}$, it is important to realise that we would like to preserve the smoothness in the transition from the linear to the nonlinear regime. In turn, this is tied to the shape of the one-halo terms on scales $0.5\lesssim k\,{\rm Mpc}\,h^{-1}\lesssim 1$, where the two-halo contribution becomes subdominant. In this regime the one-halo terms are well approximated by their large-scale limit $P_{1\rm h}(k\rightarrow 0)$, thus suggesting that

[TABLE]

is a good choice, one that only depends on the ratio $n_{\rm vir}^{\rm real}/n_{\rm vir}^{\rm pseudo}$.

In Eq. (53), the transition rate from linear to nonlinear scales is governed by the parameter $k_{\star}~{}$111111As we shall see below, $k_{\star}$ is derived from perturbation theory and is largely independent of the one-halo contribution. However, its specific value should be interpreted with caution, in that equally good alternatives to the exponential function in Eq. (53) can provide different solutions to Eq. (55).: for $k_{\star}\rightarrow 0$ the halo model reaction collapses to the ratio of one-halo terms; in the opposite case, $k_{\star}\rightarrow\infty$, the parameter $\mathcal{E}$ loses any role, and the reaction reduces to its definition in Mead (2017). We determine this scale using standard perturbation theory (SPT) (Koyama et al., 2009; Brax & Valageas, 2013; Bose & Koyama, 2016) by solving the equation

[TABLE]

where

[TABLE]

and we set $k_{0}=0.06$ $h\,{\rm Mpc}^{-1}$, which we found to be the largest wavenumber that can both ensure reliable perturbative predictions and keep the inaccuracies induced by the exponential sensitivity to $k_{\star}$ under control (see App. B and Carlson et al., 2009). Expressions for the second order corrections $P_{22}$, $P_{13}$ and $P_{13}^{\Psi}$ to the linear power spectrum are given in App. B. Note that alternative perturbation schemes can also be used in Eq. (55), such as the Lagrangian Perturbation Theory for modified gravity recently developed in Aviles & Cervantes-Cota (2017).

The role of the two-halo correction factor in Eq. (53) becomes clear in the limiting case $k_{\star}\rightarrow\infty$, where it goes to unity. Mead (2017) showed that this form of the reactions matches smooth dark energy simulations at percent level or better on scales $k\lesssim 1$ $h\,{\rm Mpc}^{-1}$for $z=0$ (see also Sec. 5 below). On quasi-linear scales this remarkable agreement can be understood in terms of standard perturbation theory. Schematically, the nonlinear matter power spectrum is a non-trivial function of the linear power spectrum obtained through some operator $\mathcal{K}$, i.e. $P(k)=\mathcal{K}[P_{\rm\scriptscriptstyle L}(k)]$ (Bernardeau et al., 2002). Provided that gravitational forces remain unchanged, then Eq. (51) enforces $P^{\rm real}\approx P^{\rm pseudo}$ on linear and quasi-linear scales. This is no longer true for modifications of gravity, since the structure of the $\mathcal{K}_{\rm MG}$ operator is altered by different mode-couplings and screening mechanisms. We correct for this fact by including the two-halo pre-factor in Eq. (53), so that finite $k_{\star}$ roughly encapsulates the extent of the mismatch between $\mathcal{K}_{\rm MG}$ and $\mathcal{K}_{\rm GR}$. For comparison, in F5 at $z=0$ we have $k_{\star}=0.3$ $h\,{\rm Mpc}^{-1}$, whereas in nDGPm it becomes $k_{\star}=0.95$ $h\,{\rm Mpc}^{-1}$, which reflects the different screening efficiency on large scales between the chameleon and Vainshtein mechanisms.

Although our choice of $k_{0}$ is commonly regarded as well within the quasi-linear regime, screening mechanisms in modified gravity induce nonlinearities on large scales that can be more important than in GR. For this, the determination of $k_{\star}$ can be complicated by inaccuracies specific to the perturbation theory employed, and to reduce their impact on the halo model reactions we take advantage of the following two facts: (i) on large scales we expect $P_{\rm NoScr}^{\rm real}\approx P^{\rm pseudo}$, where $P_{\rm NoScr}^{\rm real}$ denotes the nonlinear matter power spectrum of the real cosmology assuming there is no screening mechanism; (ii) the kernel operators $\mathcal{K}_{\rm MG}^{\rm Scr}$ and $\mathcal{K}_{\rm MG}^{\rm NoScr}$ for the screened and unscreened real cosmology, respectively, have a similar structure (see App. B). Therefore, at least in principle, the ratio $P_{\rm Scr,SPT}^{\rm real}/P_{\rm NoScr,SPT}^{\rm real}$ could give a better description of the reaction on large scales than the obvious candidate $P_{\rm Scr,SPT}^{\rm real}/P_{\rm SPT}^{\rm pseudo}$. Hereafter we will use $P_{\rm NoScr,SPT}^{\rm real}$ instead of $P_{\rm SPT}^{\rm pseudo}$ in Eq. (55), which in spite of being a sub-optimal strategy in some cases (see right panel of Figure 11) produces the most consistent behaviour across the cosmological models we have tested, as shown in Sec. 5.

In summary, halo model reactions provide a fast (we only need one-loop SPT for a single wavenumber) and general framework to map accurate nonlinear matter power spectra in $\Lambda$CDM to other non-standard cosmologies. We apply this method to $f(R)$ gravity, nDGP and evolving dark energy, and test its performance in Sec. 5 against the cosmological simulations described in the next section.

4 Simulations

The simulations of $f(R)$ gravity and DGP models used in this work were run using ecosmog (Li et al., 2012, 2013a; Li et al., 2013b), which has been developed to simulate the structure formation in various subclasses of models within the Horndeski family of theories. ecosmog is an extension of the simulation code ramses (Teyssier, 2002), which is a particle-mesh code employing adaptive mesh refinement to achieve high force resolution. The simulations are dark matter only and run in boxes with comoving size $512\,h^{-1}$Mpc using $1024^{3}$ simulation particles. Other basic information can be found in Table 2. The initial conditions of the simulations are generated using 2lptic (Crocce et al., 2006), which calculates the particle initial displacements and peculiar velocities up to second order in Lagrangian perturbations, allowing us to start from a relative low initial redshift $z_{\rm ini}=49$. To isolate the effect of nonlinearities we use identical phases for the initial density field in all cases. The linear power spectra used to generate the initial conditions are computed using camb (Lewis et al., 2000), with $\Omega_{\rm m}=0.3072$, $\Omega_{\Lambda}=0.6928$, $h=0.68$, $\Omega_{\rm b}=0.0481$ for all simulations. Importantly, following Eq. (51) the normalisation – and shape in $f(R)$ gravity – of the initial linear power spectra are different in the real and pseudo simulations. Since at early times deviations from GR are negligible, simulations in $\Lambda$CDM and modified gravity share the same initial conditions set by the $\Lambda$CDM power spectrum,

[TABLE]

with $\sigma_{8}(z=0)=0.8205$. The pseudo runs (to which we will apply the halo model reactions) have different initial conditions, generated using modified gravity linear power spectra at the final redshift, $z_{\rm f}$, and then rescaled with the $\Lambda$CDM linear growth to the starting redshift as

[TABLE]

where in this work $z_{\rm f}=0$ or 1. By evolving the initial real and pseudo power spectra, Eqs. (57) and (58), with the modified and standard laws of gravity, respectively, Eq. (51) will be automatically satisfied at $z_{\rm f}$. We extract the nonlinear matter power spectrum from our particle snapshots using the public code powmes (Colombi et al., 2009).

Simulations of dark energy models were run using a modified version of Gadget-2 that allows for the $\{w_{0},w_{a}\}$ parametrisation under the assumption that the dark energy is homogeneous. Initial conditions were generated at $z_{\rm ini}=199$ using N-GenIC (Springel, 2015), a code that calculates initial particle displacements and peculiar velocities based on the Zeldovich approximation. Our simulations take place in 200 ${\rm Mpc}\,h^{-1}$ boxes and use $512^{3}$ particles. Note that since we are concerned only with ratios of power spectra the overall resolution requirements on the simulations are less stringent than if we were interested in the absolute power spectra. We checked that our simulated reactions were insensitive to the realisation, box size, particle number and softening up to the wave numbers we show.

Differently from the modified gravity runs, we fix $\sigma_{8}=0.8$ for all the evolving dark energy models. Then, for the real cosmologies the amplitudes of the initial density field are determined by

[TABLE]

where $D_{\rm DE}(z)$ is the linear growth of a specific dark energy model, while for the pseudo counterparts one simply replaces $P_{\rm\rm\scriptscriptstyle L}^{\rm MG}$ with $P_{\rm\rm\scriptscriptstyle L}^{\rm DE}$ in Eq. (58).

5 Results

Here we test our halo model reactions (see Sec. 3.3) against the same quantities constructed from the real and pseudo cosmological simulations described in Sec. 4. To help get a better sense of the performance of our method, for each real cosmology $\mathcal{C}$ we also compute the standard ratios $P_{\mathcal{C}}^{\rm real}(k)/P_{\Lambda{\rm CDM}}(k)$, where our theoretical prediction for the real nonlinear power spectrum is obtained as

[TABLE]

where $\mathcal{R}(k)$ is given by Eq. (53). To test our modified gravity predictions we calculate $P_{\mathcal{C}}^{\rm pseudo}(k)$ from HMcode (Mead et al., 2015b, 2016) or using the measurement from the simulations directly. For evolving dark energy, instead, we use the pseudo nonlinear power spectrum given by the Coyote Universe emulator (Heitmann et al., 2014), and in so doing we illustrate how one could predict the real power spectrum for cosmologies beyond the concordance model with the aid of a carefully designed $\Lambda$CDM-like emulator (Giblin et al. in prep.).

5.1 $f(R)$ gravity

Figure 1 shows the matter power spectrum reactions calculated using Eq. (53) for the F5 cosmology at $z=0$ (left panel) and $z=1$ (right panel) in comparison to the measured reactions in the $N$-body simulations. For the virial halo mass function entering $P_{1\rm h}$ (see Eq. 49) we first adopt the approach developed in Lombriser et al. (2013), which incorporates both self- and environmental-screening in the spherical collapse (see App. A). At $z=0$ the halo concentrations, virial radii and mass functions are good enough to give percent level predictions on scales $k\lesssim 1$ $h\,{\rm Mpc}^{-1}$. A deviation of a few percent is however visible at $z=1$ starting on scales as large as $k\approx 0.2$ $h\,{\rm Mpc}^{-1}$. In App. C we show that changes in the halo profiles only affect scales $k\gtrsim 0.5$ $h\,{\rm Mpc}^{-1}$, suggesting that the observed inaccuracies could be caused by a mismatch between the predicted virial mass function ratio $n_{\rm vir}^{\rm real}/n_{\rm vir}^{\rm pseudo}$ and the same quantity measured in simulations. Indeed, Cataneo et al. (2016) found that, for halo masses defined by spherical regions with an average matter density 300 times the mean matter density of the Universe, the halo mass function of Lombriser et al. (2013) can deviate up to 10% from the simulations.

Given the complexity of measuring the virial halo mass function in $f(R)$ simulations121212Due to the nature of the chameleon screening, in $f(R)$ gravity the virial overdensity depends on both the mass of the halo and the gravitational potential in its environment. Things are much simpler in DGP, where by virtue of the Vainshtein screening both dependencies disappear., we investigate changes in the reactions induced by a more accurate description of the halo abundances with the fits provided in Cataneo et al. (2016). There, however, the calibration of the $n_{\Delta}^{f(R)}/n_{\Delta}^{\rm\Lambda CDM}$ ratios was performed for $\Delta=300$, whereas for our purposes we need $\Delta=\Delta_{\rm vir}$. We go from one mass definition to the other with the scaling relations outlined in Hu & Kravtsov (2003) (for a first application to $f(R)$ gravity see Schmidt et al., 2009). Inevitably, this transformation suffers from inaccuracies in $c_{\rm vir}$ and $\Delta_{\rm vir}$, which we attempt to compensate for by adjusting the $M_{300}(M_{\rm vir})$ relation so that the new rescaled mass function provides a present-day halo model reaction that is at least as good as the reaction obtained when using the Lombriser et al. (2013) virial mass function (dotted line in the left panel of Figure 1)131313In practice, we start with the Hu & Kravtsov (2003) relation $M_{300}=\mathcal{Q}(M_{\rm vir})M_{\rm vir}$, and make the replacement $\mathcal{Q}(M_{\rm vir})\rightarrow\mathcal{Q}^{\prime}(M_{\rm vir})=\min[{\tilde{a}}\mathcal{Q}^{f(R)}({\tilde{b}}M_{\rm vir}),\mathcal{Q}^{\rm GR}(M_{\rm vir})]$, where ${\tilde{a}}$ and ${\tilde{b}}$ are $\mathcal{O}(1)$ free parameters fine-tuned to reach the required accuracy in the halo model reaction at $z=0$. The use of the minimum operator ensures that the $f(R)$ conversion factor matches the corresponding GR value for masses large enough to fully activate the chameleon screening.. We then take the ratio of the Cataneo et al. (2016) rescaled mass function to that of Lombriser et al. (2013), and treat this quantity as a correction factor for the latter. To find the required adjustment at high redshifts we shift the $z=0$ correction by an amount $\Delta\log_{10}M_{\rm vir}$ inferred from the redshift evolution of the ratio of the two halo mass functions over the range $z\in[0,0.5]$ (see central panel of Figure 4 in Cataneo et al., 2016). A simple extrapolation to $z=1$ gives $\Delta\log_{10}M_{\rm vir}=1$. Although far from being a rigorous transformation, the resulting halo model reaction now agrees to better than 1% down to $k\approx 3$ $h\,{\rm Mpc}^{-1}$, as shown in the right panel of Figure 1. Figure 2 illustrates that similar considerations are also valid for the F6 cosmology, where we use the same mass shift $\Delta\log_{10}M_{\rm vir}=1$ to go from the $z=0$ to the $z=1$ mass function correction. In all cases, deviations in the highly nonlinear regime are most likely caused by inaccurate $c$-$M$ relations. We leave the study of $f(R)$ gravity reactions on small scales derived from proper virial quantities for future work.

Figure 3 and Figure 4 show the relative change of the matter power spectrum in $f(R)$ gravity with respect to GR for the F5 and F6 models, respectively. The left panels present the best-case scenario, that is, when “perfect” knowledge of the pseudo power spectrum is available. In this case, since the uncertainties come entirely from our halo model predictions, we can obviously reproduce the power spectrum ratios at the same level of accuracy of our reactions. For now the pseudo information comes directly from our simulations, but it is not hard to imagine a specifically designed emulator capable of generating the nonlinear matter power spectrum of $\Lambda$CDM cosmologies with non-standard initial conditions. We will analyse the requirements for such emulator in a future work (Giblin et al. in prep.). In the right panels we compute $P^{\rm pseudo}$ with HMcode to demonstrate that currently, together with publicly available codes, our method can achieve 2% accuracy on scales $k\lesssim 1$ $h\,{\rm Mpc}^{-1}$in modified gravity theories characterised by scale-dependent linear growth.

For a comparison to a range of other methods for modelling the nonlinear matter power spectrum in $f(R)$ and other chameleon gravity models, we refer to Figures 4 and 5 in Lombriser (2014), noting that the majority of these methods rely on fitting parameters in contrast to the approach discussed here.

5.2 DGP

Spherical collapse dynamics is much simpler in nDGP (Schmidt et al., 2010), with both the linear overdensity threshold for collapse, $\delta_{\rm c}$, and the corresponding average virial halo overdensity, $\Delta_{\rm vir}$, being only functions of redshift. For instance, in GR one has $\Delta_{\rm vir}(z=0)=335$ and $\Delta_{\rm vir}(z=1)=200$ for $\Omega_{\rm m}=0.3072$, whereas these values become $\Delta_{\rm vir}(z=0)=283$ and $\Delta_{\rm vir}(z=1)=178$ in our nDGPm cosmology. This fact allows us to extract the virial halo mass function directly from our simulations, and test that accurate $n_{\rm vir}^{\rm real}/n_{\rm vir}^{\rm pseudo}$ ratios do indeed produce accurate halo model reactions. Figures 5 and 6 show that, after refitting the virial Sheth-Tormen mass function to the same quantity from simulations, halo model predictions reach percent level accuracy on scales $k\lesssim 1$ $h\,{\rm Mpc}^{-1}$(see App. C for details on the halo mass function calibration to simulations). Moreover, since the Vainshtein radius (Eq. 23) for the most massive halos is of order a few megaparsecs, we expect small changes caused by the screening mechanism on large scales, i.e. $k\lesssim 0.1$ $h\,{\rm Mpc}^{-1}$. In other words, although in our calculations we keep the two-halo correction factor in Eq. (53), it contributes only marginally to improving the performance of our reaction functions. This is evident from the perturbation theory predictions shown in Figure 12. Once again, deviations on scales $k\gtrsim 1$ $h\,{\rm Mpc}^{-1}$could be primarily sourced by inaccurate real and pseudo halo concentrations, and is the subject of future investigation.

We study the ability of the halo model reactions to reproduce the relative difference of the nDGP power spectrum from that of standard gravity when combined with the pseudo matter power spectrum from either HMcode or the simulations. Figures 7 and 8 confirm that with current codes also scale-independent modifications of gravity can be predicted within 2% over the range of scales relevant for this work.

5.3 Dark energy

Figure 9 shows the reaction functions for the evolving dark energy cosmologies listed in Table 1. The left panel contains essentially the same $z=0$ information of Figure 2 in Mead (2017), with the notable difference that here we compute the spherical collapse virial overdensities including the dark energy contribution to the potential energy, and do not assume energy conservation during collapse (Schmidt et al., 2010) (expressions for the individual terms entering the virial theorem can be found in App. A). The right panel shows the same quantity at $z=1$. At both redshifts the halo model reactions based on the standard Sheth-Tormen mass function fits can capture very well the measurements from simulations down to the transition scale between the two- and one-halo terms. Also in this case, we attribute inaccuracies on small scales mainly to the inadequacy of the Dolag et al. (2004) and Bullock et al. (2001) halo concentration prescriptions for the real and pseudo cosmologies, respectively.

In Figure 10 we consider two representative dark energy models, DE2 and DE3, and compare their matter power spectra to that of $\Lambda$CDM with the same initial conditions. Their particular equations of state enhance the growth of structure in one case and suppress it in the other. Here, we employ the Coyote Universe Emulator (Heitmann et al., 2014) not only for our baseline pseudo power spectra, but also as a substitute for the real and $\Lambda$CDM cosmology simulations. This serves as an example to illustrate a straightforward application of the reaction functions: extend the cosmological parameter space of matter power spectrum emulators designed for the concordance cosmology only, without the need to run model-dependent and expensive cosmological simulations. For the evolving equation of state of DE3 we use the emulator extension code PKequal built upon the work presented in Casarini et al. (2016). Knowing that the output from the emulator is 1-2% accurate on scales $k\lesssim 1$ $h\,{\rm Mpc}^{-1}$, from the previous results in Figure 9 we can expect similar agreement between our reaction-based power spectra and those obtained from the emulator itself. This is indeed the case except in the range $0.02\lesssim k\,{\rm Mpc}\,h^{-1}\lesssim 0.5$, where the interpolation process within the emulator fails to capture the correct dependence on $\omega_{\rm b}$ and $w$ because the specific values we use sit on the edge of its domain of applicability141414The Coyote Universe emulator accepts values $0.0215<\omega_{\rm b}<0.0235$ and $-1.3<w<-0.7$, while for our two evolving dark energy models we have $\omega_{\rm b}=0.0245$, $w^{\rm(DE2)}=-1.3$ and $w_{\rm eff}^{\rm(DE3)}=-0.84$. Since our background baryon density resides outside the domain of the emulator, we set it to the maximum value allowed. This has virtually no impact on our halo model responses, in that they depend only weakly on $\omega_{\rm b}$..

6 Conclusions

The spatial distribution of matter in the Universe and its evolution with time emerge from the interplay of gravitational and astrophysical processes, and are inextricably linked to the nature of the cosmic matter-energy constituents. The power spectrum is an essential statistic describing the clustering of matter in the Universe, and lies at the heart of probes of the growth of structure such as cosmic shear and galaxy clustering. Measurements of these quantities from the next generation of large-volume surveys are expected to reach percent level uncertainty – upon careful control of systematics – on scales where nonlinearities and baryonic physics become important. It is notoriously difficult to predict the matter power spectrum in this regime to such a degree of accuracy, yet these scales contain a wealth of information on currently unanswered questions, e.g. the nature of dark energy, the sum of neutrino masses and the extent of baryonic feedback mechanisms.

In this work we focused on modelling the nonlinear matter power spectrum in modified theories of gravity and evolving dark energy cosmologies. We extended the reaction method of Mead (2017) using the halo model to predict the nonlinear effects induced by new physics on the matter power spectrum of specifically designed reference cosmologies. These fiducial – pseudo – cosmologies mimic the linear clustering of the target – real – cosmologies, yet their evolution is governed by standard gravity with $\Lambda$CDM expansion histories (which are either quick to simulate with current resources, if not already available with emulators). We showed that by applying the halo model reactions to the nonlinear matter power spectrum of the pseudo cosmologies we are able to recover the real counterpart to within 1% on scales $k\lesssim 1$ $h\,{\rm Mpc}^{-1}$for all cases under study. Remarkably, our methodology does not involve fitting the power spectra measured in simulations at any stage. Instead, having access to accurate ratios of the halo mass function in the real cosmology to that in the pseudo cosmology is crucial to achieve the observed performance. Not including this information from the simulation degrades the accuracy to $\lesssim 3\%$. The halo model reactions can also be used to predict the matter power spectrum in the highly nonlinear regime. However, this requires additional knowledge of the average structural properties of the dark matter halos as well as the inclusion of baryonic effects (see, e.g, Schneider et al., 2018). We leave these improvements for future work (Cataneo et al. in prep.).

In the case of the dark energy models we adopted the Coyote Universe emulator for the pseudo matter power spectrum (i.e. for $w=-1$), which we then combined with our halo model reaction to obtain the real expected quantity. By comparing this prediction to the real output of the emulator (i.e. for $w\neq-1$) we showed that emulators trained on pure $\Lambda$CDM models can be accurately extended to non-standard cosmologies in an analytical way, thus substantially increasing their flexibility while simultaneously reducing the computational cost for their design. However, applications of this strategy to scale-dependent modifications of gravity (such as $f(R)$ models) necessitate of a more elaborate $\Lambda$CDM emulator that takes as input also information on the linearly modified shape of the matter power spectrum (Giblin et al. in prep.). Together with suitable halo model reactions, this emulator can also be employed to predict the nonlinear total matter power spectrum in massive neutrino cosmologies (Cataneo et al. in prep.), where the presence of a free streaming scale induces a scale-dependent linear growth (Lesgourgues & Pastor, 2006).

Given that our method builds on the halo model, the halo mass function and the spherical collapse model are absolutely central for our predictions. Contrary to the standard halo model calculations, however, the accuracy of our reactions strongly depends on the precision of the pseudo and real halo mass functions. This opens up the possibility of combining in a novel way cosmic shear and cluster abundance measurements, for example. Moreover, quite general modifications of gravity – with their screening mechanisms – can be implemented in the spherical collapse calculations through the nonlinear parametrised post-Friedmannian formalism of Lombriser (2016).

In summary, halo model reactions provide a fast, accurate and versatile method to compute the real-space nonlinear matter power spectrum in non-standard cosmologies. Successful implementations in redshift-space by Mead (2017) pave the way for applications to redshift-space distortions data as well. Altogether, these features make the halo model reactions an attractive alternative in-between perturbative analytical methods and brute force emulation, and promise to be an essential tool in future combined-probe analyses in search of new physics beyond the standard paradigm.

Acknowledgements

We thank Philippe Brax, Patrick Valageas, Alejandro Aviles and Jorge Cervantes-Cota for crosschecking our perturbation theory calculations, and David Rapetti and Simon Foreman for useful discussions. MC and CH acknowledge support from the European Research Council under grant number 647112. LL acknowledges support by a Swiss National Science Foundation (SNSF) Professorship grant (No. 170547), a SNSF Advanced Postdoc.Mobility Fellowship (No. 161058), the Science and Technology Facilities Council Consolidated Grant for Astronomy and Astrophysics at the University of Edinburgh, and the Affiliate programme of the Higgs Centre for Theoretical Physics. AM has received funding from the European Union‘s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 702971 and acknowledges MDM-2014-0369 of ICCUB (Unidad de Excelencia Maria de Maeztu). SB is supported by Harvard University through the ITC Fellowship. BL is supported by the European Research Council via grant ERC-StG-716532-PUNCA, and by UK STFC Consolidated Grants ST/P000541/1 and ST/L00075X/1. This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1, ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.

Appendix A Spherical collapse in modified gravity and quintessence

We shall briefly review expressions for the force enhancement $\mathcal{F}$ in $f(R)$ and DGP gravity used in the modified spherical collapse calculation in Sec. 3.1 as well as its impact and the impact of dark energy domination on the virial theorem.

The force enhancement $\mathcal{F}$ adopted here for the spherical collapse calculation in $f(R)$ gravity is given by (Lombriser et al., 2013)

[TABLE]

which uses the thin-shell approximation (Khoury & Weltman, 2004; Li & Efstathiou, 2012) with thickness $\Delta R\ll R_{\rm\scriptscriptstyle TH}$. The expression is also adopted for the thick-shell limit and interpolates between the small-field ($\mathcal{F}=0$) and large-field ($\mathcal{F}=1/3$) regimes, which correspond to the two limiting scenarios studied in the $f(R)$ spherical collapse calculation of Schmidt et al. (2009). Furthermore, for the $f(R)$ functional form Eq. (13) one finds (Lombriser et al., 2013)

[TABLE]

where the normalised top-hat radius

[TABLE]

needs to be solved in both the halo (h) and the environment (env) using Eq. (32), that now becomes

[TABLE]

where we set $\mathcal{F}=0$ for the $\Lambda$CDM environment and primes denote derivatives with respect to $\ln a$. To solve Eq. (64) we use the initial conditions $y=1$ and $y^{\prime}=-\delta_{\rm i}/3$ obtained from the linear theory in an Einstein-de Sitter Universe, which at early times is an accurate description for all models in Sec. 2 when the contribution from radiation is ignored. We then iteratively adjust $\delta_{\rm i}$ until the condition $y(a_{\rm coll})=0$ (i.e. $R_{\rm\scriptscriptstyle TH}=0$) is satisfied at the desired time of collapse, $a_{\rm coll}$, within some small tolerance. For the environment, instead, we follow Cataneo et al. (2016). Also note that Eq. (62) can easily be generalised for the family of chameleon models (Lombriser et al., 2014).

In DGP gravity, the force modification becomes (Schmidt et al., 2010)

[TABLE]

where $x\equiv R_{\rm\scriptscriptstyle TH}/R_{\rm V}$, with the Vainshtein radius and $\beta$ function given in Eqs. (23) and (20), respectively. In smooth quintessence cosmologies $\mathcal{F}=0$, and the sole effect on the spherical collapse dynamics enters through the non-standard background expansion. In both DGP and quintessence there is no contribution from the environment, hence $y_{\rm h}$ alone describes the full evolution. Moreover, $\mathcal{F}$ can be generalised and parametrised to cover the range of different screening mechanisms (Lombriser, 2016).

The virial theorem for a general metric theory of gravity remains unchanged with respect to its formulation in GR, and reads

[TABLE]

where $W$ is the potential energy of the system and $T$ its kinetic energy. However, energy is not conserved for evolving dark energy and modified gravity scenarios, and the virial radius cannot be related to the turnaroud radius in the usual way (Lahav et al., 1991). Instead, one must find the time of virialisation, $a_{\rm vir}$, that satisfies Eq. (66) when all contributions to $W$ are considered. More specifically, the Newtonian, scalar field and background potential energies take the form (Schmidt et al., 2010)

[TABLE]

where

[TABLE]

and the kinetic energy can be written as

[TABLE]

Combining eqs. (66)-(71) together with Eqs. (36)-(37) allows us to find the virial radius $R_{\rm vir}$ as a function of halo mass, and for various theories of gravity as well as expansion histories.

Appendix B Perturbation theory

The nonlinear evolution of matter perturbations in modified gravity has been extensively studied in Koyama et al. (2009), Brax & Valageas (2012), Brax & Valageas (2013) and Bose & Koyama (2016). Here we summarise their results and provide explicit expressions for the computation of next-to-leading-order corrections to the linear power spectrum.

The continuity and Euler equations describing the evolution of the matter density perturbations, $\delta$, and peculiar velocity, ${\bf v}$, are

[TABLE]

where the gravitational potential $\Psi$ in modified gravity theories with screening mechanisms depends nonlinearly on the matter density perturbations, as in Eqs. (9) and (18). Assuming vanishing vorticity, the velocity field can be fully described by its divergence $\theta\equiv\nabla\cdot{\bf v}/(aH)$. By defining the vector field

[TABLE]

and expanding the gravitational potential up to third order in the perturbations, the fluid equations in Fourier space read

[TABLE]

where primes denote derivatives with respect to $\eta\equiv\ln a$, ${\bf k}_{1\cdots n}\equiv{\bf k}_{1}+\cdots+{\bf k}_{n}$, and repeated indices are summed over. The left hand side of Eq. (75) controls the evolution of linear perturbations, with the matrix

[TABLE]

The quantity $\epsilon$ measures the time- and scale-dependent linear departure from GR, which in $f(R)$ gravity and DGP reads

[TABLE]

where the mass

[TABLE]

with

[TABLE]

and the DGP function $\beta$ is given in Eq. (20).

The sources of nonlinearities are confined to the right hand side of Eq. (75), where the second order vertices are

[TABLE]

[TABLE]

with

[TABLE]

and the only modified gravity non-vanishing third order vertex is $\sigma_{2;111}$, with

[TABLE]

[TABLE]

The vertex $\gamma_{2;11}=0$ in standard gravity, whereas

[TABLE]

[TABLE]

The one-loop power spectrum corrections require the density and velocity fields up to third order in the perturbations, that is, we need $\varpi=\varpi^{(1)}+\varpi^{(2)}+\varpi^{(3)}$. For this, we shall first find the linear solution to Eq. (75), $\varpi^{(1)}$, and then derive recursively the higher order solutions with the help of the retarded Green function, which is obtained by setting the right hand side of Eq. (75) to a Dirac delta function (see Brax & Valageas, 2012, 2013, for more details). We can write the linear ansatz as

[TABLE]

where the growing mode, $D_{+}(k,\eta)$, satisfies the linear second-order equation

[TABLE]

with Einstein-de Sitter initial conditions $D_{+}(\eta_{\rm i})=D_{+}^{\prime}(\eta_{\rm i})=e^{\eta_{\rm i}}$ at the initial time $\eta_{\rm i}$. The decaying mode, $D_{-}$, also satisfies this equations and can be directly computed as

[TABLE]

where the Wronskian of $D_{+}$ and $D_{-}$ is given by

[TABLE]

By using the retarded Green function

[TABLE]

where $\Theta(\eta_{1}-\eta_{2})$ denotes the Heaviside function, we obtain the second and third order density perturbations as the convolutions

[TABLE]

[TABLE]

with the operators $\mathcal{K}^{(2)}$ and $\mathcal{K}^{(3)}$ built from the second and third order vertices, $\gamma_{m}$ and $\sigma_{m}$, respectively (see Brax & Valageas, 2013, for more details). In Eq. (93) the matrix $\mathcal{\tilde{G}}$ takes the explicit form

[TABLE]

where $D_{+/-i}$ is shorthand for $D_{+/-}(k_{i},\eta_{i})$.

Finally, we can write the equal-time two-point correlation – equivalent to Eq. (56) – as

[TABLE]

which produces the following expressions for the one-loop integrals

[TABLE]

[TABLE]

[TABLE]

where all indices run from 1 to 2, the linear power spectra are taken at $\eta=0$, and we have used the two-point linear correlator

[TABLE]

To deactivate the screening mechanisms we simply set $\gamma_{2;11}=\sigma_{2;111}=0$, while keeping the linear deviation $\epsilon$ given in Eqs. (77)-(78). This results in what we call the real ‘no-screening’ power spectrum, $P_{\rm NoScr}^{\rm real}$. The pseudo power spectrum up to one-loop, instead, follows from also imposing $\epsilon=0$, such that the only difference compared to the standard cosmology is in the initial conditions, i.e. in the shape and/or amplitude of the linear power spectrum. Figures 11 and 12 show the SPT reaction functions on quasi-linear scales for $f(R)$ gravity and nDGP, where the pseudo nonlinear power spectrum is computed using either the pseudo (solid lines) or ‘no-screening’ (dashed lines) one-loop correction. Although the former approach is theoretically motivated, accuracy arguments make the latter preferable for all cosmologies and redshifts investigated in this work.

Appendix C Halo mass function and concentration tests

In Sec. 3.3 we showed that the real to pseudo halo mass function ratio, $n_{\rm vir}^{\rm real}/n_{\rm vir}^{\rm pseudo}$, controls the halo model reaction on scales $k\gtrsim 0.1$ $h\,{\rm Mpc}^{-1}$. In particular, for the nDGP cosmologies we explicitly checked that fitting the semi-analytical halo mass functions directly to our simulations results in better performance of the halo model reactions for wavenumbers $0.1\lesssim k\,{\rm Mpc}\,h^{-1}\lesssim 1$. Below we briefly explain how these fits to the halo mass functions extracted from our simulations were carried out.

We measure the mean halo abundances and their uncertainties from simulations as outlined in Cataneo et al. (2016), with the important difference that, here, halos identified with the rockstar halo finder (Behroozi et al., 2013) have masses defined in spherical volumes with mean density $\Delta_{\rm vir}$ times the background comoving matter density, as in Eq. (36). The virial overdensity depends on redshift, matter content and theory of gravity through Eq. (37). Compared to Cataneo et al. (2016) we also use a finer mass binning, with a bin size $\Delta\log_{10}M=0.1$. Figures 13 and 14 present a rescaled version of the halo mass functions (i.e. the large-scale limit of the one-halo integrand Eq. 49) as well as the $n_{\rm vir}^{\rm real}/n_{\rm vir}^{\rm pseudo}$ ratios. The standard halo mass function fits for the parameters in Eq. (41) systematically underpredict the simulation measurements. Refitting these parameters to each real and pseudo simulation separately can change the halo mass function ratio $n_{\rm vir}^{\rm real}/n_{\rm vir}^{\rm pseudo}$ to an extent relevant for the halo model reactions. Note that when this ratio remains largely unaffected by the refitting, the predicted reactions also show very little variation (compare the right panels of Figure 14 and Figure 6). Table 3 summarises the refitted ST halo mass function parameters for the real and pseudo nDGP cosmologies used in this work.

Given the structure of the one-halo term Eq. (49), one might argue that inaccurate halo profiles are partially responsible for the few percent mismatch between the predicted and measured reactions in Figure 1 and Figure 5. If that was the case, then having accurate halo mass functions would not necessarily correspond to having accurate halo model reactions. However, Figure 15 shows otherwise, in that even extreme variations of the halo concentrations cause important changes only on scales $k\gtrsim 0.5$ $h\,{\rm Mpc}^{-1}$.

From these tests we conclude that accurate virial halo mass function ratios are central to determining $\lesssim 1\%$ accuracy for the halo model reactions, and that the mean concentration of massive halos is well described by Eq. (46), at least as long as this relation is employed in the reaction ratios. In future works we will investigate the implications of refitting the $c$-$M$ relation to that measured in simulations, which has the potential to ameliorate our halo model predictions in the highly nonlinear regime.