Lagrangian perturbation theory for modified gravity

TL;DR

This paper develops a formalism for calculating Lagrangian displacement fields in modified gravity theories, accounting for scale-dependent forces and screening effects, and applies it to compute power spectra in $f(R)$ gravity.

Contribution

It introduces a third-order perturbation theory formalism for modified gravity, including frame-lagging terms, and demonstrates its application to $f(R)$ gravity power spectrum calculations.

Findings

Nonlinearities cause force screening in scalar-tensor theories.

Frame-lagging terms are crucial at large scales for consistency.

The formalism accurately computes 1-loop power spectra in $f(R)$ gravity.

Abstract

We present a formalism to compute Lagrangian displacement fields for a wide range of cosmologies in the context of perturbation theory up to third order. We emphasize the case of theories with scale dependent gravitational strengths, such as chameleons, but our formalism can be accommodated to other modified gravity theories. In the non-linear regime two qualitative features arise. One, as is well known, is that nonlinearities lead to a screening of the force mediated by the scalar field. The second is a consequence of the transformation of the Klein-Gordon equation from Eulerian to Lagrangian coordinates, producing frame-lagging terms that are important especially at large scales, and if not considered, the theory does not reduce to the $\Lambda$CDM model in that limit. We apply our formalism to compute the 1-loop power spectrum and the correlation function in $f(R)$ gravity by using different resummation schemes. We further discuss the IR divergences of these formalisms.

Equations

$${\mathcal{L}}_{ST}=\frac{1}{16\pi }\sqrt{-g}\left[f(\phi )R-g(\phi ){\nabla }_{\mu }\phi {\nabla }^{\mu }\phi -V(\phi )\right]+{\mathcal{L}}_m(\Psi ,h(\phi )g_{\mu \nu }),$$

$$\mathcal{L}=\frac{1}{16\pi G}\sqrt{-g}\left[\phi R-\frac{{\omega }_{\text{BD}}(\phi )}{\phi }{\nabla }_{\mu }\phi {\nabla }^{\mu }\phi -V(\phi )\right]+{\mathcal{L}}_m(\Psi ,g_{\mu \nu }),$$

$${\mathcal{L}}_R=\sqrt{-g}(R+f(R)).$$

$$\mathcal{L}=-\frac{1}{2}{Z}^{\mu \nu }(\phi ,\partial \phi ,{\partial }^2\phi ){\partial }_{\mu }\phi {\partial }_{\nu }\phi -V(\phi )+\beta (\phi )T_{\mu }^{\mu },$$

$$\frac{1}{a^2}{\nabla }^2\psi$$

$$(3+2{\omega }_{\text{BD}})\frac{1}{a^2}{\nabla }^2\phi$$

$$(3+2{\omega }_{\text{BD}})\frac{k^2}{a^2}\phi =8\pi G\rho \delta -\mathcal{I}(\phi ).$$

$$\delta \mathcal{I}(\phi )$$

$$+\frac{1}{6}\int \frac{d^3{k}_1{d}^3{k}_2{d}^3{k}_3}{(2\pi {)}^6}{\delta }_D(\text{k}-{\text{k}}_{123})M_3({\text{k}}_1,{\text{k}}_2,{\text{k}}_3)\phi ({\text{k}}_1)\phi ({\text{k}}_2)\phi ({\text{k}}_3)+\cdots ,$$

$$G_{\mu \nu }+f_R{R}_{\mu \nu }-{\nabla }_{\mu }{\nabla }_{\nu }{f}_R-\left(\frac{f}{2}-\square f_R\right)g_{\mu \nu }=8\pi G{T}_{\mu \nu },$$

$$3\square f_R=R(1-f_R)+2f-8\pi G\rho ,$$

$$\frac{3}{a^2}{\nabla }^2\delta f_R=-8\pi G\bar{\rho }\delta +\delta R,$$

$$\delta R=\sum_{i}\frac{1}{n!}M_{n}(\delta f_{R})^{n},\qquad M_{n}\equiv\frac{d^{n}R(f_{R})}{df_{R}^{n}}\Bigg{|}_{f_{R}=\bar{f}_{R}},$$

$$m=\sqrt{\frac{M_1}{3}}.$$

$$f(R)=-M^2\frac{c_1(R/M^2{)}^n}{c_2(R/M^2{)}^n+1},$$

$$f_R\simeq f_{R0}{\left(\frac{\bar{R}}{R}\right)}^{n+1},$$

$$\bar{R}=6(\dot{H}+2H^2)=3H_0^2({\Omega }_{m0}{a}^{-3}+4{\Omega }_{\Lambda }).$$

$$M_1(a)=\frac{3}{2}\frac{H_0^2}{|f_{R0}|}\frac{({\Omega }_{m0}{a}^{-3}+4{\Omega }_{\Lambda }{)}^3}{({\Omega }_{m0}+4{\Omega }_{\Lambda }{)}^2},$$

$$M_2(a)=\frac{9}{4}\frac{H_0^2}{|f_{R0}{|}^2}\frac{({\Omega }_{m0}{a}^{-3}+4{\Omega }_{\Lambda }{)}^5}{({\Omega }_{m0}+4{\Omega }_{\Lambda }{)}^4},$$

$$M_3(a)=\frac{45}{8}\frac{H_0^2}{|f_{R0}{|}^3}\frac{({\Omega }_{m0}{a}^{-3}+4{\Omega }_{\Lambda }{)}^7}{({\Omega }_{m0}+4{\Omega }_{\Lambda }{)}^6}.$$

$$\text{x}(\text{q},t)=\text{q}+\mathcal{S}(\text{q},t),$$

$$\delta (\text{x},t)=\frac{1-J(\text{q},t)}{J(\text{q},t)},$$

$$J_{ij}:=\frac{\partial x^i}{\partial q^j}={\delta }_{ij}+\frac{\partial {\Psi }^i}{\partial q^j}$$

$${\nabla }_{\text{x}}\cdot (\ddot{\Psi }(\text{q},t)+2H\dot{\Psi })=-\frac{1}{a^2}{\nabla }_{\text{x}}^2\psi (\text{x},t)$$

$$\frac{1}{a^2}{\nabla }_{\text{x}}^2\psi =4\pi G\bar{\rho }\delta (\text{x},t)-\frac{1}{2a^2}{\nabla }_{\text{x}}^2\phi .$$

$$(3+2w_{\text{BD}})\frac{k_x^2}{a^2}\phi ({\text{k}}_x,t)=8\pi G\bar{\rho }\delta ({\text{k}}_x,t)-M_1({\text{k}}_x,t)\phi ({\text{k}}_x,t)-\delta I(\phi ).$$

$$\hat{\mathcal{T}}=\frac{d^2}{d{t}^2}+2H\frac{d}{dt},$$

$$\nabla_{{\text{\bf x}}}\cdot\mathbf{\hat{\mathcal{T}}}\Psi=-4\pi G\bar{\rho}\delta({\text{\bf x}},t)+\frac{1}{2a^{2}}\nabla^{2}\varphi+\frac{1}{2a^{2}}\big{(}\nabla^{2}_{\text{\bf x}}\varphi-\nabla^{2}\varphi\big{)},$$

$$J$$

$$(J^{-1}{)}_{ij}$$

$$\Psi$$

$${\nabla }_{\text{x}}^2\phi -{\nabla }^2\phi =(J^{-1}{)}_{im}{\nabla }_m((J^{-1}{)}_{in}{\nabla }_n\phi )-{\nabla }^2\phi \sim \mathcal{O}({\lambda }^2),$$

Full text

I Introduction

General Relativity (GR) is the most successful theory of gravity we have nowadays. Indeed, it has passed so far all experimental tests ranging from 1$\mu$m to 1 AU $\sim 10^{11}$m Will:2014kxa ; Berti:2015itd . Nevertheless, mainly since the discovery of the accelerated expansion of the universe Perlmutter:1998np ; Riess:1998cb , concerns about its validity at cosmological scales have arisen. In addition, with the recent and forthcoming advent of a lot of cosmological data, especially from modern surveys Beutler:2016ixs ; Garcia-Fernandez:2016oud ; Aghamousa:2016zmz ; Racca:2016qpi ; Abell:2009aa , there is an increasing interest to test gravity at cosmological scales.

The minimal modification to the Einstein-Hilbert action compatible with the principle of generalized covariance —and the only that introduces no additional degrees of freedom— is provided by the introduction of a cosmological constant; this leads to the standard model of Cosmology, the $\Lambda$CDM model, that can explain the cosmic background radiation anisotropies with exquisite precision Ade:2015xua , the late time background evolution probed by measurements of the baryon acoustic oscillations (BAO) and supernova distance-redshift relations Betoule:2014frx , among several other observations. But this success comes with a price, the $\Lambda$ side of the theory suffers from the smallness and coincidence problems Weinberg:1988cp ; Martin:2012bt . One alternative to the cosmological constant is that of evolving dark energy, in which exotic matter fields endowed with a negative pressure counteract the gravitational attraction and generate the cosmic acceleration Ratra:1987rm ; Tsujikawa:2010sc . From the observational point of view, a practical inconvenience of evolving dark energy is that a background expansion of the Universe that is close to $\Lambda$CDM implies that collapse features are nearly indistinguishable as well, as long as the dark energy is stress-free, and as a consequence it would be difficult to call the correct theory Bertschinger:2006aw ; Bertschinger:2008zb . The case of Modified Gravity (MG) is different because the two scalar gravitational potentials of the metric are not necessarily equal even in the absence of stress sources, and therefore a theory that mimics a dark energy model at the background level can predict a very different structure formation history.

The construction of alternatives theories of gravity has been a difficult endeavor because they are restricted to reduce to GR in some limits, while still having an impact on cosmic scales in order to provide the cosmic acceleration. For this reason, cosmologists have concentrated in gravitational theories that poses a screening mechanism Khoury:2010xi . In theories with one additional degree of freedom this can be achieved by the non-linearities of its evolution equation —the Klein-Gordon equation in the case of a scalar field. The general idea is that in regions with high density or large gradients of it, the additional (the fifth) force becomes negligible, recovering GR. These effects can be achieved by the scalar field itself or by its spatial derivatives, leading to different types of screenings: as, for example, the chameleon Khoury:2003aq ; khoury:2003rn , kinetic Babichev:2009ee , or Vainshtein Vainshtein:1972sx mechanisms.

