Theoretical Accuracy in Cosmological Growth Estimation

TL;DR

This paper compares theoretical templates based on general relativity and nDGP gravity for estimating cosmological growth from redshift space distortions, highlighting the importance of consistent non-linear modeling for future surveys.

Contribution

It demonstrates the necessity of gravity-model specific non-linear treatments in RSD analysis to accurately recover growth parameters in modified gravity scenarios.

Findings

GR template fails to recover fiducial growth rate at z=1 within 1σ.

DGP template accurately recovers the fiducial growth rate.

Bias in GR template becomes significant for models with enhanced non-linearity.

Abstract

We elucidate the importance of the consistent treatment of gravity-model specific non-linearities when estimating the growth of cosmological structures from redshift space distortions (RSD). Within the context of standard perturbation theory (SPT), we compare the predictions of two theoretical templates with redshift space data from COLA (COmoving Lagrangian Acceleration) simulations in the normal branch of DGP gravity (nDGP) and General Relativity (GR). Using COLA for these comparisons is validated using a suite of full N-body simulations for the same theories. The two theoretical templates correspond to the standard general relativistic perturbation equations and those same equations modelled within nDGP. Gravitational clustering non-linear effects are accounted for by modelling the power spectrum up to one loop order and redshift space clustering anisotropy is modelled using the Taruya, Nishimichi and Saito (TNS) RSD model. Using this approach, we attempt to recover the simulation's fiducial logarithmic growth parameter $f$. By assigning the simulation data with errors representing an idealised survey with a volume of $10\mbox{Gpc}^3/h^3$, we find the GR template is unable to recover fiducial $f$ to within 1$\sigma$ at $z=1$ when we match the data up to $k_{\rm max}=0.195h$/Mpc. On the other hand, the DGP template recovers the fiducial value within $1\sigma$. Further, we conduct the same analysis for sets of mock data generated for generalised models of modified gravity using SPT, where again we analyse the GR template's ability to recover the fiducial value. We find that for models with enhanced gravitational non-linearity, the theoretical bias of the GR template becomes significant for stage IV surveys. Thus, we show that for the future large data volume galaxy surveys, the self-consistent modelling of non-GR gravity scenarios will be crucial in constraining theory parameters.

Figures

Tables

Table 1. Table 1: Summary of template performances at z = 0.5 𝑧 0.5 z=0.5 for nDGP simulations where fiducial f = 0.783 𝑓 0.783 f=0.783 .

Template	$k_{max}\left[h/\text{Mpc}\right]$	bins	$V_s{\left[\text{Gpc}/h\right]}^3$	$f\pm 2\sigma$	${\sigma }_v\pm 2\sigma \left[\text{Mpc}/h\right]$
GR	0.110	18	10	$0.774{\pm }_{0.022}^{0.034}$	$5.92{\pm }_{0.53}^{0.43}$
DGP	0.110	18	10	$0.779{\pm }_{0.029}^{0.031}$	$5.86{\pm }_{0.44}^{0.42}$
GR	0.135	22	10	$0.769{\pm }_{0.016}^{0.029}$	$5.86{\pm }_{0.16}^{0.22}$
DGP	0.135	22	10	$0.779{\pm }_{0.022}^{0.022}$	$5.86{\pm }_{0.19}^{0.19}$
GR	0.147	24	10	$0.777{\pm }_{0.025}^{0.025}$	$5.96{\pm }_{0.15}^{0.14}$
DGP	0.147	24	10	$0.786{\pm }_{0.018}^{0.019}$	$5.89{\pm }_{0.13}^{0.12}$

Table 2. Table 2: Summary of template performances for nDGP simulations at z = 1 𝑧 1 z=1 where fiducial f = 0.909 𝑓 0.909 f=0.909 .

Template	$k_{max}\left[h/\text{Mpc}\right]$	bins	$V_s{\left[\text{Gpc}/h\right]}^3$	$f\pm 2\sigma$	${\sigma }_v\pm 2\sigma \left[\text{Mpc}/h\right]$
GR	0.147	24	10	$0.895{\pm }_{0.025}^{0.028}$	$4.08{\pm }_{0.26}^{0.22}$
DGP	0.147	24	10	$0.902{\pm }_{0.025}^{0.025}$	$4.05{\pm }_{0.20}^{0.19}$
GR	0.171	28	10	$0.902{\pm }_{0.020}^{0.021}$	$4.19{\pm }_{0.09}^{0.10}$
DGP	0.171	28	10	$0.912{\pm }_{0.019}^{0.015}$	$4.16{\pm }_{0.10}^{0.10}$
GR	0.195	32	10	$0.893{\pm }_{0.018}^{0.018}$	$4.13{\pm }_{0.07}^{0.06}$
DGP	0.195	32	10	$0.908{\pm }_{0.014}^{0.014}$	$4.10{\pm }_{0.06}^{0.05}$
GR	0.195	32	20	$0.893{\pm }_{0.013}^{0.012}$	$4.13{\pm }_{0.05}^{0.04}$
DGP	0.195	32	20	$0.908{\pm }_{0.010}^{0.010}$	$4.10{\pm }_{0.04}^{0.04}$

Table 3. Table 3: Summary of template performances for GR simulations at z = 1 𝑧 1 z=1 where fiducial f = 0.859 𝑓 0.859 f=0.859 .

Template	$k_{max}\left[h/\text{Mpc}\right]$	bins	$V_s{\left[\text{Gpc}/h\right]}^3$	$f\pm 2\sigma$	${\sigma }_v\pm 2\sigma \left[\text{Mpc}/h\right]$
GR	0.171	28	10	$0.851{\pm }_{0.020}^{0.020}$	$3.94{\pm }_{0.14}^{0.13}$
DGP	0.171	28	10	$0.869{\pm }_{0.010}^{0.019}$	$4.01{\pm }_{0.11}^{0.09}$

Equations

$$d{s}^2=-(1+2\Phi )d{t}^2+a^2(1-2\Psi ){\delta }_{ij}d{x}^{i}d{x}^j.$$

$$\delta (\mathbf{x})=\frac{{\rho }_m(\mathbf{x})-\bar{\rho }}{\bar{\rho }}$$

$$\theta({\bf x})=-\frac{\nabla\cdot v_{p}({\mbox{\boldmath$x$}})}{aH(a)f}$$

$$\displaystyle a\frac{\partial\delta(\mbox{\boldmath$k$})}{\partial a}+\theta(\mbox{\boldmath$k$})=-\int\frac{d^{3}\mbox{\boldmath$k$}_{1}d^{3}\mbox{\boldmath$k$}_{2}}{(2\pi)^{3}}\delta_{\rm D}(\mbox{\boldmath$k$}-\mbox{\boldmath$k$}_{1}-\mbox{\boldmath$k$}_{2})\alpha(\mbox{\boldmath$k$}_{1},\mbox{\boldmath$k$}_{2})\,\theta(\mbox{\boldmath$k$}_{1})\delta(\mbox{\boldmath$k$}_{2}),$$

$$\displaystyle a\frac{\partial\theta(\mbox{\boldmath$k$})}{\partial a}+\left(2+\frac{aH^{\prime}}{H}\right)\theta(\mbox{\boldmath$k$})-\left(\frac{k}{a\,H}\right)^{2}\,\Phi(\mbox{\boldmath$k$})=-\frac{1}{2}\int\frac{d^{3}\mbox{\boldmath$k$}_{1}d^{3}\mbox{\boldmath$k$}_{2}}{(2\pi)^{3}}\delta_{\rm D}(\mbox{\boldmath$k$}-\mbox{\boldmath$k$}_{1}-\mbox{\boldmath$k$}_{2})\beta(\mbox{\boldmath$k$}_{1},\mbox{\boldmath$k$}_{2})\,\theta(\mbox{\boldmath$k$}_{1})\theta(\mbox{\boldmath$k$}_{2}),$$

$$\displaystyle\alpha(\mbox{\boldmath$k$}_{1},\mbox{\boldmath$k$}_{2})=1+\frac{\mbox{\boldmath$k$}_{1}\cdot\mbox{\boldmath$k$}_{2}}{|\mbox{\boldmath$k$}_{1}|^{2}},\quad\quad\beta(\mbox{\boldmath$k$}_{1},\mbox{\boldmath$k$}_{2})=\frac{(\mbox{\boldmath$k$}_{1}\cdot\mbox{\boldmath$k$}_{2})\left|\mbox{\boldmath$k$}_{1}+\mbox{\boldmath$k$}_{2}\right|^{2}}{|\mbox{\boldmath$k$}_{1}|^{2}|\mbox{\boldmath$k$}_{2}|^{2}}.$$

