Beyond consistency test of gravity with redshift-space distortions at quasi-linear scales

# Beyond consistency test of gravity with redshift-space distortions at quasi-linear scales

Atsushi Taruya Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan Research Center for the Early Universe, Graduate School of Science, The University of Tokyo, Bunkyo-ku, Tokyo 113-0033, Japan Kavli Institute for the Physics and Mathematics of the Universe, Todai Institutes for Advanced Study, the University of Tokyo, Kashiwa, Chiba 277-8583, Japan (Kavli IPMU, WPI)    Kazuya Koyama Institute of Cosmology & Gravitation, University of Portsmouth, Dennis Sciama Building, Portsmouth, PO1 3FX, United Kingdom    Takashi Hiramatsu Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan    Akira Oka Department of Physics, The University of Tokyo, Bunkyo-ku, Tokyo 113-0033, Japan
July 11, 2019
###### Abstract

Redshift-space distortions (RSD) offer an attractive method to measure the growth of cosmic structure on large scales, and combining with the measurement of the cosmic expansion history, it can be used as cosmological tests of gravity. With the advent of future galaxy redshift surveys aiming at precisely measuring the RSD, an accurate modeling of RSD going beyond linear theory is a critical issue in order to detect or disprove small deviations from general relativity (GR). While several improved models of RSD have been recently proposed based on the perturbation theory (PT), the framework of these models heavily relies on GR. Here, we put forward a new PT prescription for RSD in general modified gravity models. As a specific application, we present theoretical predictions of the redshift-space power spectra in gravity model, and compare them with -body simulations. Using the PT template that takes into account the effects of both modifications of gravity and RSD properly, we successfully recover the fiducial model parameter in -body simulations in an unbiased way. On the other hand, we found it difficult to detect the scale dependence of the growth rate in a model-independent way based on GR templates.

cosmology, large-scale structure
###### pacs:
98.80.-k, 98.62.Py, 98.65.-r
preprint: YITP-13-93

## I Introduction

Redshift-space distortions (RSD) of galaxy clustering, which appear as systematic effects in determining the redshift of each galaxy via spectroscopic measurements and manifestly break statistical isotropy Hamilton (1997); Peebles (Princeton University Press, 1980), are now recognized as a sensitive probe of the growth of structure. Combining the distance measurement of galaxies using the baryon acoustic oscillations as standard ruler (e.g., Seo and Eisenstein (2003); Blake and Glazebrook (2003); Matsubara (2004); Glazebrook and Blake (2005)), RSD offers a unique opportunity to test the theory of gravity on cosmological scales (e.g., 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. (2013)) and help us obtain a deeper understanding of the current accelerating expansion of the Universe. This is indeed one of the main goal of on-going and up-coming galaxy surveys such as the Baryon Oscillation Spectroscopic Survey of Sloan Digital Sky Survey III, the WiggleZ survey222wigglez.swin.edu.au, the Subaru Measurement of Imaging and Redshifts, the Dark Energy Survey444www.darkenergysurvey.org, the BigBOSS project555bigboss.lbl.gov/index.html, and the ESA/Euclid survey666www.euclid-ec.org, which will provide precision measurements of the power spectrum (or correlation function). The late-time cosmic acceleration, first discovered by the observations of distant type Ia supernovae Perlmutter et al. (1999); Riess et al. (1998), may be the result of a dark energy which can be realized in the presence of dynamical scalar field, or it may indicate the breakdown of general relativity (GR) on cosmological scales. The latter case requires a consistent model of gravity that explains the accelerating expansion on large scales with the modification of gravity, while neatly evading the stringent constraints on the deviation from GR at solar system scales (e.g., Khoury and Weltman (2004); Deffayet et al. (2002); Hu and Sawicki (2007); Starobinsky (2007)). In this respect, the large-scale structure offers the best opportunity to distinguish between modified gravity and dark energy models in GR, and the measurement of RSD is a powerful tool to probe gravity.

Given the high-precision measurements of RSD in the near future, accurate theoretical templates of the redshift-space power spectrum or correlation function is highly demanded in order to detect a small deviation of gravity from GR. This is indeed now active research subject, and there are many studies to accurately model RSD. The RSD measurement is basically made at the scales close to the linear regime of gravitational evolution, but the nonlinearity arising both from the gravity and the RSD is known to play a crucial role. Moreover, due to the non-Gaussian nature of RSD Scoccimarro (2004), the applicable range of linear theory prediction is fairly narrower than that in real space. Thus, beyond the linear scales, a sophisticated treatment is required for reliable theoretical predictions with a wider applicable range. Recently, based on the perturbation theory of large-scale structure, several improved models of RSD have been proposed Taruya et al. (2010); Nishimichi and Taruya (2011); Matsubara (2008); Reid and White (2009); Carlson et al. (2013); Seljak and McDonald (2011); Wang et al. (2013); Vlah et al. (2012, 2013). These models properly account for the non-Gaussian nature of RSD, and are tested against -body simulations, successfully describing redshift-space power spectrum and/or correlation function at weakly nonlinear regime. Applying these models to real observations, constraints on the growth of structure have been also obtained (e.g., Reid et al. (2012); Blake et al. (2011)).

However, it should be noted that the proposed models of RSD have been so far tested only in the case of GR. Further, beyond linear theory, the template of RSD is computed with the perturbation theory under the assumption that gravity is described by GR. Thus, the observational constraints derived from the PT-based template can be only used as a consistency test with GR, and a care must be taken in addressing the constraint on a specific model of modified gravity.

The aim of the present paper is to examine these issues based on an improved model of RSD developed by Ref. Taruya et al. (2010). The power spectrum expression of this model is similar to the one proposed by Ref. Scoccimarro (2004) and the so-called streaming model frequently used in the literature (e.g., Peebles (Princeton University Press, 1980); Fisher (1995); Hatton and Cole (1998); Scoccimarro (2004); Jennings et al. (2011)), but it includes two important PT corrections as a result of the low- expansion. Although the model also includes a phenomenological term to account for the Finger-of-God damping arising from the small-scale physics, combining the recently developed resummed PT, it successfully describes not only the matter but also the halo power spectra in -body simulations Taruya et al. (2013); Nishimichi and Taruya (2011); Ishikawa et al. (2013). It is shown that the model can be used as a theoretical template to simultaneously constrain the parameters associated with the cosmic expansion and the structure growth in an unbiased manner, and applying it to the Luminous Red Galaxy sample of Sloan Digital Sky Survey Data Release 7, a robust contraint is obtained Oka et al. (2013).

