# Renormalized Cosmological Perturbation Theory

###### Abstract

We develop a new formalism to study nonlinear evolution in the growth of large-scale structure, by following the dynamics of gravitational clustering as it builds up in time. This approach is conveniently represented by Feynman diagrams constructed in terms of three objects: the initial conditions (e.g. perturbation spectrum), the vertex (describing non-linearities) and the propagator (describing linear evolution). We show that loop corrections to the linear power spectrum organize themselves into two classes of diagrams: one corresponding to mode-coupling effects, the other to a renormalization of the propagator. Resummation of the latter gives rise to a quantity that measures the memory of perturbations to initial conditions as a function of scale. As a result of this, we show that a well-defined (renormalized) perturbation theory follows, in the sense that each term in the remaining mode-coupling series dominates at some characteristic scale and is subdominant otherwise. This is unlike standard perturbation theory, where different loop corrections can become of the same magnitude in the nonlinear regime. In companion papers we compare the resummation of the propagator with numerical simulations, and apply these results to the calculation of the nonlinear power spectrum. Remarkably, the expressions in renormalized perturbation theory can be written in a way that closely resembles the halo model.

## I Introduction: improving perturbation theory

The growth of density and velocity perturbations is a fundamental result of cosmology that can be used to test basic properties of the universe, such as the amount of dark matter and dark energy. Although the linear and weakly non-linear growth of perturbations is well understood, giving the large-scale asymptotics (tree-level amplitudes in diagrammatic perturbation theory language [1]) of the power spectrum and higher-order statistics as a function of time, the nonlinear corrections to these results (“loop corrections” [2]) are less well understood (see [3] for a review).

The standard approach of perturbation theory (hereafter PT) is to solve the equations of motion by expanding in powers of the perturbation amplitudes. This is well justified when doing tree-level PT, as one is essentially looking for asymptotic behavior (the large-scale limit, where perturbation amplitudes become very small compared to unity), however as one approaches the nonlinear scale and loop corrections become increasingly important, the validity of such a method becomes less clear. Part of the issue is that cosmological PT is not a standard PT in the sense of having a small coupling constant; its validity depends on the scale, redshift, type of initial perturbation spectrum and statistic under consideration. For example, for initial spectra with significant large-scale power (see Fig. 13 in [3]) or at high redshift for CDM models, including just one-loop corrections to the linear power spectrum gives very good agreement with numerical simulations. On the other hand, for redshifts when the effective spectral index at the nonlinear scale becomes one-loop PT deviates up to for (see Fig. 6 in [4]). Similar results hold for the bispectrum [5].

The purpose of this paper is to cast PT in a different form, which helps to understand the physical meaning of different contributions in the infinite perturbation series, following ideas initially put forward in [6]. In particular we show that different class of diagrams organize themselves into a few characteristic quantities, the most important of which (in terms of improving PT) is the non linear propagator that measures the memory of perturbations to their initial conditions as a function of scale. When an infinite number of diagrams is resummed into this quantity, the remaining (infinite number of) diagrams form a well-defined expansion where each term with loops corresponds to the effects of -mode-coupling, and dominates in a narrow range of scales.

To illustrate these ideas and hopefully motivate the reader to go through the rest of the paper (which is rather technical for those not familiar with field theory methods), a quick calculation of the nonlinear power spectrum in the Zel’dovich approximation provides a good understanding of the difference in behavior at nonlinear scales between (standard) PT and the renormalized perturbation theory (RPT hereafter) we develop in this paper. In the Zel’dovich approximation [7], the trajectories of particles away from their initial (Lagrangian) positions q are given by where is the displacement field. The matter density perturbation, defined as , can be written in Fourier space as,

(1) |

by means of the mass conservation relation . This gives for the power spectrum, defined as ,

(2) |

where we introduced , and used that from Eq. (1). In the Zel’dovich approximation the displacement satisfies , where stands for the linear density perturbation in Fourier space. Moreover, for Gaussian initial conditions is a Gaussian random field and thus , which can be used to cast Eq. (2) as,

(3) |

where and is the variance of the displacement field (and also the one-dimensional velocity dispersion in linear theory). Equation (3) is the well-known relationship between the nonlinear and linear power spectrum in the Zel’dovich approximation [8]. Although this provides a poor description of the power spectrum (with the nonlinear power being smaller than linear), having an exact solution such as Eq. (3) serves to illustrate the differences between PT and RPT. In PT, one basically expands in the amplitude of the linear spectrum, therefore

