# Lagrangian perturbation theory for modified 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 CDM model in that limit. We apply our formalism to compute the 1-loop power spectrum and the correlation function in gravity by using different resummation schemes. We further discuss the IR divergences of these formalisms.

###### pacs:

PACS## 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 1m to 1 AU m (1); (2). Nevertheless, mainly since the discovery of the accelerated expansion of the universe (3); (4), 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 (5); (6); (7); (8); (9), 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 CDM model, that can explain the cosmic background radiation anisotropies with exquisite precision (10), the late time background evolution probed by measurements of the baryon acoustic oscillations (BAO) and supernova distance-redshift relations (11), among several other observations. But this success comes with a price, the side of the theory suffers from the smallness and coincidence problems (12); (13). 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 (14); (15). 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 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 (16); (17). 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 (18). 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 (19); (20), kinetic (21), or Vainshtein (22) 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. -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 -body simulations have been provided for a few MG theories in (23); (24); (25); (26); (27); (28), while semi-analytical--body codes that rely on screening for spherical configurations (29), and the COLA approach, introduced in (30) for CDM, were implemented recently in (31); (32). 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 (33); (34); (35) for reviews in standard cosmology, and Refs. (36); (37); (38) 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) (39); (40); (41); (42); (43); (44); (45); (46); (47); (48); (49) and Eulerian Standard Perturbation Theory (SPT) (50); (51); (52); (53); (54); (55), but so many others have been tailored in order to improve its performance, for example see Refs. (56); (57); (58); (59). 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 (60) based in the formalism of the closure theory of Ref. (58). 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 (61); (62); (63); (64); (65); (66); (67); (68); (69). 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 (70). 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 (48).

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 model (71). 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. (72); (73); (74); (75); (76). 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 (31); (32) in order to provide a framework for dealing with large scales in the COLA code. More recently, the 2LPT of (32) was extended to include massive neutrinos in (77). 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 (36); (38); (37), from which we extract the models we are interested in.

We may consider general scalar-tensor theories (78); (79); (80),

(1) |

where and are functions that can be different for various theoretical frameworks. A conformal transformation decouples matter from the scalar field, and therefore in this Jordan Frame point particles follow geodesics of the new metric. One can further redefine and to have a non-minimal coupling given by

(2) |

which reproduces the Brans-Dicke theory when is a constant and . 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 gravity (81); (82), in which one replaces the Einstein-Hilbert Lagrangian density by a generic function of the Ricci scalar, such that

(3) |

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 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 (83); (84).

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 (85); (37),

(4) |

where the function stands for possible derivative self-interactions of the scalar field, its a potential, and is the trace of the matter stress-energy tensor. The coupling yields a scalar field solution that depends on the local density; for non-relativistic matter . 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 (86)— satisfying known observational constraints but allowing smoking guns to be discriminated or confirmed by near-future cosmological data, e.g. (7).

The screening of the fifth force at local scales can be realized when considering fluctuations around background values with mass and due to the presence of the functions or . According to Ref. (85), the screening can be classified as due to: i) Weak coupling, in which the function 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 (87); (88); (89) and varying dilaton (90); (91); 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 (19); (20); and iii) Large inertia, when the kinetic function 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 -mouflage models (21); (92), and of the latter are Vainshtein models (22).

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 (60). 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
(62); (93); (94); (95), 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 (94) and numerical simulations in (95),
and in a general framework, introducing the effective sound speed () of the gravity model in (62); this later
work concludes that the quasi-static approximation depends, not on the wavenumber in units of the cosmological horizon, , but on the combination ,
which must be greater than unity.^{1}

(5) | ||||

(6) |

where is the background matter energy density, its perturbation and 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

(7) |

Here the self-interaction term may be expanded following (60) as with

(8) |

where the functions 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 theories for which the are only time dependent, as we see below.

Now we consider in more detail gravity, that was first used to explain the current cosmic acceleration in (96); (97). Variations with respect to the metric of the action constructed with the Lagrangian density of Eq. (3), lead to the field equations

(9) |

where . By taking the trace to this equation we obtain

(10) |

where we use a dust-like fluid with . In cosmological pertubative treatments one splits and , where the bar indicates background quantities and . Necessary conditions to get a background cosmology consistent with the CDM model are and (71).

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

(11) |

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 and Brans-Dicke parameter . Since , we may expand the potential as

(12) |

from which the mass associated to the perturbed field can be read as

(13) |

In this work we focus our attention to the Hu-Sawicky model (71), defined by

(14) |

where the energy scale is chosen to be . In this parametrized model, at high curvature () the function 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 for the solutions are stable and the mapping to a scalar tensor gravity is allowed (98). 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 CDM model it is also necessary that (71), thus leaving two parameters to fix the model. Noting that the background value is about , Eq.(14) simplifies and leads to

(15) |

with

(16) |

We consider the case . Using Eqs.(12) and (15) we find the functions , and :

(17) |

(18) |