$$-\left(\frac{k}{aH}\right)^{2}\Phi=\frac{3\Omega_{m}(a)}{2}\mu(k,a)\,\delta(\mbox{\boldmath$k$})+S(\mbox{\boldmath$k$}),$$

$$\displaystyle S(\mbox{\boldmath$k$})$$

$${\delta }_n(\mathbf{k},a)$$

$${\theta }_n(\mathbf{k},a)$$

$$P_{ij}^{1-\mathrm{loop}}(k)=P_0(k)+P_{ij}^{22}(k)+P_{ij}^{13}(k)$$

$$\langle\delta_{0}(\mbox{\boldmath$k$})\delta_{0}(\mbox{\boldmath$k$}^{\prime})\rangle=(2\pi)^{3}\delta_{\rm D}(\mbox{\boldmath$k$}+\mbox{\boldmath$k$}^{\prime})\,P_{0}(k).$$

$$\displaystyle\langle g_{i}^{2}(\mbox{\boldmath$k$})g_{j}^{2}(\mbox{\boldmath$k$}^{\prime})\rangle$$

$$\displaystyle\langle g_{i}^{1}(\mbox{\boldmath$k$})g_{j}^{3}(\mbox{\boldmath$k$}^{\prime})+g_{i}^{3}(\mbox{\boldmath$k$})g_{i}^{1}(\mbox{\boldmath$k$}^{\prime})\rangle$$

$$P_K^S(k,\mu )=(b+f{\mu }^2{)}^2{P}_{\delta \delta }(k)$$

$$P^S(k,\mu )=\text{mbox}{D}_{\text{mbox}FoG}(fk\mu {\sigma }_v)\{P_{\delta \delta }(k)-2{\mu }^2{P}_{\delta \theta }(k)+{\mu }^4{P}_{\theta \theta }(k)+A(k,\mu )+B(k,\mu )\}$$

$$A(k,\mu )=-(k\mu )\int d^3{\mathbf{k}}^{\prime }\frac{k_z^{\prime }}{k^{\prime 2}}\{B_{\sigma }({\mathbf{k}}^{\prime },\mathbf{k}-{\mathbf{k}}^{\prime },-\mathbf{k})-B_{\sigma }({\mathbf{k}}^{\prime },\mathbf{k},-\mathbf{k}-{\mathbf{k}}^{\prime })\}$$

$$B(k,\mu )=(k\mu {)}^2\int d^3{\mathbf{k}}^{\prime }F({\mathbf{k}}^{\prime })F(\mathbf{k}-{\mathbf{k}}^{\prime })$$

$$F(\mathbf{k})=\frac{k_z}{k^2}\left[P_{\delta \theta }(k)-\frac{k_z^2}{k^2}{P}_{\theta \theta }(k)\right]$$