(4) |

where the term with gives the linear spectrum (tree-level, zero loops, i.e. ), gives one-loop () corrections, and so on. The left panel in Fig. 1 shows the for , for a CDM model with normalization and shape parameter , with nonlinear scale . Solid lines denote positive contributions, dashed lines negative ones. Note that at large enough scales only the linear spectrum () contributes, as expected, however as the nonlinear scale is approached different loop corrections become of the same order, with significant cancellations among them. Similar behavior has been seen in the case of PT for the exact dynamics, at least for some initial spectra, see e.g. [9].

In RPT, as we discuss in detail below, one attempts (among other things) to sum an infinite class of diagrams for the propagator, which (as will be shown in section III.3.2) correspond to the factor of in the Zel’dovich approximation, Eq. (3). The remaining infinite series corresponds to the result

(5) |

which is shown in the right panel in Fig. 1 for . We see that as a result of this partial resummation, the behavior of the resulting perturbation expansion is drastically altered: successive terms dominate at increasingly smaller scales, and there are no cancellations among them (all contributions are positive). Moreover, as we discuss below, RPT provides a physical understanding of the different contributions: the sum represents the contribution of -modes coupling (the larger the generation of power moves to smaller scales), the decaying exponential factor represents how the amplitude and phases of the Fourier modes differ from linear evolution from the initial conditions. Indeed, it is easy to see that Eq. (5) can be rewritten as

(6) |

where , and represents the PT kernels [see Eqs. (10-13) for a definition] in the Zel’dovich approximation [see Eq. (36) below]. The sum corresponds to all mode-coupling diagrams (with the appropriate weighting), with all “tadpole” diagrams in the language of standard PT [10] being summed to the exponential prefactor. This propagator renormalization factor is given by the square of , where represents here a functional derivative ^{1}^{1}1In the sections below we use to denote functional derivatives when there is no possibility of confusion with the density field perturbations. and is the linear growth factor. The exponential decay is caused by the deviations from linear evolution (see section III.3.2 for a derivation), and the decay length (here ) defines a characteristic scale at which individual Fourier modes do not evolve as in linear PT. The negative sign contributions in PT (see left panel in Fig. 1) are an artifact of the power series expansion of propagator renormalization; as increases these terms become large and unless one sums up the infinite series this leads to a breakdown of PT. Note that the term in Eq. (6) gives a contribution to which is proportional to the linear spectrum, i.e. , therefore if the linear spectrum were cutoff at some scale [ for ] the nonlinear effects encoded in propagator renormalization do not lead to power generation at . It is in this sense that encodes mode-mode coupling.

In summary, the identification in RPT of physical quantities such as the propagator in the perturbation series helps to cure the lack of convergence of PT in the nonlinear regime. Notice also that the behavior of RPT in Fig. 1 is similar to the halo model description of the power spectrum, the similarities are explored below in section V.

The rest of this paper is organized as follows. In the next section we review the standard results from PT and describe an alternative description necessary to develop RPT. Section III applies the result in the previous section to the calculation of statistics, in particular the power spectrum and the nonlinear propagator. Section IV derives the main results of RPT, and section V discusses the connection between RPT and the halo model. The conclusions are presented in section VI.

## Ii Dynamics

### ii.1 In Standard Form

For completeness, we review here the standard form of the equations of motion in Fourier space and their solutions in standard PT, before going into the new approach that helps us understand better the different terms in the perturbation series and thus facilitate the process of resummation.

We study the evolution of Cold Dark Matter (CDM) in the single-stream approximation, where the relevant equations of motion are conservation of mass, momentum, and the Poisson equation (see e.g. [11]). Consistent with this, the velocity field is assumed to be irrotational, then gravitational instability can be fully described in terms of the density contrast and the peculiar velocity divergence . Defining the conformal time , with the cosmological scale factor, and the conformal expansion rate , the equations of motion in Fourier space become [12, 13, 14]

(7) |

(8) |

The non linear interactions (in the right-hand side) are responsible for the coupling between different Fourier modes, given by

