Emulators for the non-linear matter power spectrum beyond $\Lambda$CDM

TL;DR

This paper introduces a fast, accurate emulator for the non-linear matter power spectrum in Hu-Sawicki $f(R)$ gravity models, enabling efficient predictions for modified gravity effects in cosmological surveys.

Contribution

The authors develop a novel, computationally inexpensive method to create emulators for modified gravity models by leveraging existing $ ext{Lambda}$CDM emulators and simulations, applicable to a range of cosmological parameters.

Findings

The emulator accurately predicts the enhancement of the matter power spectrum for $k \,\lesssim\, 10\, h\text{Mpc}^{-1}$ and $z \lesssim 3$.

The cosmology dependence of the enhancement is weak, allowing use of a single cosmology simulation for broad predictions.

The method can be extended to other modified gravity models and used for Fisher forecasts in upcoming surveys.

Abstract

Accurate predictions for the non-linear matter power spectrum are needed to confront theory with observations in current and near future weak lensing and galaxy clustering surveys. We propose a computationally cheap method to create an emulator for modified gravity models by utilizing existing emulators for $\Lambda$CDM. Using a suite of $N$-body simulations we construct a fitting function for the enhancement of both the linear and non-linear matter power spectrum in the commonly studied Hu-Sawicki $f(R)$ gravity model valid for wave-numbers $k \lesssim 5-10\, h\text{Mpc}^{-1}$ and redshifts $z \lesssim 3$. We show that the cosmology dependence of this enhancement is relatively weak so that our fit, using simulations coming from only one cosmology, can be used to get accurate predictions for other cosmological parameters. We also show that the cosmology dependence can, if needed, be included by using linear theory, approximate $N$-body simulations (such as COLA) and semi-analytical tools like the halo model. Our final fit can easily be combined with any emulator or semi-analytical models for the non-linear $\Lambda$CDM power spectrum to accurately, and quickly, produce a non-linear power spectrum for this particular modified gravity model. The method we use can be applied to fairly cheaply construct an emulator for other modified gravity models. As an application of our fitting formula we use it to compute Fisher-forecasts for how well galaxy clustering and weak lensing in a Euclid-like survey will be at constraining modifications of gravity.

Tables

Table 1. Table 1: The best-fit parameters for the enhancement of the linear power spectrum for | f R ​ 0 fid | = 10 − 5 superscript subscript 𝑓 𝑅 0 fid superscript 10 5 |f_{R0}^{\rm fid}|=10^{-5} .

	$i=0$	$i=1$	$i=2$
$b_{0i}$	3.10	2.34466	-1.86362
$c_{0i}$	34.4951	28.8637	-13.1302
$d_{0i}$	0.14654	-0.0100	-0.14944
$e_{0i}$	1.62807	0.71291	-1.41003

Table 2. Table 2: The best fit values using the data from | f R ​ 0 | = { 10 − 5 , 5 ⋅ 10 − 6 , 10 − 6 , 10 − 7 } subscript 𝑓 𝑅 0 superscript 10 5 ⋅ 5 superscript 10 6 superscript 10 6 superscript 10 7 |f_{R0}|=\{10^{-5},5\cdot 10^{-6},10^{-6},10^{-7}\} with | f R ​ 0 | fid = 5 ⋅ 10 − 6 superscript subscript 𝑓 𝑅 0 fid ⋅ 5 superscript 10 6 |f_{R0}|^{\rm fid}=5\cdot 10^{-6} (left), from | f R ​ 0 | = { 5 ⋅ 10 − 5 , 10 − 5 , 5 ⋅ 10 − 6 } subscript 𝑓 𝑅 0 ⋅ 5 superscript 10 5 superscript 10 5 ⋅ 5 superscript 10 6 |f_{R0}|=\{5\cdot 10^{-5},10^{-5},5\cdot 10^{-6}\} with | f R ​ 0 | fid = 10 − 5 superscript subscript 𝑓 𝑅 0 fid superscript 10 5 |f_{R0}|^{\rm fid}=10^{-5} (middle) and from | f R ​ 0 | = { 10 − 4 , 5 ⋅ 10 − 5 , 10 − 5 } subscript 𝑓 𝑅 0 superscript 10 4 ⋅ 5 superscript 10 5 superscript 10 5 |f_{R0}|=\{10^{-4},5\cdot 10^{-5},10^{-5}\} with | f R ​ 0 | fid = 5 ⋅ 10 − 5 superscript subscript 𝑓 𝑅 0 fid ⋅ 5 superscript 10 5 |f_{R0}|^{\rm fid}=5\cdot 10^{-5} (right).

	$i=0$	$i=1$	$i=2$