$${\delta }_D({\mathbf{k}}_1+{\mathbf{k}}_2+{\mathbf{k}}_3)B_{\sigma }({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)=⟨\theta ({\mathbf{k}}_1)\{\delta ({\mathbf{k}}_2)-\frac{k_{2z}^2}{k_2^2}\theta ({\mathbf{k}}_2)\}\{\delta ({\mathbf{k}}_3)-\frac{k_{3z}^2}{k_3^2}\theta ({\mathbf{k}}_3)⟩\}$$

$$P_{\ell }^{(S)}(k)=\frac{2\ell +1}{2}{\int }_{-1}^{1}d\mu P_{TNS}^{(S)}(k,\mu ){\mathcal{P}}_{\ell }(\mu )$$

$${\sigma }_{\mathrm{norm}}(a)=\frac{F_1^{\mathrm{fiducial}}(a)}{F_1(a)}$$

$$-2\ln \mathcal{L}=\sum _n\sum _{\ell ,{\ell }^{\prime }=0,2}\left(P_{\ell ,\mathrm{data}}^{(S)}(k_n)-P_{\ell ,\mathrm{model}}^{(S)}(k_n)\right)\text{mbox}{Cov}_{\ell ,{\ell }^{\prime }}^{-1}(k_n)\left(P_{{\ell }^{\prime },\mathrm{data}}^{(S)}(k_n)-P_{{\ell }^{\prime },\mathrm{model}}^{(S)}(k_n)\right)$$

$${\chi }_{MG}^2(k_{\max })=\frac{1}{N_{dof}}\sum _{\ell }\sum _k^{k_{\max }}\text{mbox}{Cov}_{\ell ,\ell }^{-1}(k)[P_{\ell }^{DGP}(k)-P_{\ell }^{LCDM}(k){]}^2$$

$$ϵ\frac{H}{r_c}=H^2(1-{\Omega }_m(a))$$

$$\mu (k;a)=1+\frac{1}{3\beta (a)}$$

$$\gamma_{2}(k,\mbox{\boldmath$k$}_{1},\mbox{\boldmath$k$}_{2};a)=-\frac{H_{0}^{2}}{24H^{2}\beta(a)^{3}\Omega_{rc}}\left(\frac{\Omega_{m0}}{a^{3}}\right)^{2}(1-\mu_{1,2}^{2})$$

$$\gamma_{3}(k,\mbox{\boldmath$k$}_{1},\mbox{\boldmath$k$}_{2},\mbox{\boldmath$k$}_{3};a)=\frac{H_{0}^{2}}{144H^{2}\beta(a)^{5}\Omega_{rc}^{2}}\left(\frac{\Omega_{m0}}{a^{3}}\right)^{3}(1-\mu_{1,2}^{2})(1-\mu_{1,23}^{2})$$

$$\beta (a)=1+\frac{H}{{\Omega }_{rc}}\left(1+\frac{a{H}^{\prime }}{3H}\right)$$

Full text

I Introduction

The accelerated expansion of the universe has been the focus of strenuous research since its discovery almost 20 years ago Riess et al. (1998); Perlmutter et al. (1999) with the concordance model of cosmology performing very well in explaining observational data on large scales. Most of this success is based on a homogenous and isotropic solution within Einstein’s GR. This standard model ($\Lambda$ Cold Dark Matter -LCDM) assumes that GR is the valid theory of gravity at large scales. Additionally, GR does just as well locally, having passed all tests within the solar system with flying colours. Further, a recent remarkable measurement made by the LIGO consortium has also confirmed the existence of gravitational waves Abbott et al. (2016), an early GR predicted phenomena. On the other hand, theoretically the standard model of cosmology leaves a lot to be desired with an ubiquitous and vastly dominant dark energy component having to be introduced. A natural explanation for this energy has lead to the largest fine tuning problem in physics Weinberg (1989); Martin (2012) and the resulting search for a more natural solution has brought researchers to the tricky task of modifying GR Hu and Sawicki (2007); Starobinsky (1980); Dvali et al. (2000); Gomes and Amendola (2014); Hassan and Rosen (2012); de Rham and Heisenberg (2011); de Rham and Gabadadze (2010); de Rham et al. (2011); Maartens and Koyama (2010).

A common feature of these modifications is the introduction of an additional fifth force sourced by the modified theory’s additional degree of freedom. To deal with observational constraints, these modifications make use of the so called screening mechanisms, shielding matter from the fifth force at small scales (See Koyama (2016); Sakstein (2015); Clifton et al. (2012) for reviews). This vastly diminishes the distinguishability of these models at small scales, while the unscreened fifth force outside the screened regime makes them very different from GR. The large scale structure (LSS) of the universe has thus become an ideal laboratory for testing the phenomenology of modifications to gravity, a fact widely discussed in the literature Uzan and Bernardeau (2001); Lue et al. (2004); Ishak et al. (2006); Knox et al. (2006); Koyama (2006); Chiba and Takahashi (2007); Amendola et al. (2008); Simpson et al. (2013); Terukina and Yamamoto (2012); Terukina et al. (2014); Yamamoto et al. (2010); Jain and Zhang (2008); Zhao et al. (2009a, b); Asaba et al. (2013a).

One of the most powerful sources of cosmological information is the signal encoded in anisotropic galaxy clustering observed in the redshift space 111Redshift space is a coordinate system where the radial distance to a galaxy is obtained by its observed redshifts., an effect commonly denoted as Redshift Space Distortions (RSD) Kaiser (1987).This effect is due to the non-linear mapping between real position and redshift space, the linear part coming from the usual cosmological redshift-distance relation while the non-linear part arising from a galaxy’s peculiar velocity. In most modified gravity scenarios the fifth force couples to matter, enhancing matter’s peculiar velocity and producing an imprint on the redshift space clustering statistics Linder (2008); Guzzo et al. (2008); Yamamoto et al. (2008); Song and Percival (2009); Song and Kayo (2010); Guzik et al. (2010); Song et al. (2011); Asaba et al. (2013b); Hellwing et al. (2014).

At linear scales, the observed galaxy clustering is usually translated into a measurement of the parameter combination $f\sigma_{8}$ where $f=d\ln(F_{1})/d\ln(a)$, with $F_{1}$ being the linear growth of structure and $a$ being the scale factor, while $\sigma_{8}$ normalises the linear power spectrum (see eg. Peebles (1980)). In the last decade, the extraction of the $f\sigma_{8}$ parameter from the observed RSD signal in galaxy spectroscopic surveys has become a common practice, and a number of such estimates have provided us with current state-of-the-art cosmological constraints. Here, among others, we can mention the measurements using luminous red galaxy (LRG) sample of the Sloan Digital Sky Survey Oka et al. (2014), the BOSS survey galaxy samples Samushia et al. (2014); Sanchez et al. (2013, 2014); Gil-Mar n et al. (2016); Sanchez et al. (2016) and the VIPERS survey de la Torre et al. (2013). However, an important caution needs to be made here. So far all of these state-of-the-art constraints on $f$ were obtained by employing a standard GR approach to RSD modelling, even when the aim was to place constraints on modified models of gravity. The question naturally arises: * is the flexibility of standard templates, with the inclusion of nuisance parameters such as velocity dispersions, enough to encompass MG effects or should the analysis be done in a model dependent way?* In other words, are the systematic biases induced by using a GR-based RSD approach small enough that, when compared to the statistical errors in the data, can be safely ignored? Such a question is very timely now at the advent of precision cosmology, which will be fuelled by vast amounts of data from upcoming large observational endeavours.

Whether the GR treatment gives biased results for the constraints on modified gravity model parameters or not depends on the approach used in RSD modelling and also the precision of the measurements. Clearly, if the errors in the measurements are larger than the bias introduced by the inaccurate modelling of RSD, there is no immediate call to change the RSD modelling. This condition should be tested before the analysis pipeline is applied to real data. For the TNS model Taruya et al. (2010), which has been used to measure $f\sigma_{8}$ using the power spectrum measurements of the BOSS survey Beutler et al. (2016), this has been done only for $f(R)$ gravity. Using this model, it was shown that the standard GR template gives a biased estimation of the model parameter assuming an ideal survey with a volume of $10\mbox{Gpc}^{3}/h^{3}$ at $z=1$ Taruya et al. (2013). On the other hand, in Barreira et al. (2016) they find that with their analysis of the redshift space power spectrum for the braneworld model of gravity by Dvali, Gabadadze and Porrati (DGP) Dvali et al. (2000), negligible model bias is found. In that work the authors employ the growth rate estimation pipeline of Sanchez et al. (2016). It is clear that a general way to deal with this systematic is needed, especially in preparation for stage IV surveys such as the Dark Energy Spectroscopic Instrument (DESI) 222http://desi.lbl.gov/ and the ESA/Euclid survey333www.euclid-ec.org, which will significantly reduce observational uncertainties.

It is worth mentioning that if one is able to analyse the data in a model dependant way without cost then the answer to the above question is clear. In a previous work Bose and Koyama (2016) we presented a code which includes model specific non-linearities when constructing the TNS redshift space power spectrum, giving us access to a wide variety of theoretical templates. The non-linearities are constructed in the context of standard Eulerian perturbation theory (SPT) (See Bernardeau et al. (2002) for a comprehensive review).

In this work we present predictions generated using the framework described in Bose and Koyama (2016) as well as further investigate the issue of model bias, extending the work of Barreira et al. (2016), and providing a means of quickly testing the validity of constraining non-GR models using the standard GR template. In addition, our aim is to assess the level and scales at which theoretical modelling bias becomes important and starts to affect the results in a significant way. We do this first in terms of taking the analysis to increasing non-linear scales. With respect to non-linear scales, the errors on $f\sigma_{8}$ depend sensitively on the maximum $k$ used in the analysis. By including smaller scales the predictions of each template become more unique and so model bias becomes more pronounced. The maximum scale we can include depends on the redshift Carlson et al. (2009) and so there is also a redshift dependency on the systematic given our theoretical approach. Lastly, we stress-test the GR template to obtain a limit at which it becomes significantly biased due the enhanced small-scale dynamics induced by MG. We do this by increasing the modification to gravity at non-linear scales in MG mock data.

This paper is organised as follows: In Sec.II we review the relevant theory behind the theoretical templates used to fit the RSD. In Sec.III we test the realm of validity of the framework using results from both LCDM and DGP N-body simulations as well as MG-PICOLA simulations Winther et al. (2017). We then conduct an MCMC estimation of $f$ using both the standard GR template and the full MG template for varying inclusion of scales. This is done using a suite of 20 $1\mbox{Gpc}^{3}/h^{3}$ MG-PICOLA simulations. The same analysis is then done for mock data generated using SPT in order to mimic an ideal survey. We do this for varying levels of fifth force interactions. Finally, we summarise our results and highlight future work in Sec.IV.

II Theoretical Template for Growth Estimation

II.1 The Perturbative Treatment : SPT

Our starting point will be perturbations in a universe described by the Friedman-Robertson-Walker (FRW) metric in the Newtonian gauge:

[TABLE]

The evolution of matter and velocity perturbations within the Hubble horizon, under the quasi-static approximation are given by the Euler and continuity equations. Specifically, we consider the evolution of the density contrast

[TABLE]

where $\rho_{m}(\bf{x})$ and $\bar{\rho}$ are the matter densities at a given point and the background respectively. We also consider the evolution of the peculiar velocity field $v_{p}(\bf{x})$. At the scales of interest (ie. well above the scales of virialised, collapsed objects), the cosmic velocity filed is to a good approximation curl-free (see eg.Pichon and Bernardeau (1999); Pueblas and Scoccimarro (2009); Libeskind et al. (2013)), thus as a potential flow it can be fully characterised by the divergence only part of the vector. In our analysis we will use the velocity divergence, here defined as:

[TABLE]

where $H(a)$ is the Hubble function. As usual convention, we choose to work in the Fourier space, where the many SPT expressions become simpler. Beyond linear order, mode coupling terms are introduced and in Fourier space the equations are given by (see eg.Bernardeau et al. (2002))

[TABLE]

the prime denoting a scale factor derivative w.r.t cosmic time and the mode coupling kernels, $\alpha$ and $\beta$, are given by

[TABLE]

At linear order, $\alpha=\beta=0$. Gravity enters the perturbation’s evolution through the Poisson term including the Newtonian Potential $\Phi$. We can parametrise a wide range of modifications to gravity by writing out the Poisson term as Koyama et al. (2009); Bose and Koyama (2016)

[TABLE]

where $\Omega_{m}(a)=8\pi G\rho_{m}/3H^{2}$. The function S(\mbox{\boldmathk}) is the non-linear source term up, to the third order given by Bose and Koyama (2016)

[TABLE]

The functions $\mu(k,a)$, \gamma_{2}(\mbox{\boldmathk},\mbox{\boldmathk}_{1},\mbox{\boldmathk}_{2};a) and \gamma_{3}(\mbox{\boldmathk},\mbox{\boldmathk}_{1},\mbox{\boldmathk}_{2},\mbox{\boldmathk}_{3};a) encode the modification to gravity. Particularly, \gamma_{2}(\mbox{\boldmathk},\mbox{\boldmathk}_{1},\mbox{\boldmathk}_{2};a) and \gamma_{3}(\mbox{\boldmathk},\mbox{\boldmathk}_{1},\mbox{\boldmathk}_{2},\mbox{\boldmathk}_{3};a) contain the non-linear part of the information about the theory of gravity, in particular about the screening mechanism. These functions will completely specify the difference between theoretical templates used to estimate the growth parameter later on. In Appendix A we specify the form of these functions for the Vainshtein screened DGP model of gravity Dvali et al. (2000).

Classically, the approach of SPT is to solve eq.(4) and eq.(5) perturbatively, with n-th order solutions given by

[TABLE]

where $\bm{k}_{1...n}=\bm{k}_{1}+...+\bm{k}_{n}$ and F_{i}(\mbox{\boldmathk}_{1},\mbox{\boldmathk}_{2}...) and G_{i}(\mbox{\boldmathk}_{1},\mbox{\boldmathk}_{2}...) are the $i^{th}$ order kernels which we solve perturbatively for. Using the numerical algorithm described in Taruya (2016) we can do this for the perturbations up to any order in a general way. With the perturbative solutions up to 3rd order we can construct the 1-loop power spectrum

[TABLE]

where $P_{0}(k)$ is the linear power spectrum defined as

[TABLE]

and

[TABLE]

where $g^{i}_{1}=\delta_{i}$ and $g^{i}_{2}=\theta_{i}$.

Eq.11 describes the first order contribution to the power spectrum using the density and velocity divergence fields as our expansion variables. The inclusion of the higher order loop terms has been shown to improve the prediction of theory Jeong and Komatsu (2006), an improvement more pronounced at higher redshift Carlson et al. (2009). The regime of applicability of SPT with a required percent-level precision was assessed by comparing with N-body simulations. Presently, the best estimates are given by Jeong and Komatsu (2006), who found $\Delta^{2}(k,a)\leq 0.4$. Here $\Delta^{2}(k,a)=k^{3}P_{0}(k,a)/(2\pi^{2})$ is the dimensionless power spectrum which is a function of time, and thus so are the scales of SPT applicability. This being said, it is well known that SPT suffers from convergence problems, with higher loop contributions being of comparable size to lower ones, a problem that worsens at late times Carlson et al. (2009). This problem is evident from the 1-loop power spectrum expressions. Mode coupling introduces an integral over all scales producing a small scale dependency. Errors in the small scale regime thus sneak into our large scale predictions when we go beyond linear theory. This makes comparisons to data limited with the loss of very constraining small scale information. A number of analytical treatments have been proposed in the literature Taruya and Hiramatsu (2008); Matarrese and Pietroni (2007); Valageas (2007); Matsubara (2008); Pietroni (2008); Anselmi and Pietroni (2012); Blas et al. (2016) which all go a way to improving the match with N-body simulations in the mildly non-linear regime. The effective field theory of large scale structure (EFToLSS) Baumann et al. (2012); Carrasco et al. (2012) has also made great improvements to the range of validity of theoretical predictions albeit with the aid of N-body simulations. This makes it impractical for our general approach.

So far we have discussed theoretical predictions for clustering statistics in real space. All depth information of the universe is confined to redshift space and so next we will discuss how to construct the redshift space power spectrum.

II.2 A non-linear Model for RSD

As mentioned in the introduction, RSD arises from the non-linear mapping between real and redshift space due to contamination of the redshift distances by the contribution from the line-of-sight component of galaxy peculiar velocities. First modelled by Kaiser Kaiser (1987), the anisotropy was described as a squashing effect along the line of sight because of coherent infall velocities at large linear scales.

[TABLE]

where $\mu$ is the cosine of the angle between the line of sight and $k$ while $P_{\delta\delta}(k)$ is the linear matter power spectrum, and $b=\delta_{g}/\delta_{m}$ is the linear density bias parameter, $\delta_{g}$ being the galaxy density contrast 444The use of $\mu$ here should not be confused with the function $\mu(k;a)$ which will always include its arguments.. If we move to smaller scales we find that virialised motion causes a broadening in the velocity distribution. This results in a damping effect. This highly non-linear effect was coined the Fingers of God (FoG) and is usually modelled phenomenologicallyScoccimarro (2004); Percival and White (2009); Cole et al. (1995); Peacock and Dodds (1994); Park et al. (1994); Ballinger et al. (1996); Magira et al. (2000).

We will focus our attention on the TNS model of RSD mentioned in the introduction. This model has proven to be quite successful in reproducing simulation data Nishimichi and Taruya (2011); Taruya et al. (2013); Ishikawa et al. (2014) and is founded on SPT. Further it has enjoyed a lot of success in the context of survey data comparisons with the WiggleZ Dark Energy Survey Blake et al. (2011) and the Vipers survey Pezzotta et al. (2016). For these reasons it is the natural choice for our analysis. We quote the TNS power spectrum as

[TABLE]

where $P_{\delta\delta},P_{\delta\theta}$ and $P_{\theta\theta}$ are all at 1-loop order. The correction terms, A and B are given by

[TABLE]

[TABLE]

where

[TABLE]

The cross bispectrum $B_{\sigma}$ is given by

[TABLE]

Finally $D_{FoG}(fk\mu\sigma_{v})$, where $\sigma_{v}$ is treated as a free parameter quantifying the dispersion in velocities (expressed in RSD displacement units Mpc/$h$), is not treated perturbatively but rather phenomenologically. We choose to take an exponential form $D_{FoG}(fk\mu\sigma_{v})=\exp{(-f^{2}k^{2}\mu^{2}\sigma_{v}^{2})}$ Peacock (1992). Again, this term provides the small scale damping of the power spectrum due to random, small scale, motion.

The main feature of this model is the inclusion of the $A$ and $B$ correction terms which account for higher-order interactions between the density and velocity fields. This gives the model good predictive power at weakly nonlinear scales, as shown by $N$-body comparisons Taruya et al. (2010); Nishimichi and Taruya (2011); Taruya et al. (2013). In GR these terms have been shown to enhance the power spectrum amplitude at the BAO scale and have a non-negligible effect on the acoustic features of the power spectrum Taruya et al. (2010). In modified gravity theories the $B$ term is generally expected to be enhanced because of its linear growth dependance while the $A$ term involves the 2nd order perturbations so it is not obvious how it will change. In Taruya et al. (2014) the authors show that these correction terms can be significantly different if the model of gravity is changed, specifically between GR and the $f(R)$ model by Hu and Sawicki Hu and Sawicki (2007).

With the inclusion of galaxy bias, Eq.(16) gives a prediction for an observable that can be measured from ongoing and upcoming surveys. In its presented form, one can still perform a valid comparison with matter statistics from N-body simulations. This being said, a comment on its realm of validity should also be noted: because it relies on SPT, it also suffers from a restricted realm of validity due to divergences in the PT scheme. Within the context of GR, this expression has been computed up to 2-loop order in the resummed PT scheme Taruya et al. (2013, 2009); Okamura et al. (2011); Crocce and Scoccimarro (2008); Crocce et al. (2012); Taruya et al. (2012) which has a larger range of validity over the basic SPT treatment. The standard SPT treatment still gives us a good working range of scales in the quasi non-linear regime which will be the regime for our basic test of model bias.

In the next section we present template comparisons using different data sets, highlighting when model bias becomes an issue at the level of matter statistics.

III Results

Our main goal here is to highlight the importance of model-specific non-linearities and biases that need to be accounted for in cosmological parameters estimation from galaxy spectroscopic surveys. We focus on the logarithmic growth rate, $f$, as the main parameter of interest. Thus, inspired by the approach presented in Oka et al. (2014), we will fix the amplitude of linear density perturbations, or equivalently $\sigma_{8}$, to the fiducial value. We will be dealing with data for which all cosmological parameters (ie. $\Omega_{m},\Omega_{DE},H_{0}$ etc) are already known a priori, and therefore choose to keep all of them fixed during the analysis. Allowing them all to vary would just decrease the statistical significance of the estimates and fixing them does not introduce systematics. In this way we end up with only two free parameters $\{\sigma_{v},f\}$ where $\sigma_{v}$ is the 1-dimensional velocity dispersion in the TNS model which needs to be fit to data.

In this work we focus on GR and the normal branch of the DGP model (nDGP) with $\Omega_{rc}=0.438$ where $\Omega_{rc}=1/(2H_{0}r_{c})^{2}$ and $r_{c}$ is the cross over scale (See Appendix A). The choice of the nDGP model as our guinea pig is motivated by its phenomenology, enriched by the non-linear screening Vainshtein mechanism. We assume a LCDM background in nDGP so the difference between nDGP and LCDM appears only in the structure growth. $f$ is derived from the linear versions of Eq.4 and Eq.5 (the right hand side as well as $\gamma_{2}$ and $\gamma_{3}$ being set to 0) and because $\Omega_{m}=\Omega_{m}^{fiducial}$, the only free parameter is $\Omega_{rc}$ and so we will opt to parametrise $f$ by $\Omega_{rc}$. Since $\Omega_{rc}>0$, a lower bound for $f$ is also imposed. Otherwise the priors for both $\sigma_{v}$ and $f$ are flat. For clarity, Fig.1 shows the relationship between $f$ and $\Omega_{rc}$ at $z=0.5$ and $z=1$. We note that $f(\Omega_{rc}=0)$ corresponds to the logarithmic growth in LCDM.

To get an estimate of $f$ we consider the multipoles of Eq.16. These are modelled as

[TABLE]

where $\mathcal{P}_{\ell}(\mu)$ denote the Legendre polynomials and $P_{TNS}^{(S)}(k)$ is given by eq.(16). In order to get a robust estimation of higher order multipoles, simulations with large volumes and high mass resolution are required. Since our simulations are limited in size and resolution we will limit our modelling and analysis to the monopole and quadrupole.

In $P_{TNS}^{(S)}(k,\mu)$ we normalise the constituents to the fiducial values of the linear growth $F_{1}$ by using the following factor

[TABLE]

with $\sigma_{{\rm norm}}^{2}(a)$ applied to $P^{{\rm linear}}_{ij}(k,a)$ and $\sigma_{{\rm norm}}^{4}(a)$ applied to $P^{13}_{ij}(k,a),P^{22}_{ij}(k,a),A(k,\mu,a)$ and $B(k,\mu,a)$.

As stated earlier, $f$ and $\sigma_{v}$ are treated as free parameters and we fit the first two TNS multipoles to the data with the aim of recovering our fiducial value for $f$. To do this we perform an Markov Chain Monte Carlo (MCMC) analysis using the following likelihood function

[TABLE]

where $\mbox{Cov}_{\ell,\ell^{\prime}}$ is the covariance matrix between the different multipoles. Expressions for the covariance components can be found in Appendix C of Taruya et al. (2010). For our analysis we do not consider non-Gaussianity in the covariance (See Takahashi et al. (2009) for justification of this treatment) but include the effect of shot-noise assuming the number density of an ideal future survey $\bar{n}_{g}=4\times 10^{-3}h^{3}\mbox{Mpc}^{-3}$. Further, we use linear theory to estimate the covariance matrix components. This approximation has been checked to work well within $k\leq 0.3h$/Mpc for the LCDM simulations used in Taruya et al. (2010). Given this, we have found it sufficient for our purposes to check that the covariance of both the density and velocity divergence fields are the same in both our LCDM and nDGP simulations within the scales of interest. Using the scaled covariance (the so-called decoherence function, see Chodorowski and Ciecielag (2002)) we have estimated that they agree to sub percent levels at $k\leq 0.2h$/Mpc.

We consider two redshifts, $z=1$ and $z=0.5$. This is to give an idea of the trade off between enhanced non-linearity but decreased realm of validity of SPT at lower $z$. For both these redshifts, we assume in the MCMC analysis that the errors are those characteristic for a survey with a volume of $V_{s}=10h^{-3}\mbox{Gpc}^{3}$ which is conservatively smaller than the upcoming Dark Energy Spectroscopic Instrument (DESI) Aghamousa et al. (2016) which is further around 3 times smaller than the Euclid spectroscopic survey Laureijs et al. (2011). Our number density, $\bar{n}_{g}=4\times 10^{-3}h^{3}/\mbox{Mpc}^{3}$ and volume are comparable to that of the BOSS MGS sample Ross et al. (2015) or DESI’s BGS Zhao (2014); Aghamousa et al. (2016).

Since we are only presented with $1\mbox{Gpc}^{3}/h^{3}$ box realisation of a full N-body run, using solely this data set could hamper our analysis through uncertainties connected with cosmic variance that become severe at the box-scale. Observing this we decided to eliminate that potential caveat by using an additional 20 $\times$ $1\mbox{Gpc}^{3}/h^{3}$ box realisations as our data set. We do this by employing a rather inexpensive modified gravity COLA approach (hereafter MG-PICOLA) recently presented in Winther et al. (2017) . MG-PICOLA is based on a parallel COLA implementation (PICOLA) (see Howlett et al. (2015) for details). As we have mentioned MG-PICOLA is relatively computationally inexpensive, but this advantage comes at the price of significantly limited accuracy (when compared to a same resolution N-body run) in the non-linear regime of structure formation. We have carefully performed many tests to ascertain that MG-COLA is sufficiently accurate for the purpose of the analysis presented in this work. The reader is referred to Appendix B for the details of the MG-PICOLA tests we have performed.

III.1 Real and Redshift Space: Theory vs Simulations

Here we provide comparison of our full N-body measurements and MG-PICOLA measurements to the SPT predictions. To begin, we compare the real space N-body auto and cross real space power spectra with the SPT 1-loop predictions. This provides a measure of the non-linearity captured by SPT as well as realm of validity of the PT treatment. We then compare the TNS multipoles to MG-PICOLA measurements.

III.1.1 Real Space Comparisons: N-body

We use cosmological N-body simulations of the nDGP model and LCDM and compare them to SPT predictions to gain an idea of the range of validity of the theoretical modelling. These simulations were run using the AMR code ECOSMOG (Li et al., 2012). The background cosmology is taken from WMAP9 (Hinshaw et al., 2013): $\Omega_{m}=0.281$, $h=0.697$, and $n_{s}=0.971$. The box length is $1024\mbox{Mpc}/h$ with $1024^{3}$ dark matter particles used and a starting redshift of $49$. This design sets the resulting mass resolution at $m_{p}\cong 7.8\times 10^{10}M_{\odot}h^{-1}$ and the Nyquist fluid approximation limit of $k_{Nyq}\cong\pi h$/Mpc. The most refined AMR grid were at level 16, setting a maximal force resolution at $\epsilon=1024/2^{16}=0.015\mbox{Mpc}/h$. The initial conditions were generated using MPGrafic555Available at http://www2.iap.fr/users/pichon/mpgrafic.html and both nDGP and LCDM simulations begin with the same initial seeds. The linear theory power spectrum normalisations was set to be $\sigma_{8}=0.844$. The nDGP simulation uses $\Omega_{rc}=1/4r_{c}^{2}H_{0}^{2}=0.438$. We choose to use a high $\sigma_{8}$ with these specific combination of nDGP parameters to obtain a model that would be characterised by noticeable deviations from the GR fiducial linear growth rate and at the same time would have strong enough small scale non-linearities. In this way we can expect that our nDGP bed-test should contain a signal strong enough to be detected by RSD analysis.

We use the Delaunay Tesselation Field Estimator (DTFE) method implemented in publicly avilable DTFE code by Cautun and van de Weygaert (2011). The DTFE code employs the Delanay Tesselation Field Estimation, a method described in details in Schaap and van de Weygaert (2000); van de Weygaert and Schaap (2009), which assures that the resulting smooth fields have the highest attainable resolution, are volume weighted and have suppressed sampling noise. The fields are then smoothed using top-hat filtering and we proceed to obtain density P_{\delta\delta}(k)=\langle\delta(\mbox{\boldmathk})\delta^{*}(\mbox{\boldmathk})\rangle and velocity divergence P_{\theta\theta}(k)=\langle\theta(\mbox{\boldmathk})\theta^{*}(\mbox{\boldmathk})\rangle power spectra following the method of Li et al. (2013); Hellwing et al. (2013). It is well known that energetic processes connected to the highly non-linear physics of galaxy formation affect the cosmic density field up to scales even of tens of Megaparsecs. Thus our simplistic N-body only approach could introduce some additional scatter and bias in the analysis. However, as recently shown by Hellwing et al. (2016) with a use of the state-of-the-art galaxy formation simulation, the EAGLE suite Schaye et al. (2015), the energetic baryonic processes have a negligible impact on both velocity and density fields on the scales we consider in this study. Hence, we can be assured that our analysis will not be affected by the fact that we ignore baryons completely in our modelling.

Fig.2 shows SPT does very well in modelling non-linearities at $z=1$. For the nDGP simulation we find agreement at $1(3)\%$ level up to $k=0.175(0.18),0.18(0.2)$ and $0.22(0.26)h$/Mpc for $P_{\delta\delta}$, $P_{\delta\theta}$ and $P_{\theta\theta}$ respectively. Given this we consider scales up to $k_{\rm max}=0.195h$/Mpc at $z=1$. We have fit Poisson errors to this data assuming a volume of $1\mbox{Gpc}^{3}/h^{3}$ with a shot noise term of $\bar{n}=3\times 10^{-4}h^{3}/\mbox{Mpc}^{3}$ (See Eq.27 of Zhao (2014) for example).

On the other hand, Fig.3 shows the decline in the accuracy of the SPT approach at later times. Again for the nDGP simulation, now at $z=0.5$, we find agreement at $1(3)\%$ level up to $k=0.12(0.18),0.18(0.2)h$/Mpc for $P_{\delta\delta}$ and $P_{\delta\theta}$ respectively. $P_{\theta\theta}$ is found to be very noisy around the $1$-$3\%$ band within $0.1\leq k\leq 0.2$. Given this, we consider scales up to a $k_{\rm max}=0.147h$/Mpc at $z=0.5$. The stated ranges of applicability of our SPT modelling in nDGP are found to be very similar to the GR simulation comparisons.

In Fig.4 we have plotted the ratio of the real space spectra in nDGP to LCDM simulations for $z=0.5$ and $z=1$ along with the linear predictions as dotted and dashed lines. This figure captures the effects of modified gravity on the real space spectra and growth. We note that the density spectra remain very much the same as we proceed into the quasi non-linear regime. On the other hand we find that the DGP velocity spectra becomes more suppressed as we go to smaller scales. This scale dependence is expected as larger fifth force enhancements to the velocity field can effectively reduce the correlation between particle velocities, an effect that is larger at lower redshift where there is more non-linear structure growth. In Fig.5 we further elucidate this point where we have plotted the ratio of velocity and cross spectra to their linear predictions for both LCDM and nDGP simulations. A first point is that clearly non-linearity becomes very important in the scales considered. Secondly, the nDGP simulation shows an enhanced non-linearity (and hence suppression of velocity correlations) over the LCDM simulation.

These results are consistent with previous results found for other MG simulations such as in Fig.4 and Fig.7 of Winther et al. (2015) and Fig.7 and Fig.8 of Winther et al. (2017)). They found a pattern of a constant boost at the linear scales (reflecting enhanced growth rate) with the scale-dependance beginning to be suppressed closer to non-linear scales. At non-linear scales the 5th force in nDGP theories starts to be screened effectively by the Vainshtein mechanism, recovering the Newtonian value inside most of the virialised structures (see also Falck et al. (2015)). Our analysis is limited to only quasi-linear scales where the complicated scale-dependent patterns of the Vainshtein mechanism are still not well developed. At those scales the MG-physics can still be faithfully captured by the two-parameter TNS modelling.