Equations (7) and (8) are valid in an arbitrary homogeneous and isotropic universe, which evolves according to the Friedmann equations. For simplicity here we will restrict ourselves to an Einstein-de Sitter model for which and . As a consequence, these equations can formally be solved with the following perturbative expansion,

(9) |

where only the fastest growing mode is taken into account (below we discuss how to generalize this to the full time dependence). The equations of motion, Eqs. (7-8) determine and in terms of the linear fluctuations,

(10) |

(11) |

where and are homogeneous functions of the wave vectors {} with degree zero. They are constructed from the fundamental mode coupling functions according to the recursion relations [13]

(12) | |||||

(13) | |||||

where , , , and , with , the initial perturbations.

### ii.2 In Matrix Form

The equations of motion can be rewritten in a more symmetric form [15, 6] by defining a two-component “vector”

(14) |

where the index ; and we have introduced a new time variable,

(15) |

corresponding to the number of e-folds of the scale factor. In the cosmology we are considering (), represents also the linear growth factor of perturbations in the growing mode. Using one can rewrite Eqs. (7) and (8) together as (we henceforth use the convention that repeated Fourier arguments are integrated over),

(16) |

where

(17) |

and is the symmetrized vertex matrix given by

(18) |

and zero otherwise, . An implicit integral solution to Eq. (16) can be found by Laplace transforming in the variable ,

(19) |

where denotes the initial conditions, set when the growth factor is one, and . Multiplying by the matrix,

and performing the inversion of the Laplace transform, we finally get the formal solution

(20) |

where the linear propagator is defined as ( to pick out the standard retarded propagator [15])

(21) |

for , whereas for due to causality, and as . The linear propagator is the Green’s function of the linearized version of Eq. (16) and describes the standard linear evolution of density and velocity fields from any configuration of initial perturbations. From a physical point of view, the most interesting initial conditions are given by when and are proportional random fields (instead of independent), in which case we can write

(22) |

where is a two component “vector”. This covers the usual case of growing mode initial conditions for which , for which the second term in Eq. (21) does not contribute upon contraction of with , and decaying mode initial conditions for which . For the sake of simplicity and definiteness, we will only contemplate the case throughout this paper.

### ii.3 Diagrammatic Representation of the Solution

Equation (20) for a given k-mode has a simple interpretation. The first term corresponds to the linear propagation from the initial conditions, whereas the second term contains information on non-linear interactions (mode-mode coupling). This non-linear contribution comes from the interaction of all pairs of waves and (whose sum is k due to translational invariance), at all intermediate times (with ). The interaction is characterized by the vertex and then linearly evolved in time by . As we will see shortly, this simple interpretation will lead, when applied recursively, to a graphical representation of the solution in terms of diagrams that allows the resummation of the different statistical objects of interest.

With the help of Eq. (20), we seek for an explicit expression for in the form of a series expansion [see Eq. (9)]

(24) |

where . Replacing Eqs. (23,24) into Eq. (20), we find the recursion relation satisfied by the kernels [see Eqs. (12-13)]

(25) |

where the r.h.s. has to be symmetrized under interchange of any two wave vectors. For , . Equation (25) reduces to the standard recursion relations given by Eqs. (12-13) in the limit that initial conditions are imposed in the infinite past, i.e. by replacing the lower limit of integration in Eq. (25) by , in which case only the fastest growing mode survives and . Otherwise Eq. (25) gives the full time dependence of PT solutions, including all transients from initial conditions [15].

Although in principle the solution for the kernels gives us an analytic expression for , it quickly becomes very cumbersome with increasing , therefore in practice this is not an ideal way to proceed; this is however the route taken in standard PT. Another reason apart from technical complication is the fact that the kernels (thought as functions of ) contain redundant information, since they all share the same building blocks, the vertex and linear propagator. In order to overcome these disadvantages we will introduce a graphical representation of the basic objects of the theory, and simple rules to combine them, such as to build a set of diagrams that are in one to one correspondence with each order .