In this work we are interested in the formation of cosmological structures at large scales. In the standard scenario, this process is driven by gravitational collapse of cold dark matter (CDM). Once baryons decouple from the primordial plasma, they fall into potential wells already formed by the CDM, and ultimately forming the galaxies we observe today. $N$-body simulations are the most trustable method to study this process, especially at small scales, and it is considered in some sense to lead to the correct solution to the problem. Full $N$-body simulations have been provided for a few MG theories in Zhao:2010qy ; BLi2013 ; Puchwein:2013lza ; Llinares:2013jza ; Rizzo:2016mdr ; Winther:2015wla , while semi-analytical-$N$-body codes that rely on screening for spherical configurations Winther:2014cia , and the COLA approach, introduced in Tassev:2013pn for $\Lambda$CDM, were implemented recently in Valogiannis:2016ane ; Winther:2017jof . The main disadvantage of simulations is that they are computationally very expensive; but furthermore, the understanding of how the process of structure formation takes place stay relatively hidden. For these reasons, analytical and semi-analytical methods have been developed during the last few decades; see Sahni:1995rm ; Ber02 ; Cooray:2002dia for reviews in standard cosmology, and Refs. Clifton:2011jh ; Koyama:2015vza ; Joyce:2016vqv for MG.

Of our particular interest in this work are perturbation theories of large scale structure formation, here the relevant fields are perturbed about their background cosmological values and evolution equations are provided by taking moments of the Boltzmann equation. Perturbation theory comes in different versions, each one having its advantages and disadvantages, mainly Lagrangian Perturbation Theory (LPT) Zel70 ; Buc89 ; Mou91 ; Catelan:1994ze ; BouColHivJus95 ; TayHam96 ; Mat08a ; Carlson:2012bu ; Sugiyama:2013mpa ; Vlah:2014nta ; Matsubara:2015ipa and Eulerian Standard Perturbation Theory (SPT) Juszkiewicz01121981 ; Vishniac01061983 ; Fry:1983cj ; Goroff:1986ep ; Jain:1993jh ; Carlson:2009it , but so many others have been tailored in order to improve its performance, for example see Refs. McDonald:2006hf ; Crocce:2005xy ; Taruya:2007xy ; Pietroni:2008jx . The main drawback of perturbation theory is that the expansion of fields become meaningless when higher orders become larger than lower orders, this occurring as soon as non-linearities become important, but at the weakly non-linear scales it leads to notable improvements over linear theory.

In MG, perturbation theory for large scale structure formation was developed first in 2009PhRvD..79l3512K based in the formalism of the closure theory of Ref. Taruya:2007xy . That work introduces a Fourier expansion for the screening in which each order carries its own screening, becoming evident the effect of the screening even at large scales. Other works in perturbation theory for MG, some of them including the systematic effects of redshift space distortions (RSD), were presented in Taruya:2013quf ; Brax:2013fna ; Taruya:2014faa ; 2014PhRvD..89d3509T ; Bellini:2015oua ; Taruya:2016jdt ; Bose:2016qun ; Fasiello:2017bot ; Bose:2017dtl . The formalisms developed so far are closer to SPT than to LPT. An advantage of LPT is that it is better in modeling the BAO features of the power spectrum, and it is highly successful in reproducing the correlation function around the acoustic peak Tassev:2013rta . Formally, the power spectrum derived from SPT cannot be Fourier transformed to obtain the correlation function, but even if a damping factor is introduced to make the integral convergent a double peaked structure is observed around the BAO scale. On the other hand, SPT follows better than LPT the broadband power spectrum trend line; in the latter it decays very quickly at small scales because dark matter particles follow their trajectories almost inertially and do not clump enough to form small structures leading to a deficit in the power spectrum Vlah:2014nta .

In this work we develop an LPT theory that can be applied to a wide range of MG theories, essentially to those that can be written as a scalar tensor gravity theory. When the Klein-Gordon equation for the scalar field is expressed in Lagrangian coordinates additional terms arise that can not be neglected, otherwise GR is not recovered at large scales, we treat in detail these contributions. Throughout this work we exemplify the formalism by applying it to the Hu-Sawicky $f(R)$ model Hu:2007nk . Our method can also be applied straightforwardly to dark energy models in which we can neglect the dark energy perturbations; and, by adding a source to the Lagrangian displacement equations these perturbations can be incorporated if they remain linear. Though, we do not work out these ideas here.

Since we are interested in MG models that preserve the equivalence principle, at least at the large scales, matter and Lagrangian displacements statistics should be IR-safe as it is know from the works of Refs. 1996ApJ…456…43J ; Scoccimarro:1995if ; Peloso:2013zw ; Carrasco:2013sva ; Creminelli:2013mca . Therefore, we further discuss potentially undesirable divergences during this work.

Two different second order LPT theories (2LPT) for modified gravity has been recently presented in Valogiannis:2016ane ; Winther:2017jof in order to provide a framework for dealing with large scales in the COLA code. More recently, the 2LPT of Winther:2017jof was extended to include massive neutrinos in Wright:2017dkw . Nevertheless the formalisms of these works differ. We expect our results can clarify these discrepancies.

The rest of the paper is organized as follows: In Sect.II we give a brief introduction to MG theories; in Sect.III we present the general LPT theory; in Sect.IV we show the results for the second order 2LPT theory; Sect.V is devoted to the third order; in Sect. VI we study Lagrangian displacement statistics relevant to obtain the matter power spectrum, and discuss their IR and UV divergences; in Sect.VII we compute the power spectra and correlation function using different resummation schemes and further discuss their properties; and in Sect.VIII we present our conclusions. Some computations are delegated to two appendices.

II Theories of modified gravity

The schemes treated in this work are general modified gravity theories that i) have a screening mechanism to comply with local, gravitational constraints; ii) have a fifth force due to a scalar degree of freedom that acts at intermediate scales; and iii) at large, cosmological distances GR is recovered. General modified gravity scenarios and their cosmological consequences have been reviewed in recent years Clifton:2011jh ; Joyce:2016vqv ; Koyama:2015vza , from which we extract the models we are interested in.

We may consider general scalar-tensor theories Bergmann:1968ve ; 1970ApJ…161.1059N ; PhysRevD.1.3209 ,

[TABLE]

where $f,g,h,$ and $V$ are functions that can be different for various theoretical frameworks. A conformal transformation $h(\varphi)g_{\mu\nu}\rightarrow g_{\mu\nu}$ decouples matter from the scalar field, and therefore in this Jordan Frame point particles follow geodesics of the new metric. One can further redefine $f,g,h,$ and $V$ to have a non-minimal coupling given by

[TABLE]

which reproduces the Brans-Dicke theory when $\omega_{\text{BD}}$ is a constant and $V=0$. It is known that the Lagrangian density of Eq. (2) can be transformed to the Einstein Frame (in which the Ricci scalar is not coupled to scalar fields) by means of yet another conformal transformation. In this way, it is common to study cosmological physics in the Einstein Frame.

Another common type of theory is the $f(R)$ gravity 1970MNRAS.150….1B ; 1980PhLB…91…99S , in which one replaces the Einstein-Hilbert Lagrangian density by a generic function of the Ricci scalar, such that

[TABLE]

This Lagrangian can also be transformed to different frames by using two different formalisms: a Legendre transformation relates Eq. (3) to a scalar tensor theory with null coupling parameter $\omega_{\text{BD}}$ and with no couplings to possibly added matter fields in Eq. (3). The other way is to again apply a conformal transformation to obtain a scalar tensor theory in the Einstein frame BARROW1988515 ; PhysRevD.39.3159 .

The Lagrangian densities mentioned above are the most considered in the literature. Thus, when studying screening mechanisms for the scalar field, one can consider the Einstein frame and scalar field Lagrangian of the following general form Joyce:2014kja ; Koyama:2015vza ,

[TABLE]

where the function $Z^{\mu\nu}$ stands for possible derivative self-interactions of the scalar field, $V(\varphi)$ its a potential, and $T^{\mu}_{\;\mu}$ is the trace of the matter stress-energy tensor. The coupling $\beta$ yields a scalar field solution that depends on the local density; for non-relativistic matter $T^{\mu}_{\;\mu}=-\rho$. Eq. (4) is then a generalization to Eq. (2) in which one introduces more general scalar field kinetic terms. Equation (4) encompasses the most used modified gravity theories —including Horndeski theories Horndeski:1974wa — satisfying known observational constraints but allowing smoking guns to be discriminated or confirmed by near-future cosmological data, e.g. Aghamousa:2016zmz .

The screening of the fifth force at local scales can be realized when considering fluctuations around background values with mass $m(\varphi)$ and due to the presence of the functions $Z^{\mu\nu}$ or $\beta$. According to Ref. Joyce:2014kja , the screening can be classified as due to: i) Weak coupling, in which the function $\beta$ is small in regions of high density (local scales), thus the fifth force is suppressed. At large scales the density is small and the fifth forces acts. Examples of theories of this type are symmetron PhysRevLett.104.231301 ; PhysRevD.72.043535 ; PhysRevD.77.043524 and varying dilaton DAMOUR1994532 ; PhysRevD.83.104026 ; ii) Large mass, when mass of the fluctuation is large in regions of high density and the fifth force is suppressed, and at low densities the scalar field is small and can mediate a fifth force. Examples of this type are Chameleons Khoury:2003aq ; khoury:2003rn ; and iii) Large inertia, when the kinetic function $Z$ depends on the environment, making it large in density regions, and then the coupling to matter is suppressed. One finds two cases, when the first derivatives of the scalar field are large or when the second derivatives are large. Examples of the former are $k$-mouflage models Babichev:2009ee ; PhysRevD.84.061502 , and of the latter are Vainshtein models Vainshtein:1972sx .