III.1.2 Redshift Space Comparisons : MG-PICOLA

Next we take a look at the predictive power of the TNS multipoles, providing the realm of validity of the RSD modelling. The MG-PICOLA multipoles are measured using the distant-observer approximation 666i.e we assume that the observer is located at a distance much greater then the boxsize ($r\gg 1024\mbox{Mpc}/h$), and so we treat all the lines-of-sight as parallel to the chosen Cartesian axes of the simulation box. Next, we use an appropriate velocity component ($v_{x},v_{y}$ or $v_{z}$) to disturb the position of a matter particle. and averaged over three line of sight directions. We further average over 20 MG-PICOLA simulations each of $1\mbox{Gpc}^{3}/h^{3}$ volume thus ignoring the mode covariance at and above box-size scales. This should correspond to an ideal survey with a resulting volume of $20\mbox{Gpc}^{3}/h^{3}$.

Fig.6 shows the monopole and quadrupole predictions at $z=1$ for three different values of $\sigma_{v}$ where we have fit using 32 $k$-bins up to $k_{\rm max}=0.195h$/Mpc. The reduced $\chi^{2}_{\rm red}=-2\ln{\mathcal{L}}/{\rm N_{dof}}$ is shown in brackets, where $N_{dof}$ is the total degrees of freedom at that $k_{\rm max}$ which equals twice the number of bins minus the number of parameters. In this case $N_{dof}=62$($=32\times 2-2$), since we have 2 parameters. Similarly, Fig.7 shows the same results at $z=0.5$ where we have fit using 24 $k$-bins ($N_{dof}=46$) up to $k_{\rm max}=0.147h$/Mpc as dictated by the real space power spectra comparisons.