Here, extending these previous works in GR, we put forward a prescription to compute the redshift-space power spectrum in modified gravity models. As a specific example, we explicitly compute the redshift-space power spectrum in gravity model, as one of the representative modified theory of gravity Hu and Sawicki (2007); Starobinsky (2007). The theoretical prediction based on the standard PT calculation is compared with the results of -body simulations, and a good agreement is found. Then, we will discuss the potential impact of the precision modeling of RSD on the model-independent test of GR and/or constraint on modified gravity models. We will show that a tight and unbiased constraint on modified gravity models is achieved only with an improved PT model of RSD in which the effect of modified gravity is properly taken into account in the PT calculation. With the improved PT template, testing GR will be made possible in a model-independent way, but we argue that a quantitative characterization of the small deviation from GR generally requires a prior knowledge of modified gravity models.

The paper is organized as follows. In Sec. II, we begin by briefly reviewing the model of RSD proposed by Ref. Taruya et al. (2010). Employing the standard PT calculation, we then give a prescription on how to compute the redshift-space power spectrum in modified gravity models. In Appendix A, we summarize the basic formalism of the standard PT in a general context of modified gravity models, and explicitly give expressions for the second-order PT kernels used to compute the higher-order corrections of RSD. In Sec. III, as one of the representative models of modified gravity, we consider the gravity model, and quantitatively compare the PT predictions in redshift space with -body simulations. Based on this, in Sec. IV, a potential impact of the precision PT model of RSD is discussed, particularly focusing on a precision constraint on the model parameter of modified gravity, and model-independent analysis of detecting or characterizing a small deviation from GR. Finally, Sec. V is devoted to the summary and conclusion.

## Ii Modeling redshift-space power spectrum from perturbation theory

### ii.1 An improved model of RSD

We begin by writing the exact expression for redshift-space power spectrum. Let us denote the density and velocity fields by and . Owing to the distant-observer approximation, which is usually valid for the observation of distant galaxies of our interest, one can write (e.g., Scoccimarro (2004); Bernardeau et al. (2002); Taruya et al. (2010))

 P(S)(\boldmathk)=∫d3\boldmathxei\boldmathk⋅\boldmathx⟨eikμΔuz ×{δ(\boldmathr)−∇zuz(\boldmathr)}{δ(\boldmathr′)−∇zuz(\boldmathr′)}⟩, (1)

where denotes the separation in real space and indicates an ensemble average. In the above expression, the -axis is taken as an observer’s line-of-sight direction, and we define the directional cosine by . Further, we defined , and for the line-of-sight component of the velocity field. Note that the above expression has been derived without invoking the dynamical information for velocity and density fields, i.e., the Euler equation and/or continuity equations. Thus, Eq. (1) does hold even in modified gravity models.

Eq. (1) can be re-expressed in terms of the cumulants. Then, the term in the bracket is factorized into two terms, each of which includes the exponential factor (e.g., see Eq.(6) of Ref. Taruya et al. (2010) for explicit expression). Among these, the overall factor, expressed as with being the cumulant, is responsible for the suppression of the power spectrum arising mostly from the virialized random and coherent motion on small scales. The effect of this is known to be partly non-perturbative, and seems difficult to treat petrubatively. Since it has been shown to mainly change the broadband shape of the power spectrum, we may phenomenologically characterize it with a general functional form with being a scale-independent constant. On the other hand, the remaining factor includes the term leading to the Kaiser effect in the linear regime Kaiser (1987); Hamilton (1992), and is likely to affect the structure of power spectrum on large-scales. Although there also appears the exponential factor in each term of this factor, these contributions should be small as long as we consider the large scales, and the perturbative treatment may be applied.

With the proposition given above, Ref. Taruya et al. (2010) applied the low- expansion, keeping the overall prefactor as general functional form . The resultant power spectrum expression at one-loop order becomes

 P(S)(k,μ)=DFoG[kμσv] ×{PKaiser(k,μ)+A(k,μ)+B(k,μ)}, (2)

which we hereafter call TNS model. Owing to the single-stream approximation in which the dynamics of large-scale structure is described by the density and velocity divergence , the quantities , , and are explicitly written as

 PKaiser(k,μ)=Pδδ(k)−2μ2Pδθ(k)+μ4Pθθ(k), (3) A(k,μ)=−kμ∫d3\boldmathp(2π)3pzp2 ×{Bσ(\boldmathp,% \boldmathk−\boldmathp,−\boldmathk)−Bσ(% \boldmathp,\boldmathk,−\boldmathk−\boldmathp)}, (4) B(k,μ)=(kμ)2∫d3\boldmathp(2π)3F(\boldmathp)F(\boldmathk−\boldmathp); (5) F(\boldmathp)=pzp2{Pδθ(p)−p2zp2Pθθ(p)},

where , , and respectively denote auto-power spectra of density and velocity divergence, and their cross power spectrum. The function is cross bispectrum defined by

 ⟨θ(\boldmathk1){δ(% \boldmathk2)−k22zk22θ(\boldmathk2)}{δ(\boldmathk3)−k23zk23θ(\boldmathk3)}⟩ =(2π)3δD(\boldmathk1+% \boldmathk2+\boldmathk3)Bσ(\boldmathk1,\boldmathk2,\boldmathk3). (6)

Note that in deriving Eq. (2), we do not assume any gravity model. Although the expression (2) has been originally derived based on the consideration in GR, as long as the deviation from GR is small, Eq. (2) can apply to any model of modified gravity.