Let us now focus on perturbations. When considering quasi-nonlinear scales, it is common to use the Jordan frame, as given by Eq. (2), following 2009PhRvD..79l3512K . Since our interest will be the study of kinematics well inside the cosmological horizon, it seems plausible to take the quasi-static limit, in which time derivatives of the scalar field are much less important than spatial derivatives. This approximation is non-trivial however, as it has been shown in the context of modified gravity models Brax:2013fna ; PhysRevD.92.084061 ; Noller:2013wca ; Bose:2014zba , but its correctness depends upon details of the modified gravity free parameters which vary from model to model. Therefore, it is improper to assume its correctness in general, and in fact its accuracy depends upon model parameters and wavenumbers studied. For the concrete models we will later treat in this work (Hu-Sawicky gravity) the validity of this approximation has been checked within the linear perturbation theory in Noller:2013wca and numerical simulations in Bose:2014zba , and in a general framework, introducing the effective sound speed ($c_{s}$) of the $f(R)$ gravity model in Brax:2013fna ; this later work concludes that the quasi-static approximation depends, not on the wavenumber in units of the cosmological horizon, $k/aH$, but on the combination $c_{s}k/aH$, which must be greater than unity.111In realistic, observational applications this depends on the volume spanned by surveys, and the redshifts to be considered. The bigger and deeper the survey the closer to unity the sound speed must be. Assuming these considerations are taken into account, we consider perturbations in a standard way using the Newtonian gauge, one end up with the following set of MG equations 2009PhRvD..79l3512K :

[TABLE]

where $\bar{\rho}$ is the background matter energy density, $\delta$ its perturbation and $\psi$ the gravitational potential. NL stands for possible non-linearities that may appear in the Lagrangian. As usual, these equations are coupled to the standard continuity and Euler equations for the evolution of the matter fields.

When going to the Fourier space, the Klein-Gordon equation can be written as

[TABLE]

Here the self-interaction term ${\cal I}$ may be expanded following 2009PhRvD..79l3512K as ${\cal I}(\varphi)=M_{1}(k)\varphi+\delta\mathcal{I}(\varphi)$ with

[TABLE]

where the functions $M_{i}$ are in general scale and time dependent and are determined by the particular MG model. In the following sections we will keep this dependence explicitly, although we illustrate our results with $f(R)$ theories for which the $M_{i}$ are only time dependent, as we see below.

Now we consider in more detail $f(R)$ gravity, that was first used to explain the current cosmic acceleration in Capozziello:2002rd ; Carroll:2003wy . Variations with respect to the metric of the action constructed with the Lagrangian density of Eq. (3), lead to the field equations

[TABLE]

where $f_{R}\equiv\frac{df(R)}{dR}$. By taking the trace to this equation we obtain

[TABLE]

where we use a dust-like fluid with $T^{\mu}_{\mu}=-\rho$. In cosmological pertubative treatments one splits $f_{R}=\bar{f_{R}}+\delta f_{R}$ and $R=\bar{R}+\delta R$, where the bar indicates background quantities and $\bar{R}\equiv R(\bar{f}_{R})$. Necessary conditions to get a background cosmology consistent with the $\Lambda$CDM model are $|\bar{f}/\bar{R}|\ll 1$ and $|\bar{f}_{R}|\ll 1$ Hu:2007nk .

In the quasi-static limit, the trace equation can be written as

[TABLE]

where we substracted the background terms in Eq.(10). We recognize Eq.(11) as the Klein Gordon equation for a scalar field with a potential $\delta R$ and Brans-Dicke parameter $\omega_{\text{BD}}=0$. Since $\delta R=R(f_{R})-R(\bar{f}_{R})$, we may expand the potential as

[TABLE]

from which the mass associated to the perturbed field $\delta f_{R}$ can be read as

[TABLE]

In this work we focus our attention to the Hu-Sawicky model Hu:2007nk , defined by

[TABLE]

where the energy scale is chosen to be $M^{2}=H_{0}^{2}\Omega_{m0}$. In this parametrized model, at high curvature ($R\gg M^{2}$) the function $f(R)$ approaches a constant, recovering GR with cosmological constant, while at low curvature it goes to zero, recovering GR; the manner these two behaviors are interpolated is dictated by the free parameters. Given that $f_{RR}>0$ for $R>M^{2}$ the solutions are stable and the mapping to a scalar tensor gravity is allowed Magnano:1993bd . The latter also implies that this gravitational model introduces a single one extra degree of freedom. In order to mimic the background evolution of the $\Lambda$CDM model it is also necessary that $c_{1}/c_{2}=6\Omega_{\Lambda}/\Omega_{m0}$ Hu:2007nk , thus leaving two parameters to fix the model. Noting that the background value is about $\bar{R}\gtrsim 40M^{2}\gg M^{2}$, Eq.(14) simplifies and leads to

[TABLE]

with

[TABLE]

We consider the case $n=1$. Using Eqs.(12) and (15) we find the functions $M_{1}$, $M_{2}$ and $M_{3}$:

[TABLE]

[TABLE]

[TABLE]

These functions depend only on the background evolution since they are the coefficients of the expansion of a scalar field potential about its background value. In other theories, e.g. with non-canonical kinetic terms as Galilieons or DGP, the “potential” includes spatial derivatives, thus the $M_{i}$ functions may carry scale dependences; see, for example, 2009PhRvD..79l3512K for the $M_{i}$ functions in DGP gravity.

With $n$ fixed to unity, the only free parameter is $f_{R0}$. The models constructed with it are named F4, F5, F6 and so on, corresponding to the choices $f_{R0}=-10^{-4},-10^{-5},-10^{-6}$ respectively. From these, F4 represents the largest deviation to GR, and though it is ruled out using linear theory, it is worthy to analyze its effects. These theories are indistinguishable from GR at the background, homogeneous and isotropic cosmic evolution Hu:2007nk as already was assumed in Eq. (16). Therefore throughout this work, we shall use these three models to exemplify the method with an expansion history given by a $\Lambda$CDM model with $\Omega_{m}=0.281$ and $h=0.697$.

In this work we specialize to $f(R)$ theories. Even though, by dealing consistently with the scale dependence of the $M_{i}$ functions, our results cover a wide range of models, essentially the whole Horndeski sector. Indeed, in appendix B of Ref. Bose:2016qun it is shown how to relate the $M_{i}$ functions to the Horndeski free functions.222More precisely, in that work the authors use the so-called $\gamma_{1}$ and $\gamma_{2}$ kernels for screenings, these are obtained by expanding $\delta I$ in terms of overdensities, instead of scalar fields as in Eq. (II). That is, the $\gamma_{i}$ functions are the kernels of the expansion of Eq. (III) below.

III Lagrangian perturbation theory in modified gravity

In a Lagrangian description we follow the trajectories of individual particles with initial position q and current position x. Lagrangian and Eulerian coordinates are related by the transformation rule

[TABLE]

with the Lagrangian map $\mathbf{\mathcal{S}}({\text{\bf q}},t)$ vanishing at the initial time where Lagrangian coordinates are defined. In kinetic theory we define the Lagrangian displacement velocity $\dot{\Psi}({\text{\bf q}},t)$ as a momentum average of particles’ velocities over the ensemble. For CDM scenarios, which we posit here, we can safely neglect velocity dispersions and the identity $\mathcal{S}=\mathbf{\Psi}$ follows Aviles:2015osc . Furthermore, mass conservation implies the relation between Lagrangian displacements and overdensities

[TABLE]

where the Jacobian of the transformation is

[TABLE]

and $J$ is its determinant. We make note that in writing Eq. (21) we identify the Lagrangian displacement field as a random field. We further assume this field to be Gaussian at first order in perturbation theory.

We consider the cases for which the MG theory can be written as a Brans-Dicke like model at linear level. For longitudinal modes we have to solve the system

[TABLE]

[TABLE]

In $x$-Fourier space, and in the quasi-static limit, the Klein Gordon equation is

[TABLE]

The notation ${\text{\bf k}}_{x}$ emphasizes that the wavenumbers correspond to Eulerian coordinates, that are not the correct coordinates to be transformed to the Fourier space, since equations of motion will be given in the $q$-space, therefore we shall work in $q$-Fourier space, where wavenumbers correspond to Lagrangian coordinates. For the sake of brevity, in the following we will omit the time argument $t$ when this not lead to confusions. The term $I(\varphi)=M_{1}\varphi+\delta I(\varphi)$, introduced in 2009PhRvD..79l3512K , is the expansion of the scalar field potential about its background value $\varphi$, being $M_{1}\varphi$ the linear piece and $\delta I(\varphi)$ the non-linear which should give rise to the screening of the fifth force. If the latter is not present, the solutions lead to large violations of local experiments, unless the Brans-Dicke parameter $\omega_{\text{BD}}$ is very large. On the other hand, theories which poses a screening mechanism shield the scalar field in nonlinear regions such that $\nabla^{2}_{\text{\bf x}}\varphi\approx 0$, recovering GR. Following Matsubara:2015ipa , we define the operator

[TABLE]

thus, Eqs. (23) and (24) can be written as

[TABLE]

where $\nabla_{i}={}_{,i}$ is partial derivative with respect to $q^{i}$. With this splitting of the Laplacian we recognize the scalar field as a function of Lagrangian coordinates, the term $\nabla^{2}_{\text{\bf x}}\varphi-\nabla^{2}\varphi$ has a geometrical nature that arises when transforming the Klein-Gordon equation from Eulerian to Lagrangian coordinates. We shall call it the frame-lagging term, and we show below that is non-negligible and it is necessary for recovering $\Lambda$CDM at sufficiently large scales, where the fifth force mediated by the scalar field is essentially zero.

To transform to Lagrangian coordinates we use the relations

[TABLE]

with $\epsilon_{ijk}$ the fully antisymmetric Levi-Civita symbol. In perturbation theory the relevant fields are formally expanded as

[TABLE]

and analogously for $\delta$ and $\varphi$. From now on the control parameter $\lambda$ is absorbed in the definitions of the perturbed fields. Since spatial derivatives in Lagrangian and Eulerian coordinates are related by $\nabla_{{\text{\bf x}}\,i}=(J^{-1})_{im}\nabla_{m}$, the frame-lagging term of Eq. (27) can be written as

[TABLE]

implying it is a non-linear term. The Fourier transform of $\varphi({\text{\bf q}})$ leads to

[TABLE]

where $[\,(\,\cdots\,)\,]({\text{\bf k}})$ means Fourier transform of $(\,\cdots\,)({\text{\bf q}})$ and we defined