We find that the fit of PT is comparable between GR and nDGP simulations at both redshifts, although GR does slightly better at small scales at $z=0.5$. The quadrupole also does worse than the monopole over both models at $z=0.5$ which we can attribute to the increased dependency on the velocity auto power spectrum for which the theoretical template does worse in reproducing (See Fig.3).

As was done with the power spectra, we identify the scale at which deviations from theory are within the $1(3)\%$ region. We do this using the values of $\sigma_{v}$ with the lowest $\chi^{2}$ shown in Fig.6 and Fig.7. This will be used to set a maximum Fourier mode for our statistical analysis. For the nDGP simulation at $z=1$ we find agreement at $1(3)\%$ level up to $k=0.24(0.25)h$/Mpc for $P_{0}$. $P_{2}$ is significantly noisier around the $1\%$ deviation region but matches PT up to $k=0.16h$/Mpc within $3\%$. Similarly, for $z=0.5$ we find a theory-data agreement of $1(3)\%$ up to $k=0.147(0.2)h$/Mpc for $P_{0}$ while a much worse $k=0.09(0.1)h$/Mpc for $P_{2}$. Since $P_{2}$’s contribution to Eq.23 is significantly smaller than $P_{0}$ because of their respective errors, we decide to use $k_{\rm max}=0.147h$/Mpc at $z=0.5$ and $k_{\rm max}=0.195h$/Mpc at $z=1$ for the MCMC analysis, despite the poor quadrupole fit. In the case of GR the range of validity is found to be similar.