As discussed [see Eq. (20)], there are only three basic building blocks, the initial field , the linear propagator , and the vertex . Their graphical representation are shown in Fig. 2. Since we resort to a graphical representation of the basic building blocks rather than derived quantities such as the kernels (as in done in the diagrammatic approach in standard PT [12, 13, 10]), the vertices in our diagrams are simpler than in standard PT, they always involve three lines (a consequence of quadratic nonlinearities) instead of a variable number of lines depending on the order in PT. The complication, at first sight, is that the diagrams carry information about time variables when the nonlinear interactions occur, and one integrates over these intermediate times (effectively summing over all possible interactions that happen between the initial conditions and the present time). This information is hidden in standard PT, where time evolution has been already “integrated out” (standard PT diagrams can be indeed thought as “collapsing” the diagrams shown in Fig. 3 below in the time direction). We will see here that this information can be taken advantage of, in fact making possible the process of resummation.

To obtain the set of diagrams corresponding to the th order term in Eq. (23) one has to draw all topologically different trees with vertices (branchings) and initial conditions. Each tree is obtained as follows: from the final time variable we draw a line backwards until it reaches a vertex, where the line bifurcates into two branches, with each of them continuing until they reach another vertex or an initial field at . If the branching is asymmetric it carries a factor of . This process is repeated at each vertex until all the branches end up in initial fields.

Each diagram with vertices represents an integral contributing to . The rules to obtain these integrals are as follows. Each of the initial fields is characterized by a momentum . Each branching corresponds to a vertex and represents the interaction of two incoming waves and that couple to form one outgoing k (due to the quadratic nonlinearities in the equations of motion). Each interaction happens at a given time () and conserves momentum (). The branches correspond to linear propagators representing the linear evolution of a given mode from time variables until . Notice that all the arrows throughout the diagram point away from the initial conditions, indicating the direction of forward time evolution. Finally, all the intermediate wave vectors are integrated over as well as all the time variables , each between . The diagrams up to are shown in Fig. 3. Note that it is the existence of these simple rules for finding the th order term, which are independent of , that will allow us to considerably simplify the resummation of diagrams.

As an example of these rules let us write explicitly the integrals corresponding to the and diagrams in Fig. 3,

(26) | |||

(27) |

The diagrams we construct here are very similar to any field theory with quadratic nonlinearities, perhaps the closest example to ours is that of turbulence where similar methods are well known [16, 18, 17]. There is also a recent paper [19] that applies path-integral methods to gravitational clustering in cosmology.

## Iii Statistics

### iii.1 Initial Conditions

As it stands, the integral Equation (20) can be thought as an equation for in the presence of an “external source” or forcing given by the initial conditions (i.e. ), with prescribed statistics. Here we assume that the initial conditions are Gaussian; the statistical properties of are then completely characterized by its two-point correlator

(28) |

where denotes the initial power spectrum of density fluctuations. According to Wick’s theorem, all higher order correlations of an odd number of fields vanish, whereas for an even number there are contributions corresponding to all different pairings of the fields,

(29) |

### iii.2 Power Spectrum

The goal is to calculate different correlation functions of the final field , the simplest being its two-point correlator or power spectrum,

(30) |

Equations (28-29) above give ensemble averages of a product of initial fields , with the result either vanishing or depending only on the initial power spectrum. This will allow us to develop a perturbative expansion for the final power spectrum in terms of integrals of the initial one, similar to the expansion for in terms of described in Eqs. (23) and (24). Replacing the series expansion of Eq. (23) into Eq. (30), we obtain where the -loop term is given by (assuming Gaussian initial conditions)

(31) |

We can now extend the diagrammatic representation of Section II.3 to describe each term in the above expansion. To draw each diagram that contributes to , we put one of the tree-like diagrams for against one for with their initial fields facing each other. Next, we pair the initial fields in all possible ways and glue the pairs. Each pair of initial fields with characteristic wave numbers and and indices and is then converted to an initial power spectrum . Finally, we do the same with all combinations of one tree of with one of (i.e. taking into account that after order there is more than one tree at each order). Going through this process one often obtains the same diagram. Therefore, the final weight of a diagram is given by the numerical factor for each “independent” diagram due to its trees, as described in Fig. (2), and another factor that takes into consideration the counting of the equal diagrams generated by the pairing process. All diagrams up to one-loop and some of the two-loop diagrams are shown in Fig. 5 (see [16] for a full account in the similar case of turbulence).

Two sets of diagrams have not been included in Fig. 5 since they vanish identically. One set is composed by diagrams that contain a sub-diagram linked to the main part by only one linear propagator. It can be shown that the sum of these fragments to all orders gives the ensemble average of the field, , which is zero. For a similar reason the set of unconnected diagrams vanishes.