$b_{0i}$	0.76878	0.22638	0.00759
$b_{1i}$	-0.40537	-0.10711	-0.00102
$b_{2i}$	0.00752	0.04846	0.01180
$c_{0i}$	0.02886	-0.02438	-0.02963
$c_{1i}$	-0.06382	-0.05196	0.02597
$c_{2i}$	-0.40121	-0.03518	0.07688
$d_{0i}$	1.00000	0.10901	0.12027
$d_{1i}$	0.00000	0.08189	0.02492
$d_{2i}$	0.00000	-0.05682	-0.02985
$e_{0i}$	0.36951	0.14719	0.03127
$e_{1i}$	0.10939	0.06176	0.02933
$e_{2i}$	-0.34209	-0.13138	-0.01419
$f_{0i}$	1.03544	-0.13912	-0.05656
$f_{1i}$	-0.26277	-0.13231	-0.03389
$f_{2i}$	0.23028	0.13132	0.07127
$g_{0i}$	0.20246	0.06323	0.05229
$g_{1i}$	-0.11611	0.06943	0.07807
$g_{2i}$	0.10245	0.00296	-0.09770

Table 3. Table 3: Constraints on the parameters in the Hu-Sawicki model, for the fiducial value | f R ​ 0 | = 10 − 5 subscript 𝑓 𝑅 0 superscript 10 5 |f_{R0}|=10^{-5} , coming from galaxy clustering (above), weak-lensing (middle) and combined (below) for two different values of k max subscript 𝑘 max k_{\rm max} and ℓ max subscript ℓ max \ell_{\rm max} .

GC $k_{\max }=0.15h{\mathrm{Mpc}}^{-1}$
${\sigma }_{{\Omega }_b}=4.55\%$	${\sigma }_h=0.45\%$	${\sigma }_{{\Omega }_m}=2.78\%$	${\sigma }_{n_s}=1.04\%$	${\sigma }_{{10}^9{A}_s}=23.81\%$	${\sigma }_{\log \left|f_{R0}\right|}=8.00\%$
GC $k_{\max }=0.25h{\mathrm{Mpc}}^{-1}$
${\sigma }_{{\Omega }_b}=0.23\%$	${\sigma }_h=0.15\%$	${\sigma }_{{\Omega }_m}=1.39\%$	${\sigma }_{n_s}=0.83\%$	${\sigma }_{{10}^9{A}_s}=9.52\%$	${\sigma }_{\log \left|f_{R0}\right|}=4.00\%$
WL ${\ell }_{\max }=1000$
${\sigma }_{{\Omega }_b}=9.09\%$	${\sigma }_h=2.99\%$	${\sigma }_{{\Omega }_m}=5.56\%$	${\sigma }_{n_s}=4.17\%$	${\sigma }_{{10}^9{A}_s}=9.52\%$	${\sigma }_{\log \left|f_{R0}\right|}=18.00\%$
WL ${\ell }_{\max }=3000$
${\sigma }_{{\Omega }_b}=4.55\%$	${\sigma }_h=1.49\%$	${\sigma }_{{\Omega }_m}=3.47\%$	${\sigma }_{n_s}=3.12\%$	${\sigma }_{{10}^9{A}_s}=3.76\%$	${\sigma }_{\log \left|f_{R0}\right|}=16.00\%$
GC+WL $k_{\max }=0.15h{\mathrm{Mpc}}^{-1}$, ${\ell }_{\max }=1000$
${\sigma }_{{\Omega }_b}=3.12\%$	${\sigma }_h=0.3\%$	${\sigma }_{{\Omega }_m}=1.39\%$	${\sigma }_{n_s}=0.62\%$	${\sigma }_{{10}^9{A}_s}=1.88\%$	${\sigma }_{\log \left|f_{R0}\right|}=2.00\%$
GC+WL $k_{\max }=0.25h{\mathrm{Mpc}}^{-1}$, ${\ell }_{\max }=3000$
${\sigma }_{{\Omega }_b}=1.79\%$	${\sigma }_h=0.15\%$	${\sigma }_{{\Omega }_m}=0.69\%$	${\sigma }_{n_s}=0.31\%$	${\sigma }_{{10}^9{A}_s}=0.94\%$	${\sigma }_{\log \left|f_{R0}\right|}=0.80\%$

Equations

$$B_{{\Omega }_m\mathrm{corr}}\simeq 1-a\frac{\Delta {\Omega }_m}{{\Omega }_m}\tanh {\left(\frac{k}{k_{\ast }}\right)}^b$$

$$B_{{\sigma }_8\mathrm{corr}}\simeq 1+\frac{\Delta {\sigma }_8}{{\sigma }_8}\frac{k}{(1+(k/k_{\ast }){)}^2}$$

$$\ddot{\delta }+2H\dot{\delta }=\frac{3}{2}{\Omega }_m(a)H^2\delta \left(1+\frac{1}{3}\frac{k^2}{k^2+m^2{a}^2}\right)$$

$$m^2(a)=\frac{H_0^2({\Omega }_m+4{\Omega }_{\Lambda })}{2|f_{R0}|}{\left(\frac{{\Omega }_m{a}^{-3}+4{\Omega }_{\Lambda }}{{\Omega }_m+4{\Omega }_{\Lambda }}\right)}^3$$

$$\frac{P^{\mathrm{linear}\mathrm{fid}}(k,z)}{P_{\Lambda }^{\mathrm{linear}}\mathrm{CDM}(k,z)}$$