[TABLE]

$A(k)$ is the gravitational strength in the MG cases, while $A_{0}$ is for GR.

In general we do not write a tilde over Fourier space functions, since it does not lead to confusion, but we do write a tilde over the overdensity $\delta$ to make clear that the Fourier transform of $\delta({\text{\bf x}})$ is taken with Lagrangian coordinates: it is neither the $q$-Fourier transform of $\delta({\text{\bf q}})$ nor the $x$-Fourier transform of $\delta({\text{\bf x}})$. It is given instead by $\tilde{\delta}({\text{\bf k}})=\int d^{3}qe^{-i{\text{\bf k}}\cdot{\text{\bf q}}}\delta({\text{\bf x}})$, or using Eq. (21) this is

[TABLE]

Only for the linear overdensities we have $\tilde{\delta}^{(1)}({\text{\bf k}})=\delta^{(1)}({\text{\bf k}}_{x})=\int d^{3}xe^{-i{\text{\bf k}}\cdot{\text{\bf x}}}\delta^{(1)}({\text{\bf x}})$.

The function $M_{1}$ is related to the mass of the scalar field by $M_{1}=(3+2\omega_{\text{BD}})m^{2}(a)$, which gives the range of the interaction $1/m(a)$, if it is finite we recover GR at large scales. At sufficiently small scales, such that $M_{1}\ll k^{2}/a^{2}$, we note $A(k)\rightarrow(4+2\omega_{\text{BD}})/(3+2\omega_{\text{BD}})A_{0}$, for which solar system observations restrict $\omega_{\text{BD}}>40000$ Will:2014kxa , making the theory effectively indistinguishable from GR in the absence of non-linear terms. For example, in $f(R)$ gravity the scalar field perturbation is identified with $\delta f_{R}$ as it is discussed in Sec.III, obtaining $\omega_{\text{BD}}=0$ for the Brans-Dicke parameter as can be seen directly from Eq.(11). Therefore, in this situation $A(k\rightarrow\infty)\rightarrow\frac{4}{3}A_{0}$ which would rule out the theory; nevertheless, the nonlinearities of the potential (encoded in $\delta I$ in perturbation theory) are responsible to drive the theory to GR in that limit Hu:2007nk . On the other hand, the large scale behavior is dictated by the mass $m$, related to $M_{1}$ in Eq. (13), and if it is not zero, GR is recovered at large scales.

Now, using Eqs.(22), (27), and (32), we obtain the equation of motion

[TABLE]

and up to third order the frame-lagging term is

[TABLE]

which is obtained from Eqs. (29) and (31). We see below that this term has a key role to understand the correct contributions to LPT. We write it as a Fourier expansion as follows333Throughout this paper we adopt the shorthand notations

$\underset{{\text{\bf k}}_{12\dots n}={\text{\bf k}}}{\int}=\int\prod_{i=1}^{n}\frac{d^{3}k_{i}}{(2\pi)^{3}}(2\pi)^{3}\delta_{D}({\text{\bf k}}-{\text{\bf k}}_{12\cdots n}).$

(39)

and ${\text{\bf k}}_{12\dots n}={\text{\bf k}}_{1}+{\text{\bf k}}_{2}+\cdots{\text{\bf k}}_{n}$.

[TABLE]

Below we give expressions for kernels $\mathcal{K}_{\text{FL}}$. By using Eq. (32) iteratively order by order we can write Eq.(II) as

[TABLE]

For compactness we introduced the function

[TABLE]

In LPT it is usual to multiply the equation of motion (Eq. (27)) by the determinant $J$ before performing the Fourier transform, leading to a closed equation of motion cubic in $\mathbf{\Psi}$ Matsubara:2015ipa . Since in our case the gravitational strength is scale-dependent, that approach leads to further complications. We instead use Eq. (37) and expand quantities as

[TABLE]

[TABLE]

where

[TABLE]

We are interested in finding an equation valid up to third order because their solutions lead to the first corrections to the power spectrum (the 1-loop) and correlation function. The inverse of the Jacobian matrix was expanded up to second order because it is already multiplied by the Lagrangian displacement in Eq. (37). Now, using Eq. (21), the matter perturbation is

[TABLE]

Noting from Eq. (43) that $(J^{-1})_{ij}\mathbf{\hat{\mathcal{T}}}\Psi_{i,j}=\mathbf{\hat{\mathcal{T}}}\Psi_{i,i}-\Psi_{i,j}\mathbf{\hat{\mathcal{T}}}\Psi_{i,j}+\Psi_{i,k}\Psi_{k,j}\mathbf{\hat{\mathcal{T}}}\Psi_{i,j}$, and using Eqs. (37) and (46) we find the Lagrangian displacement equation in Fourier space for third order perturbation theory:

[TABLE]

At linear order the right hand side of Eq. (III) is set to zero, leading to the Zel’dovich solution

[TABLE]

where we choose $t_{0}$ to be the present time and $D_{+}(k)$ is the fastest growing solution to

[TABLE]

we further normalize $D_{+}(k,t=t_{0})=1$ at $k=0$. The linear growth function $D_{+}$ depends only on time and the magnitude of k because $A({\text{\bf k}})=A(k)$, which is consistent with the rotational invariance of the underlying theory. This does not happen with higher order growth functions since they depend on the manner modes with different wavelengths interact.

The other solution to Eq. (49), in general decaying, is given by $D_{-}(k,t)$. It will be used in the Green function associated to the linear operator $\mathbf{\hat{\mathcal{T}}}-A(k)$:

[TABLE]

The initial time here is set to $t_{in}$ and a dot means derivative with respect to the time argument. In using Eq. (50) a numerical error may arise if the function at which the operator is applied is not zero at $t_{in}$, but can be compensated if the solution is known at that time, which is the case of models with an early Einstein-de Sitter (EdS) phase, as the one we use below as an example, as well as many others found in the literature. The alternative is set $t_{in}$ to be very small. For $\Lambda$CDM models there is no $k$-dependence and the growth functions become $D_{-}\propto H(t)$ and $D_{+}\propto H^{3}(t)\int^{t}a^{-3}(t^{\prime})H^{-3}(t^{\prime})dt^{\prime}$ 1977MNRAS.179..351H ; while for the special case of EdS, $D_{+}\propto a(t)$ and $D_{-}\propto a^{-3/2}(t)$. In Fig. 1 we plot the ratio of MG to $\Lambda$CDM linear growth functions, this is done for the F4, F5 and F6 models, and for redshifts $z=0$ and $z=0.5$; the background cosmology is fixed with $\Omega_{m}=0.281$ and $h=0.697$, as given by WMAP Nine-year results 2013ApJS..208…19H .

We end this section by writing the scalar field to first order,

[TABLE]

Since the factor $\frac{2A_{0}}{3\Pi(k)}=2a^{2}\frac{(A(k)-A_{0})}{k^{2}}$ is smaller than unity for typical values of $M_{1}$, it follows that the scalar field perturbations evolve at a slower pace than overdensities.

IV Second order: 2LPT

In this section we compute second order quantities and develop the 2LPT theory. The frame-lagging term kernel is obtained from Eq. (38) and the first order fields given in Eqs. (48) and (51),

[TABLE]

while the self interacting term in Eq. (III) becomes

[TABLE]

Hereafter $\delta_{1}$ and $\delta_{2}$ denote the linear density contrasts with wavenumbers ${\text{\bf k}}_{1}$ and ${\text{\bf k}}_{2}$ evaluated at present time.

A straightforward computation of Eq. (III), see Appendix B.1, gives the Lagrangian displacement to second order:

[TABLE]

Momentum conservation implies ${\text{\bf k}}={\text{\bf k}}_{12}\equiv{\text{\bf k}}_{1}+{\text{\bf k}}_{2}$, as it is explicit in the Dirac delta function. The growth functions are given by

[TABLE]

and the normalized growth functions are defined as

[TABLE]

The growth functions depend only on three numbers, these can be chosen as $k_{1}$, $k_{2}$ and $k_{12}$, such that alternatively we can write

[TABLE]

Both notations will be used interchangeably.

The linear operator $(\mathbf{\hat{\mathcal{T}}}-A(k))$ is not invertible over its whole domain, and $(\mathbf{\hat{\mathcal{T}}}-A(k))^{-1}h$ is only a particular solution to the differential equation $(\mathbf{\hat{\mathcal{T}}}-A(k))f=h$. The growth functions obtained from Eqs. (55)-(58) project out the linear order solution to Eq. (49), as it is required for being pure second order functions. On the other hand, numerical computations using differential equations instead, present no disadvantages if the initial conditions are carefully chosen.

For dark energy models in which we can neglect dark energy perturbations or MG scale independent models, we have $I_{2}(t)\equiv\bar{D}^{(2)}_{a}(t)=\bar{D}^{(2)}_{b}(t)$ (following the notation of Lee:2014maa ), and $\bar{D}^{(2)}_{\delta I}=\bar{D}^{(2)}_{\text{FL}}=0$. For EdS, these reduce to $I_{2}(t)=1$, and for $\Lambda$CDM $I_{2}(t)\simeq 1.01$ at the present time for $\Omega_{m0}\simeq 0.3$.

We can write the divergence of the Lagrangian displacement as

[TABLE]

with

[TABLE]

from which we can read the extension to the LPT second order kernel as

[TABLE]