As we mentioned in Sec. I, the main characteristic of the model given in Eq. (2) is the two additional terms and , which represent the higher-order coupling between velocity and density fields. It has been shown in previous studies in GR that these two terms enhance the power spectrum amplitude over the scales of baryon acoustic oscillations, and moderately but notably change the acoustic structure imprinted in the power spectrum Taruya et al. (2010). As a result, the model (2) successfully describes both the matter and halo power spectra of -body simulations at weakly nonlinear scales Taruya et al. (2010); Nishimichi and Taruya (2011); Taruya et al. (2013). These features are expected to hold qualitatively even in the modified theory of gravity, but quantitative aspect of the RSD would generally differ from that of GR, which we will study in detail.

### ii.2 Perturbation theory treatment

To compute the redshift-space power spectrum beyond linear theory, we apply the PT treatment of gravitational evolution, and calculate each term in Eq. (2) in the quasilinear regime. While the power spectrum calculation in the case of GR has been made possible with resummed PT scheme up to the two-loop order (e.g., Taruya et al. (2013) in redshift space, and Taruya et al. (2009); Okamura et al. (2011); Crocce and Scoccimarro (2008); Crocce et al. (2012); Taruya et al. (2012) in real space) and the applicable range of the PT prediction has become wider, we here work with the standard PT calculation at one-loop order for the predictions in modified theory of gravity. Although the standard PT treatment in GR is known to produce ill-behaved PT expansion that lacks good convergence properties (e.g., Crocce and Scoccimarro (2006); Carlson et al. (2009); Taruya et al. (2009); Bernardeau et al. (2012)), using the standard PT as theoretical template, we can still get a fruitful cosmological constraint at the quasilinear scales (e.g., Saito et al. (2011); Zhao et al. (2013)). In the present paper, we use the standard PT formalism developed by Ref. Koyama et al. (2009), which is suited to deal with a wide class of modified gravity models. In what follows, we separately give a prescription on how to compute the power spectrum corrections in Eq. (2).

#### ii.2.1 Non-linear Kaiser term

The term in Eq. (2) is the leading-order contribution to the redshift-space power spectrum. In the large-scale limit where the linear theory prediction is safely applied, we have , and Eq. (3) is reduced to the Kaiser formula, Kaiser (1987), where is the linear growth rate defined by with being the linear growth factor. Beyond the linear theory, a simple relation between density and velocity divergence fields no longer hold, and we need to separately evaluate the three power spectra, , and , especially in the modified gravity models Scoccimarro (2004).

In contrast to the GR, one crucial point in the modified gravity models is that a new scalar degree of freedom, sometimes referred to as the scalaron, arises and modifies the force law. In the presence of an extra scalar field, even though the conservation law of energy momentum tensor remains unchanged, the Poisson equation is inevitably modified and is coupled to the field equation for scalaron. In particular, in successful modified gravity models that have a mechanism to recover GR on small scales, the scalaron generally acquires nonlinear interaction terms, and they play an important role to recover GR on small scales. Thus, we need to properly take into account such a nonlinear interaction of the scalaron, and consistently solve the evolution equations of density and velocity divergence.

In Ref. Koyama et al. (2009), we have developed a formalism to calculate the nonlinear power spectrum in a wide class of modified gravity models, including gravity and Dvali-Gabadadze-Porratti (DGP) braneworld Dvali et al. (2000) models. The formalism perturbatively treats the effect of nonlinear scalaron, and employing the standard PT technique, we have explicitly computed the power spectrum of density field, , at one-loop order, which reproduces the -body results at quasi-linear scales. In what follows, we adopt this formalism to perturbatively compute auto- and cross-power spectra of density and velocity divergence. The basic equations for perturbations are briefly summarized in Appendix A.

To compute the one-loop power spectra, specifically in the gravity model presented below (Sec. III), the analytic calculations starting naively with the basic equations in Appendix A is technically difficult in practice, because the perturbation equations cannot be separately treated in time and scales. Instead of solving the equations for and in Appendix A, we will numerically solve the evolution equations for the power spectra called the closure equation in Ref. Koyama et al. (2009) [Eqs. (4.3)(4.4)(4.5)], which has been derived by truncating an infinite chain of the moment equations at one-loop order. Implementation and technical detail of the numerical scheme to solve the closure equations are presented in Ref. Hiramatsu and Taruya (2009) (see also Appendix A of Ref. Koyama et al. (2009)).

#### ii.2.2 A term

Next consider the term. The expression given in Eq. (4) can be rewritten with a more convenient form suited for numerical integration. Introducing the doublet777This is somewhat different from the frequently-used definition in GR, with being the linear growth rate, (e.g., Bernardeau et al. (2008); Taruya et al. (2009))., , we define the bispectrum :

 ⟨Φa(\boldmathk1)Φb(\boldmathk2)Φc(\boldmathk3)⟩ =(2π)3δD(\boldmathk1+% \boldmathk2+\boldmathk3)Babc(\boldmathk1,\boldmathk2,\boldmathk3). (7)

In terms of this, the three-dimensional integral is reduced to the sum of the two-dimensional integrals, and the final form of the term becomes Taruya et al. (2013)

 A(k,μ)=3∑n=12∑a,bμ2n(−1)a+b−1k3(2π)2 ×∫∞0dr∫1−1dx{Anab(r,x)B2ab(\boldmathp,\boldmathk−\boldmathp,−% \boldmathk) +˜Anab(r,x)B2ab(\boldmath% k−\boldmathp,\boldmathp,−\boldmathk)} (8)

with and . The non-vanishing components of and are exactly the same as those presented in Ref. Taruya et al. (2013) (see Sec. III-B2).

Since the term appears as a next-to-leading order correction, the tree-level calculation of the bispectrum is sufficient for a consistent calculation of redshift-space power spectrum at one-loop order. Expanding the doublet as , and assuming the Gaussian initial condition, the bispectrum at the tree-level order becomes

 Babc(\boldmathk1,\boldmathk2,% \boldmathk3;t) =2{F(2)a(\boldmathk2,% \boldmathk3;t)F(1)b(k2;t)F(1)c(k3;t) ×P0(k2)P0(k3)+(cyc.perm.)}, (9)

where the functions are the symmetrized standard PT kernel of the -th order perturbative solutions888Since we are interested in the late-time evolution of cosmic structure, we only consider the fastest growing term.:

 Φ(n)a(\boldmathk;t)=∫d3% \boldmathk1⋯d3\boldmathkn(2π)3(n−1)δD(\boldmathk−\boldmathk1⋯n) ×F(n)a(\boldmathk1,⋯,\boldmathkn;t)δ0(\boldmathk1)⋯δ0(\boldmathkn) (10)