By having assessed and compared the range of validities of SPT for both a LCDM and nDGP cosmologies, we have gained a handle on how enhanced dynamics due to fifth force interactions degrade SPT’s performance. Our comparisons indicate that the non-linearity generated by $\Omega_{rc}=0.438$ is mild enough not to significantly effect the range of validity. We caution however, this may not be the case for other models of gravity characterised by a larger degree of deviation from GR dynamics. Next we will test the capabilities of the theoretical templates when matching to the MG-PICOLA multipoles.

III.2 Template Performance

We begin by explicitly defining the theoretical templates. The GR template here means that we use the TNS model for the RSD where we set $\mu=1,\gamma_{2}=\gamma_{3}=0$ but treat $f$ as a free parameter parametrised by $\Omega_{rc}$. On the other hand, for the DGP template, we use $\mu,\gamma_{2},\gamma_{3}$ as given in Appendix A to model the non-linearity in the DGP model properly. It is a common practice to use a GR template to estimate $f$ (eg. Beutler et al. (2016)).

Our main intention in this work is to robustly assess whether using GR-only based RSD modelling (in other words, a model ignorant of any possible deviations from the GR picture of structure formation) on a non-GR universe is an accurate enough procedure to allow for a recovery of the fiducial value of the growth rate parameter $f$ from the data. Note that the DGP template encompasses a pure-GR scenario (by setting $\Omega_{rc}=0$) and so we should observe no biasing in using the DGP template for the GR data. Further, we want to determine on what scales and what amplitudes biasing becomes significant. As we fix the linear growth of density perturbations, $F_{1}$, both theoretical templates are equivalent on linear scales but as we include increasingly non-linear scales in the analysis the templates deviate producing the bias. This bias may be masked with the freedom of the TNS damping parameter $\sigma_{v}$, although the model specific non-linearities act beyond pure damping of small scale power and so $\sigma_{v}$ cannot perfectly substitute for incorrect perturbation modelling. In other words, model bias should still be expected at some level.

Our analysis follows the techniques described in Lewis and Bridle (2002), with the MCMC algorithm walking in $\{\sigma_{v},\Omega_{rc}\}$ parameter space. In order to measure when non-linearities become an issue in terms of theoretical modelling we complete the analysis for varying values of $k_{\rm max}$ up to the upper bounds found in the previous subsection. By going to higher $k$-modes we suppress the statistical errors. Fig.8 shows the results at $z=0.5$. On the left hand side we present the $1\sigma$(68$\%$ C.L) and $2\sigma$(95$\%$ C.L) contours for the DGP template (green) and GR template (blue) for a matching up to $k_{\rm max}=0.147h$/Mpc using 24 bins. The marginalised statistics are shown in the side panels and the fiducial value, $\Omega_{rc}=0.438$, is marked as a dashed line. We note a slight offset of the best fit values of $\Omega_{rc}$, although the fiducial value remains in the $1\sigma$ region for both templates. There is also an offset in the best fit $\sigma_{v}$ which can be interpreted as the GR template’s use of this parameter to compensate for non-linear effects in the DGP-PICOLA simulations. On the right of Fig.8, we see the marginalised best fit values for $f(\Omega_{rc})$ with their $2\sigma$ errors for varying scale inclusion ($k_{\rm max}=0.110,0.135,0.147h$/Mpc with GR template’s values shifted slightly for better visualisation). Note that as we push to higher $k_{\max}$ the errors become smaller, as expected. We here remind the reader that there is a lower bound on $f$ imposed by $\Omega_{rc}\geq 0$. As can be seen in both plots, the GR template comfortably accommodates the data within $68\%$ (for $k_{\rm max}=0.147h$/Mpc) within the templates’ validity regime.