Consider now the collapse of parallel sheets of matter. In GR it is well known that the equation of motion for the Lagrangian displacement becomes linear and the Zel’dovich approximation is the exact solution. This is because the Newtonian gravitational force, which decays as the inverse of the squared distance, becomes a constant regardless the separation of the sheets McQuinn:2015tva . Therefore, the second and higher order Lagrangian displacement kernels are zero. This can be deduced from Eq. (62) by setting $D^{(2)}_{a}=D^{(2)}_{b}$, $D^{(2)}_{\delta I}=D^{(2)}_{\text{FL}}=0$, and choosing parallel wavevectors ${\text{\bf k}}_{1}\parallel{\text{\bf k}}_{2}$, in this situation $D^{(2)}_{\text{GR}}({\text{\bf k}}_{1},{\text{\bf k}}_{2})=0$, as required. In modified theories of gravity the situation is rather different, the fifth force between the sheets is not longer a constant, and the growth function $D^{(2)}({\text{\bf k}}_{1},{\text{\bf k}}_{2})$ does not vanish since $D^{(2)}_{a}({\text{\bf k}}_{1},{\text{\bf k}}_{2})\neq D^{(2)}_{b}({\text{\bf k}}_{1},{\text{\bf k}}_{2})$. Nevertheless, at the scales the fifth force can be neglected we have to recover the GR result, this should be the case in the large scales limit of $k\rightarrow 0$. Since $D_{a}^{(2)}(k=0,k_{1},k_{2})\neq D_{b}^{(2)}(k=0,k_{1},k_{2})$ and $D^{(2)}_{\delta I}(k=0,k_{1},k_{2})=0$, the cancellation should be accomplished by the frame-lagging. Indeed, the squeezed configuration of triangles formed as ${\text{\bf k}}={\text{\bf k}}_{1}+{\text{\bf k}}_{2}$ with ${\text{\bf k}}_{12}\simeq 0$, ${\text{\bf k}}_{1}\simeq{\text{\bf p}}$ and ${\text{\bf k}}_{2}\simeq-{\text{\bf p}}$ is a planar collapse; and, by using Eqs. (52),(55)-(58),(62), and the identity $M_{1}(0)=3\Pi(0)$, it is easy to see that $D^{(2)}({\text{\bf p}},-{\text{\bf p}})=0$, recovering GR when $k^{2}/a^{2}\ll M_{1}$.

As it was mentioned in the Introduction, 2LPT equations in MG have been derived in two (to our knowledge) separate works Valogiannis:2016ane ; Winther:2017jof . In Valogiannis:2016ane the results differ with the given here, since the authors find $D_{a}^{(2)}=D_{b}^{(2)}$ and, furthermore, they do not seem to consider the frame-lagging. In Winther:2017jof their computations to second order coincide with ours. In those works the authors use 2LPT for implementations of MG into the COLA code, obtaining very satisfactory results. These works share the $\Lambda$CDM as a limit at large scales, where the 2LPT in COLA implementations is more important. Indeed, to speed-up computations, in Ref. Valogiannis:2016ane the 2LPT of $\Lambda$CDM is additionally used into the MG runs, finding also a good agreement with $N$-body simulations.

In Fig.2 we plot second order growth functions for the F4 case444The effects of F5, and F6 are qualitatively very similar but are shifted toward higher values of $k$. and we choose three different triangular configurations: (a) equilateral for which $k_{12}=k_{1}=k_{2}$; (b) orthogonal with $k_{12}=\sqrt{2}k_{1}=\sqrt{2}k_{2}$; and (c) squeezed, $k_{12}\simeq 0$, $k_{1}\simeq k_{2}$. We can read from the equations for the growth functions that in the case of equilateral configuration we have $\bar{D}_{a}^{(2)}=\bar{D}_{b}^{(2)}$. The squeezed case corresponds to the limit of large scales where the theory reduces to $\Lambda$CDM, we can see from panel (c) in Fig. 2 that actually by considering the combination $\bar{D}^{(2)}_{a}-\bar{D}^{(2)}_{b}+\bar{D}^{(2)}_{\text{FL}}$ (dotted red line in panel (c)) in this configuration we get $D^{(2)}=0$, as was shown analytically above. The frame-lagging second order growth function for equilateral configuration is exactly zero and for the orthogonal is very similar in shape to the squeezed one, but reaching smaller values. In panel (d) we show the screening normalized growth functions for the three considered triangular configurations. We make note that the functions $\bar{D}_{a}^{(2)}$ and $\bar{D}_{b}^{(2)}$ do not tend to 1 as $k$ goes to zero, the EdS value, instead their limits are $\sim 1.009$, the $I_{2}(t_{0})$ value for the chosen cosmological parameters in the $\Lambda$CDM model.

Now, we calculate the matter and scalar field perturbation to second order. From Eq. (46), $\tilde{\delta}^{(2)}({\text{\bf k}})=[-\Psi_{i,i}^{(2)}+\frac{1}{2}(\Psi_{i,i}^{(1)}\Psi_{j,j}^{(1)}+\Psi_{i,j}^{(1)}\Psi_{j,i}^{(1)})]({\text{\bf k}})$, then

[TABLE]

and the scalar field becomes

[TABLE]

with

[TABLE]

We shall use the above perturbations in the following section to find the third order Lagrangian displacement field.

V Third order

In this section we compute the expansion of the relevant fields up to third order. Using Eqs. (38), (61) and (65), the frame-lagging kernel is

[TABLE]

The third order screening gives

[TABLE]

with kernel

[TABLE]

A straightforward but lengthy computation, see Appendix B.2, leads to

[TABLE]

with normalized growth functions

[TABLE]

and growth functions

[TABLE]

The growth functions depend on three wavevectors, but these are constrained to form a quadrilateral ${\text{\bf k}}={\text{\bf k}}_{1}+{\text{\bf k}}_{2}+{\text{\bf k}}_{3}$; thus, they depend only on 6 numbers for a given k. Furthermore, the relevant configurations for these quadrilaterals in 2-point statistics are the so-called double squeezed in which ${\text{\bf k}}=-{\text{\bf k}}_{1}$ and ${\text{\bf p}}\equiv{\text{\bf k}}_{3}=-{\text{\bf k}}_{2}$, which left us with only 3 degrees of freedom, highly simplifying the analysis.

For the $\Lambda$CDM model we have $D^{(3)}_{A}=D^{(3)}_{B}=(\mathbf{\hat{\mathcal{T}}}-A_{0})^{-1}(\frac{9}{7}D_{+}\mathbf{\hat{\mathcal{T}}}D_{+}^{2})$, and since $\mathbf{\hat{\mathcal{T}}}D_{+}^{2}=2D_{+}\mathbf{\hat{\mathcal{T}}}D_{+}+2\dot{D}_{+}^{2}$, we get $D^{(3)}_{A,B}=\frac{5}{7}6(\mathbf{\hat{\mathcal{T}}}-A_{0})^{-1}(\frac{3}{2}\Omega_{m}H^{2}D_{+}^{3}(\frac{3}{5}+\frac{2}{5}f^{2}/\Omega_{m}))$, with $f=d\ln\delta_{L}/d\ln a$ the growth factor of linear matter perturbations. For $f=\Omega_{m}^{1/2}$ we obtain the result for EdS, $D^{(3)}_{A,B}=\frac{5}{7}D_{+}^{3}$, meaning that the normalized growth functions $\bar{D}^{(3)}_{A,B}$ are exactly 1. While for $\Lambda$CDM we get $\bar{D}^{(3)}_{A,B}(t_{0})\simeq 1.02$ for $\Omega_{m0}\simeq 0.3$. Analogously, we obtain $D^{(3)}_{C,D,E}=6(\mathbf{\hat{\mathcal{T}}}-A_{0})^{-1}(\frac{3}{2}\Omega_{m}H^{2}D_{+}^{3})$, which reduces to $D^{(3)}_{C,D,E}=D_{+}^{3}$ for EdS, while for $\Lambda$CDM we find $\bar{D}^{(3)}_{C,D,E}\simeq 1.02$ at the present time. This analysis shows that the third order displacement field given by Eq. (V) reduces to the standard case when the theory is scale independent.

We can write the third order Lagrangian displacement kernel as

[TABLE]

with

[TABLE]

In Eq.(82) we omit to write a transverse part because it cancels when it is contracted with k.

In 1-loop statistics, the third order growth function should be symmetrized by summing over all permutations —contrary to cyclic permutations which are sufficient in $\Lambda$CDM. By doing this, we obtain

[TABLE]

where in writing “$(\,{\text{\bf p}}\rightarrow-{\text{\bf p}}\,)$”, we assume that $M_{2}({\text{\bf p}},{\text{\bf k}})$ is symmetric in k and p, otherwise it should be symmetrized. It is far from obvious that Eq. (V) goes to its $\Lambda$CDM form as k goes to zero. But by taking the limit case $k=0$ to Eq. (V), several cancellations lead to $D^{(3)symm}(0,-{\text{\bf p}},{\text{\bf p}})=0$, as it should be the case since double squeezed configurations reduce to the 1-dimensional collapse in the limit of $k\rightarrow 0$. The expression for $\mathcal{K}_{\delta I}^{(3)symm}$ is somewhat large and not necessary for recovering GR, thus we do not write it here.

We finish this section by making some comments on the success of the Zel’dovich approximation at the large scales; for further discussion see Tassev:2013pn ; Whi14 ; McQuinn:2015tva ; Cusin:2016zvu . The second order Lagrangian displacement trivially reduces to a planar (one dimensional) collapse when $k\rightarrow 0$ simply because only two modes with ${\text{\bf k}}_{1}$ and ${\text{\bf k}}_{2}$ interact and by conservation of momentum ${\text{\bf k}}_{1}\rightarrow-{\text{\bf k}}_{2}$. For the third order Lagrangian displacement the wavelength vectors should form a quadrilateral ${\text{\bf k}}={\text{\bf k}}_{1}+{\text{\bf k}}_{2}+{\text{\bf k}}_{3}$, again due to momentum conservation; and moreover, statistics as the power spectrum or the correlation function rely on double-squeezed configurations which also reduce to planar collapse in the case $k\rightarrow 0$. Thus, 2-point statistics at large scales essentially probe planar collapses, for which Zel’dovich approximation is exact.

VI Lagrangian displacements 2- and 3-point functions

Lagrangian displacement power spectra and bispectra are 2- and 3-rank tensors, of the form $\langle\Psi_{i}(k_{1})\Psi_{j}(k_{2})\rangle$ and $\langle\Psi_{i}(k_{1})\Psi_{j}(k_{2})\Psi_{k}(k_{3})\rangle$, that we contract with related momenta to construct scalars. In this section we are interested in those combinations that are necessary for matter statistics at 1-loop. These special combinations are defined in Appendix A, Eqs.(118)-(122). Straightforward calculations using the Lagrangian displacements up to third order yield

[TABLE]

with $r=p/k$ and $x=\hat{{\text{\bf k}}}\cdot\hat{{\text{\bf p}}}$. The normalized growth functions in Eqs. (85) and (86) are evaluated as (using the notation of Eq (60))

[TABLE]