with . The function is the initial density field, for which we assume Gaussian statistics. The statistical property of is characterized by the power spectrum:

 ⟨δ0(\boldmathk)δ0(\boldmathk′)⟩=(2π)3δD(\boldmathk1+% \boldmathk′)P0(k) (11)

The remaining task in computing the term is to evaluate the PT kernels up to the second order, which can be done analytically. This is also the case with the gravity model given below, although some numerical works are involved. In Appendix A, based on the basic equations, we derive the explicit functional form of the PT kernels and , and summarize the procedure to compute these kernels. The explicit expression for the PT kernels in gravity and DGP models is also presented.

#### ii.2.3 B term

Similar to the term, the term given in Eq. (5) is reduced to the sum of two-dimensional integrals. From Refs. Taruya et al. (2010, 2013), the resultant expression becomes

 B(k,μ)=4∑n=12∑a,b=1μ2nk3(2π)2∫∞0dr∫1−1dx ×Bnab(r,x)Pa2(k√1+r2−2rx)Pb2(kr)(1+r2−2rx)a, (12)

Here, the coefficients are the same as presented in Appendix B of Ref. Taruya et al. (2010) 999In the case of GR, Eq. (12) exactly coincides with Eq. (A4) of Ref. Taruya et al. (2010), but the definition of power spectra is somewhat different. . For the one-loop order of the redshift-space power spectrum, the linear-order power spectra are sufficient to evaluate Eq. (12), and we have

 Pab(k;t)=F(1)a(k;t)F(1)b(k;t)P0(k) (13)

with the linear PT kernel .

## Iii Redshift-space distortions in f(R) gravity model

In this section, as an illustrative example showing the RSD beyond linear scales in modified gravity model, we compute the redshift-space power spectrum in the gravity model, and compare the PT prediction with results of -body simulations.

### iii.1 f(R) gravity

The gravity is one of the representative gravity models, and it has a mechanism to recover GR on small scales. Generalizing the Einstein-Hilbert action to include an arbitrary function of the scalar curvature , the model is given by:

 S=∫d4x√−g[R+f(R)2κ2+Lm], (14)

where and is the Lagrangian of the ordinary matter. This theory is equivalent to the Brans-Dicke theory with parameter , but there is a non-trivial potential. This can be seen from the trace of modified Einstein equations:

 3□fR−R+fRR−2f=−κ2ρ, (15)

where and is a Laplacian operator and we assumed the matter dominated universe. We can identify as the scalaron, i.e., extra scalar field, and its perturbations are defined as

 φ=δfR≡fR−¯¯¯fR, (16)

where the bar indicates that the quantity is evaluated on the background universe. Here, we consider the cases with and . These conditions are necessary to have the background close to cold dark matter (CDM) cosmology. Then the perturbations for scalaron satisfy

 31a2∇2φ=−κ2ρmδ+δR,δR≡R(fR)−R(¯¯¯fR). (17)

Note that this is nothing but the equation for the BD scalar perturbations with .

In what follows, we consider the specific function of the form

 f(R)∝RAR+1, (18)

where is a constant with dimensions of length squared Hu and Sawicki (2007). If we take the limit , we obtain and cosmological constant does not appear. For high curvature , on the other hand, can be expanded as

 f(R)≃−2κ2ρΛ+|fR0|¯R20R, (19)

where is determined by . The quantity is the background curvature today, and we defined (see e.g., Marchini and Salvatelli (2013); Lombriser et al. (2012); Yamamoto et al. (2010); Okada et al. (2013); Schmidt et al. (2009a) for recent cosmological constraints on ).

In the setup of -body simulations and PT calculation below, we will take , and assume that the background expansion just follows the CDM history with the same .

### iii.2 N-body simulations

We use the subset of cosmological -body simulations presented in Ref. Li et al. (2013); Jennings et al. (2012). The data set of -body simulations were created by the -body code for modified gravity models, ECOSMOG Li et al. (2012a), which is a modification of the mesh-based N-body code, RAMSES Teyssier (2002). With cubic boxes of side length Gpc and particles, the initial conditions were generated at redshift using MPgrafic, according to the linear matter power spectrum determined by the cosmological parameters: , , , , , 111111The value of indicated in Ref. Jennings et al. (2012) is a typo..

In the analysis presented below, we mainly consider GR and gravity with (labeled as CDM and F4 in Ref. Jennings et al. (2012)), and focus on the output results at . In Ref. Jennings et al. (2012), with independent realizations, the matter power spectra are measured in redshift space, applying the distant-observer approximation. Adopting the line-of-sight direction perpendicular to each side of the simulation box, the density field is assigned on grids with cloud-in-cell interpolation, and the power spectra are measured for three different line-of-sight directions. Here, in comparison with the PT model, we use the power spectra averaged over all line-of-sight directions and realizations. With this treatment, the measured power spectra tend to be smooth and the outliers disappear, while the estimation of the error covariance is rather complex because of the duplicated power spectrum measurement with same realizations. In what follows, for the analysis of the fitting and parameter estimation, we assume a hypothetical galaxy survey of the volume Gpc, and consider the statistical error limited by the cosmic variance. Unless otherwise stated, the error bars of the -body results indicate the 1- error of this hypothetical survey computed with the linear power spectrum (see Appendix C of Ref. Taruya et al. (2010))121212Eq. (C4) of Ref. Taruya et al. (2010) includes typos. In the parenthesis of the third line, it should be correctly replaced with .

Finally, in addition to the redshift-space power spectra, we also compare the real-space power spectra of the density and velocity fields in Ref. Jennings et al. (2012), which are used to estimate the valid range of PT predictions.

### iii.3 Comparison with N-body simulations