$$+d(a)\left|\frac{\log (k)k}{k-1}\right|\arctan (e(a)k)$$

$$X(z,f_{R0})=X_0(r)+X_1(r)(a-1)+X_2(r)(a-1{)}^2$$

$$\frac{P_{f(R)}(k,z)}{P_{\Lambda \mathrm{CDM}}(k,z)}=1+b(z,f_{R0})\frac{(1+c(z,f_{R0})\cdot k)}{(1+d(z,f_{R0})\cdot k)}\times$$

$$\times \arctan {\left(e(z,f_{R0})\cdot k)\right)}^{f(z,f_{R0})+g(z,f_{R0})\cdot k}$$

$$X(z,f_{R0})=X_0(r)+X_1(r)(a-1)+X_2(r)(a-1{)}^2$$

$$X_i(r)=X_{i0}+X_{i1}r+X_{i2}{r}^2$$

$$\log \mathcal{L}=-\frac{1}{2}\sum _{k,z}(P(k,z)-P^{\mathrm{fid}}(k,z))C^{-1}(P(k,z)-P^{\mathrm{fid}}(k,z))$$

$$P_g(k,\mu ,z)=\frac{D_{A,f}^2(z)/H_f(z)}{D_A^2(z)/H(z)}(b(z)+f(z){\mu }^2{)}^2\times$$

$$\times e^{-k^2{\mu }^2{\sigma }^2}{P}_{\mathrm{DM}}(k,z)$$

$$F_{ij}=\frac{V_{\mathrm{survey}}}{8{\pi }^2}{\int }_{-1}^1\mathrm{d}\mu {\int }_{k_{\min }}^{k_{\max }}\mathrm{d}k\left[\frac{\partial D}{\partial {\theta }_i}{D}^{-1}\frac{\partial D}{\partial {\theta }_j}{D}^{-1}\right]$$

$$C_{ij}(\ell )=\frac{9}{4}{\int }_0^{\infty }\mathrm{d}z\frac{W_i(z)W_j(z)H^3(z){\Omega }_m^2(z)}{(1+z{)}^4}\times$$

$$\times P_{\mathrm{DM}}(k=\ell /r(z),z)$$

$$F_{ij}=f_{\mathrm{sky}}\sum _{\ell }^{{\ell }_{\max }}\frac{(2\ell +1)\Delta \ell }{2}\text{tr}\left[\frac{\partial C}{\partial {\theta }_i}{\text{Cov}}^{-1}\frac{\partial C}{\partial {\theta }_j}{\text{Cov}}^{-1}\right]$$

Full text

I Introduction

One of the objectives of many next generation surveys such as Euclid Scaramella et al. (2014) and LSST Ivezić et al. (2008) is to look for and constrain any deviations from the predictions of general relativity (GR). Modifications of gravity have been studied quite extensively over the last decade (see e.g. Koyama (2016); Clifton et al. (2012) and references within). Such modifications, when they are in agreement with local and astrophysical tests of gravity, usually reduce to having most of their interesting effects in the non-linear regime of structure formation.

Testing such models, and extracting the maximum information that is contained in the data gathered from current and future galaxy and weak-lensing surveys, require us to include non-linear scales. This requires theoretical predictions for the matter power spectrum on these scales.

Currently this is either done using semi-analytical predictions and/or fits like HALOFIT Smith et al. (2003); Takahashi et al. (2012) and HMCode Mead et al. (2015, 2016) (which is implemented in commonly used Boltzmann codes such as CAMB Lewis and Bridle (2002) and CLASS Lesgourgues (2011)) or using an emulator such as CosmicEmu Lawrence et al. (2010); Heitmann et al. (2014); Casarini et al. (2016) and EuclidEmulator Euclid Collaboration et al. (2018). An emulator is constructed by performing a large number of N-body simulations in the parameter-space and then performing an interpolation to obtain the power spectrum for any parameter combination of interest. This can be quite expensive to make, but once it’s created it can provide non-linear spectra (usually from linear spectra) almost for free.

For $\Lambda{\rm CDM}$ both of the approaches above have been adopted and used to provide constraints from observations. For the case of modified gravity models the only fit provided so far is MGHALOFIT Zhao (2014); Hojjati et al. (2011) which has a modified HALOFIT that was calibrated using N-body simulation data to the Hu-Sawicki model Hu and Sawicki (2007). For a coupled Dark Energy model, a fitting function was developed in Casas et al. (2016) and was used to forecast the constraints on the coupling parameter $\beta$. To be able to derive constraints, or to provide forecasts for how well future experiments will constrain deviations from GR, we need to have an accurate model for the non-linear matter power spectrum. This generally has to be derived on a model by model basis.

However very recently an interesting semi-analytical method based on the halo model and nonlinear perturbation theory was proposed in Cataneo et al. (2018) and shown to be able to get $1-3\%$ accuracy out to $k\sim 1\,h{\rm Mpc}^{-1}$.