Fig.9 shows the results at $z=1$ this time matching up to a $k_{\rm max}=0.195h$/Mpc. It is in this case that we find the GR template struggling to capture the full shape of the multipoles and in the left hand plot we see that the fiducial parameter lies outside the $1\sigma$ region and just within the $2\sigma$ region. On the other hand, the DGP template is centred around the fiducial value. On the right hand side of Fig.9 we show the results for an analysis similar to that of the left hand side but for a larger survey with a volume of $V_{s}=20\mbox{Gpc}^{3}/h^{3}$, which is the estimated volume of the Dark Energy Spectroscopic Instrument (DESI) Aghamousa et al. (2016). In this case the GR template fails to capture the fiducial at the $2\sigma$ level implying that inadequate theoretical modelling for such a large survey could introduce a significant error in parameter inference.

Fig.10 shows the improvement in constraints as we increase $k_{\rm max}$. We have included an annotation with the same analysis done for a larger survey with a volume of $V_{s}=20\mbox{Gpc}^{3}/h^{3}$ fitting up to $k_{\rm max}=0.195h$/Mpc. We see that the DGP template does consistently well in reproducing the fiducial while the GR begins to fail at around $k_{\rm max}=0.171h$/Mpc.

As a consistency check of our analysis, we make use of the GR simulations at $z=1$. Again, by assuming a survey volume of $V=10\mbox{Gpc}^{3}/h^{3}$, we repeated the analysis and obtained constraints on $\Omega_{rc}$ using the DGP template, as well as a constraint on $f$ using the GR template without the induced prior on $f$ coming from our parametrisation using $\Omega_{rc}$. Fig.11 shows the results. Both contours recover the fiducial parameters within $1\sigma$.

As seen from the left panel of Fig.11, the nDGP model with $\Omega_{rc}=0.438$ can be excluded with high significance ($\gg 2\sigma$) by a survey with a volume $V_{s}=10\mbox{Gpc}^{3}/h^{3}$ fitting up to $k_{max}=0.171h$/Mpc, if our universe is described by GR. To quantify this, we computed the following quantity using the N-body measurements

[TABLE]

where $\mbox{Cov}^{-1}_{\ell,\ell}$ is the covariance matrix between the multipoles assuming an ideal survey of $V_{s}=10\mbox{Gpc}^{3}/h^{3}$. We also computed the same quantity using SPT. Fig.12 shows the results up to $k_{\rm max}=0.2h$/Mpc at $z=0.5$ and $z=1$. The results clearly shows that our ability to distinguish between LCDM and nDGP increases with $k_{max}$. Also we find that the $\chi^{2}$ from SPT is very similar to those obtained from simulation. This indicates that the TNS model is capable of providing the same level of significance of deviations from LCDM as the full non-linear treatment, making it a good estimator of structure formation in this range of scales.

We have compiled the template results in Table.1, Table.2 and Table.3.

III.3 An Ideal Survey: SPT Mock Data

Being model dependant, we can expect that both the scales affected as well as the magnitude of theoretical bias will in general depend on the specific phenomenology of a given gravity model. Therefore, if one wants to precisely estimate the importance of such theoretical bias for a given set of real galaxy spectroscopic data, one would need to run N-body simulations for each model under consideration and then perform a similar analysis as in the previous section. This is obviously not practical and in this section we provide a means of getting a first indication of whether or not model bias is an issue for a given model, in other words, whether or not the model can be safely encompassed by the GR template within the relevant range of scales.

We proceed as follows. First, multipole data is produced for a given model of gravity using SPT up to some valid $k_{\rm max}$. Then the covariance matrix for this data is computed as was done for the N-body data using the parameters of an ideal survey. Finally, this data is given Gaussian errors using the covariance matrix. This provides an easily produced, idealistic, simulated mock data set which can be done for any model of gravity described by the framework discussed in the first section. A statistical analysis as done in the previous section can then be performed on this data. Here we do this for the nDGP model with the same fiducial model parameter previously used, $\Omega_{rc}=0.438$ but with $\sigma_{8}=0.87$. All other cosmological parameters are the same as the nDGP N-body simulation. We choose a fiducial $\sigma_{v}=5.5\mbox{Mpc}/h$ and use the ideal survey parameters $V_{s}=10\mbox{Gpc}^{3}/h^{3}$ and $\bar{n}_{g}=4\times 10^{-3}h^{3}/\mbox{Mpc}^{3}$. Only $z=1$ is considered in this section and because SPT underestimates the non-linear effects we extend our statistical analysis to $k_{\rm max}=0.2h$/Mpc to include more non-linear scales.

In addition to including more scales in the analysis one can also ask for which models does the model-bias systematic become important, or how much enhanced dynamics induced by modified gravity is needed to detect a significant deviation from the GR template. One could investigate this by creating mock data for larger values of $\Omega_{rc}$ although this becomes quite unrealistic. In fact the value of $\Omega_{rc}=0.438$ is already ruled out by BOSS LOWZ and CMASS data to within $2\sigma$ (See Barreira et al. (2016)), with the authors placing an upper bound of 0.36. What we choose to do instead is rescale the non-linear mode mixing which governs the change of the scalar field’s non-linear derivative interactions - the source of screening. The rescaling is done by introducing the parameter $\alpha$. We will scale $\gamma_{2}$ by $\alpha$ and $\gamma_{3}$ by $\alpha^{2}$ in the Euler equation’s non-linear source term S(\mbox{\boldmathk}) (Eq.8). This is equivalent to enhancing the nDGP second order mode mixing term by $\alpha$ (See Koyama et al. (2009) for details). Then, by setting $\alpha=1$ we obtain the usual nDGP model but for values larger than unity the model changes to one with enhanced non-linearities. By tuning $\alpha$ we will be able to test the capabilities of the GR template to cover model non-linearities given an idealistic survey. Note this approach is not meant to represent a realistic or viable model but rather to be illustrative.

Fig.13 show the results of the MCMC analysis for $\alpha=1$. The left hand plot shows the $1\sigma$ and $2\sigma$ contours for a matching done up to $k_{\rm max}=0.2h$/Mpc. We see that in this case both templates well recover the fiducial $\Omega_{rc}$ although the GR template fails to recover the fiducial $\sigma_{v}$ as expected. The right hand plot shows that the templates are comparable in their best fit value for $f$ as well as their $2\sigma$ constraints at all values of $k_{\rm max}$. We direct the reader again to Fig.9 which shows both best fits for $\Omega_{rc}$ when comparing to the N-body data. In contrast to results presented in Fig.13, we find the values as being substantially smaller than the fiducial value. As mentioned this may reflect a combined effect of insufficient statistics due to box size and also a use of too high a $k_{\rm max}$, where the PT-based templates begin to fail. The lower-than-fiducial preferences of both templates in Fig.9 would be consistent with a break down of SPT. This is because the loop corrections begin to over predict the power spectrum as SPT breaks down and so to compensate theory prefers a smaller value of $\Omega_{rc}$. Further, the breakdown of SPT for GR cosmologies at $z=1$ has been shown to be around $k_{\rm max}=0.147h$/Mpc Taruya et al. (2010); Carlson et al. (2009). This is curious as it would be noticed in Fig.2 and Fig.6. Fig.11 suggests that if the offset is indeed due to a theoretical breakdown, it is not to a sufficient extent to push the fiducial value out of the $1\sigma$ region. In any case, we expect more light to be shed on this issue with more data.

Finally, Fig.14 illustrates that if we set $\alpha=15$ the GR template completely fails to recover the fiducial value even at large (more linear) scales, with $\Omega_{rc}=0.438$ still lying outside the $2\sigma$ errors for $k_{\rm max}=0.1h$/Mpc (right hand side plot). This is an extreme case with the non-linearities becoming far too important at linear scales. Fig.15 shows the deviation between the templates for $\sigma_{v}$ and $\Omega_{rc}$ set to fiducial values, clearly indicating when the enhancement enters the linear regime and away from $\sigma_{v}$’s ability to suppress it. With $\alpha=13$ we note 1-2$\%$ deviations at scales as large as $k=0.05h$/Mpc. The bias starts to become important out of $\sigma_{v}$’s reach at $\alpha=10$ (See Fig.6 for an indication of $\sigma_{v}$’s impact range).

IV Summary and Conclusions