Before presenting the redshift-space power spectrum, we first separately compute the contribution of each term in the power spectrum expression (2), and compare it with -body simulation. Fig. 1 shows the results of standard PT calculation at one-loop order at in the GR (dotted) and F4 ( gravity with , solid) cases, with the same cosmological parameters as adopted in -body simulations. The (left) and (right) terms are plotted together with the auto- and cross-power spectra of density and velocity divergence fields, multiplied by . According to Eqs. (8) and (12), the and terms are expanded as and , and we here plot the scale-dependent coefficients and (: magenta, : cyan, : green, : yellow). Overall, the resultant amplitude of power spectrum corrections in F4 is rather larger than that in GR. The acoustic signature are clearly seen in both GR and F4 cases not only for the real-space quantities but also the RSD correction (i.e., term). The reason for a larger amplitude in gravity is mainly attributed to the scale-dependent enhancement of the linear growth factor on small scales in gravity, and with a large value of , the mechanism to recover GR is still inefficient at quasi-linear scales. This implies that the redshift-space power spectrum can be quite different between GR and F4, and the RSD corrections (i.e., and terms) would play an important role.

To see the domain of applicability of the standard PT calculation, we next plot in Fig. 2 the real-space power spectra , , and from -body simulations, and compare those with the PT results. According to a phenomenological rule calibrated with -body simulations Nishimichi et al. (2009); Taruya et al. (2009), the standard PT power spectrum at one-loop order is expected to agree with -body simulations at Mpc with an accuracy of level. Although the simulation results show somewhat noisy behavior, the PT predictions seem to work well at least at the scales indicated by the empirical rule, where the deviation from linear theory is around in both GR and F4. The result suggests that even in the presence of a substantial difference in the linear growth, the nonlinear gravitational growth itself does not change so much between GR and modified gravity models. This would be probably true as long as the mechanism to recover GR is still inefficient at quasi-linear scales.

Keeping in mind the applicability of PT calculation, we now focus on the redshift-space power spectrum, and compare the PT calculations with -body simulations. In Fig. 3, top panels show the monopole () and quadrupole () moments of power spectrum multiplied by , while bottom panels present the ratio of monopole and quadrupole spectra to the linear theory prediction taking only account of the Kaiser effect. The multipole power spectrum is defined by

 P(S)ℓ(k)=2ℓ+12∫1−1dμP(S)(k,μ)Pℓ(μ), (20)

with being the Legendre polynomials. The PT results based on an improved model of RSD [i.e., TNS model, Eq. (2)] are depicted as solid lines, while the results without correction terms are also shown in dashed lines. In both cases, we adopt the Gaussian damping function in computing PT predictions:

 DFoG(kμσv)=exp[−(kμσv)2]. (21)

Here, the velocity dispersion is a free parameter and is determined by fitting the model predictions to the -body results of monopole and quadrupole spectra up to Mpc (indicated by vertical arrows), corresponding to the valid range of PT. Note that we also examined the Lorentzian form, but the choice of the damping function did not change the results as long as we consider the applicable range of standard PT one-loop.

Fig. 3 shows that the model (2) successfully describes the -body results of RSD in both GR and model. Although the applicable range of standard PT one-loop is limited, the and terms still play an important role. In the presence of these terms, the acoustic signature of redshift-space power spectrum tends to be smeared compared to the real-space power spectrum, and this indeed improves the agreement with -body simulations. In the panels of Fig. 3, we show the reduced chi-squared statistic defined by131313Strictly speaking, the non-vanishing monopole and quadrupole moments of redshift-space power spectra yield a non-zero covariance between them. This is true even in the Gaussian statistics. However, the magnitude of covariance is shown to be fairly small at large scales Taruya et al. (2010, 2011), and the impact of covariance is ignorable in our analysis.

 χ2red=1ν∑ℓ=0,2∑i[P(S)ℓ,N-body(ki)−P(S)ℓ,PT(ki)]2[ΔP(S)ℓ(ki)]2, (22)

with the quantity being the number of degrees of freedom. Here, the statistical error is estimated from the cosmic variance error assuming the survey volume Gpc. The number of Fourier bins in the above summation can be inferred from the maximum wavenumber shown in Fig. 3, depicted as vertical arrows. The resultant taking account of the and terms (i.e., TNS model) are clearly lower than those ignoring the corrections.

To show the quantitative difference of RSD between GR and gravity, Fig. 4 shows the ratio of quadrupole-to-monopole ratio in F4 to that in GR, i.e., . Note that the errorbars of the -body simulation shown in the panel are not the cosmic variance error, but are estimated from the -body data of the realizations for a particular line-of-sight direction. The linear theory predicts a slight enhancement of the ratio, while the actual -body result rather shows a noticeable reduction at small scales. This basically comes from a stronger suppression of the power spectra in gravity, as shown in Fig. 3 (see bottom panel). Fig. 5 summarizes the fitting results of the parameter together with the resultant reduced chi-squared. At , the fitted value of the velocity dispersion is relatively large in F4 by . Neglecting the correction terms, the relative difference of between model and GR is more prominent (), although the values themselves are even smaller than those taking account of the and terms. As a result, the PT prediction ignoring the corrections exhibits a strong damping behavior in Fig. 4, and tends to deviate from -body simulations at small scales. By contrast, the prediction with and terms (i.e., TNS model) faithfully traces the -body trend beyond the applicable range of the standard PT.

Finally, while we mainly presented the results at , we briefly comment on other cases at , where we also examined the F5 case ( gravity with ). All the results are summarized in Fig. 5. At , the nonlinear clustering is strongly developed, and the applicable range of standard PT one-loop is quite limited. Nevertheless, with a limited fitting range of Mpc, the PT results show an excellent performance with , and the prediction including the and terms gives a better agreement with -body simulations. Fig. 5 shows that the fitted value of velocity dispersion in gravity is generally larger than that in GR, roughly consistent with the one estimated with linear theory (solid lines). This suggests that a stronger damping of the power spectrum amplitude may be a good indicator for modified gravity, as pointed out by Ref. Jennings et al. (2012) (see also Ref. Lam et al. (2012, 2013)). Note, however, that the actual value of depends on the underlying model of RSD. Further, our observable is not dark matter but galaxy distribution, which does not faithfully trace the dark matter distribution. A careful study is needed, and we leave this issue to future work.

## Iv Implications

Having confirmed that the PT model of RSD works well at quasi-linear scales, in this section, we discuss the potential impact of the PT template on the RSD measurement at quasilinear scales.