In this paper we will consider the Hu-Sawicki model as this is a representative model when it comes to using cosmology and astrophysics to constrain deviations from GR. Our aim is to provide the community with a precise fitting function for this model that can be used for example for making forecasts for future surveys. Instead of providing a fit for $P_{f(R)}$ directly, we present a fitting function for the enhancement $P_{f(R)}/P_{\Lambda{\rm CDM}}$ (both for the linear and non-linear power spectrum) as function of scale, redshift and the model parameter $f_{R0}$ which controls the size of the deviations from GR (GR is recovered as $f_{R0}\to 0$). This enhancement, as we will show, has a fairly weak cosmology dependence and we can therefore fit it using simulations from only one cosmology saving a lot of computational time. Our function have been fitted using a large suite of available N-body simulation data and can easily be incorporated in a Boltzmann code like CAMB or CLASS to scale from a non-linear $P(k)$ for $\Lambda{\rm CDM}$ (created for example using an emulator like EuclidEmulator ) to a non-linear $P(k)$ for the Hu-Sawicki model.

The requirements for the accuracy of the matter power spectrum is dictated by large upcoming surveys like LSST Ivezić et al. (2008) and Euclid Scaramella et al. (2014). Estimates for how accurate the power spectrum needs to be to take full advantage of the statistical power of such surveys varies from $1-2\%$ Huterer and Takada (2005) down as small as $0.5\%$ Hearin et al. (2012) for scales $k\lesssim 10\,h{\rm Mpc}^{-1}$. However this is ignoring111Also note that common N-body algorithms in state of the art codes disagree at the $\sim 1\%$ level already at $k=1\,h{\rm Mpc}^{-1}$ and at the $\sim 3\%$ level at $k=10\,h{\rm Mpc}^{-1}$ Schneider et al. (2016). model uncertainties in the way baryonic feedback affects the matter power spectrum. Baryonic effects are expected (from simulations and observations) to be as large as $10-30\%$ for scales $1\lesssim k\lesssim 10\,h{\rm Mpc}^{-1}$ Chisari et al. (2018); Schneider and Teyssier (2015). The accuracy of the newly released EuclidEmulator is quoted as being $\sim 1\%$ accurate for $k\lesssim 1\,h{\rm Mpc}^{-1}$ and for redshifts $z\lesssim 3.5$. Based on these considerations our aim is to produce a fit that is $\sim 1\%$ accurate for scales $k\lesssim 1\,h{\rm Mpc}^{-1}$ and $<5\%$ accurate for scales $1<k<10\,h{\rm Mpc}^{-1}$ and covering redshifts $z\lesssim 3.5$.

The setup of this paper is as follows: in Sec. II we describe the simulations we have used, in Sec. III we discuss the cosmology dependence of the enhancement $P_{f(R)}/P_{\Lambda{\rm CDM}}$, in Sec. IV we describe the fitting function we have created together with some tests, in Sec. V we show an application of the fitting formula by computing forecasts for how well galaxy clustering and weak-lensing in a Euclid-like survey will be at constraining $f(R)$ gravity before we conclude in Sec. VI.

II Simulations

We take advantage of a large set of simulations to make the fit.222All the power spectrum data that we used are available at https://github.com/HAWinther/FofrFittingFunction For more about how the modified gravity simulations are performed see e.g. Winther et al. (2015) and references within.

The main simulation suite we use is ELEPHANT Cautun et al. (2018) (WMAP9 cosmology) which has $N=1024^{3}$ particles, $L=1024\,{\rm Mpc}/h$, $\Omega_{m}=0.281$, $\Omega_{b}=0.046$, $h=0.6974$, $n_{s}=0.971$ and $\sigma_{8}=0.820$ ($A_{s}=2.3\cdot 10^{-9}$) with $|f_{R0}|=\{10^{-4},10^{-5},10^{-6},0\}$. These simulations were run with the ECOSMOG code Li et al. (2012).

For the same cosmology as above we have also run extra simulations (also with the ECOSMOG code) using $N=256^{3}$ particles in a $L=200\,{\rm Mpc}/h$ box with $|f_{R0}|=\{10^{-5},5\cdot 10^{-6},2\cdot 10^{-5},5\cdot 10^{-5},0\}$. This simulation suite contain simulations of $\Lambda{\rm CDM}$ and $f(R)$ gravity with $|f_{R0}|=10^{-5}$ for $\sigma_{8}=0.88$ and $0.72$ that allows us to test the $\sigma_{8}$ dependence of the modified gravity power spectrum enhancement.

We also use, mainly for testing and validation, the DUSTGRAIN simulations suite Giocoli et al. (2018) (Planck 2015 cosmology) which has $N=768^{3}$ particles, $L=750\,{\rm Mpc}/h$, $\Omega_{m}=0.31345$, $\Omega_{b}=0.0481$, $h=0.6731$, $n_{s}=0.9658$ and $\sigma_{8}=0.842$ ($A_{s}=2.199\cdot 10^{-9}$) with $|f_{R0}|=\{10^{-4},5\cdot 10^{-5},10^{-5},0\}$. This simulation suite contain simulations of $\Lambda{\rm CDM}$ and $f(R)$ gravity with $|f_{R0}|=10^{-5}$ and with $\Omega_{m}=0.2$ and $0.4$ that allows us to test the $\Omega_{m}$ dependence on the power spectrum. All these simulations have the same value of $\sigma_{8}$ and were run with the MG-GADGET code Puchwein et al. (2013).