For future use in the process of resummation, we define here the “principal cross section” of a diagram as the line that splits the diagram into two cutting only through initial power spectra. It represents the points where the two trees have been glued together. It is shown as a dashed line for the different diagrams in Fig. 5. Each diagram in Fig. 5 is in one-to-one correspondence with a certain integral. The rules for obtaining these integrals are the same to those described in Section II.3. For example, the diagram in Fig. 5, leads to the well known power spectrum for growing mode initial conditions ()

(32) |

and the integral corresponding to diagram (one loop) in the same figure gives,

(33) |

where we introduced the shorthand notations ( for growing mode initial conditions) and

where we made use of the symmetry with respect to the dashed line in Fig (5) diagram 8, and defined,

### iii.3 Non Linear Propagator

#### iii.3.1 Definition and Diagrammatic Representation

Nonlinearities can be thought as modifying the linear propagator defined in Eq. (21), leading to propagator renormalization. The non-linear propagator is defined by

(34) |

Note that translation invariance requires the condition , familiar from two-point statistics. Indeed, translation by some vector q, , gives . Strictly speaking, this object is not a propagator in the usual sense of a Green’s function, since nonlinearities are involved, i.e. it is not true that in the nonlinear case . But is the generalization of the linear propagator in the sense that its diagrammatic representation also has one entry and one exit arrow. It also generates all diagrams for the two-point function that can be divided into two subdiagrams by cutting a single initial power spectrum , similar to what happens with in the diagram for the linear power spectrum (diag. 1 in Fig. 5)

(35) |

where we have explicitly separated the linear part from the non-linear contributions. The fully non-linear propagator represents the ensemble averaged response of the final density and velocity divergence fields to variations in the initial conditions. In more physical terms it quantifies how much information of the initial distribution of a k-mode remains in the final state at the same k, in an ensemble average sense. In paper II [20], we show that for Gaussian initial conditions , therefore can also be thought as a measure of the cross-correlation between final and initial conditions, or indeed as a genuine propagator (or Green’s function) in two-point sense. At very large scales, we expect linear evolution to scale initial conditions by the growth factor; consequently, the non-linear corrections in Eq. (35) vanish as . As one approaches the non-linear regime interactions gradually “erase” the initial distributions; therefore, we expect to vanish as . See section III.3.2 below for the behavior of in the Zel’dovich approximation, the case of the exact dynamics is presented in paper II [20] together with comparisons against numerical simulations.

The series defining the nonlinear propagator can be represented in a diagrammatic fashion. For Gaussian initial conditions, only the odd terms in the expansion for in Eq. (35) will contribute to . The th non-linear correction corresponding to loops is obtained by a functional derivative of . This can be done diagrammatically in a straightforward way. From each tree characterizing in Fig. 3 there will be contributions, each obtained after dropping one initial condition represented by a blank circle (this accounts for the functional derivative). The linear propagator left at that branch will carry a wave number k. Finally the ensemble average of the remaining ’s is done by pairing the fields according to Gaussianity.

Figure 6 shows all diagrams up to two loops for the non-linear propagator. Notice that each diagram has one entering arrow and one exiting, connected by a chain of linear propagators that runs through the diagram without intersecting initial conditions. This path is called the “principal path”, and will be used later in the process of resummation.

#### iii.3.2 A simple example: Zel’dovich Approximation

To illustrate the behavior of the nonlinear propagator we go back to the discussion in the introduction about the Zel’dovich Approximation. Let us work out the density propagator, which assuming growing mode initial conditions follows directly from taking the derivative of the final density with respect to the initial one. For simplicity we consider only the fastest growing mode at each order, so we can use the standard PT kernel representation, Eqs. (9), (10) and (12), and take advantage of the simple form the density field kernel takes in this approximation [21]

(36) |

where . Due to Gaussian initial conditions, only odd terms in the series expansion for contribute, therefore

(37) |

where the term gives just the scale (and growth) factor. Using Eq. (36) this series can be summed up right away,

(38) |

where . The extension of this result to the exact dynamics for density and velocity fields and comparison with measurements in numerical simulations is given in [20].

### iii.4 Vertex

The vertex defined in Eq. (18) is the third basic object of the theory, and is also renormalized due to higher-order nonlinear interactions. We define the symmetric full vertex through the relation

(39) |

thus