In testing gravity with RSD, a primary goal would be to clarify whether GR really holds on cosmological scales or not. In this respect, the measurement of the linear growth rate, , provides an important clue, and with the GR-based PT template, we may look for a possible deviation of from GR prediction. One important property in a large class of modified gravity models, including gravity, is the scale dependence of . Thus, a detection of scale-dependent immediately implies the deviation of gravity from GR. The crucial question is how well one can detect or characterize such scale dependence in a model-independent manner.

Alternatively, we may consider some specific gravity models, and try to directly constrain the models themselves. In this case, the linear growth rate might not be an appropriate indicator to characterize a possible deviation from GR. Rather, one tries to directly constrain the model parameter of modified gravity (e.g., in the case of gravity). Then, with the prior assumption of the specific gravity models, the question is how well we can accurately constrain the model parameter based on the PT template of RSD in an un-biased way.

Below, we will separately consider these two issues, and examine the parameter estimation analysis. Note that we will adopt below of Eq. (22) to estimate the goodness-of-fit. Strictly speaking, this is not entirely correct, because the nonlinear gravitational evolution induces the non-Gaussian contribution, which produces non-vanishing power spectrum covariances between diffrent Fourier modes. However, it is shown in the GR case that as long as we consider the quasi-linear scales at moderately high redshift, the off-diagonal componets are small enough, and the diagonal components can be approximately described by the simple Gaussian contribution, leading to a negligible influence on the parameter etstimation (e.g., Takahashi et al. (2009, 2011)). We thus expect that the same would be true in our case of the gravity model that is close to GR, and our simple treatment with Gaussian error contribution would be validated at quasi-linear scales.

### iv.1 Constraining model parameters of modified gravity

Let us first consider the model-dependent analysis to constrain the model parameter of modified gravity, assuming the gravity with as our fiducial gravity model. For specific functional form with Eq. (19) [or Eq. (18)], the parameter is the only parameter characterizing a deviation of gravity from GR. Thus, the test of gravity is made possible with constraining the model parameter by fitting the theoretical template to the data set of redshift-space power spectrum. Here, as a simple demonstration, we ignore the effect of galaxy bias, and allowing to flow, we fit the PT template to the -body data at .

Fig. 6 shows the results of parameter estimation based on the Markov chain Monte Carlo (MCMC) technique. Assuming the hypothetical survey limited by the cosmic variance error with the survey volume Gpc, the best-fit value of and the 1- statistical error are derived, and are plotted (top) as function of maximum wavenumber, , together with the reduced chi-squared statistic (bottom), where represents the range of the wavenumber used for parameter estimation. Note here that the number of free parameters is two, i.e., and . Accordingly, the derived constraint is rather tight, and a slight discrepancy between the template and data can lead to a biased estimation of the . Fig. 6 shows that only the improved model of RSD computed in gravity (filled circles) recovers the fiducial out to  Mpc, corresponding to the applicable range of standard PT one-loop. A slight change of the PT template, depicted as open circles and filled triangles, leads to a biased estimation of the model parameter. Ignoring the damping function (crosses) further adds a large systematic error. This is even true at  Mpc.

Left panel of Fig. 7 shows the representative result of the two-dimensional constraints on and taken from Fig. 6, where we fix the maximum wavenumber to  Mpc. The meaning of color types are the same as in Fig. 6, and in each error contour, inner and outer contours respectively represent the - ( C.L.) and - ( C.L.) constraints. Overall, the degeneracy between and is weak, and the result suggests that at the scales accessible by PT template, the model parameter can be constrained down to from future RSD measurements.

Note, however, that this is only true when we properly take account of the effect of modified gravity in computing the PT template. Most of the analysis in the literature considered the effect of modified gravity only in the linear growth rate and incorporated it into the GR-based template to constrain the model parameter using the measurements of RSD (e.g., Yamamoto et al. (2010); Okada et al. (2013) for recent works). The right panel of Fig. 7 indeed demonstrates such a case. That is, we adopt the GR-based PT template in which the effect of modified gravity is only incorporated in the linear growth rate . In GR, the velocity-divergence field is known to be factorized as , where is perturbatively expanded as . As a result, at a given redshift, the PT template of the redshift-space power spectrum is described as the function of , and , i.e., . Since the growth rate controlls the strength of RSD, we naively expect that simply incorporating the scale-dependent in modified gravity into the PT template allows us to faithfully constrain the model parameter .

However, this actually leads to a biased estimation of the model parameter , as shown in the contour with orange color of Fig. 7. The reason for the large systematic bias is ascribed to the fact that the modification of gravity not only alters the linear growth rate but also affects the shape of the real-space power spectra because of the scale-dependent growth, as clearly shown in Fig. 2. Thus, for an unbiased estimation of , we need to additionally incorporate the effect of gravity bias, that accounts for the relative difference of the clustering amplitude between GR and gravity, into the PT template. The contour with magenta color is the results taking account of this gravity bias, simply assuming the following relation:

 δn-body,F4(\boldmathk)=b(k)δPT,GR(\boldmathk);b(k)=1+A2k21+A1k, (23)

where is the density field in -body simulation, whilst is the density field for the PT calculation. The function characterizes the scale-dependent growth relative to the GR prediction, and we adopt here the functional form similar to those frequently used to model the galaxy bias (e.g., Cole et al. (2005); Sanchez and Cole (2008)). Allowing the parameters and to float, the result marginally reproduces the fiducial value of , and the goodness-of-fit quantified by is improved. With the increased number of free parameters, however, constraining power is significantly reduced, and the size of error contour indeed becomes large (c.f. left panel of Fig. 7). This proves that the heterogeneous PT template is insufficient to tightly constrain the model parameter of modified gravity, and a full PT modeling taking proper account of the modified gravity is required for unlocking the full power of precision RSD measurement.

### iv.2 Model-independent detection of a small deviation from GR