All the $f(R)$ simulations mentioned above have corresponding $\Lambda{\rm CDM}$ simulations run with the same initial conditions that allows us to extract ratios $P_{f(R)}/P_{\Lambda{\rm CDM}}$ that (on the largest scales) are free of cosmic variance. For each simulation we have about $\sim 30$ redshifts between $z=0$ and $z=3$ that we use to compute the fit.

The power spectra used for the fitting functions have been estimated using POWMES Colombi et al. (2009) and (for DUSTGRAIN) by codes made by the authors. These codes have been tested and shown to give accurate results so we don’t expect any bias due to different power-spectrum evaluation codes. What could give rise to a bias is that some of the data we use comes from different N-body codes, but as shown in Winther et al. (2015) even though the actual power spectrum varies between different code-types Schneider et al. (2016) the enhancement $P_{f(R)}(k)/P_{\Lambda{\rm CDM}}(k)$ (with both spectra computed by the same code) does not and this is all that goes into our fit below.

The effect of the size of the box, mass-resolution and cosmic variance on the matter power spectrum was investigated in Li et al. (2013). We find that the spectra we extract can be trusted down to scales of $k\sim 5-10\,h{\rm Mpc}^{-1}$ depending on redshift.

III Variation of enhancement with cosmological parameters

The fit for $P_{f(R)}/P_{\Lambda{\rm CDM}}$ we perform in this paper is using N-body data from one single cosmology. The reason is that the cosmology dependence of this ratio is expected to be weak for almost all of the standard parameters with the possible exception of $\Omega_{m}$ and $\sigma_{8}$ as these correlate with the efficiency of screening and with the growth-rate of the matter density perturbation. There is also one additional effect that is potentially significant, which is how degenerate the enhancement is with baryonic feedback (see e.g. Puchwein et al. (2013); Arnold et al. (2018) for a discussion on the size of these effects compared to the modified gravity enhancement).

In this section we will go through the different cosmological parameters and check how much the modified gravity enhancement changes. We will use linear perturbation theory, the halo model, fast approximative N-body simulations Winther et al. (2017) and full N-body simulations (for the cases we have this available) to investigate this.

By the halo model we mean the prediction of Lombriser et al. (2014), which combines the modified linear power spectrum with a modified 1-halo contribution and a quasilinear correction motivated by higher-order perturbation theory Koyama et al. (2009). It incorporates the chameleon mechanism through an implementation of the thin-shell approximation in the spherical collapse model Lombriser et al. (2013). This generates a mass and environment dependent spherical collapse density, from which an environmentally averaged modified peak threshold is determined that is used to compute the $f(R)$ modification and chameleon screening effects on the halo mass function and concentration determining the 1-halo contribution. A comparison to other modelling techniques of the modified nonlinear matter power spectrum in $f(R)$ gravity can be found in Lombriser (2014).

The approximate N-body simulations we use is the COLA implementation of $f(R)$ gravity. COLA simulations are $\mathcal{O}(100-1000)$ faster than high-resolution N-body simulations, but can reproduce the enhancement of the power spectrum to % -level accuracy down to fairly non-linear scales $k\lesssim 1-5\,h{\rm Mpc}^{-1}$.

III.0.1 Massive neutrinos

The effects of massive neutrinos are highly degenerate with a modified gravity signal since massive neutrinos decrease the growth of structure on small scales, while modifications of gravity usualy enhances the growth. However we don’t expect a big change in the enhancement (i.e. for fixed cosmological parameters) and this is what we see in Fig. 1 for N-body simulations with $|f_{R0}|=10^{-4}$ Baldi et al. (2014). The variation is seen to be at the sub-percent level for $k\lesssim 1\,h{\rm Mpc}^{-1}$ and $\sim 4\%$ around $k=1\,h{\rm Mpc}^{-1}$ for the large value $m_{\nu}=0.6$ eV.

III.0.2 $\Omega_{m}$

Varying $\Omega_{m}$ changes the growth of linear perturbations so could potentially have a big effect on the enhancement. In Fig. 2 we show how much the enhancement changes for a fairly large variation of $\Omega_{m}$ (from $\Omega_{m}=0.3$ to $\Omega_{m}=0.2$ and $\Omega_{m}=0.4$) based on N-body simulations. This figure also shows that the correction is accurately captured by linear theory and/or the halo model.

A simple fit to results computed using linear theory at $z=0$ shows that the enhancement of the power-spectrum $B_{\Omega_{m}\rm~{}corr}\equiv P(k)/P_{\Omega_{m}=0.3}(k)$ given by

[TABLE]