and in Eq. (89):

[TABLE]

We emphasize that in the notation of the left hand side of the above equation, the growth functions are symmetric in their arguments. The evaluation of the third order growth function was already given above, in Eq. (V). In this section we consider the F4 model, since their MG effects are more evident than in F5 and F6, but the qualitative features are not affected.

In Fig. 3 we show plots for the $Q_{1}$ and $Q_{2}$ functions considering their different contributions. Terms containing only frame-lagging are shown with black curves, while the red dot-dashed have no frame-lagging contributions, the sum of both gives the full $Q_{1}$ or $Q_{2}$ functions (blue curves), and the orange curves are for the $\Lambda$CDM model. For $Q_{1}$, the cancellation due to the frame-lagging is about one part in a ten thousand at $k=0.001$, showing the importance of these terms at large scales. For $Q_{2}$ the cancellation is about one order of magnitude at the same scale. Analogously, in Fig. 4 we plot the $R_{1}$ and $R_{2}$ functions. In these cases the cancellations are not sufficient to drive these functions to their $\Lambda$CDM counterparts. More detailed inspections of their large scales behavior reveals the reason, as we show below.

To emphasize the importance of frame-lagging terms we consider first the $Q_{1}$ function. By letting the external momentum $k$ goes to zero in Eq. (85) we have

[TABLE]

with $\mathcal{A}=\bar{D}^{(2)}_{a}({\text{\bf p}},-{\text{\bf p}})+\bar{D}^{(2)}_{\text{FL}}({\text{\bf p}},-{\text{\bf p}})$ and $\mathcal{B}=\bar{D}^{(2)}_{b}({\text{\bf p}},-{\text{\bf p}})$. The $\bar{D}^{(2)}_{\delta I}$ growth function vanishes for this configuration, otherwise it would contribute to the function $\mathcal{A}$. Expanding the integrand about $k=0$ and discarding terms with odd powers in $x$ since they become zero after the angular integration, we get

[TABLE]

These squeezed configurations survey large scales and then should reduce to GR. We explicitly showed in Sect. IV that $\mathcal{A}-\mathcal{B}=\bar{D}^{(2)}_{a}({\text{\bf p}},-{\text{\bf p}})-\bar{D}^{(2)}_{b}({\text{\bf p}},-{\text{\bf p}})+\bar{D}^{(2)}_{\text{FL}}({\text{\bf p}},-{\text{\bf p}})=0$, thanks to the frame-lagging. By this reason the leading term that survives in Eq. (VI) is of order $\mathcal{O}((k/p)^{4})$, yielding

[TABLE]

In panel (D) of Fig. 2 we plotted $\bar{D}^{(2)}_{b}({\text{\bf p}},-{\text{\bf p}})$. It starts to depart from GR value early, at yet linear scales, but the integrand above is quickly damped by the $p^{-2}$ factor, thus the net effect is small. For F4 we find $C_{Q_{1}}\simeq 1.03$, for $\Lambda$CDM $C_{Q_{1}}=I^{2}_{2}(t=t_{0})\simeq 1.02$, and for EdS $C_{Q_{1}}$ is exactly one. We note that there is a residual $k$ dependence in $C_{Q_{2}}$ since the external momentum takes small, but non-zero values. On the other hand, if we neglect the frame-lagging, then $\mathcal{A}-\mathcal{B}\neq 0$ and the leading term in the expansion of Eq. (VI) is of order zero and we get

[TABLE]

This explains that in Fig. 3 the function $Q_{1}$ without frame-lagging approaches a constant for low-$k$, and due to the $p^{2}$ term the function takes large values.

Analogous calculations for $Q_{2}$ leads to

[TABLE]

while without the frame-lagging one would obtain $Q_{2}^{\text{FL}}(k\rightarrow 0)\propto k^{2}\int dpP^{2}_{L}(p)$. The other loop functions yield

[TABLE]

The difference in the derivation of $R_{2}$ is that the second order growth functions are evaluated as $\bar{D}^{(2)}(0,-{\text{\bf p}})$, for which frame-lagging terms give small negative contributions, as observed in Fig. 4. The case of $R_{1}$ was guessed and then tested numerically. We have checked the goodness of these approximations numerically, showing that all the $C_{X}$ functions are close to unity, but more relevant for this discussion is that they depend very weakly on $k$. Therefore, the $Q$ and $R$ functions have the same $k$-dependence at large scales in MG and GR.

From Eqs. (94) and (96), we now notice that the large scales modes of $Q_{1}$ and $Q_{2}$ functions receive negligible contributions from small scales, making them indistinguishable for MG models that reduce to $\Lambda$CDM at those scales, as depicted in Fig. 3. On the other hand, the $R_{1}$ and $R_{2}$ functions at low-$k$ (see Eqs. (97) and (98)) reduce to the linear power spectrum times the variance of the Lagrangian displacements, and as a result they take larger values than to those of $\Lambda$CDM, as shown in Fig. 4.

The analysis of the $Q_{3}$ function is simpler because it derives from products of linear displacement fields, thus it has the same form in MG and GR. This function is not necessary for LPT matter statistics, but it is required for computing the SPT power spectrum, as explained in Appendix A. Its large scale limit is

[TABLE]

Now, for scale invariant universes with $P_{L}(p)\propto p^{n}$ as we consider in the following, we fix the external $k$ momentum and let the internal momentum $p$ goes to infinity. The equations derived above in the limit $k\rightarrow 0$ can do this job for power law power spectra, and in this way we observe the $Q$ functions have UV divergences for $n\geq 1/2$ and the $R$ functions for $n\geq-1$; these are the same UV divergences that have the $P_{22}$ and $P_{13}$ pieces of the SPT power spectrum, respectively. We notice that if frame-lagging is not considered, UV divergences for $Q_{1}$ are present for $n\geq-3/2$, while for $Q_{1}$ when $n>-1/2$.

On the other hand, IR divergences may show up both at ${\text{\bf p}}\rightarrow 0$ and ${\text{\bf p}}\rightarrow{\text{\bf k}}$. In Ref. Carrasco:2013sva $P_{22}$ is written in such a way both divergences appear at momentum zero. We follow that work and define

[TABLE]

From Eqs. (114) and (120), $Q_{3}$ can be written as

[TABLE]

In the second equality we have split the region of integration in two pieces separated by the $p\sim k$ divergence. In the second integral of the third equality we redefined the variable ${\text{\bf p}}={\text{\bf k}}-\mathbf{\tilde{p}}$. In the last we use the symmetry $q_{3}({\text{\bf k}},{\text{\bf p}})=q_{3}({\text{\bf k}},{\text{\bf k}}-{\text{\bf p}})$. In this way, $Q_{3}(k)$ has IR divergences only for ${\text{\bf p}}\rightarrow 0$. Fixing the external momentum, a standard computation leads to write, with a little abuse of notation,

[TABLE]

and the leading IR divergence appears for spectral index $n\leq-1$, while sub-leadings divergences reveal for $n\leq-3$.

For $Q_{2}(k)$, we find an IR divergence that does not appear in $\Lambda$CDM. Expanding the angular integrand of Eq. (86) about $r=0$ (or equivalently $p=0$) with power law potential and discarding odd powers of $x$, we get

[TABLE]

Here $\mathcal{A}=\bar{D}^{(2)}_{a}(0,{\text{\bf k}})-\bar{D}^{(2)}_{\delta I}(0,{\text{\bf k}})$ and $\mathcal{B}=\bar{D}^{(2)}_{a}(0,{\text{\bf k}})$, while $\bar{D}^{(2)}_{\text{FL}}(0,{\text{\bf k}})=0$ for this configuration. This is not the same squeezed configuration found above in Eq. (VI), and then $\mathcal{A}-\mathcal{B}\neq 0$. As a result we get

[TABLE]

Although this IR divergence does not appear in $\Lambda$CDM, it is safe since is revealed for $n\leq-3$. For $R_{2}$ and $Q_{1}$ functions the IR divergences when the internal momentum is sent to zero appear for $n<-3$, as in $\Lambda$CDM.

On the other hand for ${\text{\bf p}}\rightarrow{\text{\bf k}}$, functions $Q_{1}$ and $Q_{2}$ present divergences that goes as $(p-k)^{n+2}$ and $(p-k)^{n+1}$, respectively, which are of the same form as in $\Lambda$CDM. These results rely in the approximation of letting fixed the growth functions as $\bar{D}^{(2)}({\text{\bf k}},0)$ as follows from Eq. (90).

For $R_{2}$ as ${\text{\bf p}}\rightarrow{\text{\bf k}}$ we have the evaluation $\bar{D}^{(2)}(-{\text{\bf k}},{\text{\bf k}})$, which is the same squeezed configuration found in Eq. (VI) for which $\mathcal{A}-\mathcal{B}=0$. Thus, we can substitute $\mathcal{A}$ for $\mathcal{B}$ in the integrand of Eq. (89). Since this is a function of $k$ only it can be factorized and pulled out of the integral. Thus $R_{2}$ has no IR divergence for ${\text{\bf p}}\rightarrow{\text{\bf k}}$ as in $\Lambda$CDM. In the absence of frame-lagging $\mathcal{A}\neq\mathcal{B}$ and we obtain a $\log|p-k|$ divergence after expanding Eq. (89) about $p=k$ and performing the angular integration.

The intuition developed in this section suggest that IR and UV divergences for $Q$ and $R$ functions are the same in MG and in $\Lambda$CDM. Though, a rigorous proof is still lacking; in particular we did not analyze the divergences of $R_{1}$ function. We stress that in our derivations the role played by the frame-lagging terms is determinant, and further, we assumed MG models that converge to GR at large scales.

We finalize this section by showing the power spectrum $a(k)\equiv P_{\nabla\cdot\Psi\,\nabla\cdot\Psi}$ of the divergence of Lagrangian displacement fields, $\nabla_{i}\Psi^{i}$, at 1-loop. This is given by $a(k)=k^{i}_{1}k^{j}_{2}\langle\Psi_{i}^{(1)}(k_{1})\Psi_{j}^{(1)}(k_{2})\rangle_{c}^{\prime}+k^{i}_{1}k^{j}_{2}\langle\Psi_{i}^{(2)}(k_{1})\Psi_{j}^{(2)}(k_{2})\rangle_{c}^{\prime}+2k^{i}_{1}k^{j}_{2}\langle\Psi_{i}^{(1)}(k_{1})\Psi_{j}^{(3)}(k_{2})\rangle_{c}^{\prime}$, where the notation $\langle(\cdots)\rangle^{\prime}$ means that we omit a Dirac delta function. Using Eqs. (118) and (121) this is