Consider next the model-independent test of GR, and discuss how well we can characterize or detect the scale dependence of the linear growth rate, . Here, for illustrative purpose, we examine the two simple cases. One is to divide the power spectrum data into several wavenumber bins, and in each bin, we try to estimate to see a possible deviation from spatially homogeneous . The other case is to assume a specific functional form of , and to constrain its parameters. In both cases, similar to the analysis shown in right panel of Fig. 7, we adopt the GR-based PT template with an improved model of RSD (i.e., TNS model), and take account of the gravity bias in Eq. (23). We then fit the template to the monopole and quadrupole power spectra at measured from -body simulations of gravity with .

Fig. 8 shows the result of MCMC analysis for the binned linear growth rate, where we set  Mpc, and divide the power spectrum data into three equal bins. Dotted and solid lines represents the linear growth rate of the gravity with and without binning, while the vertical errorbars of the binned results indicate the - statistical uncertainty derived from the MCMC analysis, marginalized over other nuisance parameters. Note that number of free parameters is . The best-fit value of in each bin is close to the fiducial value, but slightly away from linear theory prediction except for the central bin. As a result, the errorbars share almost the same value of , and no notable trend of the scale-dependent growth is found from the binned estimate of .

Fig. 9 examines the other case, in which we assume a specific functional form of given below:

 fapprox(k)=f0[1+ϵtanh(k/kc)]. (24)

Allowing the parameters , and to float, we perform the MCMC analysis. Again, the number of free parameters is , and we set  Mpc. Note that for the scales of our interest, Eq. (24) is shown to accurately describe the scale-dependent linear growth rate of gravity, and fitting directly Eq. (24) to the linear theory prediction at , we obtain , , and  Mpc.

The left panel of Fig. 9 shows the two-dimensional projected errors on the parameters, , , and . The MCMC analysis of the RSD measurement favors non-zero values of these parameters, strongly indicating a deviation from spatially constant . However, a closer look at two-dimensional contours reveals a substantial difference in between the best-fit result and the directly fitted value, . As a result, the MCMC result is unable to reproduce the underlying scale-dependent linear growth rate. Right panel of Fig. 9 shows the constraint on the scale-dependence of . Based on the best-fit values and the associated - errors shown in left panel of Fig. 9, the best-fit curve is plotted in red solid line, and its - statistical uncertainty is shown in red shaded region. The scale-dependence inferred from the MCMC result is rather stronger than that of the linear theory prediction (black dashed).

These two examples imply that the model-independent detection and characterization of the scale-dependent are generally difficult, and the results are rather sensitive to the choice of parameterized form of the linear growth rate. This is presumably because each of the parameters characterizing the scale-dependent cannot be determined locally, but rather it must be estimated with a wide range of wavenumber. Then the parameters tend to be highly correlated with each other, leading to a biased estimation. In this respect, a sophisticated treatment with principal component analysis may provide a way to robustly detect a scale-dependent .

## V Conclusion

In this paper, we studied how well we can clarify the nature of gravity at large scales with redshift-space distortions (RSD), especially focusing on the quasilinear regime of the gravitational evolution. While most of previous works have been done with the theoretical template assuming GR as underlying gravity theory, we here developed a new perturbation theory (PT) prescription for RSD in the general context of the modified gravity models. Extending our previous works on the improved model of RSD proposed by Ref. Taruya et al. (2010), we applied the standard PT framework by Ref. Koyama et al. (2009), which has been formulated to deal with a wide class of modified gravity models, to the computation of the redshift-space power spectrum. As a specific application, in this paper, we consider the gravity model, and compared the PT prediction of RSD with results of -body simulations. Albeit the limited applicable range of the standard PT, the PT results successfully describe the -body simulations, and the predicted monopole and quadrupole spectra quantitatively agree with -body results.

Then, we next considered how well we can characterize and/or constrain the deviation of gravity from GR. One obvious approach is to first assume a specific modified theory of gravity as an underlying gravity model and constrain their model parameters. Using the PT as a theoretical template, we performed the parameter estimation analysis, and checked if the theoretical template correctly recovers the fiducial value of the model parameter in the -body simulations. Adopting the improved model of RSD [TNS model, Eq. (2)], a full PT template calculated in the modified gravity model was found to reproduce the correct model parameter, while a slight deficit in the PT template led to a biased parameter estimation. As another approach, we have also examined the model-independent analysis, and based on the PT template calculated in GR, we tried to characterize the scale-dependent linear growth rate from monopole and quadrupole power spectra. Without assuming any modified gravity model, the parameterization of the scale-dependent linear growth rate is necessary, and the parameters characterizing are highly correlated in general. Our simple two examples suggest that the results are highly sensitive to the choice of parameterization, and it is generally difficult to characterize the scale-dependence of in an unbiased manner unless employing some sophisticated methods such as principal component analysis.

Throughout the paper, we have worked with the standard PT, but the standard PT is known to have a bad convergence property. While we can still get a fruitful constraint on modified gravity models, resummed PT schemes with a wide applicable range are highly desirable to improve the observational constraint. A development of improved PT template in redshift space is an important future direction (see Brax and Valageas (2012, 2013) for recent attempt). Another important issue is the application of the present prescription to the real measurement of RSD. With full PT implementation of the theoretical template, a tight cosmological constraint is expected to be obtained in a robust and unbiased way. In doing this, however, a proper account of the galaxy bias would be crucial. Although the present paper mainly focused on the matter power spectrum, the galaxy bias would be also affected by the modification of gravity, and this may produce a non-trivial scale-dependent shape of the observed power spectrum. In fact, the -body study of the halo clustering properties has revealed that the halo bias in gravity is systematically lower than that in GR (e.g., Schmidt et al. (2009b)). Since even the velocity dispersion and clustering amplitude of the dark matter distribution in gravity differ from those in GR (see Figs. 2 and 3 in Sec. III), coupled with the nonlinear gravity, this could impose a non-trivial trend in the halo/galaxy bias (see e.g., Lombriser et al. (2013); Ferraro et al. (2011); Li et al. (2012b); Fontanot et al. (2013) for recent study on the abundance and clustering of halos and galaxies). Hence, a careful study of the halo/galaxy bias is necessary together with extensive tests with -body mock catalogs toward an unbiased test of gravity.