with $\frac{\Delta\Omega_{m}}{\Omega_{m}}=\frac{\Omega_{m}-0.3}{0.3}$, $a=0.105$, $b=1.4$ and $k_{*}=0.16(10^{-5}/|f_{R0}|)^{1/2}\,h{\rm Mpc}^{-1}$ (in general these parameters will depend on redshift). This gives us, for example, that a $10\%$ change in $\Omega_{m}$ leads only to a $\sim 1\%$ change in the enhancement so this is not a large effect, but it’s straightforward to compute the correction using linear theory or the halo model if needed.

III.0.3 Clustering amplitude

In linear theory there is no variation with $\sigma_{8}$ (or more technically speaking with $A_{s}$, the primordial amplitude). However the amount of screening on non-linear scales depends crucially on how clustered matter is so this parameter could have a significant impact on non-linear scales. Even in $\Lambda{\rm CDM}$ an enhancement of $\sigma_{8}$ leads to a greater enhancement of the clustering on non-linear scales than as predicted by linear theory. In Fig. 3 we show how large this variation is based on N-body simulations together with predictions from both linear theory and the halo model.

A fit to our simulations with varying $\sigma_{8}$ ($0.72$, $0.8$ and $0.88$) shows that the effect of varying $\sigma_{8}$ can be described by a multiplicative correction to the enhancement which for $|f_{R0}|=10^{-5}$ at $z=0$ and for $k\lesssim 10\,h{\rm Mpc}^{-1}$ is approximately given by

[TABLE]

where $k_{*}=1.2\,h{\rm Mpc}^{-1}$ and $\frac{\Delta\sigma_{8}}{\sigma_{8}}\equiv\frac{\sigma_{8}-0.8}{0.8}$. A $10\%$ deviation of $\sigma_{8}$ from it’s fiducial value ($0.8$) leads to a $\sim 2-3\%$ level deviation in the enhancement for medium wave-numbers.

III.0.4 Other parameters

Linear perturbation theory predicts zero variation with other cosmological parameters such as the Hubble constant $H_{0}$ (when the fitting function is expressed in terms of $k$ in units of $\,h{\rm Mpc}^{-1}$) and the spectral index $n_{s}$. Changing the spectral index does modify the amplitude of clustering on small scales and could also influence screening, however we have checked using the halo model and COLA simulations that the expected variation is $\ll 1\%$ for all scales of interest within a reasonable variation in these parameters (here defined to be $3\sigma$ of the Planck 2018 cosmological constraints).

IV Fitting functions

IV.1 Fitting function for the linear enhancement

At the level of linear perturbations the growth of matter perturbations are determined by

[TABLE]

where

[TABLE]

from which it follows that the enhancement of the linear matter power spectrum for a general model is simply the enhancement of a fiducial model evaluated at $k^{*}=k(f_{R0}/f_{R0}^{\rm fid})^{1/2}$. We perform a fit using the fitting function

[TABLE]

which is constructed to interpolate between the expected low and high $k$ limits ($\sim 1+bk^{2}$ for small $k$ and $\sim\log(k)$ for large $k$). Similar functional forms for the fitting functions of the enhancement within a coupled Dark Energy model were used in Casas et al. (2016). The functions, $X=b,c,d,e$, above are written as

[TABLE]

where $r=\log(f_{R0}/f_{R0}^{\rm fid})$ leaving us with $12$ free parameters to fit. The best-fit from taking $|f_{R0}|^{\rm fid}=10^{-5}$ using $\Omega_{m}=0.281$ is shown in Table 1 and in Fig. 4 we show a comparison to the true result.

IV.2 Fitting function for the non-linear enhancement

The fitting function we use for our fit is given by

[TABLE]

where (for $X=b,c,d,e,f,g$)

[TABLE]

with

[TABLE]

where $r=\log(f_{R0}/f_{R0}^{\rm fid})$. This makes $54$ free parameters for the full $f_{R0}$, scale and redshift dependence.

We choose to make three different fits: one using $|f_{R0}|=\{10^{-4},5\cdot 10^{-5},10^{-5}\}$ (high), one using $|f_{R0}|=\{10^{-5},5\cdot 10^{-6},10^{-6}\}$ (medium) and one using $|f_{R0}|=\{10^{-5},5\cdot 10^{-6},10^{-6},10^{-7}\}$ (low) and then interpolate between these (overlapping) fits. The best-fit parameters we find are given in Table 2 and the agreement as a function of scale and redshift can be seen in Fig. 5. This figure shows that the fit is good to $\lesssim 1$% for most scales and redshifts, with the exception of the smallest scales for the largest values of $|f_{R0}|$.

In Fig. 6 we perform a test of our fitting function by predicting the enhancement for $|f_{R0}|=2\cdot 10^{-5}$, a value that was not used to generate the fit. This value is in the middle of $|f_{R0}|=10^{-5}$ and $|f_{R0}|=5\cdot 10^{-5}$ and should therefor give us a good estimate for the accuracy of our fit. The agreement is within $\sim 1-2$% for $k\lesssim 1\,h{\rm Mpc}^{-1}$ and $\sim 2-3$% for $k\lesssim 10\,h{\rm Mpc}^{-1}$.