[TABLE]

where $P_{L}(k)$ coincides with the linear power spectrum of matter fluctuations. In Fig. 5 we show a plot of $a(k)$ with their different contributions. We note that the piece without frame-lagging has negative values.

VII CDM 2-point statistics

In this section we compute the power spectrum and correlation function for F4, F5 and F6 models using the schemes described in Appendix A. The building blocks are the $Q_{i}$ and $R_{i}$ functions already computed in the previous section. As discussed in the introduction, LPT statistics are better for modeling the BAO feature while SPT is better in the broadband power spectrum, and for that reason we explore both approaches here. Specifically, we use the Convolution Lagrangian Perturbation Theory (CLPT) of Ref. Carlson:2012bu , the Lagrangian Resummation Theory (LRT) Mat08a , and the SPT* explained in Appendix A.

For our numerical computations we fix cosmological parameters to $\Omega_{m}=0.281$, $\Omega_{b}=0.046$, $h=0.697$, $n_{s}=0.971$, and $\sigma_{8}=0.82$, corresponding to the best fit of the WMAP Nine-year results 2013ApJS..208…19H . The linear $\Lambda$CDM matter power spectrum is computed with the CAMB code Lewis:1999bs , and the linear MG is obtained by multiplying it by the square of linear growth factor,

[TABLE]

We compute the CLPT power spectrum in Eq. (129) from the Hankel transformations of Eqs. (131), (132) and (133) and using FFTLog numerical methods, as described in Hamilton:1999uv .555In practice, the sums in Eqs. (137) are performed up to some $\ell_{max}$. In the code we use for this work, this maximum is adaptive in the interval $k=0.001$ – $3\,h/\text{Mpc}$ ranging from $\ell_{\text{max}}=10$ to 30. We specifically make use of the public code released in Vlah:2016bcl .666https://github.com/martinjameswhite/CLEFT_GSM., that is supplied with precomputed $\Xi_{\ell}$ functions, defined in Vlah:2015sea , which are combinations of the $X$, $Y$, $V$ and $T$ functions defined in Eqs. (125)-(128).

The LRT power spectrum, written in this section as $P_{\text{LRT}}(k)=e^{-k^{2}\sigma^{2}_{L}}P_{L}(k)+\text{NL}$, with NL denoting the loop integrals, and

[TABLE]

the total variance of the divergence of linear Lagrangian displacements, is obtained directly by providing Eq. (142) with the $Q$ and $R$ functions.

Since LRT and CLPT power spectra decay quickly (in $k$), we expand them to obtain the 1-loop SPT* power spectrum —which in EdS coincides with the 1-loop SPT power spectrum— as explained in the Appendix A. Thus, we analyze also the power spectrum $P_{\text{SPT}^{*}}(k)=(1-\sigma^{2}_{L}k^{2})P_{L}(k)+\text{NL}$ of Eq. (139).

In Figs. 6 and 7 we show the ratio of the different power spectra and the linear power spectrum of the considered model with the BAO removed, this is done for CLPT, LRT and SPT*, and for F4, F5, F6 and $\Lambda$CDM models at redshift $z=0$. Dashed blue curves correspond to the CLPT scheme, solid red are for the SPT*, brown dashed for LRT, and long-dashed black for the linear theory. As a reference, for MG models we plot the linear $\Lambda$CDM model as well (depicted by dotted green curves). The vertical lines denote the expected maximum wavelengths $k_{max}$ of validity of the perturbative Lagrangian formalisms. For these, we adopt the prescription in Ref. Mat08a relying on the damping scale, such that $k_{max}=\sigma_{\text{L}}^{-1}/2$; other prescriptions are possible, but do not differ substantially. We note that in MG models non-linearities starts to be more relevant at larger scales since clustering is more efficient due to the extra force.

LPT formulations perform better for the correlation function; it was indeed noted that even the Zel’dovich approximation provides very good accuracy about the BAO peak; see for example Tassev:2013rta ; Whi14 for recent discussions. We compute the CLPT correlation functions given by Eq. (138). To do it we feed the code of Vlah:2015sea 777https://github.com/alejandroaviles/CLEFT. with precomputed $Q_{i}$ and $R_{i}$ functions. We show the results in Fig. 8, with the lower panel showing the relative difference when comparing to the $\Lambda$CDM Zel’dovich approximation in the $\Lambda$CDM, which is computed from Eq. (131). At the BAO scale, the differences in the correlation among the three MG and $\Lambda$CDM models are smaller than $2\%$ since at these scales the fifth force is almost negligible. Soon below the BAO scales the differences depart considerably for F4, which is not surprising since this model is ruled out even at the linear level Hojjati:2015ojt .

In Figs. 9 and 10 we show the SPT* power spectra for two redshifts $z=0$ and $z=0.5$ with and without screenings. As a reference we use $N$-body simulations data Liaknowled performed with the ECOSMOG code of Ref. Li:2011 . The presented data are the average of two realizations in boxes of size $1024\,\text{Mpc}/h$ and with $1024^{3}$ dark matter particles. It is difficult to tell what is the scale at which SPT* is valid since the data is noisy at large scales,888We stress that $N$-body codes for MG are numerically quite expensive and are aimed mainly to test the small scales. Thus they use relatively small boxes and large scales turn out to be highly affected by cosmic variance. but we expect to have a good accuracy up to the $k_{max}$ given above. Nevertheless, we can note that the trend line of the data is properly followed by the analytical method.

Since small scales structures are subject to negligible tidal forces from bulk flows provided by long wavelength fluctuations, these are only affected by an overall translation which has no net effect in their matter $N$-point statistics. The flaw of LRT is that bulk flows are split in $q$ dependent and independent pieces [see Eq. (141)] and does not cancel at zero Lagrangian separation. Indeed, in the exponential prefactor of Eq. (142) we find the term $\sigma^{2}_{L}$, that has an IR divergence for power law power spectrum with spectral index $n\leq-1$, as a consequence the bulk flows largely affect the small scales statistics, this is a manifestation of the breaking the Galilean invariance (see for example Sugiyama:2013mpa for recent discussion in LRT). This effect is observed in Fig. 6 and 7 where a deficit in the power spectrum with respect to the CLPT results is quite evident. Instead, in CLPT all linear contributions are kept in the exponential [see Eq. (129)], including the term

[TABLE]

which is IR-safe for $n>-3$.999For power spectrum with $n\leq-3$, tidal forces from long modes can not be neglected at any scale Carrasco:2013sva .

The SPT* power spectrum is usually written as $P_{\text{SPT*}}(k)=P_{L}(k)+P_{22}(k)+P_{13}(k)$, in terms of $Q$ and $R$ functions these are

[TABLE]

The loop integrals have in general IR divergences for power law power spectrum with spectral index $n<-1$, but they cancel out exactly when the sum is considered. This is a consequence of the discussion above but is more general due to the equivalence principle Carrasco:2013sva ; Creminelli:2013mca . We have already found in Eq. (102) IR divergences in $Q_{3}$ when $p\rightarrow 0$ and $p\sim k$, which are twice that of $\sigma^{2}_{L}P_{L}$, but in Eq. (109) $Q_{3}$ is multiplied by a $1/2$ factor; thus, when summing up the two SPT leading divergences cancel. This is an old known result since the work of 1983MNRAS.203..345V , and afterwards noted in 1996ApJ…456…43J ; Scoccimarro:1995if that the same should happen for sub-leading divergences; we indeed have found other sub-leading divergences for $R$ and $Q$ functions. Our reasoning of Sect. VI suggests that the additional, sub-leading IR divergences with power spectral index $n>-3$ are the same as those of $\Lambda$CDM, making the theory IR-safe.

Even for MG theories with universal couplings, the principle of equivalence can be violated since some bodies may develop screening against the additional force while some others may not Hui:2009kc , but for long wavelength fluctuations this is irrelevant as long as they essentially exert only the Newtonian force upon the small structures. A brute force approach is elusive, but the absence of IR divergences (for $n>-3$) in MG theories that reduce to GR at large scales is in principle ensured by the equivalence principle Valageas:2013cma .

VIII Conclusions

In this work we found a generalization to alternative theories of gravity for the Lagrangian displacement field up to third order in perturbation theory. In doing this, an LPT theory for the study of large scale structure formation was developed. Our formalism is suitable for Horndeski models Bose:2016qun , this is the same range of validity that share several investigations in SPT that came after to the pioneering work of Ref. 2009PhRvD..79l3512K . The basic requirement is that the theory can be recasted as a scalar tensor gravity at linear order, all non linear contributions to the Klein-Gordon equation, including possible kinetic terms, can be written as sources and treated perturbatively. In theories that have a screening mechanism these contributions are responsible to the recover GR at small scales, by making the Laplacian of the scalar field negligible.

Since in LPT is standard to work in $q$-Fourier space, one needs to transform the Poisson equation accordingly. In these theories the Poisson equation has two sources: the standard source of matter perturbations can be written in terms of Lagrangian displacement straightforwardly due to mass conservation, but the terms arising from the scalar field must be properly transformed. By doing it, terms to compensate spatial derivatives in Eulerian space appear into the theory, we call these geometrical contributions frame-lagging terms, and have a crucial role in the LPT for MG, which is more evident at large scales where the theory should reduce to GR. We have discussed this limit from different perspectives. A special one is the 1-dimensional collapse for which the Zel’dovich approximation is exact (up to shell-crossing) in $\Lambda$CDM, while in MG this is not the case since the force mediated by the scalar field does not decay as the square of the inverse of the distance. But at scales larger than the range of the fifth force the only force in play is the standard Newtonian force and we have to recover Zel’dovich as the exact solution, or in other words, the Lagrangian displacement growth functions at orders higher than 1 should vanish. The formalism presented here is capable to capture this physical fact, mathematically through cancellations provided by the frame-lagging contributions.