###### Acknowledgements.
We are grateful to Elise Jennings for providing us the power spectrum data of -body simulations. We also thank Baojiu Li and Gong-bo Zhao for allowing us to use the simulation data. This work is supported in part by a Grant-in-Aid for Scientific Research from the Japan Society for the Promotion of Science (No. 23740186 for T.H and No. 24540257 for A.T). K.K is supported by STFC grant ST/K0090X/1, the European Research Council and the Leverhulme trust. T.H. acknowledges a support from MEXT HPCI Strategic Program.

## Appendix A Basic equations for perturbations and second-order kernels

In this appendix, after briefly reviewing the formalism developed in Ref. Koyama et al. (2009), we derive general expressions for the second-order PT kernel , as well as the linear growth factor in modified gravity models. In Sec. A.1, we begin by reviewing the framework to treat the evolution of matter fluctuations in modified gravity models. We then develop the perturbation theory and derive the PT kernels up to the second-order in Sec. A.2. The explicit expressions for PT kernels are given in specific modified gravity models, i.e., gravity and Dvali-Gabadadze-Porratti (DGP) models.

### a.1 Evolution equations

Let us first consider the matter sector. Apart from the force law of gravity, the basic equations governing the evolution of matter sector is basically described by the conservation of energy momentum tensor, which would remain unchanged even in the modified gravity model. Hence, under the single-stream approximation, the matter fluctuations are treated as pressureless fluid flow, whose evolution equations are the continuity and Euler equations:

 ∂δ∂t+1a∇⋅[(1+δ)\boldmathv]=0, (25) ∂\boldmathv∂t+H\boldmath% v+1a(\boldmathv⋅∇)⋅\boldmathv=−1a∇ψ, (26)

where is the Newton potential.

On the other hand, for the gravity sector, there may appear a new scalar degree of freedom referred to as the scalaron, which results in a large-distance modification to the gravity. On large scales, the scalaron mediates the scalar force, and behaves like the Brans-Dicke scalar field without potential and self-interactions, while it should acquire some interaction terms on small scales, which play an important role to recover GR. Indeed, for several known mechanisms such as chameleon and Vainshtain mechanisms (e.g., Khoury and Weltman (2004); Deffayet et al. (2002)), the nonlinear interaction terms naturally arise and eventually become dominant, leading to a recovery of GR. As a result, even on subhorizon scales, the Poisson equation is modified, and is coupled to the field equation for scalaron with self-interaction term:

 1a∇2ψ=κ22ρmδ−12a2∇2φ, (27) (3+2ωBD)1a2∇2φ=−2κ2ρmδ−I(φ) (28)

with and being the Brans-Dicke parameter. Here, we have used the quasi-static approximation and neglected the time derivatives of the perturbed quantities compared with the spatial derivatives. This treatment is always valid as long as we consider the evolution of matter fluctuations inside the Hubble horizon. The function represents the nonlinear self-interaction, and it can be expanded as

 I(φ)=M1(k)+12∫d3% \boldmathk1d3\boldmathk2(2π)3δD(%\boldmath$k$−\boldmathk12) ×M2(\boldmathk1,\boldmath% k2)φ(\boldmathk1)φ(\boldmathk2)+⋯ (29)

In Fourier space, Eqs. (25)–(28) can be reduced to a more compact form. Assuming the irrotationality of fluid quantities, the velocity field is expressed in terms of velocity divergence, . Then, we have Koyama et al. (2009),

 H−1∂δ(\boldmathk)∂t+θ(\boldmathk)=−∫d3\boldmathk1d3% \boldmathk2(2π)3δD(\boldmathk−\boldmath% k12) ×α(\boldmathk1,\boldmathk2)θ(\boldmathk1)δ(% \boldmathk2), (30) H−1∂θ(\boldmathk)∂t+{2+˙HH2}θ(\boldmathk)+κ2ρm2H2{1+13(k/a)2Π(% \boldmathk)}δ(\boldmathk) =−12∫d3\boldmathk1d3\boldmathk2(2π)3δD(\boldmathk−% \boldmathk12)β(\boldmathk1,\boldmathk2)θ(\boldmathk1)θ(\boldmathk2) −12(ka)2S(\boldmathk)H2, (31)

with the mode-coupling kernels, and given by

 α(\boldmathk1,\boldmathk2)=1+\boldmathk1⋅\boldmathk2|\boldmathk1|2, β(\boldmathk1,\boldmathk2)=(\boldmathk1⋅\boldmathk2)|\boldmathk1+\boldmathk2|2|\boldmathk1|2|\boldmathk% 2|2.

In the above, the function characterizes the deviation of the Newton constant from GR, while the quantity is originated from the non-linear interactions of the scalaron, which is responsible for the recovery of GR at small scales. The functional form of these are obtained from the Poisson equation and field equation for scalaron, and the expressions relevant for the second-order perturbations are respectively given by Koyama et al. (2009):

 Π(k)=13{(3+2ωBD)k2a2+M1(k)}, S(k)=−16Π(k)(κ2ρm3)2∫d3\boldmathk1d3\boldmathk2(2π)3δD(\boldmathk−\boldmathk12) ×M2(\boldmathk1,\boldmathk2)δ(\boldmathk1)δ(\boldmathk2)Π(k1)Π(k2)+⋯. (32)

Here, in deriving the last expression, we perturbatively solve the scalaron field in terms of using Eqs. (28) and (29). The explicit functional form of or , and depends on actual modified gravity models, which will be specified later.

### a.2 PT kernels

Let us now perturbatively solve the evolution equations (30) and (31). Consider first the linear-order solutions. Ignoring all the non-linear terms in Eqs. (30) and (31), we obtain

 ¨δ(1)(\boldmathk)+2H˙δ(1)(\boldmathk)−κ2ρm2{1+13(k/a)2Π(k)}δ(1)(\boldmathk)=0, θ(1)(\boldmathk)=−1H˙δ(1)(\boldmathk) (33)

The solutions and can be formally expressed as

 δ(1)(\boldmathk;t)=D(k;t)δ0(% \boldmathk),θ(1)(\boldmathk;t)=−˙D(k;t)Hδ0(\boldmathk), (34)

where the function is the initial density field [see Eq. (11)]. The function is the linear growth factor, and satisfies the following evolution equation:

 ¨D+2H˙D−κ2ρm2{1+13(k/a)