(19) |

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 functions may carry scale dependences; see, for example, (60) for the functions in DGP gravity.

With fixed to unity, the only free parameter is . The models constructed with it are named F4, F5, F6 and so on, corresponding to the choices 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 (71) 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 CDM model with and .

In this work we specialize to theories. Even though, by dealing consistently with the scale dependence of the
functions, our results cover a wide range of models, essentially the whole Horndeski sector. Indeed, in appendix B of Ref. (67) it is
shown how to relate the functions to the Horndeski free functions.^{2}

## 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

(20) |

with the Lagrangian map vanishing at the initial time where Lagrangian coordinates are defined. In kinetic theory we define the Lagrangian displacement velocity 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 follows (99). Furthermore, mass conservation implies the relation between Lagrangian displacements and overdensities

(21) |

where the Jacobian of the transformation is

(22) |

and 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

(23) |

(24) |

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

(25) |

The notation 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 -space, therefore we shall work in -Fourier space, where wavenumbers correspond to Lagrangian coordinates. For the sake of brevity, in the following we will omit the time argument when this not lead to confusions. The term , introduced in (60), is the expansion of the scalar field potential about its background value , being the linear piece and 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 is very large. On the other hand, theories which poses a screening mechanism shield the scalar field in nonlinear regions such that , recovering GR. Following (49), we define the operator

(26) |

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

(27) |

where is partial derivative with respect to . With this splitting of the Laplacian we recognize the scalar field as a function of Lagrangian coordinates, the term 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 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

(28) | ||||

(29) |

with the fully antisymmetric Levi-Civita symbol. In perturbation theory the relevant fields are formally expanded as

(30) |

and analogously for and . From now on the control parameter is absorbed in the definitions of the perturbed fields. Since spatial derivatives in Lagrangian and Eulerian coordinates are related by , the frame-lagging term of Eq. (27) can be written as

(31) |

implying it is a non-linear term. The Fourier transform of leads to

(32) |

where means Fourier transform of and we defined

(33) | ||||

(34) | ||||

(35) |

is the gravitational strength in the MG cases, while 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 to make clear that the Fourier transform of is taken with Lagrangian coordinates: it is neither the -Fourier transform of nor the -Fourier transform of . It is given instead by , or using Eq. (21) this is

(36) |

Only for the linear overdensities we have .

The function is related to the mass of the scalar field by , which gives the range of the interaction , if it is finite we recover GR at large scales. At sufficiently small scales, such that , we note , for which solar system observations restrict (1), making the theory effectively indistinguishable from GR in the absence of non-linear terms. For example, in gravity the scalar field perturbation is identified with as it is discussed in Sec.III, obtaining for the Brans-Dicke parameter as can be seen directly from Eq.(11). Therefore, in this situation which would rule out the theory; nevertheless, the nonlinearities of the potential (encoded in in perturbation theory) are responsible to drive the theory to GR in that limit (71). On the other hand, the large scale behavior is dictated by the mass , related to 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

(37) |

and up to third order the frame-lagging term is

(38) |

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
follows^{3}

(40) |

Below we give expressions for kernels . By using Eq. (32) iteratively order by order we can write Eq.(II) as

(41) |

For compactness we introduced the function

(42) |

In LPT it is usual to multiply the equation of motion (Eq. (27)) by the determinant before performing the Fourier transform, leading to a closed equation of motion cubic in (49). 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

(43) |

(44) |

where

(45) |

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

(46) |

Noting from Eq. (43) that , and using Eqs. (37) and (46) we find the Lagrangian displacement equation in Fourier space for third order perturbation theory:

(47) |

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

(48) |

where we choose to be the present time and is the fastest growing solution to

(49) |

we further normalize at . The linear growth function depends only on time and the magnitude of k because , 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 . It will be used in the Green function associated to the linear operator :

(50) |

The initial time here is set to 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 , 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 to be very small. For CDM models there is no -dependence and the growth functions become and (100); while for the special case of EdS, and . In Fig. 1 we plot the ratio of MG to CDM linear growth functions, this is done for the F4, F5 and F6 models, and for redshifts and ; the background cosmology is fixed with and , as given by WMAP Nine-year results (101).

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

(51) |

Since the factor is smaller than unity for typical values of , 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),

(52) |

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

(53) |

Hereafter and denote the linear density contrasts with wavenumbers and evaluated at present time.

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

(54) |

Momentum conservation implies , as it is explicit in the Dirac delta function. The growth functions are given by

(55) | ||||

(56) | ||||

(57) | ||||

(58) |

and the normalized growth functions are defined as

(59) |

The growth functions depend only on three numbers, these can be chosen as , and , such that alternatively we can write

(60) |

Both notations will be used interchangeably.

The linear operator is not invertible over its whole domain, and is only a particular solution to the differential equation . 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 (following the notation of (102)), and . For EdS, these reduce to , and for CDM at the present time for .

We can write the divergence of the Lagrangian displacement as

(61) |

with