Throughout this article we apply the formalism by using the Hu-Sawicky $f(R)$ theory Hu:2007nk , more precisely to the so-called F4, F5, and F6 models. We compute 2-point statistics both for the Lagrangian displacement and for the matter perturbations. For the matter power spectrum we use recent $N$-body simulations Liaknowled as a reference. LPT is very successful in reproducing the acoustic peak in the correlation function, and we use two different resummation schemes: CLPT and LRT. Though, LPT is not satisfactory in following the broadband power spectrum. For this reason, we also use a scheme that reduces to SPT for the EdS case Mat08a —in MG it is not clear if one recovers SPT exactly since frame-lagging terms do not seem to cancel out— showing good agreement between the analytic theory and simulations.

We further discussed UV and IR divergences in matter and Lagrangian displacement statistics. Our analysis is not exhaustive, but it suggests that the schemes of SPT* and CLPT are IR-safe in MG.

Acknowledgements.

The authors would like to thank Baojiu Li for useful discussions and for sharing the $N$-body simulation data. We also thank to Zvonimir Vlah and Matteo Cataneo for suggestions and discussions. A.A. is supported by Cátedras CONACyT project No.574. The authors acknowledge financial support from CONACyT Project 269652 and Fronteras Project 281.

Appendix A Lagrangian displacement polyspectra and the matter power spectrum

This appendix briefly reviews the methods we use to compute the power spectrum and the correlation function. For completeness we write the relevant equations but we refer the reader to the original papers for detailed derivations.

The matter power spectrum in LPT is given by TayHam96

[TABLE]

where $\Delta_{i}=\Psi_{i}({\text{\bf q}}_{2})-\Psi_{i}({\text{\bf q}}_{1})$ are the Lagrangian displacement differences, and ${\text{\bf q}}={\text{\bf q}}_{2}-{\text{\bf q}}_{1}$ the Lagrangian coordinate difference. By using the cumulant expansion theorem we can recast the power spectrum as Carlson:2012bu

[TABLE]

where $A_{ij}=\langle\Delta_{i}\Delta_{j}\rangle_{c}$ and $W_{ijk}=\langle\Delta_{i}\Delta_{j}\Delta_{k}\rangle_{c}$. Here we show terms up to third order in the Lagrangian displacement since we are interested in the first, 1-loop corrections to the matter two point functions. The products $\Delta_{i}\Delta_{j}$ have contributions evaluated at the same point, the so-called zero-lag, and contributions given at a separation q.

The direct computation of the LPT power spectrum of Eq. (112) is difficult since it involves highly oscillatory integrands. This has been done in Vlah:2014nta by using expansions in spherical Bessel functions (see below) and in Sugiyama:2013mpa by similar methods involving Legendre expansions of the Lagrangian displacement differences. In part because of these numerical complications, alternative resummation schemes that Taylor expand some terms in the exponential of Eq. (112) exist in the literature. We below review those we use in this work.

We first define the polyspectra $C^{(n_{1}\cdots n_{N})}_{i_{1}\cdots i_{N}}({\text{\bf k}}_{1},\dots,{\text{\bf k}}_{N})$ at order $(n_{1}+\cdots+n_{N})$ as in Mat08a

[TABLE]

which in terms of the Lagrangian displacement kernels can be written as

[TABLE]

Higher order polyspectra do not enter in 1-loop calculations. The following scalar functions are constructed

[TABLE]

Notice these functions are all of order $\mathcal{O}(P_{L}^{2})$. By using decomposition in vectors $\hat{q}_{i}$ the spectra and bispectra of Lagrangian displacement differences take the form Carlson:2012bu ; Vlah:2014nta

[TABLE]

with

[TABLE]

where $P_{L}(k)$ is the linear matter power spectrum. (In general “$L$” denotes a linear piece and “loop” a pure 1-loop piece.) These functions vanish as Lagrangian coordinates separation goes to zero, implying that the power spectrum does not have contribution for arbitrary small scales. This is consistent with Eq. (111).

The formalism of Convolution Lagrangian Perturbation Theory (CLPT) as presented in Vlah:2015sea relies on expanding the loop contributions of Eq. (112) while keeping the linear terms in the exponential,101010CLPT was introduced in Carlson:2012bu where also the $A_{ij}^{\text{loop}}$ was kept exponentiated. leading to

[TABLE]

This can be separated in Zel’dovich and loops contributions as

[TABLE]

with (for ${\text{\bf k}}\neq 0$)

[TABLE]

By using the expansions Vlah:2014nta ; Vlah:2015sea

[TABLE]

the integrals in Eqs.(131), (132) and (133) reduce to one-dimensional Hankel transformations.

In the CLPT scheme, the correlation function can be written compactly by Fourier transforming the power spectrum of Eq. (129) and using Gaussian integrations, resulting in

[TABLE]

where $G_{ij}=A^{-1}_{ij}-g_{i}g_{j}$, $\Gamma_{ijk}=A^{-1}_{ij}g_{k}+A^{-1}_{jk}g_{i}+A^{-1}_{ki}g_{j}-g_{i}g_{j}g_{k}$, and $g_{i}=A^{-1}_{ij}(r_{j}-q_{j})$.

Now, a straightforward expansion of the exponential in Eq. (129) leads to Vlah:2014nta ; Vlah:2015sea

[TABLE]

where

[TABLE]

is the 1-dimensional Lagrangian displacement variance. The power spectrum of Eq. (139) coincides with the SPT power spectrum in the $\Lambda$CDM case. Note that this may not be the case for MG, since it is not clear that the frame-lagging terms cancel out, and SPT should be free of these terms. Nevertheless, the behavior the power spectrum in this scheme is similar to what is expected from an SPT theory, then we use it in Sect. VII and denote it as SPT*.

Other scheme widely used in the literature is the Lagrangian Resummation Theory (LRT) of Matsubara Mat08a , here the $A_{ij}$ matrix splits as

[TABLE]

and the zero-lag term $X(q\rightarrow\infty)=2\sigma^{2}$ is kept in the exponential while the rest is expanded, leading to

[TABLE]

(This equation coincides with Eq. (35) of Mat08a .) By a further expansion of the exponential we arrive back at Eq. (139).

The exponential prefactor in Eq. (142) is responsible for the smearing of the BAO peak, and leads to the suppression of the power spectrum at small scales. Note also that the loop variance $\sigma^{2}_{\text{loop}}$ is neglected in the exponential since it is already a $\mathcal{O}(P_{L}^{2})$ quantity. If it is considered, the corrections become large and the BAO peak is over-suppressed Vlah:2014nta .

We note that LRT is inconsistent since $A_{ij}$ should vanish at small separations. This is a manifestation of the breaking of Galilean invariance in the Matsubara formalism, and it is discussed in Sect.VII. Nevertheless, this scheme has the advantages to be computationally simpler, it is straightforward to add biased tracers and RSD Matsubara:2008wx , and furthermore, its correlation function shows very good agreements with $N$-body simulations. For these reasons, we also consider it in this work.

Appendix B Computation of growth functions

B.1 Second order

To second order, the equation of motion (Eq. (III)) can be written as

[TABLE]

We compute the result for each source $S_{i}$:

Source $S_{1}$:

[TABLE]

In the fourth equality we used the linear order equation $\mathbf{\hat{\mathcal{T}}}D_{+}(k_{2})=A(k_{2})D_{+}(k_{2})$, and in the fifth equality we symmetrized the integrand. We get for the displacement field

[TABLE]

with

[TABLE]

Source $S_{2}$:

[TABLE]

and the displacement field becomes

[TABLE]

with

[TABLE]

Source $S_{3}$ corresponds to the screening and it leads to

[TABLE]

with

[TABLE]

Source $S_{4}$ corresponds to the frame-lagging source, leading to

[TABLE]

with

[TABLE]

Using Eqs. (145),(146),(148),(149),(150),(151),(152), and (153), and rearranging terms we arrive to Eq. (54).

B.2 Third order

The equation of motion to third order is given by

[TABLE]

We compute first for source $S_{1}$:

[TABLE]

then

[TABLE]

with

[TABLE]

Similarly for the other sources:

$S_{2}=[\Psi^{(1)}_{i,j}\mathbf{\hat{\mathcal{T}}}\Psi^{(2)}_{j,i}]({\text{\bf k}})$:

[TABLE]

with

[TABLE]

$S_{3}=-A(k)[\Psi^{(2)}_{i,j}\Psi^{(1)}_{j,i}]({\text{\bf k}})$:

[TABLE]

with

[TABLE]

$S_{4}=-A(k)[\Psi^{(1)}_{i,j}\Psi^{(2)}_{j,i}]({\text{\bf k}})$:

[TABLE]

with

[TABLE]

Source $S_{5}=-[\Psi_{i,k}^{(1)}\Psi_{k,j}^{(1)}\mathbf{\hat{\mathcal{T}}}\Psi_{j,i}^{(1)}]({\text{\bf k}})$:

[TABLE]

with

[TABLE]

$S_{6}=\frac{A(k)}{3}[\Psi_{i,k}^{(1)}\Psi_{k,j}^{(1)}\Psi_{j,i}^{(1)}]({\text{\bf k}})$:

[TABLE]

with

[TABLE]

$S_{7}=\frac{A(k)}{6}[\Psi_{i,i}^{(1)}\Psi_{j,j}^{(1)}\Psi_{k,k}^{(1)}]({\text{\bf k}})$:

[TABLE]

with

[TABLE]

$S_{8}=\frac{A(k)}{2}[\Psi_{l,l}^{(1)}\Psi_{i,j}^{(1)}\Psi_{j,i}^{(1)}]({\text{\bf k}})$:

[TABLE]

with

[TABLE]

$S_{9}=\frac{k^{2}/a^{2}}{6\Pi(k)}\delta I^{(3)}({\text{\bf k}})$:

[TABLE]

with

[TABLE]

$S_{10}=\frac{M_{1}}{3\Pi(k)}\frac{1}{2a^{2}}[(\nabla^{2}_{\text{\bf x}}\varphi-\nabla^{2}\varphi)]^{(3)}({\text{\bf k}})$:

[TABLE]

with

[TABLE]

Using Eqs. (156)-(175) and rearranging terms we arrive to Eq. (V).