In Fig. 7 we test our fitting function by comparing it to simulations with a different background cosmology from that is used to make the fit. As expected from the discussion in the previous section the agreement is very good.

These numbers should be contrasted to the enhancement of the matter power spectrum itself relative to $\Lambda{\rm CDM}$ which is typically $10-40\%$.

To test the fitting function we create mock data and try to perform a fit to $P(k,z)$ for a Euclid-like survey with $V=50(\text{Gpc}/h)^{3}$ with $n_{\rm gal}=10^{-3}(\,{\rm Mpc}/h)^{-3}$ using the diagonal likelihood

[TABLE]

where $C^{-1}=\frac{Vk^{2}\Delta k}{4\pi^{2}(P^{\rm fid}(k,z)+1/n)^{2}}$. The sum is over $6$ evenly spaced $z$-bins from $z=0$ to $z=2$ and $30$ logarithmically spaced $k$-bins between $k=10^{-4}\,h{\rm Mpc}^{-1}$ and $k=5\,h{\rm Mpc}^{-1}$. This assumes Gaussian fluctuations on all scales which significantly underestimate the errors on non-linear scales and also does not take into account uncertainty of unknown baryonic physics. Thus any fit based on this would be completely dominated by the smallest scales. To get more realistic errors we try to take this into account by imposing a minimum $1\%$ error on non-linear scales starting at $k=0.5\,h{\rm Mpc}^{-1}$ and growing to $10\%$ at $k=10\,h{\rm Mpc}^{-1}$. The fiducial power spectrum is generated using the Eisenstein-Hu fitting function for $\Lambda{\rm CDM}$ Eisenstein and Hu (1998, 1999), converted to a non-linear power spectrum using HALOFIT and finally multiplied by the $f(R)$ enhancement found in simulations with $|f_{R0}|=2\cdot 10^{-5}$. The result of fits to $P(k,z)$ around $z=1$ (where our fit is seen to deviate a bit from the N-body result) can be found in Fig. 9. We also performed this test for other redshifts with similar results.

A fitting function, based on HALOFIT , for the Hu-Sawicki model already exists in the literature, namely MGHALOFIT . In Fig. 8 we show a comparison of our fitting formula to MGHALOFIT which shows that our fit performs better.

V Forecast

In this section we will give an example application of the fitting formula presented in this paper by using it to compute forecasts for how well the Hu-Sawicki model can be constrained in future surveys (see Alonso et al. (2017) for forecasts for general scalar-tensor theories). We adopt the Fisher-matrix formalism to the main cosmological observables for next-generation galaxy surveys, namely Galaxy Clustering (GC) and Weak Lensing (WL). The matter power spectrum is computed using CAMB with the HALOFIT prescription to get the non-linear power spectrum for $\Lambda{\rm CDM}$ and then use our fitting-formula to go from $\Lambda{\rm CDM}$ to $f(R)$. We neglect cross-correlations between GC and WL. To perform our forecasts we use the survey parameters for a Euclid-like mission Laureijs et al. (2011).

For GC the main observable is the observed (redshift-space) galaxy power spectrum which we model as

[TABLE]

where $D_{A}$ is the angular diameter distance, $H$ is the Hubble function, $\mu$ is the cosine of the angle between the line of sight and the Fourier vector $\vec{k}$, $f(z)$ is the growth function and $\sigma^{2}\equiv\sigma_{r}^{2}+\sigma_{v}^{2}$ is a parameter parametrising errors induced by spectroscopic redshift measurements ($\sigma_{r}$) and the fingers-of-god effect ($\sigma_{v}$) that we marginalise over. A subscript f denotes the value in the fiducial cosmology which we take to be a $\Lambda{\rm CDM}$ cosmology with parameters $n_{s}=0.96$, $10^{9}A_{s}=2.126$, $h=0.67$, $\Omega_{m}=0.32$, $\Omega_{b}=0.05$, $\Omega_{\Lambda}=0.68$, $m_{\nu}=0.06$ eV and $\sigma_{v}=300\text{ km/s}$. The fiducial $f(R)$ parameter is taken to be $|f_{R0}|=10^{-5}$.

The Fisher matrix for GC is taken to be

[TABLE]

where the data vector $D=P_{g}(k,\mu,z)+1/n(z)$ with $V_{\rm survey}$ being the survey volume and $n(z)$ is the galaxy number density. We compute the constraints for two different values of $k_{\rm max}=0.15$ and $0.25$ both with $k_{\rm min}=0.008\,h/$Mpc. For computing the Fisher matrix we used $9$ $z$-bins in the range $z=0.95$ and $z=1.75$.

The second probe we include is WL cosmic shear: the distortions in the ellipticities of galaxies due to bending of light around large cosmic structures. The cosmic shear is computed using $10$ redshift bins in redshift range $0<z<2.5$. The cosmic shear at a redshift bin $i$ is correlated with the cosmic shear at another redshift bin $j$ since light coming from each bin will propagate through some of the same structures on the way to us. The cross power spectrum of cosmic shear in bin $i$ and $j$ is determined by the underlying dark matter power spectrum via