The work presented here was motivated by the question of whether constraints of the growth rate derived with a GR PT-template give an unbiased measurement of growth for the case of a universe described by a modified gravity. This question is very relevant in the context of the next stage of cosmological surveys. We have here provided a first level analysis of the theoretical model bias systematic and have presented a method for quick assessment of the model specific significance of this systematic.

Firstly, we compared nDGP MG-PICOLA simulation data with the GR and nDGP theoretical predictions for the redshift space power spectrum. We used the TNS model of RSD which has been validated against both GR and modified gravity simulations Taruya et al. (2010, 2014). This was done at the level of dark matter clustering and only the first two multipoles were considered. Idealised future survey parameters were adopted in the analysis. We found out that the small-scale velocity damping term $\sigma_{v}$ included in the TNS model provides a flexibility through which the template can, to some extent, accommodate the enhanced small-scale clustering of the nDGP model. This was clearly indicated by higher values of $\sigma_{v}$ attained by the GR-template fit.

Both templates perform well in recovering the simulation’s fiducial parameter at low redshift. We point out that the real space analysis done in Barreira et al. (2016) concluded that no nDGP model bias is evident at redshifts up to $z=0.57$ which is consistent with our results at $z=0.5$. That being said, a full comparison of our results is difficult as there are a number of differences between their analysis and the one done here. In particular the use of different RSD models, the modelling of survey errors and their inclusion of galaxy bias which provides more fitting freedom to the GR template. We do find that at high redshift the GR prediction becomes increasingly biased and the difference between the two templates is greater. Using $V_{s}=20\mbox{Gpc}^{3}/h^{3}$, which will be realisable with stage IV surveys, we find systematically biased estimates of the GR template, with it failing to recover the fiducial parameter to within $2\sigma$ at $z=1$. This apparent bias might be due to specific limitations intrinsic to the SPT approach. The current analysis is left suggestive with robustness sought in additional theoretical modelling (for example including galaxy bias) or fuller treatment of non-linear scales such as using the EFT approach Baumann et al. (2012); Carrasco et al. (2012). We leave for future work such analysis tailored for the detailed specifications of future generation galaxy surveys. Nonetheless, what is clear is that for the $V_{s}=20\mbox{Gpc}^{3}/h^{3}$ case the $1\sigma$ regions of both templates do not intersect which implies that the template’s predictions are inconsistent at that level. Again, the inclusion of galaxy bias may relieve this.

In our second analysis we created mock data from SPT predictions by adding Gaussian noise generated using the errors derived from an idealised survey. Two data sets were created at $z=1$ using varying levels of model dependent non-linearity. This was done to simulate modified gravity models which have an enhanced non-linear source term. We find that by increasing the non-linear contributions to the higher order density and velocity perturbations, the GR template fails to recover the fiducial $\Omega_{rc}$, with the model bias being unimportant up to a non-linear contribution of around 10 times the base value. Above this the GR predictions become very biased and at 15 times the base non-linearity the GR template fails to recover the fiducial value even at scales of $k=0.1h$/Mpc. This exercise provides an indication on what scales and at what level of enhanced small-scale clustering a modified gravity model has to be consistently treated in RSD modelling in order to avoid significant theoretical biases that otherwise would diminish the desired accuracy of growth rate estimates. On this note, the creation of mock data can be done for any model of gravity within the framework discussed in Bose and Koyama (2016) giving an avenue for assessing the importance of theoretical model bias in growth rate estimation from a given data set. The data quality of stage IV surveys indicates that this test will be important and is essential if we wish to put trusted constraints on modified gravity parameter space.

In future analyses, the inclusion of higher multipoles would reduce the statistical errors and so further highlight model bias. On the other hand, we have only dealt with dark matter statistics and in a real survey galaxy bias needs to be included. This increases the degrees of freedom of the theoretical templates and gives the GR prediction more ability to cover the model bias systematic. The effect galaxy bias inclusion is left to a future work. We conclude that if the model exhibits sufficient non-linearity then model bias becomes an issue for modified gravity constraints using upcoming spectroscopic surveys such as DESI. This effect is more prominent at higher redshifts where we benefit from an increased spatial range of applicability of our theoretical template. Further, the improvement of theoretical predictions, such as through the EFToLSS Baumann et al. (2012); Carrasco et al. (2012), will more loudly call for proper theoretical treatment of modified gravity models before attempts to obtain robust constraints are made.

Finally, we comment on the analysis contained in this work which validate the use of MG-PICOLA Winther et al. (2017) to check theoretical predictions. In Appendix B we perform a number of tests which show that the COLA approach to modified gravity consistently gives the correct non-linear predictions within SPT’s realm of validity, up to an additional damping. This can be well captured by the TNS models extra degree of freedom, and with loose constraints on $\sigma_{v}$ one may even push to smaller scales, with EFT approaches say, and perform similar data-theory analyses.

Acknowledgments

We thank the anonymous referee for useful and insightful comments and suggestions which benefited the scientific quality of this paper. BB is supported by the University of Portsmouth. KK, HAW and WAH 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. WAH also acknowledges the support from the Polish National Science Center under contract $\#$UMO- 2012/07/D/ST9/02785.

Appendix A Appendix A: DGP gravity

In the DGP model of gravity Dvali et al. (2000) we live on a 4 dimensional brane embedded in 5 dimensional Minkowski spacetime, giving this theory a crossover scale $r_{c}$, which is the ratio between the 5D Newton gravitational constant and the 4D Newton gravitational constant. $r_{c}$ is the only free parameter of the theory and we parametrise it as $\Omega_{rc}=1/(4r_{c}^{2}H_{0}^{2})$, $H_{0}$ being the Hubble parameter today. The modified Friedman equation in this theory is given by

[TABLE]

where $\epsilon=\pm 1$. The solution for $\epsilon=1$ is known to be ghostly and so we consider what is called the normal branch with $\epsilon=-1$ (nDGP). In this branch acceleration is achieved through a dark energy constant as in GR. We also impose a background history following LCDM, done by tuning the dark energy equation of state. The non-linear interaction terms for this theory are found to be (for details on the Poisson equation form in this theory see Koyama et al. (2009) for example)

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

$\mu_{i,j}$ is the cosine of the angle between \mbox{\boldmathk}_{i} and \mbox{\boldmathk}_{j}, \mbox{\boldmathk}_{ij}=\mbox{\boldmathk}_{i}+\mbox{\boldmathk}_{j} and $H^{\prime}=dH/da$.

Appendix B Appendix B: N-body vs MG-PICOLA

In this Appendix we validate our use of MG-PICOLA in Sec.3B. We will refer to a single PICOLA simulation with the same initial conditions as N-body as COLA1 and we will refer to the averaged measurements from the 20 PICOLA runs as COLA20. We perform a number of tests listed below.

We compare the real space spectra from the full N-body simulation to COLA1. This serves to test the accuracy of PICOLA’s evolution of structure. Fig.16 shows that the COLA method reproduces the full non-linear real space spectra to within $2\%$ up to $k=0.2h$/Mpc at $z=0.5$ and $z=1$. This asserts that the full non-linear dynamics and evolution is sufficiently captured by MG-PICOLA at the scales of interest. 2. 2.

We then compare the multipoles from the full N-body simulation to COLA1 and test for fiducial parameter recovery using the nDGP template for both measurements. We use a direct FFTW estimation of the multipoles from the N-body data. Fig.17 shows the multipole comparisons. The redshift space multipoles show less damping in MG-PICOLA simulations compared to the full N-body measurements due to less non-linear structures in these simulations, which give less fingers of god effects. Since the TNS model has the free parameter $\sigma_{v}$, that models this non-linear effect, the reduced damping can be accounted for by a smaller value of $\sigma_{v}$. Fig.18 shows the results from an MCMC analysis using multipole data sets from COLA1 and N-body. It shows that we get a very good match in the marginalised posterior distribution of our parameter of interest $\Omega_{rc}$ and the contours are only shifted along $\sigma_{v}$. We use $1\mbox{Gpc}^{3}/h^{3}$ survey errors in accordance with the size of the simulations. 3. 3.

Finally, we compare the redshift space multipoles from COLA1 to COLA20. This checks to see if the initial phases used in the full N-body simulation are outside the variance of the 20 runs. Fig.19 shows that COLA1 is within the variance of the 20 runs and nothing is unusual about the N-body’s initial seeds. We also note that the initial condition for N-body was generated by MPGrafic which uses the Zeldovich approximation while MG-PICOLA employs 2nd order Lagrangian Perturbation theory to generate initial conditions.