[TABLE]

where $r(z)$ is the co-moving distance and $W$ is a a window function given by the photometric redshift distribution function and the galaxy number density distribution.

The Fisher matrix for WL is given by

[TABLE]

where $C$ is a matrix with elements $C_{ij}$, $f_{\rm sky}=0.36$ ($15000\text{ deg}^{2}$) is the sky-fraction covered by the survey, $\Delta\ell$ is the size of the $\ell$-bins, $\ell_{\rm max}$ is the maximum multipole number and Cov is the WL covariance matrix. In this paper we consider the two values $\ell_{\rm max}=1000$ and $3000$. We used $100$ logarithmically spaced $\ell$-bins between $\ell_{\rm min}=100$ and $\ell_{\rm max}$ and we used $10$ evenly spaced $z$-bins between $z=0.15$ and $z=2.5$. Apart from the particular numbers quoted above we use the same setup as Casas et al. (2017) so see this paper for more details.

In Fig. 10 we show the Fisher forecast constraints we obtained using the fiducial value $|f_{R0}|=10^{-5}$, which is slightly below the best constraints from current cosmological data (see e.g. Dossett et al. (2014)). We show how the results change when going from only using fairly linear scales to including more and more non-linear scales ($k_{\rm max}=0.15\,h{\rm Mpc}^{-1},\ell_{\rm max}=1000$ versus $k_{\rm max}=0.25\,h{\rm Mpc}^{-1},\ell_{\rm max}=3000$) in the forecast. In Table 3 we show the marginalized constraints on the cosmological parameters for the different cases we have considered.

GC and WL are individually able to constrain $\log_{10}|f_{R0}|$ to $\sim 5\%$ and $\sim 15\%$ respectively depending on how many non-linear modes we include in the forecast. The modified gravity parameter $f_{R0}$ is seen to be mostly degenerate with the clustering amplitude $A_{s}$. Combining GC and WL we are able to break this degeneracy and bring the potential constraints down to $\sim 1-2\%$ ($\Delta f_{R0}\lesssim 2\cdot 10^{-6}$). However we caution that this is a simplified forecast not taking baryonic effects on the matter power spectrum or including massive neutrinos which both are known to be degenerate with a potential modified gravity signal Baldi et al. (2014); Arnold et al. (2019).

VI Conclusions

In this paper we have presented a fitting function for the linear and non-linear matter power spectrum of the Hu-Sawicki $f(R)$ model using power spectra computed from N-body simulations for several different values of the model parameters. This is one of the most studied modified gravity models and often the fiducial choice when trying to constrain modified gravity effects in observational data.

We have shown that the enhancement has a weak cosmology dependence which allows us to make the fit for a fixed cosmology. Any cosmology dependence can be accurately included by using inexpensive tools like linear theory, the halo model or COLA simulations.

We have demonstrated that the fitting function is accurate over a large range of scales and redshifts. We also tested it against the existing MGHALOFIT code and found that our fitting function generally performs better.

One can easily integrate our fitting function in any approach that produces a non-linear matter power spectrum for $\Lambda{\rm CDM}$. With this paper we provide the fitting function already implemented333This can be found at https://github.com/HAWinther/FofrFittingFunction in many common programming languages like C, Fortran and Python plus an implementation in both CAMB and CLASS .

Finally, as an application, we have used the fitting functions to compute Fisher forecasts for how well a Euclid-like survey will be at constraining the Hu-Sawicki model. We find that the potential constraints from combining GC and WL from a Euclid-like survey, when including a reasonable amount of non-linear scales in the forecast, are at the $\Delta f_{R0}\sim 2\cdot 10^{-6}$ level. This is without taking into account baryonic effects in the power spectrum and including massive neutrinos, but nevertheless shows the potential constraining power in future survey when including non-linear scales.

The same approach as used in this paper can likely be applied to cheaply create emulators for other non-standard gravity models, possibly in conjunction with the semi-analytical method of Cataneo et al. (2018).

Acknowledgment

KK and HAW are supported by the European Research Council through 646702 (CosTesGrav). KK is also supported by the UK Science and Technologies Facilities Council grants ST/N000668/1. GBZ is supported by NSFC Grant No. 11673025, and by a Royal Society-Newton Advanced Fellowship. LL is supported by a Swiss National Science Foundation Professorship grant (No. 170547). SC acknowledges support from CNRS and CNES grants. BL acknowledge the support of the UK STFC Consolidated Grants (ST/P000541/1 and ST/L00075X/1) and Durham University. BL is also supported in part by the European Union’s Horizon 2020 research and innovation program. MB acknowledges support from the Italian Ministry for Education, University and Research (MIUR) through the SIR individual grant SIMCODE (project number RBSI14P4IH), from the grant MIUR PRIN 2015 ”Cosmology and Fundamental Physics: illuminating the Dark Universe with Euclid”, and from the agreement ASI n.I/023/12/0 “Attività relativealla fase B2/C per la missione Euclid”.