A&A 465, 725-747 (2007)
DOI: 10.1051/0004-6361:20066832

Large-N expansions applied to gravitational clustering

P. Valageas

Service de Physique Théorique, CEA Saclay, 91191 Gif-sur-Yvette, France

Received 28 November 2006 / Accepted 26 January 2007

Abstract
We develop a path-integral formalism to study the formation of large-scale structures in the universe. Starting from the equations of motion of hydrodynamics (single-stream approximation) we derive the action which describes the statistical properties of the density and velocity fields for Gaussian initial conditions. Then, we present large-N expansions (associated with a generalization to N fields or with a semi-classical expansion) of the path-integral defined by this action. This provides a systematic expansion for two-point functions such as the response function and the usual two-point correlation. We present the results of two such expansions (and related variants) at one-loop order for a SCDM and a $\Lambda$CDM cosmology. We find that the response function exhibits fast oscillations in the non-linear regime with an amplitude which either follows the linear prediction (for the direct steepest-descent scheme) or decays (for the 2PI effective action scheme). On the other hand, the correlation function agrees with the standard one-loop result in the quasi-linear regime and remains well-behaved in the highly non-linear regime. This suggests that these large-N expansions could provide a good framework to study the dynamics of gravitational clustering in the non-linear regime. Moreover, the use of various expansion schemes allows one to estimate their range of validity without the need of N-body simulations and could provide a better accuracy in the weakly non-linear regime.

Key words: gravitation - cosmology: theory - cosmology: large-scale structure of Universe

   
1 Introduction

The large-scale structures we observe in the present universe (such as galaxies and clusters of galaxies) have formed thanks to gravitational instability (Pebbles 1980) which amplified the small density perturbations created in the early universe (e.g. through quantum fluctuations generated during an inflationary phase, Liddle & Lyth 1993). Moreover, the power increases at small scales as in the CDM model (Peebles 1982) which leads to a hierarchical scenario where small scales become non-linear first. Then, small structures merge in the course of time to build increasingly large objects. Therefore, at large scales or at early times one can use a perturbative approach to describe the density field (e.g. Bernardeau et al. 2002, for a review). This is of great practical interest for observations such as weak-lensing surveys (e.g. Bernardeau et al. 1997) and CMB studies (e.g. Seljak & Zaldarriaga 1996). At smaller scales one uses N-body numerical simulations to obtain fits to the density power-spectrum (e.g. Smith et al. 2003) which can be compared with observations.

In the weakly non-linear regime one may apply the standard perturbative expansion over powers of the initial density fluctuations (e.g., Fry 1984; Goroff et al. 1986). This procedure is used within the hydrodynamical framework where the system is described by density and velocity fields rather than by the phase-space distribution $f({\vec x},{\vec v},t)$, that is one starts from the Euler equation of motion rather than the Vlasov equation. This is quite successful at tree-order where one can obtain the leading-order contribution to all p-point correlations and reconstruct the density probability distribution (Bernardeau 1994), assuming Gaussian initial fluctuations as in usual inflationary scenarios (Liddle & Lyth 1993). However, several problems show up when one tries to compute higher-order corrections. First,ultraviolet divergences appear for linear power-spectra $P_{\rm L}(k) \propto k^n$ with $n\geq -1$ (Scoccimarro & Frieman 1996b) which would break the self-similarity seen in numerical simulations. This may also be interpreted as a consequence of the breakdown of the hydrodynamical description beyond shell-crossing (Valageas 2002). However, for n<-1 the one-loop correction yields a good agreement with numerical simulation (Scoccimarro & Frieman 1996b). Secondly, being an expansion over powers of $P_{\rm L}(k)$ this standard perturbative approach yields a series of terms which grow increasingly fast into the non-linear regime at higher orders and cannot relax at small scales. This makes the series badly behaved and it seems difficult to go deeper into the non-linear regime in this manner.

On the other hand, a good description of the weakly non-linear regime becomes of great practical interest as cosmological probes now aim to constrain cosmological parameters with an accuracy of the order of $1\%$ to keep pace with CMB observations. In particular, weak-lensing surveys mainly probe this transition regime to non-linearity and a good accuracy for the matter power-spectrum is required to make the most of observations and constrain the dark energy equation of state (Munshi et al. 2007). Another probe of the expansion history of the Universe is provided by the measure of the baryon acoustic oscillations which also requires a good description of weakly non-linear effects (Seo & Eisenstein 2007; Koehler et al. 2007).

As a consequence, it is worth investigating other approaches which may partly cure some of these problems and provide a better description of the non-linear regime. For instance, McDonald (2006) proposed a renormalization group method to improve the standard perturbative approach. On the other hand, Crocce & Scoccimarro (2006a,b) presented a diagrammatic technique to organize the various terms which arise in the standard perturbative expansion and to perform some infinite partial resummations. Besides, using the asymptotic behavior of the vertices of the theory they managed to complete the resummation in the small-scale limit for the response function R (also called the propagator).

In this article, we apply to the hydrodynamical framework the method presented in Valageas (2004) for the Vlasov equation. This approach introduces a path-integral formalism to derive the statistical properties of the system from its action S. Then, one applies a "large-N expansion'' (which can be seen as similar to the semi-classical expansion over powers of $\hbar$ or can be derived from a generalization to N fields) to compute the quantities of interest such as the two-point correlation. This offers a different perspective to address the non-linear regime of gravitational clustering which complements other approaches and may also serve as a basis for other approximation schemes than the ones described in this article. For instance, previous approaches can be recovered by expanding the path-integral obtained in this paper in different manners (e.g. the usual perturbative expansion corresponds to the expansion over the non-Gaussian part of the action S). The two large-N expansions discussed in this article are two other means of performing systematic expansions to compute the required correlation functions. This also amounts to reorganize the standard perturbative expansion by performing some partial infinite resummations. In this fashion, one can hope to obtain new expansion schemes which are better suited to describe the non-linear regime of gravitational clustering.

First, we recall in Sect. 2 the equations of motion obtained within the hydrodynamical framework. Next, we derive the path-integral formulation of the problem in Sect. 3 and we describe two possible large-N expansions in Sect. 4. Finally, we present our numerical results in Sects. 5-7, focusing on a critical-density universe. Since our formalism applies equally well to any cosmology we describe our results for a $\Lambda$CDM universe in Sect. 8. Finally, we conclude in Sect. 9.

   
2 Equations of motion

   
2.1 Fluid approach

At scales much larger than the Jeans length both the cold dark matter and the baryons can be described as a pressureless dust. Moreover, before orbit crossing one can use a hydrodynamical approach (in the continuum limit where the mass m of particles goes to zero). Then, at scales much smaller than the horizon where the Newtonian approximation is valid, the equations of motion read (Peebles 1980):

$$\frac{\partial\delta}{\partial\tau} + \nabla.[(1+\delta) {\vec v}] = 0,$$ (1)

$$\frac{\partial{\vec v}}{\partial\tau} + {\cal H}{\vec v}+ ({\vec v}.\nabla) {\vec v}= - \nabla \phi,$$ (2)

$$\Delta \phi = \frac{3}{2} \Omega_{\rm m}{\cal H}^2 \delta ,$$ (3)

where $\tau=\int {\rm d}t/a$ is the conformal time (and a the scale factor), ${\cal H}={\rm d}\ln a/{\rm d}\tau$ the conformal expansion rate and $\Omega_{\rm m}$ the matter density cosmological parameter. Here $\delta$ is the matter density contrast and ${\vec v}$ the peculiar velocity. Since the vorticity field decays within linear theory (Peebles 1980) we take the velocity to be a potential field so that ${\vec v}$ is fully specified by its divergence  $\theta=\nabla.{\vec v}$. Then, in Fourier space Eqs. (1)-(3) read (Goroff et al. 1986):

 

$$\frac{\partial\delta({\vec k},\tau)}{\partial\tau} + \theta({\vec... ...\alpha({\vec k}_1,{\vec k}_2) \theta({\vec k}_1,\tau) \delta({\vec k}_2,\tau) ,$$ (4)

 

$$\frac{\partial\theta({\vec k},\tau)}{\partial\tau} + {\cal H}\the... ... \beta({\vec k}_1,{\vec k}_2) \theta({\vec k}_1,\tau) \theta({\vec k}_2,\tau) ,$$ (5)

where $\delta_{\rm D}$ is the Dirac distribution, the coupling functions $\alpha$ and $\beta$ are given by:

$$\alpha({\vec k}_1,{\vec k}_2)= \frac{({\vec k}_1+{\vec k}_2).... ..._1+{\vec k}_2\vert^2 ({\vec k}_1.{\vec k}_2)}{2 k_1^2 k_2^2} ,$$ (6)

and we defined the Fourier transforms as:

$$\delta({\vec k}) = \int\frac{{\rm d}{\vec x}}{(2\pi)^3} {\rm e}^{-{\rm i}{\vec k}.{\vec x}} \delta({\vec x}) .$$ (7)

As in Crocce & Scoccimarro (2006a,b) let us define the two-component vector $\psi$ as:

$$\psi({\vec k},\eta) = \left(\begin{array}{c} \psi_1({\vec k},... ...eta) \\ -\theta({\vec k},\eta)/{\cal H}f \end{array}\right) ,$$ (8)

where we introduced the time coordinate $\eta$ defined from the linear growing rate D₊ of the density contrast:

$$\eta= \ln \frac{D_+(\tau)}{D_{+0}} , \;\;\; f=\frac{{\rm d}\!... ... a} = \frac{1}{{\cal H}}\frac{{\rm d}\!\ln D_+}{{\rm d}\tau} ,$$ (9)

with D₊₀ is the value of D₊ today at redshift z=0 and:

$$\frac{{\rm d}^2 D_+}{{\rm d}\tau^2}+{\cal H}\frac{{\rm d}D_+}{{\rm d}\tau} = \frac{3}{2}\Omega_{\rm m}{\cal H}^2 D_+ .$$ (10)

Then, the equations of motion (4)-(5) can be written as:

$${\cal O}(x,x').\psi(x') = K_{\rm s}(x;x_1,x_2) . \psi(x_1) \psi(x_2) ,$$ (11)

where we introduced the coordinate $x=({\vec k},\eta,i)$ where i=1,2 is the discrete index of the two-component vectors. In Eq. (11) and in the following we use the convention that repeated coordinates are integrated over as:

$${\cal O}(x,x').\psi(x') = \!\! \int\!{\rm d}{\vec k}'{\rm d}\... ...'}({\vec k},\eta;{\vec k}',\eta') \psi_{i'}({\vec k}',\eta') .$$ (12)

The matrix ${\cal O}$ reads:

$$\! {\cal O}(x,x') \! = \! \left( \begin{array}{cc} \frac{\par... ... \delta_{\rm D}({\vec k}-{\vec k}') \delta_{\rm D}(\eta-\eta')$$ (13)

whereas the symmetric vertex $K_{\rm s}(x;x_1,x_2)=K_{\rm s}(x;x_2,x_1)$ writes:

 

$$K_{\rm s}(x;x_1,x_2) = \delta_{\rm D}({\vec k}_1+{\vec k}_2-{\vec... ...) \delta_{\rm D}(\eta_2-\eta) \gamma^{\rm s}_{i;i_1,i_2}({\vec k}_1,{\vec k}_2)$$ (14)

with:

$$\gamma^{\rm s}_{1;1,2}({\vec k}_1,{\vec k}_2)= \frac{\alpha({... ...ec k}_1,{\vec k}_2)= \frac{\alpha({\vec k}_1,{\vec k}_2)}{2} ,$$ (15)

$$\gamma^{\rm s}_{2;2,2}({\vec k}_1,{\vec k}_2)= \beta({\vec k}_1,{\vec k}_2) ,$$ (16)

and zero otherwise (Crocce & Scoccimarro 2006a).

   
2.2 Linear regime

At large scales or at early times where the density and velocity fluctuations are small one can linearize the equation of motion (11) which yields ${\cal O}.\psi_{\rm L}=0$. Then, the linear growing mode is merely:

$$\psi_{\rm L}(x) = {\rm e}^{\eta} \delta_{\rm L0}({\vec k}) \left(\begin{array}{c} 1 \\ 1 \end{array}\right) ,$$ (17)

where $\delta_{\rm L0}({\vec k})$ is the linear density contrast today at redshift z=0. We shall define the initial conditions by the linear growing mode (17) (which is consistent with the fact that we only kept the curl-free part of the peculiar velocity in Eqs. (4), (5)). Moreover, we assume Gaussian homogeneous and isotropic initial conditions defined by the linear power-spectrum $P_{\rm L0}(k)$:

$$\langle\delta_{\rm L0}({\vec k}_1) \delta_{\rm L0}({\vec k}_2... ...ngle= P_{\rm L0}(k_1) \delta_{\rm D}({\vec k}_1+{\vec k}_2) .$$ (18)

It is convenient to define the power per logarithmic wavenumber $\Delta ^2(k)$by:

$$\Delta^2(k) = 4\pi k^3 P(k) ,$$ (19)

whence:

$$\langle\delta({\vec x}_1)\delta({\vec x}_2)\rangle= \int_0^{\... ... x}_1-{\vec x}_2\vert}{k\vert{\vec x}_1-{\vec x}_2\vert} \cdot$$ (20)

Thus, for a CDM cosmology the linear power $\Delta^2_{\rm L0}(k)$grows as k⁴ at small k and as $\ln k$ at high k. Then, the linear two-point correlation function $G_{\rm L}(x_1,x_2)$ reads:

 

$$G_{\rm L}(x_1,x_2) \langle\psi_{\rm L}(x_1) \psi_{\rm L}(x_2) \rangle$$ =

$${\rm e}^{\eta_1+\eta_2} P_{\rm L0}(k_1) \delta_{\rm D}({\vec k}_1+{\vec k}_2) \left(\begin{array}{cc} 1 & 1 \\ 1 & 1 \end{array}\right) .$$ = (21)

As in Valageas (2004) it is convenient to introduce the response function R(x₁,x₂) (which is also called the propagator in Crocce & Scoccimarro 2006a,b) defined by the functional derivative:

$$R(x_1,x_2) = \langle\frac{\delta \psi(x_1)}{\delta\zeta(x_2)} \rangle_{\zeta=0} ,$$ (22)

where $\zeta(x)$ is a "noise'' added to the r.h.s. in Eq. (11). Thus, R(x₁,x₂) measures the linear response of the system to an external source of noise. Because of causality it contains an Heaviside factor $\theta(\eta_1-\eta_2)$ since the field $\psi(x_1)$ can only depend on the values of the "noise'' at earlier times $\eta_2\leq\eta_1$. Moreover, it satisfies the initial condition:

$$\eta_1\rightarrow\eta_2^+\!\!: \; R(x_1,x_2) \rightarrow \delta_{\rm D}({\vec k}_1-{\vec k}_2) \delta_{i_1,i_2} .$$ (23)

In the linear regime, the response $R_{\rm L}$ may be obtained by first calculating the linear solution $\psi_{\rm L}(x\vert\zeta)$ of the linearized equations of motion ${\cal O}.\psi_{\rm L}(x\vert\zeta) = \zeta$ and next taking the functional derivative (22). Alternatively, one can directly obtain the matrix $R_{\rm L}$ from the initial condition (23) and the linear dynamics ${\cal O}.R_{\rm L}=0$ (as implied by the definition (22) and ${\cal O}.\psi_{\rm L}=0$). For the Einstein-de-Sitter cosmology ( $\Omega_{\rm m}=1,\Omega_{\Lambda}=0$) where the factors $\Omega_{\rm m}/f^2$ are constant and equal to unity, this can be easily solved by looking for eigenmodes or performing a Laplace transform as in Crocce & Scoccimarro (2006a). This yields:

 

$${R_{\rm L}(x_1,x_2) = \theta(\eta_1-\eta_2) \delta_{\rm D}({\vec ... ...}{5} \left(\begin{array}{cc} 2 & -2 \\ -3 & 3 \end{array}\right) \right\}\cdot$$ (24)

For other cosmological parameters one must solve numerically the linearized equation of motion for $R_{\rm L}$. Note that in all cases it exhibits the factor $\theta(\eta_1-\eta_2) \delta_{\rm D}({\vec k}_1-{\vec k}_2)$ times a term which is independent of wavenumbers ${\vec k}_1,{\vec k}_2$.

   
3 Path-integral formalism

As in Valageas (2004) we can apply a path-integral approach to the hydrodynamical system (11) since we are only interested in the statistical properties of the density and velocity fields (we do not look for peculiar solutions of the equations of motion). Let us briefly recall how this can be done (also Martin et al. 1973; Phythian 1977). In order to include explicitly the initial conditions we rewrite Eq. (11) as:

$${\cal O}. \psi = K_{\rm s} . \psi \psi + \mu_i$$ (25)

with $\psi=0$ for $\eta<\eta_i$ and:

$$\mu_i(x) = \delta_{\rm D}(\eta-\eta_i) {\rm e}^{\eta_i} \delt... ...}({\vec k}) \left(\begin{array}{c} 1 \\ 1 \end{array}\right) .$$ (26)

Thus, the source $\mu_i$ merely provides the initial conditions at time $\eta_i$, obtained from the linear growing mode (17). We shall eventually take the limit $\eta_i\rightarrow -\infty$. Next, we define the generating functional Z[j] by:

$$Z[j] = \langle{\rm e}^{j.\psi} \rangle= \int [{\rm d}\mu_i] \... ...m e}^{j.\psi[\mu_i]-\frac{1}{2} \mu_i . \Delta_i^{-1} . \mu_i}$$ (27)

where we took the average over the Gaussian initial conditions:

$$\langle\mu_i\rangle=0, \;\;\; \langle\mu_i(x_1)\mu_i(x_2)\rangle= \Delta_i(x_1,x_2) .$$ (28)

All statistical properties of the field $\psi$ may be obtained from Z[j]. It is convenient to write Eq. (27) as:

 

$$Z[j] =\int [{\rm d}\mu_i] [{\rm d}\psi] \; \vert\det M\vert \; \d... ..._{\rm s}.\psi\psi) {\rm e}^{j.\psi-\frac{1}{2} \mu_i . \Delta_i^{-1} . \mu_i} ,$$ (29)

where the Jacobian $\vert\det M\vert$ is defined by the functional derivative $M=\delta\mu_i/\delta\psi$. We can compute the latter by factorizing the operator $\partial/\partial\eta$:

$$\! \det M \! = \det~({\cal O}-2 K_{\rm s}.\psi) = \det~(\part... ...tial\eta) \det~(1+{\tilde{\cal O}}-2 {\tilde{K}_{\rm s}}.\psi)$$ (30)

with:

$$\! {\tilde{\cal O}}(x,x') \! = \! \delta_{\rm D}({\vec k}-{\v... ...) & \frac{3\Omega_{\rm m}}{2f^2}(\eta') -1 \end{array}\right)$$ (31)

(note that we separated the identity part out of ${\tilde{\cal O}}$) and:

 

$${\tilde{K}_{\rm s}}(x;x_1,x_2) = \delta_{\rm D}({\vec k}_1+{\vec ... ...elta_{\rm D}(\eta_2-\eta_1) \gamma^{\rm s}_{i;i_1,i_2}({\vec k}_1,{\vec k}_2) .$$ (32)

Then, we have:

$$\det M = \det~(\partial/\partial\eta) \; {\rm e}^{{\rm Tr}\ln~(1+{\tilde{\cal O}}-2 {\tilde{K}_{\rm s}}.\psi)}$$ (33)

and:

 

$${\rm Tr}\ln~(1+{\tilde{\cal O}}-2 {\tilde{K}_{\rm s}}.\psi) \sum_{q=1}^{\infty} \frac{(-1)^{q-1}}{q} {\rm Tr}({\tilde{\cal O}}-2 {\tilde{K}_{\rm s}}.\psi)^q$$ =

$${\rm Tr}({\tilde{\cal O}}-2 {\tilde{K}_{\rm s}}.\psi) ,$$ = (34)

where we used the fact that the trace of powers $q\geq 2$ vanishes because of the Heaviside factors $\theta(\eta-\eta')$ (Zinn-Justin 1989; Valageas 2004). Next, the trace ${\rm Tr}{\tilde{\cal O}}$ is an irrelevant constant whereas ${\rm Tr}(2 {\tilde{K}_{\rm s}}.\psi)=0$ (we regularize as usual the gravitational interaction at large distances to remove ambiguities, e.g. perform angular integrals first, see also Valageas 2004). Therefore, the Jacobian in the functional integral (29) is a constant which can be absorbed into the normalization. Then, introducing an imaginary ghost field $\lambda$ to express the Dirac as an exponential and performing the Gaussian integration over $\mu_i$ we obtain:

 

$$\int [{\rm d}\mu_i] [{\rm d}\psi] [{\rm d}\lambda] \; {\rm e}^{j.... ...\mu_i-{\cal O}.\psi+K_{\rm s}.\psi\psi) - \frac{1}{2}\mu_i.\Delta_i^{-1}.\mu_i}$$                    Z[j] =

$$\int [{\rm d}\psi] [{\rm d}\lambda] \; {\rm e}^{j.\psi + \lambda.(-{\cal O}.\psi+K_{\rm s}.\psi\psi) + \frac{1}{2}\lambda.\Delta_i.\lambda} .$$ = (35)

Thus, the statistical properties of the system (11) are described by the action $S[\psi,\lambda]$ defined by:

$$S[\psi,\lambda] = \lambda.({\cal O}.\psi-K_{\rm s}.\psi\psi) - \frac{1}{2}\lambda.\Delta_i.\lambda.$$ (36)

Moreover, we can note that adding a "noise'' $\zeta$ to the r.h.s. of Eq. (11) amounts to change $\mu_i\rightarrow\mu_i+\zeta$ which translates into $S\rightarrow S-\lambda.\zeta$. Therefore, functional derivatives with respect to $\zeta$ are equivalent to insertions of the ghost field $\lambda$. In particular, we have:

$$R(x_1,x_2) = \langle\psi(x_1) \lambda(x_2) \rangle, \;\;\; \langle\lambda\rangle=0 , \;\;\; \langle\lambda\lambda\rangle=0 .$$ (37)

   
4 Large-N expansions

The path-integral (35) can be computed by expanding over powers of its non-Gaussian part (i.e. over powers of $K_{\rm s}$). This actually yields the usual perturbative expansion over powers of the linear power-spectrum $P_{\rm L}$ (see also Valageas 2001, 2004, for the case of the Vlasov equation of motion). On the other hand, the path-integral (35) may also be studied through a large-N expansion as in Valageas (2004). Thus, one considers the generating functional Z_N[j,h] defined by:

$$Z_N[j,h] = \int [{\rm d}\psi] [{\rm d}\lambda] \; {\rm e}^{N[j.\psi+h.\lambda-S[\psi,\lambda]]}$$ (38)

and one looks for an expansion over powers of 1/N, taking eventually N=1into the results. As discussed in Valageas (2004) the large-N expansions may also be derived from a generalization of the dynamics to N fields $\psi^{(\alpha)}$. This yields the same results once we take care of the long-wavelength divergences which constrain which subsets of diagrams need to be gathered.

   
4.1 Steepest-descent method

A first approach to handle the large-N limit of Eq. (38) is to use a steepest-descent method (also called a semi-classical or loopwise expansion in the case of usual Quantum field theory with $\hbar=1/N$). This yields for auxiliary correlation and response functions G₀ and R₀ the equations (Valageas 2004):

   

$${\cal O}(x,z).G_0(z,y) = 0$$ (39)

$${\cal O}(x,z).R_0(z,y) = \delta_{\rm D}(x-y)$$ (40)

$$R_0(x,z).{\cal O}(z,y) = \delta_{\rm D}(x-y)$$ (41)

whereas the actualcorrelation and response functions obey:

   

$${\cal O}(x,z).G(z,y) = \Sigma(x,z).G(z,y) + \Pi(x,z).R^T(z,y)$$ (42)

$${\cal O}(x,z).R(z,y) = \delta_{\rm D}(x-y) + \Sigma(x,z).R(z,y)$$ (43)

$$R(x,z).{\cal O}(z,y) = \delta_{\rm D}(x-y) + R(x,z).\Sigma(z,y)$$ (44)

where the self-energy terms ${\Sigma }$ and $\Pi$ are given at one-loop order by:

  

$$\Sigma(x,y) = 4 K_{\rm s}(x;x_1,x_2) K_{\rm s}(z;y,z_2) R_0(x_1,z) G_0(x_2,z_2)$$ (45)

$$\Pi(x,y) = 2 K_{\rm s}(x;x_1,x_2) K_{\rm s}(y;y_1,y_2) G_0(x_1,y_1) G_0(x_2,y_2).$$ (46)

Note that Eqs. (39)-(44) are exact and that the expansion over powers of 1/N only enters the expression of the self-energy (45)-(46). Here we only kept the lowest-order terms (see Valageas 2004) for the next-order terms). We also took the limit $\eta_i\rightarrow -\infty$ so that terms involving $\Delta_i$ vanish. From Sect. 2.2 and Eqs. (39)-(41) we find that the auxiliary matrices G₀ and R₀ are actually equal to their linear counterparts:

$$G_0= G_{\rm L}, \;\;\; R_0 = R_{\rm L} .$$ (47)

   
4.2 2PI effective action method

As described in Valageas (2004) a second approach is to first introduce the double Legendre transform $\Gamma[\psi,G]$ of the functional $W=\ln Z$ (with respect to both the field $\psi$ and its two-point correlation G) and next to apply the 1/N expansion to $\Gamma$. This "2PI effective action''method yields the same Eqs. (42)-(44) and the self-energy shows the same structure as (45)-(46) where G₀ and R₀ are replaced by G and R. Thus, the direct steepest-descent method yields a series of linear equations which can be solved directly whereas the 2PI effective action method gives a system of non-linear equations (through the dependence on G and R of ${\Sigma }$ and $\Pi$) which must usually be solved numerically by an iterative scheme. However, thanks to the Heaviside factors appearing in the response R and the self-energy ${\Sigma }$ these equations can be solved directly by integrating forward over time $\eta _1$.

   
4.3 General properties

As discussed in Valageas (2004) both the steepest-descent and 2PI effective action methods agree with the standard perturbative analysis over powers of $P_{\rm L0}(k)$ up to the order used for the self-energy (e.g. up to order $P_{\rm L0}^2$ if we only consider the one-loop terms (45)-(46)). As compared with the standard perturbative approach, the two schemes described above also include two different infinite partial resummations, as can be seen from Eqs. (42)-(44) which clearly generate terms at all orders over  $P_{\rm L0}$ for G and R even if ${\Sigma }$ and $\Pi$ are only linear or quadratic over  $P_{\rm L0}$.

We can note that the equations we obtain for the hydrodynamical system (1)-(3) are simpler than for the collisionless system studied in Valageas (2004) which is described by the Vlasov-Poisson system. Indeed, here the mean $\langle\psi\rangle$ vanishes at all orders. This can be explicitly checked from Eqs. (1)-(3) and in the derivation of Eqs. (39)-(46). This is not the case for the Vlasov dynamics where Eq. (47) does not hold.

Using a diagrammatic technique Crocce & Scoccimarro (2006a,b) derived Eqs. (42)-(46) in an integral form (i.e. without the differential operator ${\cal O}$ in the l.h.s.) by first integrating the equation of motion (11). Of course, our path-integral approach can also be applied to the integral form of Eq. (11). This amounts to write Eq. (11) as $\psi=\psi_{\rm L}+{\tilde K}_{\rm s} \psi\psi$ where $\psi_{\rm L}$ is the linear growing mode and ${\tilde K}_{\rm s}$ the integral form of the vertex $K_{\rm s}$. Then, the procedure presented in Sect. 3 can be applied to this integral form of the equation of motion, as described for instance in Valageas (2001). In this article we preferred to keep Eq. (11) in its differential form so that the r.h.s. does not contain an integral over time. Then the self-energy terms $\Sigma(\eta_1,\eta_2)$and $\Pi(\eta_1,\eta_2)$ only depend on the response and correlation at the same times $(\eta_1,\eta_2)$, see Eqs. (45), (46). By contrast, the integrated vertex ${\tilde K}_{\rm s}$ would lead to self-energies which depend on the values of the response and correlation at all past times which would entail less efficient numerical computations.

   
4.4 Integral expression for the two-point correlation

Thanks to statistical homogeneity and isotropy the matrices G₀,G and $\Pi$are symmetric and of the form:

$$G(x_1,x_2) = \delta_{\rm D}({\vec k}_1+{\vec k}_2) G_{i_1,i_2}(k_1;\eta_1,\eta_2) ,$$ (48)

with:

$$G_{i_1,i_2}(k;\eta_1,\eta_2) = G_{i_2,i_1}(k;\eta_2,\eta_1) ,$$ (49)

whereas the matrices R₀,R and ${\Sigma }$ are of the form:

$$R(x_1,x_2) = \delta_{\rm D}({\vec k}_1-{\vec k}_2) \theta(\eta_1-\eta_2) R_{i_1,i_2}(k_1;\eta_1,\eta_2) .$$ (50)

Besides, from Eqs. (42)-(44) we see that the correlation Gcan also be written in terms of the response R as:

$$G(x_1,x_2) = R \times G_0(\eta_i) \times R^T + R.\Pi.R^T$$ (51)

where the second product is defined as in (12) whereas the first product does not contain any integration over time:

 

$$R \times G_0(\eta_i) \times R^T = \delta_{\rm D}({\vec k}_1+{\vec... ..._1,\eta_i) G_{0;i_1',i_2'}(k_1;\eta_i,\eta_i) R_{i_2,i_2'}(k_1;\eta_2,\eta_i) ,$$ (52)

and we let $\eta_i\rightarrow -\infty$. The physical meaning of Eq. (51) is clear. The first term in the r.h.s. means that the fluctuations at the initial time $\eta_i$ are merely transported forward in time through the matrix R. This is the only non-zero term in the linear regime (with $R=R_{\rm L}$ and $G_0=G_{\rm L}$). The effect of the non-linear dynamics is to modify the transport matrix R and to add a second term to the r.h.s. of Eq. (51). The latter has the meaning of a source term which produces fluctuations with two-point correlation $\Pi(\eta_1',\eta_2')$ at the times $(\eta_1',\eta_2')$ which are next transported forward to later times $(\eta_1,\eta_2)$ by the matrices $R(\eta_1,\eta_1')$ and $R^T(\eta_2',\eta_2)$. Thus, the "dot product'' means that we need to integrate over all fluctuations produced at all earlier times $\eta_1'<\eta_1$ and $\eta_2'<\eta_2$. This interpretation considers $\Pi$ as an external source for an "effective'' linear dynamics. For the system which we study here the "source'' $\Pi$ is not external but corresponds to the non-linear interaction of the field $\psi$, as seen from the explicit expression (46). Moreover, the linear operator which would define the transport with time also depends on G and R through ${\Sigma }$, see Eqs. (45) and (44).

Alternatively, Crocce & Scoccimarro (2006a) pointed out that the expression (51) is somewhat similar to the result obtained within a phenomenological halo model (Seljak 2000). Indeed, the latter model splits the non-linear power-spectrum into two terms, one that dominates in the linear regime (2-halo term) and the other that dominates in the highly non-linear regime (1-halo term). They could be identified with the first and second terms of Eq. (51). However, it is not clear whether the analogy can be pursued much further. In particular, Eq. (51) is explicitly dynamical (i.e. it explicitly involves integrations over past events) whereas the halo model gives a static expression (it writes the two-point correlation as integrals over the current halo distribution, which is the basic time-dependent quantity).

Note that within the 2PI effective action method, even if R would be treated independently of G the expression (51) does not give a quadratic dependence of G on R since $\Pi$ depends quadratically on G. Besides, solving Eq. (51) perturbatively over G₀ at fixed R actually generates terms at all orders q over powers G₀^q. By contrast, in the steepest-descent approach where $\Pi$ only depends on G₀the expression (51) is fully explicit and quadratic over both $G_0(=G_{\rm L})$ and R if the latter are treated independently. This is a signature of the additional resummations involved in the 2PI effective action method as compared with the direct steepest-descent approach.

An advantage of Eq. (51) is that once $\Pi$ and R are known, we obtain an explicit expression for G. Hence we do not need to solve Eq. (42) (by solving for R we have already performed the "inversion'' of ${\cal O}-\Sigma$). A second advantage of Eq. (51) is that it provides an expression for G which is clearly symmetric. Moreover, since we start from a two-point correlation $G_0(=G_{\rm L})$ which is positive, as shown by Eq. (21), we see from Eqs. (46) and (51) that both $\Pi$ and G are positive. This also holds for the 2PI effective action approach which may be obtained by iterating the system (42)-(46) (substituting G and R into the self-energy).

   
4.5 Integration over wavenumber for the self-energy

Using the symmetry (49) for $\Pi$ and the form (50) for ${\Sigma }$we see that $\Pi(x_1,x_2)$ and $\Sigma(x_1,x_2)$ only need to be computed at times $\eta_1\geq\eta_2$. Besides, thanks to the various Dirac factors the expressions (45)-(46) only involve an integral over one wavenumber ${\vec k}'$. Thus, using the expression (14) of the kernel $K_{\rm s}$ we can write from Eqs. (45), (46) the self-energy as:

 

$$\Sigma_{i_1 i_2}(k;\eta_1,\eta_2) = 4 \gamma^{\rm s}_{i_1;j_1,j_2... ... l_1}(k';\eta_1,\eta_2) G_{j_2 l_2}(\vert{\vec k}-{\vec k}'\vert;\eta_1,\eta_2)$$ (53)

and:

 

$$\Pi_{i_1 i_2}(k;\eta_1,\eta_2) = 2 \gamma^{\rm s}_{i_1;j_1,j_2}({... ..._1}(k';\eta_1,\eta_2) G_{j_2 l_2}(\vert{\vec k}-{\vec k}'\vert;\eta_1,\eta_2) ,$$ (54)

where we integrate over ${\vec k}'$ and we sum over j₁,j₂,l₁,l₂. In order to keep the symmetry ${\vec k}'\leftrightarrow{\vec k}-{\vec k}'$ in the numerical integration we use symmetric elliptic coordinates:

$${\vec k}' = \frac{{\vec k}}{2} + {\vec q}\;\;\; \mbox{with} \... ...h\zeta\sin\mu\sin\phi \\ \cosh\zeta\cos\mu \end{array}\right.$$ (55)

where ${\vec k}$ is taken along the third axis and $\zeta\geq0, 0\leq\mu\leq\pi, 0\leq\phi<2\pi$.

At equal times $\eta_1=\eta_2$ we can obtain from Eq. (53):

$$\Sigma(k;\eta,\eta) = - \frac{4\pi}{3} k^2 \int_0^{\infty} {\... ...\; \left(\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right) .$$ (56)

The explicit expression (56) shows that within the direct steepest-descent method we need a linear power-spectrum $P_{\rm L0}(k)$ with a logarithmic slope $n={\rm d}\!\ln P_{\rm L0}/{\rm d}\!\ln k$ such that n>-1 at large scales and n<-1 at small scales to obtain finite results. These constraints are the same as those which appear in the standard perturbative approach. The divergence at long wavelengths eventually cancels out if we only consider the equal-time correlations of the density field, as seen in Jain & Bertschinger (1996) (also Vishniac 1983). As pointed out by these authors, this can be understood from the fact that contributions to the velocity field from long-wave modes correspond to an almost uniform translation of the fluid and therefore should not affect the growth of density structures on small scales. However, these infrared contributions can no longer be ignored when one studies different-times correlations or the velocity field itself. The same problem appears when one considers the collisionless dynamics as described by the Vlasov equation (Valageas 2004). The small-scale divergences are more harmful as they do not cancel in such a fashion. In fact, they may be understood as the signature of anomalous terms (which do not scale as integer powers of $P_{\rm L0}$) which cannot be recovered by the standard perturbative scheme as discussed in Valageas (2002). It is probably necessary to go beyond the hydrodynamical approximation (1)-(3) to cure this small-scale problem (e.g. use the Vlasov equation which does not break down at shell-crossing). On the other hand, note that within the 2PI effective action scheme these constraints apply to the non-linear power G₂₂(k) which may have a different asymptotic behavior than the linear power-spectrum at small scales. Therefore, the range of initial conditions where the integral in Eq. (56) converges at high k may be different within this scheme.

   
5 Direct steepest-descent method

We present in the following sections our numerical results which we compare to fits obtained from direct numerical simulations and to the standard perturbative analysis. We consider an Einstein-de-Sitter cosmology with $\Omega_{\rm m}=1,\Omega_{\Lambda}=0$, a shape parameter $\Gamma=0.5$ and a normalization $\sigma_8=0.51$ for the linear power-spectrum $P_{\rm L0}(k)$ and a Hubble constant H₀=50 km s^-1 Mpc^-1 (i.e. the reduced Hubble parameter is h=0.5) as in Smith et al. (2003) so as to compare our results with N-body simulations. Since we mainly wish to investigate in this article the properties of the large-N expansions described in Sect. 4 we first focus on the Einstein-de-Sitter cosmology where both $G_{\rm L}$ and $R_{\rm L}$ have simple analytical forms (21)-(24). We shall discuss the case of a $\Lambda$CDM cosmology in Sect. 8 below.

We first study in this section the results obtained from the direct steepest-descent method of Sect. 4.1 which involves the auxiliary two-point functions G₀ and R₀. Since the latter are equal to their linear counterparts (47) the correlation G₀ is given by Eq. (21) whereas R₀ is obtained by solving Eq. (40) forward in time, in agreement with the discussion below Eq. (23). For the Einstein-de-Sitter cosmology which we consider here R₀ is actually given by the explicit expression (24).

   
5.1 Self-energy

From G₀ and R₀ we need to compute the self-energy from Eqs. (45), (46). In fact, for the Einstein-de-Sitter cosmology the time-dependence of the self-energy can be factorized as follows using Eqs. (21), (24). First, the auxiliary response R₀ can be written from Eq. (24) as:

 

$$R_0(x_1,x_2) = \theta(\eta_1-\eta_2) \delta_{\rm D}({\vec k}_1-{\... ...t[ {\rm e}^{\eta_1-\eta_2} R_0^+ + {\rm e}^{-3(\eta_1-\eta_2)/2} R_0^- \right],$$ (57)

with:

$$R_0^+= \left(\begin{array}{cc} 3/5 & 2/5 \\ 3/5 & 2/5 \end{ar... ...begin{array}{cc} 2/5 & -2/5 \\ -3/5 & 3/5 \end{array}\right) ,$$ (58)

where we separated the contributions from the growing mode R₀⁺ and the decaying mode R₀^-. Then, from Eq. (53) the self-energy term $\Sigma_0$ which is linear over R₀ also splits into two terms as:

 

$$\Sigma_0(x_1,x_2) = \theta(\eta_1-\eta_2) \delta_{\rm D}({\vec k}... ...eta_1} \Sigma_0^+(k_1) + {\rm e}^{-\eta_1/2+5\eta_2/2} \Sigma_0^-(k_1) \right],$$ (59)

with:

 

$$\Sigma_{0;i_1 i_2}^{\pm}(k) = 4 \sum_{j_1 j_2 l_1 l_2}\int{\rm d}... ...^{\rm s}_{l_1;i_2,l_2}({\vec k},-{\vec q}) \; P_{\rm L0}(q) R_{0;j_1 l_1}^{\pm}$$ (60)

where we used Eq. (21). Besides, from the expressions (15)-(16) of the vertices $\gamma^{\rm s}$ we obtain:

$$\Sigma_0^-(k) = \left(\begin{array}{cc} \Sigma_{0;22}^+ -\Sig... ...;12}^+ \\ -\Sigma_{0;21}^+ \Sigma_{0;11}^+ \end{array}\right)$$ (61)

which is consistent with Eq. (56). In particular, we have:

$$\Sigma_{0;11}^+ + \Sigma_{0;22}^+ = - \frac{4\pi}{3} k^2 \int_0^{\infty} {\rm d}q P_{\rm L0}(q) .$$ (62)

Moreover, at high k we have the asymptotic behaviors:

$$k\rightarrow\infty\!\!: \;\;\; \Sigma_{0;11}^+ \simeq \Sigma_... ...}{5} \frac{4\pi}{3} k^2 \int_0^{\infty} {\rm d}q P_{\rm L0}(q)$$ (63)

$$k\rightarrow\infty\!\!: \;\;\; \Sigma_{0;22}^+ \simeq \Sigma_... ...{5} \frac{4\pi}{3} k^2 \int_0^{\infty} {\rm d}q P_{\rm L0}(q).$$ (64)

Of course, we can check that this is consistent with Eq. (62).

Figure 1: The self-energy terms $\Sigma _{0;ij}^+(k)$ of Eq. (60) as a function of wavenumber k. We display the four components $\Sigma _{0;11}^+$ (solid line), $\Sigma _{0;21}^+$ (dot-dashed line), $\Sigma _{0;12}^+$ (dotted line) and $\Sigma _{0;22}^+$ (dashed line). They are all negative except for $\Sigma _{0;11}^+$ which is positive at k < 0.2 h Mpc-1 and negative at higher k. All terms are of the same magnitude and grow roughly as k2 in agreement with Eq. (62).

Open with DEXTER

We display in Fig. 1 the self-energy term $\Sigma_0^+(k)$as a function of wavenumber k, from which $\Sigma_0$ can be obtained using Eqs. (59), (61). All components $\Sigma_{0;ij}^+$ are of the same magnitude and negative except for $\Sigma _{0;11}^+$ which is positive at large scales $k\la 0.2~h$ Mpc^-1. All terms grow as k² at high wavenumbers and we recover the asymptotic behaviors (63)-(64). Note that the self-energy decreases at large scales so that the two-point functions obtained from Eqs. (42)-(44) will converge to the linear asymptotics

Figure 2: The self-energy $\Sigma _0(k;\eta _1,\eta _2)$ as a function of forward time $\eta _1$, for $\eta _2=-1.4$ (i.e. z2=3) and wavenumbers k=0.1,1 and $10 \times h$ Mpc-1. The line styles are as in Fig. 1. The diagonal components vanish for $\eta _1<\eta _2$ whereas the off-diagonal components vanish for $\eta _1\leq \eta _2$. All terms are negative except for $\Sigma _{0;11}$ which becomes positive at $\eta _1>-0.8$ for k=0.1 h Mpc-1.

Open with DEXTER

We show in Fig. 2 the evolution forward over time $\eta _1$ of the self-energy $\Sigma _0(k;\eta _1,\eta _2)$ for time $\eta _2=-1.4$ (i.e. z₂=3) and wavenumbers k=0.1,1 and $10 \times h$ Mpc^-1 (from bottom to top). The absolute value of the self-energy is larger for higher k in agreement with Fig. 1. All components are negative (except for  $\Sigma _{0;11}$ at $\eta _1>-0.8$ in the case k=0.1 h Mpc^-1) and exhibit a smooth dependence with time following Eq. (59) which goes to  ${\rm e}^{2\eta_1}$ at large $\eta _1$.

	Figure 3: The self-energy $\Sigma _0(k;\eta _1,\eta _2)$ as a function of backward time $\eta _2$, for $\eta _1=0$ (i.e. z1=0) and wavenumbers k=0.1,1 and $10 \times h$ Mpc-1. The line styles are as in Fig. 1. All terms are negative except for $\Sigma _{0;11}$ which becomes positive at $\eta _2<-0.6$ for k=0.1 h Mpc-1.
Open with DEXTER

We show in Fig. 3 the evolution backward over time $\eta _2$ of the self-energy $\Sigma _0(k;\eta _1,\eta _2)$ for time $\eta _1=0$ (i.e. z₁=0) and wavenumbers k=0.1,1 and $10 \times h$ Mpc^-1 (from bottom to top). The behavior agrees with Figs. 1, 2. Following Eq. (59) $\Sigma_0$ converges to a constant for $\eta_2\rightarrow-\infty$.

	Figure 4: The self-energy $\Pi _0(k)$ as a function of wavenumber k. At high k all components are very close and positive. At $k\leq 0.06~h$ Mpc-1 the off-diagonal components become negative.
Open with DEXTER

Finally, using the factorized form (21) of the auxiliary two-point correlation the self-energy $\Pi _0$ can be written as:

$$\Pi_0(x_1,x_2) = \delta_{\rm D}({\vec k}_1+{\vec k}_2) {\rm e}^{2\eta_1+2\eta_2} \Pi_0(k_1) ,$$ (65)

with:

 

$$\Pi_{0;i_1 i_2}(k) = 2 \sum_{j_1 j_2 l_1 l_2} \int{\rm d}{\vec q}... ... q},-{\vec k}+{\vec q}) P_{\rm L0}(q) P_{\rm L0}(\vert{\vec k}-{\vec q}\vert) .$$ (66)

In agreement with the symmetry (49) we can check that $\Pi_{0;21}(k)=\Pi_{0;12}(k)$.

Since the time-dependence (65) is quite simple we only display the self-energy term $\Pi _0(k)$ of Eq. (66) as a function of wavenumber k in Fig. 4. At high k all components are very close and positive. The dependence on wavenumber k follows the smooth growth at small scales of the logarithmic power $\Delta^2_{\rm L0}(k)=k^3 P_{\rm L0}(k)$. As for $\Sigma_0$ the self-energy term $\Pi _0$ decreases very fast at low k so that the two-point functions obtained from Eqs. (42)-(44) will converge to the linear asymptotics at large scales.

   
5.2 Response function

From the self-energy obtained in Sect. 5.1 we compute the response function from Eq. (43). Thanks to the Heaviside factor within the term $\Sigma_0(x,z)$ the r.h.s. only involves earlier times hence Eq. (43) can be integrated forward over time $\eta _1$ at fixed $\eta _2$ to give $R(k;\eta _1,\eta _2)$. Note that each wavenumber k evolves independently thanks to the Dirac factors $\delta_{\rm D}({\vec k}_1-{\vec k}_2)$ within both $\Sigma_0$ and R as in Eq. (50). Thus all processes associated with mode-coupling between different wavenumbers ${\vec k}$ are contained in the calculation of the self-energy $\Sigma_0$ (and $\Pi _0$), as in Eqs. (60), (66).

Equation (43) also writes for $\eta_1>\eta_2$:

  

$$\frac{\partial R_1}{\partial\eta_1} - R_2 = \Sigma_{0;11} . R_1 + \Sigma_{0;12} . R_2$$ (67)

$$\frac{\partial R_2}{\partial\eta_1} - \frac{3}{2} R_1 + \frac{1}{2} R_2 = \Sigma_{0;21} . R_1 + \Sigma_{0;22} . R_2$$ (68)

where (R₁,R₂)=(R₁₁,R₂₁) or (R₁₂,R₂₂), indeed in this case the two columns of the matrix R are not coupled. Using Eqs. (59) and (61) the first Eq. (67) also reads:

 

$$\frac{\partial R_1}{\partial\eta_1}(\eta_1,\eta_2) - R_2(\eta_1,\eta_2) \int_{\eta_2}^{\eta_1} {\rm d}\eta \; \biggl\lbrace \Bigl( {\rm e... ..._{0;11}^+ + {\rm e}^{-\eta_1/2+5\eta/2} \Sigma_{0;22}^+ \Bigr) R_1(\eta,\eta_2)$$ =

$$+ \Bigl( {\rm e}^{2\eta_1} \Sigma_{0;12}^+ - {\rm e}^{-\eta_1/2+5\eta/2} \Sigma_{0;12}^+ \Bigr) R_2(\eta,\eta_2) \biggr\rbrace.$$ (69)

Taking advantage of the simple dependence on time $\eta _1$ of the r.h.s. it is possible to eliminate the integrals over $\eta$ by taking two derivatives of Eq. (69) with respect to $\eta _1$. This yields:

 

$$\frac{\partial^3 R_1}{\partial\eta_1^3} - \frac{3}{2} \frac{\part... ...} + \frac{5}{2} \Sigma_{0;11}^+ R_1 + \frac{5}{2} \Sigma_{0;12}^+ R_2 \right] .$$ (70)

In a similar fashion, one obtains from Eq. (68):

 

$$\frac{\partial^3 R_2}{\partial\eta_1^3} - \frac{3}{2} \frac{\part... ...} + \frac{5}{2} \Sigma_{0;21}^+ R_1 + \frac{5}{2} \Sigma_{0;22}^+ R_2 \right] .$$ (71)

These equations apply to $\eta_1>\eta_2$ and need to be supplemented by initial conditions obtained from Eqs. (67), (68) and Eq. (23) which yield at point $\eta_1=\eta_2^+$:

  

$$\left( \begin{array}{cccccc} R_1 & R_2 & \frac{\partial R_1}{\par... ...{\partial\eta_1^2} & \frac{\partial^2 R_2}{\partial\eta_1^2} \end{array}\right)$$                         X^T = (72)

$$\left( \begin{array}{cccccc} 1 & 0 & 0 & \frac{3}{2} & \frac{3}{2... ...{2\eta_1}(\Sigma_{0;11}^+ + \Sigma_{0;22}^+) & - \frac{3}{4} \end{array}\right)$$ = (73)

for the pair (R₁₁,R₂₁) and:

$$X^T = \left( \begin{array}{cccccc} 0 & 1 & 1 & -\frac{1}{2} &... ...2\eta_1}(\Sigma_{0;11}^+ + \Sigma_{0;22}^+) \end{array}\right)$$ (74)

for the pair (R₁₂,R₂₂). One can easily check that the linear response $R_{\rm L}$ of Eq. (24) is a solution of Eqs. (70), (71) with $\Sigma_0=0$. The advantage of these differential equations over the integro-differential Eqs. (67), (68) is that they no longer require the computation of integrals over past values of the response R. Therefore, one does not need to save the values of R on a grid with a very fine mesh and one can advance over large time steps $\eta _1$. Moreover, each time-step runs faster since one does not need to compute lengthy integrals as in Eqs. (67), (68). We use an adaptative Runge-Kutta algorithm to solve numerically between two grid points the Eqs. (70), (71), written as a system of coupled first-order linear equations over the components of the vector X of Eq. (72).

	Figure 5: The non-linear response $R(k;\eta _1,\eta _2)$ as a function of forward time $\eta _1$, for $\eta _2=-1.4$ (i.e. z2=3) and wavenumbers k=1 ( left panel) and $10 \times h$ Mpc-1 ( right panel). We also plot the linear response $R_{\rm L}$ which shows a simple exponential growth.
Open with DEXTER

We can note that Eqs. (70), (71) are somewhat similar to Bessel functions in terms of the scale factor $a_1={\rm e}^{\eta_1}$ but are of degree three rather than two. In the small-scale limit $\vert\Sigma_0^+\vert\rightarrow\infty$ we can look for an asymptotic solution of the form:

$$R_{1,2}(\eta_1) = \rho_{1,2}(\eta_1) {\rm e}^{{\rm i} \omega {\rm e}^{\eta_1}}$$ (75)

where the functions $\rho_i(\eta_1)$ may be expanded as:

$$\rho_i(\eta_1)= \rho_i^{(0)} + \frac{1}{\omega} \rho_i^{(1)} + \frac{1}{\omega^2} \rho_i^{(2)} + ...$$ (76)

and $\omega^2 \sim \vert\Sigma_0^+\vert$. We obtain at order $\omega^3$ after substituting into Eqs. (70), (71):

$$\omega= \pm \sqrt{-\Sigma_{0;11}^+ - \Sigma_{0;22}^+} = \pm k \left(\frac{4\pi}{3} \int {\rm d}q P_{\rm L0}(q)\right)^{1/2}$$ (77)

where we used Eq. (62). At order $\omega^2$ we obtain:

  

$$\frac{\partial\rho_1^{(0)}}{\partial\eta_1} - \rho_2^{(0)} = 0$$ (78)

$$\frac{\partial\rho_2^{(0)}}{\partial\eta_1} - \frac{3}{2} \rho_1^{(0)} + \frac{1}{2} \rho_2^{(0)} = 0$$ (79)

where we used the asymptotic behavior (63), (64) to simplify factors of the form $\Sigma_{0;ij}^+/\omega^2$. By comparison with Eqs. (67), (68), we see that these equations are precisely the equations obeyed by the linear response $R_{\rm L}=R_0$. Moreover, the initial conditions at $\eta_1=\eta_2$ are precisely the same if we redefine at order 0:

 

$$R_i^{(0)}(\eta_1) \rho_i^{(0)}(\eta_1) \cos~[\omega({\rm e}^{\eta_1}-{\rm e}^{\eta_2})]$$ =

$$\rho_i^{(0)}(a_1) \cos~[\omega(a_1-a_2)] .$$ = (80)

Therefore, at high k the "envelope'' $\rho_i^{(0)}(\eta_1)$ of the non-linear response $R_{ij}(\eta_1)$ follows exactly the linear response $R_{\rm L}$. In addition, the response R exhibits periodic oscillations over the time variable a₁ with a frequency which increases at high wavenumbers as $\omega\propto k$. Therefore, at one-loop order within the direct steepest-descent expansion the memory of initial conditions is not really erased at small scales or late times since the amplitude of the response R does not decay but only exhibits strong oscillations.

	Figure 6: The non-linear response $R(k;\eta _1,\eta _2)$ as a function of backward time $\eta _2$, for $\eta _1=0$ (i.e. z1=0) and wavenumbers k=1 ( left panel) and $10 \times h$ Mpc-1 ( right panel). We also plot the linear response $R_{\rm L}$ which shows a simple exponential growth.
Open with DEXTER

	Figure 7: The non-linear response function $R(k;\eta _1,\eta _2)$ as a function of wavenumber k, at times $\eta _1=0,\eta _2=-1.4$, i.e. z1=0,z2=3. The horizontal lines are the linear response $R_{\rm L}$ of Eq. (24) which is independent of k.
Open with DEXTER

Figure 8: The density two-point correlation $G_{11}(k;\eta _1,\eta _2)$ as a function of time $\eta _1$, for $\eta _2=-1.4$ (i.e. z2=3) and wavenumbers k=1 ( left panel) and $10 \times h$ Mpc-1 ( right panel). We plot the logarithmic power $\Delta ^2=4\pi k^3G_{11}$. We display the linear power $\Delta ^2_{\rm L}$ of Eq. (21) (dashed line), the usual one-loop perturbative result $\Delta ^2_{\rm 1 loop}$ of Eq. (81) (dotted line), the full non-linear power $\Delta ^2$ from Eq. (51) (solid line) and the contribution  $\Delta ^2_i$ from Eq. (52) (dot-dashed line). Note that G11 whence $\Delta ^2$ can be negative, as shown by the oscillations displayed in the figure.

Open with DEXTER

We first display in Fig. 5 the evolution forward over time $\eta _1$ of the non-linear response $R(k;\eta _1,\eta _2)$ for wavenumbers k=1 (left panel) and $10 \times h$ Mpc^-1 (right panel). We do not show the case k=0.1 h Mpc^-1 where the response R closely follows the linear response $R_{\rm L}$ given in Eq. (24). In agreement with the analysis performed through Eqs. (75)-(80) we find that the non-linear response R exhibits oscillations $\sim$ $\cos(\omega a_1)$ with a frequency $\omega$which grows at higher k and an amplitude which follows the linear response $R_{\rm L}$. Moreover, all components R_ij oscillate in phase.

Next, we show in Fig. 6 the evolution backward over time $\eta _2$of the non-linear response R. It shows a behavior similar to the forward evolution displayed in Fig. 5 with respect to $\eta _1$. At large scales we recover the smooth increase (in absolute values) at early times $\eta _2$ of linear theory (24) whereas at small scales we obtain fast oscillations $\sim$ $\cos(\omega a_2)$ related to Eqs. (75)-(80) with an amplitude which follows again the linear response $R_{\rm L}$.

Finally, since at equal times the response function obeys Eq. (23) we only display in Fig. 7 the response at unequal times $(\eta _1=0,\eta _2=-1.4)$. At low k we recover the linear response $R_{\rm L}$. At higher k above 0.2 h Mpc^-1 the non-linear response obtained at one-loop order within this steepest-descent approach departs from the linear prediction and shows oscillations with increasingly high frequencies while their amplitude follows the linear response $R_{\rm L}$ at high k. This behavior is a result of the oscillatory behavior with time $\eta _1$displayed in Fig. 5 above and analyzed in Eqs. (75)-(80).

   
5.3 Two-point correlation

Finally, from Eq. (51) we compute the non-linear two-point correlation function G obtained through the steepest-descent method at one-loop order. We compare our results in Figs. 8-11 at redshifts z=0,3 with a fit $G_{\rm nb}$ from numerical simulations (Smith et al. 2003), the linear prediction $G_{\rm L}$ (Eq. (21)) and the usual one-loop result $G_{\rm 1 loop}$ obtained from the standard perturbative expansion over powers of the linear density field. The standard one-loop power-spectrum can be written (e.g., Jain & Bertschinger 1994; Scoccimarro & Frieman 1996a,b) as:

$$P_{\rm 1 loop} = P_{\rm L} + P_{22} + P_{13} ,$$ (81)

where $P_{\rm L}$ is the linear contribution and P₂₂ and P₁₃ are of order  $P_{\rm L}^2$ with:

 

$$P_{22}(k;\eta_1,\eta_2) = {\rm e}^{2\eta_1+2\eta_2} \int{\rm d}{\... ...0}(q) P_{\rm L0}({\vec k}-{\vec q}) 2 F_2^{(s)}({\vec q},{\vec k}-{\vec q})^2 ,$$ (82)

and:

 

$$P_{13}(k;\eta_1,\eta_2) = \frac{1}{2} \left( {\rm e}^{3\eta_1+\et... ...(k) \int{\rm d}{\vec q}P_{\rm L0}(q) 6 F_3^{(s)}({\vec q},-{\vec q},{\vec k}) .$$ (83)

The angular kernels F₂^(s) and F₃^(s) are obtained from the couplings $\alpha$ and $\beta$ of Eq. (6). In particular we have (Goroff et al. 1986):

$$F_2^{(s)}({\vec k}_1,{\vec k}_2) = \frac{5}{7}+\frac{{\vec k}... ...{7} \left( \frac{{\vec k}_1.{\vec k}_2}{k_1k_2}\right)^2 \cdot$$ (84)

On the other hand, after integration over angles in Eq. (83) one obtains (Makino et al. 1992):

 

$${P_{13}(k;\eta_1,\eta_2) = \frac{1}{2} \left( {\rm e}^{3\eta_1+\eta_2} + {\rm e}^{\eta_1+3\eta_2} \right) P_{\rm L0}(k) }$$

$$\times \pi \int{\rm d}q P_{\rm L0}(q) q^2 \left[ \frac{6k^6-79k^4q^2+50k^2q^4-21q^6}{63k^2q^4} \right.$$

$$+ \left. \frac{(q^2-k^2)^3(7q^2+2k^2)}{42k^3q^5} \ln\left\vert\frac{k+q}{k-q}\right\vert \right] \cdot$$ (85)

From Eqs. (82)-(85) we can see that at large wavenumbers $k\rightarrow\infty$ the integrals over q in Eqs. (82), (85) are dominated by the region where the slope n of the linear power-spectrum is of order $n \simeq -1$. This yields:

 

$$k\rightarrow\infty\!\!: \;\; P_{\rm 1 loop} \simeq - {\rm e}^{2\e... ...\eta_2)-1\right) k^2 P_{\rm L0}(k) \frac{4\pi}{3} \int {\rm d}q P_{\rm L0}(q) .$$ (86)

Therefore, at unequal times $\eta_1\neq\eta_2$ the standard one-loop prediction becomes negative at high k and the logarithmic power (as defined below in Eq. (87)) grows as $\Delta_{\rm 1 loop}^2 \propto - {\rm e}^{3\eta_1} k^2$ with k and $\eta _1$, in agreement with Figs. 8, 11. Note that the leading term (86) vanishes at equal times $\eta_1=\eta_2$. This is related to the cancellation of the infrared divergences which appear for linear power-spectra $P_{\rm L}\propto k^n$ with $n \leq -1$ when one considers equal-times statistics (Vishniac 1983; Jain & Bertschinger 1996).

The dependence ${\rm e}^{3\eta_1} k^2$ merely reflects the order $P_{\rm L}^2$ of the perturbative expansion. In order to go deeper into the non-linear regime one may consider higher order terms but the latter grow increasingly fast and one would need at least partial resummations to obtain a well-behaved series (e.g. Crocce & Scoccimarro 2006a). This is precisely what the large-N expansions discussed in this article attempt to perform.

Figure 9: The logarithmic power $\Delta ^2(k)$ of Eq. (19) at redshift z=0, that is for equal times $\eta _1=\eta _2=0$ (z1=z2=0). We display the linear power $\Delta ^2_{\rm L}$ (dashed line), the usual one-loop perturbative result $\Delta ^2_{\rm 1 loop}$ (dotted line), the full non-linear power $\Delta ^2$ from Eq. (51) (solid line) and the contribution $\Delta ^2_i$ from Eq. (52) (dot-dashed line). We also show the fit $\Delta ^2_{\rm nb}$ (upper solid line) from numerical simulations (Smith et al. 2003). All quantities (including $\Delta ^2_i$ which is the square of an oscillating function) are positive.

Open with DEXTER

	Figure 10: The logarithmic power $\Delta ^2(k)$ as in Fig. 9 but at redshift z=3.
Open with DEXTER

	Figure 11: The logarithmic power $\Delta ^2(k)$ as in Fig. 9 but at unequal times $(\eta _1=0,\eta _2=-1.4)$ (i.e. z1=0,z2=3).
Open with DEXTER

We first show in Fig. 8 the evolution over time $\eta _1$ of the density two-point correlation $G_{11}(k;\eta _1,\eta _2)$, which is symmetric with respect to $\eta_1\leftrightarrow\eta_2$. We actually display the "logarithmic power'' $\Delta^2(k;\eta_1,\eta_2)$ defined as in Eq. (19) by:

$$\Delta^2(k;\eta_1,\eta_2) = 4\pi k^3 G_{11}(k;\eta_1,\eta_2) .$$ (87)

At equal times $\eta_1=\eta_2$ it is equal to the usual logarithmic power and we have $\Delta^2>0$. However, at different times $\eta_1\neq\eta_2$we can have G₁₁<0 whence $\Delta^2<0$. We present the cases k=1 (left panel) and $k=10\times h$ Mpc^-1 (right panel) in Fig. 8 for $\eta _2=-1.4$. For clarity we only display the density-density correlation G₁₁. Other components G_ij show a similar behavior.

The oscillations show that G₁₁ can indeed be negative at unequal times. At these small scales the non-linear power $\Delta ^2$ is already significantly different from the linear power  $\Delta ^2_{\rm L}$ but we can check in the left panel that $\Delta ^2$ matches the usual one-loop power  $\Delta ^2_{\rm 1 loop}$of Eq. (81) at small $\Delta ^2$. As recalled in Sect. 4.3 the one-loop large-N expansions indeed match the usual perturbative expansion up to one-loop order (i.e. $P_{\rm L}^2$). At later times while  $\Delta ^2_{\rm 1 loop}$ grows very fast as ${\rm e}^{3\eta_1}$ from Eq. (83) the prediction $\Delta ^2$ of the one-loop steepest-descent approach remains well-behaved of order $\Delta ^2_{\rm L}$. Thus, we can see in the right panel that at very small scales where $\Delta ^2_{\rm 1 loop}$ has become exceedingly large $\Delta ^2$ remains well-controlled. In fact, following the behavior of the response function analyzed in Sect. 5.2 it exhibits fast oscillations with an amplitude close to $\Delta ^2_{\rm L}$. We also show in Fig. 8 the power $\Delta ^2_i$ obtained from the first term alone of Eq. (51) which is detailed in Eq. (52). This corresponds to the "non-linear transport'' of the initial linear fluctuations without taking into account the fluctuations $\Pi$ generated at all times by the non-linear dynamics. We can see that this term is sub-dominant in the mildly non-linear regime (left panel) but happens to be again of same order as $\Delta ^2$in the highly non-linear regime.

Next, we show in Fig. 9 the logarithmic power $\Delta ^2(k)$today at equal redshifts z₁=z₂=0 (i.e. equal times $\eta _1=\eta _2=0$). In agreement with Eqs. (51), (52) and the subsequent discussion both $\Delta ^2$ and $\Delta ^2_i$ are positive ( $\Delta ^2_i$ is the square of an oscillating function). At large scales $\Delta ^2$ converges to the linear regime $\Delta ^2_{\rm L}$. At smaller scales we can check again that $\Delta ^2$ first follows the usual one-loop power $\Delta ^2_{\rm 1 loop}$ and next goes back to $\Delta ^2_{\rm L}$ whereas  $\Delta ^2_{\rm 1 loop}$ keeps increasing very fast in magnitude. Unfortunately, $\Delta ^2$ does not agree with the fit $\Delta ^2_{\rm nb}$ from numerical simulations (Smith et al. 2003) on a larger range of k as compared with  $\Delta ^2_{\rm 1 loop}$ (although the rough agreement between $\Delta ^2_{\rm nb}$ and $\Delta ^2_{\rm 1 loop}$ at small scale is purely coincidental). On the other hand, whereas adding higher-order terms in the standard perturbative expansion may not improve much the agreement with N-body simulations (especially since it would lead to increasingly steep terms at high k) the series obtained in the large-N approach is likely to be well-behaved as various terms do not "explode'' at small scales. However, this would require intricate calculations. We can note that although $\Delta ^2_i$ is of order $\Delta ^2$ one cannot neglect the fluctuations $\Pi$ generated by the non-linear dynamics and the sum yields a smooth power $\Delta ^2$.

We also display in Fig. 10 the logarithmic power $\Delta ^2(k)$at equal redshifts z₁=z₂=3. We recover the same behaviors as those obtained at redshift z=0.

Finally, we show in Fig. 11 the logarithmic power $\Delta^2(k;\eta_1,\eta_2)$ at unequal times $\eta _1=0,\eta _2=-1.4$, that is at redshifts z₁=0,z₂=3. We find again that $\Delta ^2$ matches the usual one-loop power $\Delta ^2_{\rm 1 loop}$ at large scales and follows its change of sign at $k\sim 0.8~h$ Mpc^-1. Note indeed that at unequal times $\Delta ^2$ need not be positive. At smaller scales, the steepest-descent large-N result $\Delta ^2$ departs from the fast growing $\Delta ^2_{\rm 1 loop}$ and shows a series of oscillations of moderate amplitude, which again follow $\Delta ^2_{\rm L}$. These features suggest that the first change of sign in Fig. 11, which is common to both $\Delta ^2_{\rm L}$ and $\Delta ^2_{\rm 1 loop}$, may be real and not a mere artefact of the one-loop expansion. Note that in spite of the oscillations with time seen in Figs. 8, 11 the large-N approach automatically ensures that $\Delta^2>0$ at equal times thanks to the structure of Eq. (51). This property holds at all orders over 1/N since Eq. (51) is independent of the expressions of the self-energies $\Sigma,\Pi$.

   
6 2PI effective action approach

We now present the results obtained at one-loop order from the 2PI effective action approach recalled in Sect. 4.2. Thus, we need to solve the system of coupled Eqs. (42)-(46) where in the expression (45)-(46) of the self-energy the auxiliary two-point functions G₀ and R₀ are replaced by G and R. Thanks to causality, which leads to the Heaviside factor $\theta(\eta_1-\eta_2)$ of Eq. (50) within both R and ${\Sigma }$, we solve the system (42)-(46) by moving forward over time.

Thus, we use a grid $\eta^{(n)}$ with $n=1,..,N_{\eta}$for the time variables. At the earliest time-step $\eta^{(1)}$ we initialize all matrices by their linear value at $\eta_1=\eta_2=\eta^{(1)}$. Next, once we have obtained all two-point functions up to time $\eta^{(n)}$ (i.e. over $\eta_1\leq\eta^{(n)}$ and $\eta_2\leq\eta^{(n)}$, initially n=1) we advance to the next time-step $\eta^{(n+1)}$ as follows. The response R at equal times $\eta_1=\eta_2=\eta^{(n+1)}$ is first obtained from Eq. (23). Next, we move backward over time $\eta_2=\eta^{(n)},...,\eta^{(1)}$ at fixed $\eta_1=\eta^{(n+1)}$ by using Eqs. (43) and (45). That is, to compute $R(\eta^{(n+1)},\eta^{(i)})$ with i=n,..,1 we use for each value of i the integro-differential Eq. (43) to move over $\eta _1$ from $\eta_1=\eta^{(n)}$ up to $\eta_1=\eta^{(n+1)}$ at fixed $\eta_2=\eta^{(i)}$. Thanks to the Heaviside factors the r.h.s. of Eq. (43) only involves $\Sigma(\eta_1,\eta)$ and $R(\eta,\eta_2)$ with $\eta_1\geq\eta\geq\eta_2$ which are already known. Besides, once $R(\eta^{(n+1)},\eta^{(i)})$ has been obtained we compute $\Sigma(\eta^{(n+1)},\eta^{(i)})$ from Eq. (45) before moving downward to step i-1. Since R and ${\Sigma }$ vanish for $\eta_2>\eta_1$ we have actually obtained in this fashion R and ${\Sigma }$ for all $\eta_1,\eta_2\leq\eta^{(n+1)}$ (indeed $R(\eta^{(i)},\eta^{(n+1)})=0$ for $i\leq n$).

Next, we compute $\Pi$ and G from Eqs. (46) and (51). At fixed $\eta_1=\eta^{(n+1)}$ we now move forward over time $\eta_2=\eta^{(1)},...,\eta^{(n)}$. Indeed, the r.h.s of Eq. (51) only involves $\Pi(\eta_1',\eta_2')$ with $\eta_1'\leq\eta_1$ and $\eta_2'\leq\eta_2$ which is already known. Besides, once $G(\eta^{(n+1)},\eta^{(i)})$ has been obtained we compute $\Pi(\eta^{(n+1)},\eta^{(i)})$ from Eq. (46) before moving forward to step i+1 and we use the symmetry (49) to derive $G(\eta^{(i)},\eta^{(n+1)})$ as well as $\Pi(\eta^{(i)},\eta^{(n+1)})$.

In fact, at each step i for $(\eta^{(n+1)},\eta^{(i)})$ the integrals in the r.h.s. of Eqs. (43), (51) also involve the values $R(\eta^{(n+1)},\eta^{(i)})$ and $G(\eta^{(n+1)},\eta^{(i)})$ which are being computed if we use the boundary points in the numerical computation of the time-integrals. Rather than using open formulae for the integrals we perform a few iterations at each time-step $\eta^{(n)}$. This procedure converges over a few loops.

Finally, to speed-up the numerical computation we use finite elements over Fourier space ${\vec k}$ to store the structure of the self-energy terms ${\Sigma }$and $\Pi$. Thus, since all matrices only depend on the wavenumber modulus $k=\vert{\vec k}\vert$ we use a grid k⁽ⁿ⁾ with n=1,...,N_k. Next, in order to interpolate for all values of k we can write for instance any matrix such as G(k) as:

$$G(k) = \sum_{n=1}^{N_k} G^{(n)} \chi^{(n)}(k) ,$$ (88)

with:

$$G^{(n)}=G(k^{(n)}) \hspace{0.3cm} \mbox{and} \hspace{0.3cm} \chi^{(i)}(k^{(j)})= \delta_{ij} ,$$ (89)

where $\delta_{ij}$ is the Kronecker symbol. The basis functions $\chi^{(n)}(k)$ define our interpolation scheme. We could use basis functions which all extend over the full range $0<k<\infty$ (i.e. a spectral method built from cardinal functions) but we prefer here to use simple finite elements with $\chi^{(n)}(k)=0$ if k<k^(n-1) or k>k⁽ⁿ⁺¹⁾. For convenience we actually use triangular functions over the variable $\kappa=(k-1)/(k+1)$ which maps $]0,\infty[$ to ]-1,1[:

$$\chi^{(n)}(\kappa)= \left\{ \begin{array}{ccc} {\displaystyl... ...kappa^{(n)} \leq \kappa \leq \kappa^{(n+1)} \end{array}\right.$$ (90)

which is equivalent to use linear interpolation over $\kappa$ for G(k). Then, we can take advantage of the fact that the kernels $\gamma^{\rm s}_{i;i_1,i_2}({\vec k}_1,{\vec k}_2)$ in Eqs. (53), (54) do not depend on time to compute once for all the integrals over wavenumber needed for the self-energy ${\Sigma }$ and $\Pi$. Thus, we define the quantities:

 

$$\Sigma_{\nu_1 \nu_2}^{(n;n_1,n_2)} = 4 \int{\rm d}{\vec q}\; \gam... ...{\vec k}^{(n)}) \chi^{(n_1)}(q) \chi^{(n_2)}(\vert{\vec k}^{(n)}-{\vec q}\vert)$$ (91)

and:

 

$$\Pi_{\nu_1 \nu_2}^{(n;n_1,n_2)} = 2 \int{\rm d}{\vec q}\; \gamma^... ...{(n)}+{\vec q}) \chi^{(n_1)}(q) \chi^{(n_2)}(\vert{\vec k}^{(n)}-{\vec q}\vert)$$ (92)

where we noted $\nu=(i;i_1,i_2)$ the indices of the vertices $\gamma^{\rm s}_{i;i_1,i_2}$. These quantities $\Sigma_{\nu_1 \nu_2}^{(n;n_1,n_2)}$and $\Pi_{\nu_1 \nu_2}^{(n;n_1,n_2)}$ only depend on the vertices $\gamma^{\rm s}$given in Eqs. (15), (16) and on the grid which we use over wavenumbers hence they are computed at the beginning of the code within initialization subroutines. Next, the self-energy at a given time $\eta_1,\eta_2$ is obtained as:

 

$$\Sigma_{i_1 i_2}(k^{(n)};\eta_1,\eta_2) = \sum_{j_1 j_2 l_1 l_2} ... ...n_1,n_2)} R_{j_1 l_1}^{(n_1)}(\eta_1,\eta_2) G_{j_2 l_2}^{(n_2)}(\eta_1,\eta_2)$$ (93)

with $\nu_1=(i_1;j_1,j_2)$ and $\nu_2=(l_1;i_2,l_2)$, and:

 

$$\Pi_{i_1 i_2}(k^{(n)};\eta_1,\eta_2) = \sum_{j_1 j_2 l_1 l_2} \su... ...n_1,n_2)} G_{j_1 l_1}^{(n_1)}(\eta_1,\eta_2) G_{j_2 l_2}^{(n_2)}(\eta_1,\eta_2)$$ (94)

with $\nu_1=(i_1;j_1,j_2)$ and $\nu_2=(i_2;l_1,l_2)$. Therefore, using Eqs. (93), (94) we obtain the self-energy by advancing over times $\eta_1,\eta_2$ as described above without performing new integrations over wavenumbers at each step (they are replaced by the discrete sums over n₁,n₂ with weights $\Sigma_{\nu_1 \nu_2}^{(n;n_1,n_2)}$or $\Pi_{\nu_1 \nu_2}^{(n;n_1,n_2)}$).

   
6.1 Response R and self-energy ${\Sigma }$

Figure 12: The non-linear response $R(k;\eta _1,\eta _2)$ as a function of forward time $\eta _1$, for $\eta _2=-1.4$ (i.e. z2=3) and wavenumbers k=1 ( left panel) and $10 \times h$ Mpc-1 ( right panel), as in Fig. 5 but for the 2PI effective action method. For comparison we also plot the first half-oscillation of the response obtained in Sect. 5 from the direct steepest-descent method (curves labeled "sd'') which was shown in Fig. 5.

Open with DEXTER

We first display in Fig. 12 the evolution forward over time $\eta _1$of the response $R(k;\eta _1,\eta _2)$. We can see that the non-linear response exhibits oscillations as for the steepest-descent result of Fig. 5 but its amplitude now decays at large times $\eta _1$ instead of following the linear envelope. This can be understood analytically from the following simple model. Let us consider the equation:

$$\frac{\partial R}{\partial\eta_1} = \sigma \int_{\eta_2}^{\et... ...d}\eta \; {\rm e}^{\eta_1+\eta} R(\eta_1,\eta) R(\eta,\eta_2)$$ (95)

for a one-component response $R(\eta_1,\eta_2)$. This is a simplified form of the system (43), (45) where we used the linear time-dependence (21) for G to obtain a closed equation for R. The parameter $\sigma\leq 0$ represents the amplitude of the self-energy ${\Sigma }$ at the wavenumber of interest. Then, the linear regime $R_{\rm L}$ corresponds to $\sigma=0$ and yields the solution:

$$\eta_1 \geq \eta_2\!\!: \;\;\; R_{\rm L}(\eta_1,\eta_2) = 1 ,$$ (96)

where we used the initial condition $R(\eta,\eta)=1$ at equal times. Next, the steepest-descent method presented in Sects. 4.1 and 5 amounts to replace $R(\eta_1,\eta)$ by $R_{\rm L}(\eta_1,\eta)$ into the integral in the r.h.s. of Eq. (95). This yields a linear equation for Rwhich can be transformed into a differential equation as in Sect. 5.2 by taking the derivative of Eq. (95) with respect to $\eta _1$. Since we only have one component one differentiation is sufficient to remove the integral and we obtain:

$$\frac{\partial^2 R}{\partial\eta_1^2} - \frac{\partial R}{\pa... ...hence} \;\;\; \frac{\partial^2 R}{\partial a_1^2} = \sigma R ,$$ (97)

where we introduced the scale-factor $a_1={\rm e}^{\eta_1}$ as the new time-variable. Therefore, we obtain the solution:

$$R(a_1,a_2) = \cos~[\omega(a_1-a_2)] \;\;\; \mbox{with} \;\;\; \omega= \sqrt{-\sigma} .$$ (98)

Thus, as in Sect. 5.2 we obtain within the steepest-descent method an oscillating response function with an amplitude given by the linear response $R_{\rm L}$ and a frequency over the time-variable a₁ given by $\sqrt{-\sigma}$. Finally, in the 2PI effective action method we need to solve the non-linear Eq. (95) which also reads:

$$\frac{\partial R}{\partial a_1} = \sigma \int_{a_2}^{a_1} {\rm d}a \; R(a_1,a) R(a,a_2) ,$$ (99)

where we changed variables from $\eta$ to $a={\rm e}^{\eta}$. We can look for a solution of Eq. (99) of the form R(a₁,a₂)=R(a₁-a₂)so that the r.h.s. of Eq. (99) becomes a simple convolution:

$$\frac{{\rm d}R}{{\rm d}a} = \sigma \int_0^a {\rm d}a' \; R(a-a') R(a') .$$ (100)

Introducing the Laplace transform ${\tilde R}(s)$ defined by:

$${\tilde R}(s)= \int_0^{\infty} {\rm d}a \; {\rm e}^{-s a} R(a) ,$$ (101)

and taking the Laplace transform of Eq. (100) we obtain:

$$s {\tilde R}(s) - R(0) = \sigma {\tilde R}(s)^2 .$$ (102)

Using R(0)=1 and requiring that ${\tilde R}(s)$ vanishes for $s\rightarrow +\infty$ yields:

$${\tilde R}(s) = \frac{\sqrt{s^2+4\omega^2}-s}{2\omega^2} ,$$ (103)

where we introduced $\omega$ as defined in Eq. (98). This is a well-known Laplace transform and Eq. (103) directly gives:

$$R(a_1,a_2)= \frac{J_1[2\omega(a_1-a_2)]}{\omega(a_1-a_2)} ,$$ (104)

where J₁ is a Bessel function of the first kind. Thus, we see that within the 2PI effective action method the response function exhibits oscillations as in the direct steepest-descent method. However, their amplitude is no longer given by the linear response and becomes much smaller at late times along with a higher frequency. For the simple model (95) the amplitude of the response function decays as a₁^-3/2 instead of being constant while the frequency becomes twice larger. This explains the qualitative behavior seen in Fig. 12. In particular, note that in agreement with the simple model (95)-(104) the comparison of the right panel of Fig. 12 with Fig. 5 shows that the frequency of the oscillations is indeed about twice larger in the non-linear regime for the 2PI effective action result than for the steepest-descent prediction. Of course, the result shown in Fig. 12 is more complicated than the simple model (95) since the two-point correlation Gwhich enters ${\Sigma }$ whence R also depends non-linearly on R through Eq. (42) or the more explicit expression (51). However, as we shall check in Sect. 7.3 the dependence of ${\Sigma }$ on G does not change much the behavior of the response R(see Fig. 29).

On the other hand, we can see in Fig. 12 that the first quarter of oscillation, from 1 to 0, of the response R follows rather closely the curves labeled "sd'' obtained from the steepest-descent method presented in Sect. 5 which were also shown in Fig. 5. Besides this agreement remains valid at quite large k (right panel) where two-point functions obtained from both methods are generically very different as shown by the figures below. This property can be understood from Eq. (56) which shows that the dependence on wavenumber $\Sigma(k;\eta,\eta)\propto k^2$of the self-energy ${\Sigma }$ at equal times is identical for both approaches (at one-loop order). Moreover, Eq. (56) shows that the normalization is governed by the value of the two-point correlation $G_{22}(k';\eta,\eta)$ at the point where the logarithmic slope is n=-1(this integral actually corresponds to the mean square velocity $\langle v^2\rangle$). At redshift z=0 this wavenumber is still within the linear regime (for the SCDM cosmology which we consider here) therefore the self-energy ${\Sigma }$ at equal times $\eta_1=\eta_2$ obtained within the steepest-descent method and the 2PI effective action approach are very close until z=0. Then, the early time-evolutions at $\eta_1\simeq\eta_2$ of the response R obtained within both methods from Eq. (43) are very close. We can see from Fig. 12 that this agreement holds until the response R first vanishes. Beyond this point the steepest-descent method yields increasingly large oscillations (Fig. 5) whereas the 2PI effective action yields small oscillations which decay to zero (Fig. 12).

Next, we show in Fig. 13 the evolution backward over time $\eta _2$of the response R. In a manner consistent with Fig. 12 it exhibits oscillations which are strongly damped at large time separations |a₁-a₂|, whence at early times $\eta _2$, while the frequency is also larger than for the steepest-descent result shown in Fig. 6. Again the behavior at nearly equal times $\eta_2\simeq\eta_1$ is close to the prediction of the steepest-descent method presented in Sect. 5 until the response first vanishes. The response at unequal times $(\eta _1=0,\eta _2=-1.4)$ is shown as a function of wavenumber k in Fig. 14. At large scales we recover the linear regime whereas at small scales we obtain oscillations which show a fast decay into the non-linear regime. This is consistent with the time-dependence displayed in Figs. 12, 13.

Therefore, in contrast with the steepest-descent approach we find that within the 2PI effective action method the memory of initial conditions is in some sense erased as the response function decays for large time separation. This agrees with the intuitive expectation that within the real collisionless gravitational dynamics the details of the initial conditions are erased at small scales. Indeed, after shell-crossing one can expect for instance that virialization processes build halos which mainly depend on a few integrated quantities which characterize the collapsed region (such as the mass, initial overdensity, angular momentum, ...) and exhibit relaxed cores (e.g. isotropic velocity distributions), as suggested by numerical simulations. However, we must note that the system studied in this paper is defined from the hydrodynamical equations of motion (1)-(3) which break down at shell-crossing. Therefore, there is no guarantee a priori that such small-scale relaxation processes should come out from Eqs. (1)-(3). It appears through the 2PI effective action method presented in this paper and Figs. 12-14 that this is actually the case. In fact, studying the large-k behavior of the kernel $K_{\rm s}$ Crocce & Scoccimarro (2006b) managed to resum all dominant diagrams in this large-k limit and obtained a Gaussian decay $R(k) \sim {\rm e}^{-k^2}$. In our case, the simple model (95) suggests that the one-loop 2PI effective action method only yields a power-law decay as in Eq. (104). Since at the one-loop order considered in this article we do not resum all diagrams it is not really surprising that we do not recover the "exact'' Gaussian decay. However, note that the large-N expansion does not give a mere Taylor series expansion of such a Gaussian decay (which would blow up at large k or large times) but already captures at a qualitative level the damping of the response R at small scale or at large time separation thanks to the non-linearity of the evolution equation obeyed by R.

	Figure 13: The non-linear response $R(k;\eta _1,\eta _2)$ as a function of backward time $\eta _2$, for $\eta _1=0$ (i.e. z1=0) and wavenumbers k=1 ( left panel) and $10 \times h$ Mpc-1 ( right panel). We also plot in the left panel the first half-oscillation of the response obtained in Sect. 5 from the direct steepest-descent method (curves labeled "sd'') which was shown in Fig. 6.
Open with DEXTER

	Figure 14: The non-linear response function $R(k;\eta _1,\eta _2)$ as a function of wavenumber k, at times $\eta _1=0,\eta _2=-1.4$. We also plot the first half-oscillation of the response obtained in Sect. 5 from the direct steepest-descent method (curves labeled "sd'') which was shown in Fig. 7.
Open with DEXTER

	Figure 15: The self-energy $\Sigma (k;\eta _1,\eta _2)$ as a function of forward time $\eta _1$, for $\eta _2=-1.4$ (i.e. z2=3) and wavenumbers k=1 ( left panel) and $10 \times h$ Mpc-1 ( right panel). The curves labeled "sd'' are the prediction of the steepest-descent method as in Fig. 2.
Open with DEXTER

	Figure 16: The self-energy $\Sigma (k;\eta _1,\eta _2)$ as a function of backward time $\eta _2$, for $\eta _1=0$ (i.e. z1=0) and wavenumbers k=1 ( left panel) and $10 \times h$ Mpc-1 ( right panel). The curves labeled "sd'' are the prediction of the steepest-descent method as in Fig. 3.
Open with DEXTER

	Figure 17: The self-energy $\Sigma (k)$ as a function of wavenumber k at unequal redshifts z1=0,z2=3. The curves labeled "sd'' are the prediction of the steepest-descent method.
Open with DEXTER

We display in Figs. 15 and 16 the evolution with times $\eta _1$ and $\eta _2$ of the self-energy ${\Sigma }$. At large scales (k<0.1 h Mpc^-1) and small time separations we recover the results obtained within the direct steepest-descent method but at smaller scales we obtain a series of oscillations with an amplitude which decreases at large time separations. This follows from the behavior of the non-linear response R analyzed in Figs. 12-14 above. Indeed, the self-energy ${\Sigma }$ depends linearly on the response R, see Eq. (93). Again the behavior at nearly equal times $\eta_1\simeq\eta_2$ is similar for both approaches (see in particular the left panel of Fig. 15). The right panel of Fig. 15 shows that at high wavenumbers the envelope of the oscillations of ${\Sigma }$ decays as a power-law of $a_1={\rm e}^{\eta_1}$ at late times. This agrees with the analysis of the simple model (95). In a similar fashion, Fig. 17 which displays the dependence on wavenumber k of $\Sigma (k;\eta _1,\eta _2)$, at fixed times $\eta _1=0,\eta _2=-1.4$, matches the results obtained from the direct steepest-descent method at low k (except for  $\Sigma_{12}$) and exhibits decaying oscillations at high k. We can note in Fig. 17 that the term $\Sigma_{12}$ (dotted line) does not match the steepest-descent result in the limit $k \rightarrow 0$. This is due to the behavior of the associated kernels $\Sigma_{\nu_1 \nu_2}^{(n;n_1,n_2)}$ of Eq. (91). More precisely, using Eqs. (15)-(15), Eq. (53) reads in this case:

$$\Sigma_{12}(k) \int{\rm d}{\vec q}\biggl \lbrace R_{11}(q) G_{21}(\vert{\vec k}-... ... k}.({\vec k}-{\vec q})]({\vec k}.{\vec q})}{k^2 \vert{\vec k}-{\vec q}\vert^2}$$ =

$$- R_{12}(q) G_{22}(\vert{\vec k}-{\vec q}\vert) \frac{[{\vec k}.({\vec k}-{\vec q})]^2 q^2}{k^2 \vert{\vec k}-{\vec q}\vert^4}$$

$$+ R_{21}(q) G_{11}(\vert{\vec k}-{\vec q}\vert) \frac{({\vec k}.{\vec q})^2}{k^2q^2}$$

$$- R_{22}(q) G_{12}(\vert{\vec k}-{\vec q}\vert) \frac{({\vec k}.{... ...\vec k}.({\vec k}-{\vec q})]}{k^2 \vert{\vec k}-{\vec q}\vert^2} \biggl \rbrace$$ (105)

where we did not write the time-dependence. At large wavenumbers q the integrand $S_{12}({\vec q})$ simplifies to:

$$S_{12}({\vec q}) \sim \frac{({\vec k}.{\vec q})^2}{k^2q^2} [ R_{21} G_{11} + R_{22} G_{12} - R_{11} G_{21} - R_{12} G_{22}] .$$ (106)

Therefore, in the linear regime (or for the steepest-descent self-energy $\Sigma_0$) where $G_{ij}(q;\eta_1,\eta_2) = {\rm e}^{\eta_1+\eta_2} P_{\rm L0}(q)$the two-point correlations G_ij factorize and we have:

$$q \gg \! k \! : S_{0;12}({\vec q}) \sim P_{\rm L0}(q) [R_{\rm L;21}+R_{\rm L;22}-R_{\rm L;11}-R_{\rm L;12}].$$ (107)

For a CDM-like power-spectrum $P_{\rm L0}(q) \sim q^{-3}\ln q$ at large wavenumbers while the linear response $R_{\rm L}(q)$ is constant with respect to q. Hence the integration over ${\vec q}$ yields at first sight a logarithmic divergence for $q\rightarrow\infty$. In fact, we can check from Eq. (24) that for the linear response the terms in the bracket in Eq. (107) actually cancel so that the integral is well-defined for the linear two-point functions (the integrand decreases as 1/q⁴ at large q) and it is dominated by wavenumbers $q \sim k$, in agreement with the results obtained in Sect. 5.1. However, for the 2PI effective action method where we must use the non-linear two-point functions in Eq. (106) this cancellation no longer holds at weakly non-linear scales $q_{\rm NL}$ where they depart from the linear predictions. The integral still converges because of the decay of the non-linear response at high q but this implies that for $k \rightarrow 0$ the self-energy term $\Sigma_{12}$ is dominated by the contribution of wavenumbers $q \sim q_{\rm NL}$ associated with the non-linear transition and not by wavenumbers $q \sim k$. Therefore, $\Sigma_{12}$ does not converge to the steepest-descent prediction $\Sigma_{0;12}$ at large scales $k \rightarrow 0$. By contrast, for other terms the integrand behaves as:

$$q \! \gg \! k \! : \; S_{11}({\vec q}) \! \sim \! S_{22}({\ve... ...} R G , \;\; S_{21}({\vec q}) \! \sim \! \frac{k^2}{q^2} R G .$$ (108)

The additional power of 1/q ensures that the logarithmic divergence of Eq. (107) does not appear and at large scales the self-energy $\Sigma_{ij}$ is dominated by the contribution from wavenumbers $q \sim k$hence it matches the results obtained within the steepest-descent approach.

   
6.2 Correlation G and self-energy $\Pi$

Figure 18: The density two-point correlation $G_{11}(k;\eta _1,\eta _2)$ as a function of time $\eta _1$, for $\eta _2=-1.4$ (i.e. z2=3) and wavenumbers k=1 ( left panel) and $10 \times h$ Mpc-1 ( right panel). We plot the logarithmic power $\Delta ^2=4\pi k^3G_{11}$. We display the full non-linear power $\Delta ^2$ from Eq. (51) (solid line) and the contribution $\Delta ^2_i$ from Eq. (52) (dot-dashed line), as well as the usual one-loop perturbative result $\Delta ^2_{\rm 1 loop}$ of Eq. (81) (dotted line) and the prediction of the steepest-descent method ("sd'', dotted line) also shown in Fig. 8.

Open with DEXTER

Figure 19: The logarithmic power $\Delta ^2(k)$ at redshift z=0, that is at equal times z1=z2=0. We display the full non-linear power $\Delta ^2$ from Eq. (51) (solid line) and the contribution $\Delta ^2_i$ from Eq. (52) (dot-dashed line), as well as the linear power $\Delta ^2_{\rm L}$(dashed line), the usual one-loop perturbative result $\Delta ^2_{\rm 1 loop}$ of Eq. (81) (dotted line), the prediction of the steepest-descent method ("sd'', dotted line) also shown in Fig. 9 and the fit $\Delta ^2_{\rm nb}$ from numerical simulations (Smith et al. 2003). All quantities (including $\Delta ^2_i$) are positive.

Open with DEXTER

	Figure 20: The logarithmic power $\Delta ^2(k)$ at redshift z=3, that is at equal times z1=z2=3.
Open with DEXTER

Finally, we present in this section our results for the two-point correlation G and the self-energy $\Pi$. We first display in Fig. 18 the evolution with time $\eta _1$ of the density "logarithmic power'' $\Delta^2(k;\eta_1,\eta_2)$ defined in Eq. (87). As shown in the left panel at early times and large scales we recover the steepest-descent result since we must also match the usual one-loop expansion (81). At later times or at smaller scales we obtained oscillations which are strongly damped at large time separations in contrast with Fig. 8. This difference of behaviors of the correlations G between the direct steepest-descent method and the 2PI effective action scheme follows from the difference already seen in terms of the response R analyzed in Sect. 6.1. The right panel of Fig. 18 is particularly interesting as it shows that in the highly non-linear regime the two-point correlation exhibits a peak at equal times $\eta_1=\eta_2$ but whereas both the standard one-loop expansion and the steepest-descent method yield large oscillations at larger time separations the 2PI effective action approach leads to a strong damping. Thus, this is the only method (among those presented in this paper) which yields (at one-loop order) a damping of correlations at small scales and large time separations which is expected to reflect the qualitative behavior of the exact gravitational dynamics.

Next, we show in Figs. 19, 20 the logarithmic power $\Delta^2(k;\eta_1,\eta_2)$ at equal redshifts z₁=z₂=0 and z₁=z₂=3as a function of wavenumber k. We find again that our results match the steepest-descent prediction as well as the usual one-loop power (81) at large scales. At small scales, despite the damping at unequal times shown in Fig. 14 for the response function and in Fig. 18 for the correlation function we obtain a steady growth of the power $\Delta ^2(k)$in between the linear prediction $\Delta_{\rm L}^2$ and the usual one-loop prediction $\Delta ^2_{\rm 1 loop}$. Besides, the agreement with the results from numerical simulations is better than for the steepest-descent prediction. Of course, because of the damping of the response R at large time separations analyzed above we can check that the contribution $\Delta ^2_i$ from Eq. (52) associated with the transport of the initial fluctuations becomes negligible in the non-linear regime (contrary to what happens within the steepest-descent approach, see Figs. 9, 10). However, the continuous generation of fluctuations, described by the self-energy $\Pi$ in Eq. (51) and associated with the non-linearity of the dynamics, sustains a high level of density fluctuations $\Delta ^2(k)$into the non-linear regime. The latter seen at equal times $\eta_1=\eta_2=\eta$ are related to the values of $\Pi(k;\eta_1',\eta_2')$at nearly equal times $\eta_1'\simeq\eta_2'\simeq\eta$.

We show in Fig. 21 the logarithmic power $\Delta^2(k;\eta_1,\eta_2)$ at unequal times as a function of wavenumber k. We find again that our result matches the steepest-descent prediction as well as the usual one-loop power (81) at large scales. However, the correlation now quickly decays at small scales in the non-linear regime. This follows from the damping at large time separations analyzed above for the response function which leads to a decorrelation at small scales and large time separations.

	Figure 21: The logarithmic power $\Delta ^2(k)$ at unequal times $(\eta _1=0,\eta _2=-1.4)$ (i.e. z1=0,z2=3).
Open with DEXTER

	Figure 22: The self-energy $\Pi (k;\eta _1,\eta _2)$ as a function of forward time $\eta _1$, for $\eta _2=-1.4$ (i.e. z2=3) and wavenumbers k=1 and $10 \times h$ Mpc-1.
Open with DEXTER

	Figure 23: The self-energy $\Pi (k)$ as a function of wavenumber k at equal redshifts z1=z2=0. We also plot the prediction of the steepest-descent method (curves labeled "sd'').
Open with DEXTER

	Figure 24: The self-energy $\Pi (k)$ as a function of wavenumber k at unequal redshifts z1=0,z2=3. We also plot the prediction of the steepest-descent method (curves labeled "sd'').
Open with DEXTER

We display in Fig. 22 the self-energy $\Pi$ as a function of forward time $\eta _1$. We do not show the results of the steepest-descent method which follow a mere exponential growth from Eq. (65) At weakly non-linear scales ($k\sim 1~h$ Mpc^-1) we recover at early times the exponential growth of the steepest-descent prediction whereas at late times we obtain a fast decay. At highly non-linear scales (k=10 h Mpc^-1) we obtain a peak at equal times $\eta_1=\eta_2$ and strongly damped oscillations at large time separations. This is clearly consistent with the results obtained for the correlation G₁₁ displayed in Fig. 18.

Next, we show in Fig. 23 the self-energy $\Pi$ at equal redshifts z₁=z₂=0 as a function of wavenumber k. We again recover the steepest-descent prediction at large scales. However, at small scales we now obtain a steady growth whereas the steepest-descent method yields a smooth decline. Therefore, in agreement with Figs. 19, 20 we find that the coupled system of Eqs. (46) (with G₀ replaced by G) and (51) entails at equal times a stronger growth of the correlation G and of the self-energy $\Pi$ into the non-linear regime as compared with the non-coupled equations associated with the steepest-descent method.

Finally, we display in Fig. 24 the self-energy $\Pi$ at unequal redshifts z₁=0,z₂=3 as a function of wavenumber k. In agreement with our results for other two-point functions we recover the steepest-descent prediction at large scales and decaying oscillations at small scales.

   
7 Simplified response approximations

The results obtained in the previous sections suggest that a convenient approximation would be to use a simple analytical form for the response function R and to compute the correlation G and the self-energy $\Pi$from Eqs. (46), (51). This avoids the computation of the self-energy term ${\Sigma }$ and the numerical integration of the integro-differential Eq. (43) which requires rather small time-steps. By contrast, Eq. (51) is a Volterra integral equation of the second kind which is usually better behaved. We shall investigate two such approximations associated with the direct steepest-descent approach and the 2PI effective action scheme.

   
7.1 Simplified response for the steepest-descent scheme

Within the framework of the direct steepest-descent scheme studied in Sect. 5 we found that the response function exhibits in the non-linear regime a series of oscillations with an envelope given by the linear theory prediction. Therefore, following Eq. (80) we consider the approximate response $R_{\rm app}$ defined as:

$$R_{\rm app}(x_1,x_2)= R_{\rm L}(x_1,x_2) \cos~[\omega(k)(a_1-a_2)] ,$$ (109)

where $R_{\rm L}(=R_0)$ is the linear response given in Eq. (24), $a={\rm e}^{\eta}$is the scale factor and $\omega(k)$ is given by Eq. (77).

Here it is interesting to see how expression (109) may be interpreted. In the linear regime the density and velocity fluctuations $\psi_{\rm L}$ are merely amplified by a time-dependent factor (proportional to the scale-factor) which implies in particular that the physics is exactly local: $\psi_{\rm L}({\vec x},\eta)$ only depends on the initial conditions (and a possible external noise) at the same location ${\vec x}$. Thus in real space we have $R_{\rm L}({\vec x}_1,{\vec x}_2)\propto \delta_{\rm D}({\vec x}_1-{\vec x}_2)$ and in Fourier space, using the factorization (50) associated with statistical homogeneity, $R_{\rm L}(k)$ is constant with respect to k, in agreement with Eq. (24). On the other hand, the cosine dependence of Eq. (109) yields in real space:

$$R_{\rm app}(x_1,x_2) = R_{\rm L}(a_1,a_2) {\cal P}({\vec x}_1-{\vec x}_2;a_1,a_2)$$ (110)

where $R_{\rm L}(a_1,a_2)$ is the time-dependent factor of the linear response (24) and we introduced the distribution ${\cal P}({\vec r};a_1,a_2)$ obtained from Eq. (109) as:

$${\cal P}({\vec r};a_1,a_2) = \frac{-1}{4\pi r} \delta_{\rm D}'(r-d) \;\; \mbox{with} \;\; d= (a_1-a_2) r_0$$ (111)

where $r=\vert{\vec r}\vert$ and:

$$r_0= \left(\frac{4\pi}{3} \int {\rm d}q P_{\rm L0}(q)\right)^{1/2} .$$ (112)

Thus, Eqs. (109)-(112) may be interpreted as the displacement of the linear fluctuations over a random distance ${\vec r}$ with a probability distribution ${\cal P}(r)$ (whereas the linear regime obviously corresponds to the distribution $\delta_{\rm D}({\vec r})$). This non-linear operation corresponds to some diffusion of the initial conditions (over angles of the shift ${\vec r}$) but the singular distribution over the length $\vert{\vec r}\vert$ leads to the mere cosine tail in Fourier space instead of a strong decay. This can be contrasted to the Zeldovich approximation (Zeldovich 1970) where the trajectory of particles with initial comoving location ${\vec q}$ is written as ${\vec x}({\vec q},a)= {\vec q}+{\vec r}$ with $\nabla_{{\vec q}}.{\vec r}=-\delta_{\rm L}({\vec q},a)$. Thus the displacement field is Gaussian which leads to a Gaussian decay at high k(Crocce & Scoccimarro 2006a). However, we must emphasize that we only write Eq. (110) for comparison but this is not necessarily the best way to understand the result of the direct-steepest approach (at this stage expression (110) is only a mere re-writing of Eq. (109)). In particular, there is no reason why the non-linear response should be factorized as in Eq. (110) as $R_{\rm L}$ times a displacement term.

Figure 25: The logarithmic power $\Delta ^2(k)$ at redshift z=0. We display the ratio $\Delta ^2(k)/\Delta ^2_{\rm L}(k)$, where $\Delta ^2_{\rm L}$ is the linear power, for the usual one-loop perturbative result $\Delta ^2_{\rm 1 loop}$ (dotted line), a fit $\Delta ^2_{\rm nb}$ (upper solid line) from numerical simulations (Smith et al. 2003), the result of the steepest-descent method ("sd'', dotted line) shown in Fig. 9 and the approximate power $\Delta ^2_{\rm app}$ (solid line) obtained from the approximation (109).

Open with DEXTER

Next, we use the steepest-descent prescription (65)-(66) for the self-energy $\Pi$ and we compute numerically the correlation Gfrom the explicit expression (51). In fact, using Eq. (65) the integrals over time may be performed analytically which yields:

 

$$G_{\rm app}(k;\eta_1,\eta_2) = G_{\rm L}(k;\eta_1,\eta_2) \cos~(\omega a_1) \cos~(\omega a_2) + a_1^2 a_2^2 R(\omega a_1).\Pi_0(k).R(\omega a_2)^T$$ (113)

with:

$$R(\omega a)= f_+(\omega a) R_0^+ + f_-(\omega a) R_0^- .$$ (114)

The matrices $R_0^{\pm}$ are given in Eqs. (58) and we introduced the functions:

$$f_+(x)= \int_0^1 {\rm d}t \; \cos~[x(1-t)] = \frac{\sin~(x)}{x}$$ (115)

and:

 

$$f_-(x) = \int_0^1 {\rm d}t \; t^{5/2} \cos~[x(1-t)] = \frac{5}{2x... .../2}} \left[ \cos~(x) S(\!\sqrt{2x/\pi}) - \sin~(x) C(\!\sqrt{2x/\pi}) \right] ,$$ (116)

where C(x) and S(x) are the Fresnel integrals. At early times or large scales the first term dominates in Eq. (113) and we recover the linear regime $G\rightarrow G_{\rm L}$. At late times and small scales we must take into account both terms and for large $\omega a_i$ we obtain:

 

$$\omega a_i \gg 1\!\!: \; G(k;a_1,a_2) \simeq a_1 a_2 \biggl \lbra... ...{+T} \; \frac{\sin~(\omega a_1)\sin~(\omega a_2)}{\omega^2} \biggl \rbrace\cdot$$ (117)

Thus, we can see that at late times the non-linear correlation grows as $G(k;a_1,a_2) \sim a_1 a_2$ with oscillations for $a_1\neq a_2$. This agrees with the results presented in Sect. 5.3.

	Figure 26: The ratio $\Delta ^2(k)/\Delta ^2_{\rm L}(k)$ as in Fig. 25 but at redshift z=3.
Open with DEXTER

We show our results for the equal-time density correlation G₁₁ (i.e. the power $\Delta ^2(k)$) in Figs. 25, 26 at redshifts z=0 and z=3. We display the ratio to the linear power $\Delta ^2(k)/\Delta ^2_{\rm L}(k)$ to zoom over the weakly non-linear region and we also plot the exact one-loop steepest-descent prediction ("sd'') shown in Fig. 9 for comparison. Figures 25, 26 clearly show that the exact one-loop steepest-descent prediction matches the standard one-loop perturbative result at large scales. Moreover, Fig. 25 suggests that at z=0 the fit from numerical simulations (Smith et al. 2003) may be off by $10\%$ at $k \sim 0.3~h$ Mpc^-1. The power $\Delta ^2_{\rm app}$ follows roughly the behavior of the "exact'' one-loop steepest-descent prediction but it is smaller by up to $15\%$ and does not match very well the one-loop perturbative result. This failure at intermediate scales is not really surprising since Eq. (80) was derived in the large-k limit using the asymptotic behaviors (63)-(64). However, as can be seen in Fig. 1 they are only reached for k>10 h Mpc^-1. Therefore, the approximation (109) is not sufficient to obtain accurate results over weakly non-linear scales. This shows that the results obtained for the two-point correlation are rather sensitive to the approximations used for the response R. In fact, since the direct steepest-descent method is quite simple and easy to compute the approximation (109) is not very useful for precise quantitative predictions and it is best to use the rigorous method presented in Sect. 5.

   
7.2 Simplified response for the 2PI effective action scheme

Within the 2PI effective action method we found in Sect. 6 that the response function exhibits oscillations as for the steepest-descent prediction but their amplitude decays at large times. Following the analysis of the simple model (95) and Eq. (104) we consider the approximate response $R_{\rm app}$ defined as:

$$R_{\rm app}(x_1,x_2)= R_{\rm L}(x_1,x_2) \frac{J_1[2\omega(k)(a_1-a_2)]} {\omega(k)(a_1-a_2)}\cdot$$ (118)

Of course, this expression can again be written in the form (110) in real space, with a distribution ${\cal P}(r)$ which is now given by:

$${\cal P}(\vert{\vec r}\vert;a_1,a_2) = \frac{\theta(r<2d)}{8\pi^2 d^3 \sqrt{1-(r/2d)^2}} ,$$ (119)

where the top-hat factor $\theta(r<2d)$ obeys $\theta=1$ for r<2d and $\theta=0$ for r>2d. The decay of the response function at high k corresponds to the smoother distribution (119) as compared with Eq. (111), but the discontinuity at r=2d leads to a slow power-law damping instead of an exponential decay. Again, Eq. (119) is a mere rewriting of Eq. (118) and we do not claim here that the distribution ${\cal P}(r)$ is anything more than a phenomenological tool.

Then, we use Eq. (46) (with G₀ replaced by G) and Eq. (51) to compute the self-energy $\Pi$ and the correlation G, following the procedure described in Sect. 6.

Figure 27: The logarithmic power $\Delta ^2(k)$ at redshift z=0. We display the ratio $\Delta ^2(k)/\Delta ^2_{\rm L}(k)$, where $\Delta ^2_{\rm L}$ is the linear power, for the usual one-loop perturbative result $\Delta ^2_{\rm 1 loop}$ (dotted line), a fit $\Delta ^2_{\rm nb}$ (upper solid line) from numerical simulations (Smith et al. 2003), the result of the 2PI effective action method ("2PI'', dotted line) shown in Fig. 19 and the approximate power $\Delta ^2_{\rm app}$ (solid line) obtained from the approximation (118).

Open with DEXTER

	Figure 28: The ratio $\Delta ^2(k)/\Delta ^2_{\rm L}(k)$ as in Fig. 27 but at redshift z=3.
Open with DEXTER

We show our results for the power $\Delta ^2(k)$) in Figs. 27, 28 at redshifts z=0 and z=3. We again display the ratio to the linear power $\Delta ^2(k)/\Delta ^2_{\rm L}(k)$ and we also plot the exact 2PI effective action prediction ("2PI'') shown in Fig. 19 for comparison. Figures 27, 28 again clearly show that the exact 2PI effective action prediction matches the standard one-loop perturbative result at large scales. In a fashion similar to the approximation (109) we find that the power $\Delta ^2_{\rm app}$ follows roughly the behavior of the exact one-loop 2PI effective action prediction but it is smaller by up to $15\%$ and does not match very well the one-loop perturbative result. This shows again that the two-point correlation is rather sensitive to the value of the response function.

   
7.3 Mixed approaches

Finally, we consider in this section mixed approaches which use elements from both the steepest-descent and 2PI effective action schemes. The idea is to simplify the computation involved in the 2PI effective action method by separating the computations of the pairs $(\Sigma,R)$ and $(\Pi,G)$. That is, we first compute the self-energy ${\Sigma }$and the response R from the coupled Eqs. (43), (45), where we use G₀ and R in the expression (45) for ${\Sigma }$(the steepest-descent scheme uses G₀ and R₀ whereas the 2PI effective action scheme uses G and R). Thus, the pair $(\Sigma,R)$ becomes independent of $(\Pi,G)$ and can be computed in a first step. Since we keep the explicit dependence on R in the self-energy ${\Sigma }$the evolution Eq. (43) for the response function is no longer linear but quadratic, hence we expect to recover the damping discussed in Sect. 6.1 (see the analysis of Eq. (95)). We display our results in Fig. 29 which shows that the non-linear response exhibits indeed a damping close to the one obtained in Fig. 14 for the full 2PI effective action scheme.

	Figure 29: The non-linear response function $R(k;\eta _1,\eta _2)$ as a function of wavenumber k, at times $\eta _1=0,\eta _2=-1.4$, for the mixed approaches where G0 is used in the computation of ${\Sigma }$.
Open with DEXTER

Next, we compute the pair $(\Pi,G)$ from the Eqs. (46) and (51). Here we can consider two procedures as we may use either G₀ or G in the computation of $\Pi$. Clearly, these mixed schemes still agree with the usual one-loop perturbative results. These two-step approaches are slightly faster than the full 2PI effective action method for numerical computations. Indeed, as explained in Sect. 6, within the 2PI effective action approach the equations obtained for two-point functions and self-energies form a coupled system of non-linear equations. This leads to an iterative scheme for their numerical computation (which converges within $\sim$7 steps at worst for the scales we consider here). On the other hand, at each step the computation of the correlation G takes most of the computer time (especially at late times) because of the double integral over times in Eq. (51) but it converges faster (at fixed R) than the response R (at fixed G). Therefore, it saves time to separate the computations of R (many fast iteration steps) and G (few long iteration steps).

Figure 30: The logarithmic power $\Delta ^2(k)$ at redshift z=0, that is at equal times z1=z2=0. We display the linear power $\Delta ^2_{\rm L}$ (dashed line), the usual one-loop perturbative result $\Delta ^2_{\rm 1 loop}$ (dotted line), a fit $\Delta ^2_{\rm nb}$ (lower solid line) from numerical simulations (Smith et al. 2003), the mixed approach $\Delta ^2_{\Pi _0}$ based on $\Pi _0$ (dashed line) and the mixed approach $\Delta ^2_{\Pi }$ where $\Pi$ is coupled to G (upper solid line).

Open with DEXTER

	Figure 31: The logarithmic power $\Delta ^2(k)$ at redshift z=3, that is at equal times z1=z2=3.
Open with DEXTER

Figure 32: The ratio $\Delta ^2(k)/\Delta ^2_{\rm L}(k)$ at redshift z=0 of the logarithmic power $\Delta ^2$ to the linear power $\Delta ^2_{\rm L}$. We display the ratio $\Delta ^2/\Delta ^2_{\rm L}$ for the usual one-loop perturbative result $\Delta ^2_{\rm 1 loop}$ (dotted line), a fit $\Delta ^2_{\rm nb}$ (lower solid line) from numerical simulations (Smith et al. 2003), the mixed approach $\Delta ^2_{\Pi _0}$ based on $\Pi _0$ (dashed line) and the mixed approach $\Delta ^2_{\Pi }$ where $\Pi$ is coupled to G (upper solid line).

Open with DEXTER

	Figure 33: The ratio $\Delta ^2(k)/\Delta ^2_{\rm L}(k)$ as in Fig. 32 but at redshift z=3.
Open with DEXTER

We display in Figs. 30-33 our results for the power $\Delta ^2(k)$ at redshifts z=0 and z=3 for both schemes. The dashed curve $\Delta ^2_{\Pi _0}$ shows the power obtained when the self-energy term $\Pi _0$ is used (as in the steepest-descent approach) so that the correlation G is directly given by the explicit expression (51). The solid curve  $\Delta ^2_{\Pi }$ corresponds to the case where we use the non-linear correlation G in the expression (46) of the self-energy, so that Eqs. (46) and (51) are a coupled non-linear system which we solve iteratively (as in the 2PI effective action scheme). We can see that both methods exhibit a similar behavior until $\Delta^2(k)/\Delta^2_{\rm L}(k)\sim 1.4$. Farther into the non-linear regime the "$\Pi _0$-scheme'' breaks from the non-linear growth of the power-spectrum and shows oscillations which converge close to the linear power, in a manner reminiscent of the steepest-descent approximation, whereas the coupled "$\Pi$-scheme'' follows this non-linear growth. On the other hand, in the weakly non-linear regime displayed in Figs. 32-33 the predicted power $\Delta ^2(k)$is suppressed as compared with the usual one-loop result. Depending on the redshift this can improve or worsen the global agreement with numerical simulations. However, Fig. 33 suggests that when the usual one-loop result shows a better global agreement with simulations it is a mere coincidence as the simulation result separates last from the large-N predictions (and the closer agreement with the usual one-loop result at smaller scales has no theoretical basis). Nevertheless, the limited accuracy of the numerical simulations prevents one from drawing definite conclusions. We can note in Fig. 30 that in the highly non-linear regime the approximation  $\Delta ^2_{\Pi }$ grows larger than for the full 2PI effective action prediction displayed in Fig. 19. Indeed, in the latter case the non-linear response is further damped by the strong growth of the two-point correlation G (this would correspond to a larger $\sigma$ whence $\omega$ in the simple model (104)) and this "negative feedback'' within the coupled system R,G leads to a smaller G as compared with the approximation  $\Delta ^2_{\Pi }$ where this feedback is neglected. In fact, at very high k (beyond the range shown in Fig. 30) the power $\Delta ^2_{\Pi }$ appears to diverge. It is not clear whether this is due to the finite resolution of the numerical computation or to a true divergence. However, for practical purposes this is not a serious shortcoming since it occurs at small scales beyond the range of validity of this approach (one would need to go beyond one-loop order and probably beyond the hydrodynamical approach as shell-crossing can no longer be neglected).

The good agreement in the weakly non-linear regime of these various schemes enables one to evaluate their range of validity without the need to compare their prediction with N-body simulations, which could be of great practical interest. In fact, Fig. 32 suggests that this procedure can provide a better accuracy than numerical simulations in this regime, where the latter may still be off by $10\%$.

   
8 Application to a $\Lambda$CDM universe

We finally describe in this section how the results obtained in the previous sections can be applied to more general cosmologies, such as a $\Lambda$CDM universe. As seen in Sect. 2.1 our formalism applies equally well to any cosmology. We only need to use the relevant matrix ${\cal O}$, that is to take into account the time dependence of the ratio  $\Omega_{\rm m}/f^2$ in Eq. (13), and to recall that $\eta$ is the logarithm of the linear growth factor, see Eq. (9). In practice, a widely used approximation is to take $\Omega_{\rm m}/f^2 \simeq 1$ so that the results obtained for the critical-density universe directly apply to the case  $\Lambda\neq 0$ and the only change is the time-redshift relation  $\eta \leftrightarrow z$. We shall use this simple approximation here and compare our results to numerical simulations. Thus, we consider the cosmological parameters $\Omega_{\rm m}=0.3, \Omega_{\Lambda}=0.7, \sigma_8=0.9$ and H₀=70 km s^-1 Mpc^-1as in Smith et al. (2003).

	Figure 34: The non-linear response function $R(k;\eta _1,\eta _2)$ given by the direct steepest-descent method as a function of wavenumber k for a $\Lambda$CDM universe, at redshifts z1=0,z2=3. The horizontal lines are the linear response $R_{\rm L}$ of Eq. (24) which is independent of k.
Open with DEXTER

	Figure 35: The non-linear response function $R(k;\eta _1,\eta _2)$ given by the 2PI effective action method as a function of wavenumber k for a $\Lambda$CDM universe, at redshifts z1=0,z2=3. We also plot the first half-oscillation of the response obtained from the direct steepest-descent method (curves labeled "sd'') which was shown in Fig. 34.
Open with DEXTER

We first plot in Figs. 34, 35 the non-linear response functions we obtain within the one-loop steepest-descent and 2PI effective action methods. We can check that we recover the behavior obtained for the critical-density universe displayed in Figs. 7, 14. The steepest-descent approach yields again oscillations in the non-linear regime with an amplitude given by the linear response at high k whereas the 2PI effective action method gives damped oscillations. The transition to the non-linear regime occurs at slightly smaller wavenumber because of the larger normalization $\sigma_8$ of the linear power-spectrum.

We show in Figs. 36-39 our results for the two-point correlation G₁₁, that is the non-linear density power spectrum. We again recover the behavior analyzed in previous sections for the critical-density universe. Both large-N expansion schemes agree with the standard one-loop expansion at large scales. At small scales the steepest-descent result goes back to values close to the linear power $\Delta ^2_{\rm L}$ whereas the 2PI effective action result keeps growing in a manner similar to the usual perturbative result and numerical simulations. The degree of agreement between the N-body simulations and either one of the standard and 2PI results depends on redshift, that is on the slope of the linear power-spectrum at the relevant scales.

Thus, the large-N expansions developed in this article apply in the same manner for different cosmologies and exhibit the same properties.

Figure 36: The logarithmic power $\Delta ^2(k)$ at redshift z=0 for a $\Lambda$CDM universe. We display the linear power $\Delta ^2_{\rm L}$ (dashed line), the usual one-loop perturbative result $\Delta ^2_{\rm 1 loop}$ (dotted line), a fit $\Delta ^2_{\rm nb}$ (lower solid line) from numerical simulations (Smith et al. 2003), the steepest-descent prediction (dotted line "sd'') and the 2PI effective action result (solid line "2PI'').

Open with DEXTER

	Figure 37: The logarithmic power $\Delta ^2(k)$ at redshift z=3 for a $\Lambda$CDM universe.
Open with DEXTER

Figure 38: The ratio $\Delta ^2(k)/\Delta ^2_{\rm L}(k)$ at redshift z=0 of the logarithmic power $\Delta ^2$ to the linear power $\Delta ^2_{\rm L}$ for a $\Lambda$CDM universe. We display the ratio $\Delta ^2/\Delta ^2_{\rm L}$ for the usual one-loop perturbative result $\Delta ^2_{\rm 1 loop}$ (dotted line), a fit $\Delta ^2_{\rm nb}$ (lower solid line) from numerical simulations (Smith et al. 2003), the steepest-descent prediction (dashed line "sd'') and the 2PI effective action result (solid line "2PI'').

Open with DEXTER

	Figure 39: The ratio $\Delta ^2(k)/\Delta ^2_{\rm L}(k)$ as in Fig. 38 but at redshift z=3.
Open with DEXTER

   
9 Conclusion

In this article we have studied the predictions at one-loop order of large-Nexpansions for the two-point correlations associated with the formation of large-scale structures in the expanding universe. Focussing on weakly non-linear scales we have used the hydrodynamical equations of motion as a starting point. Then, we have recalled the path-integral formalism which describes the statistical properties of the system, assuming Gaussian initial conditions, and we have presented two possible large-N expansions. Next, we have described in details the numerical predictions obtained for a SCDM cosmology.

First, we have presented a direct steepest-descent approach where the self-energy is expressed in terms of the linear response and correlation. This allows simple computations as the time-dependence of the self-energy can be factorized (for a critical density universe). Moreover, the non-linear response R can be computed directly from this self-energy as it is not coupled to the non-linear correlation G and it obeys a simple linear equation. This also allows a detailed analysis of its behavior. We have found that at one-loop order this approach yields a non-linear response which exhibits strong oscillations in the non-linear regime (small scales or late times) with an amplitude which is given by the linear response. Therefore, the response does not decay but only exhibits increasingly fast oscillations of the form $R_{\rm L} \cos~[(a_1-a_2)k/k_*]$. However, once the response is integrated over time such oscillations lead to some damping as compared with the linear response. Then, we have computed the non-linear correlation which can be explicitly written in terms of this non-linear response which transports forward over time the initial density and velocity fluctuations as well as the fluctuations generated at all previous times by non-linear couplings. We have checked that in the quasi-linear regime the non-linear correlation agrees with the usual one-loop prediction (obtained from the standard perturbative expansion over powers of the initial fluctuations). In the highly non-linear regime the correlation separates from the steep growth displayed by the standard one-loop result and converges back to the linear amplitude. Unfortunately, this behavior does not improve the agreement with N-body numerical simulations. Nevertheless, it suggests that this large-N expansion is well-behaved as the two-point functions do not "explode'' in the highly non-linear regime. Contrary to the standard perturbative expansion, where higher-order terms grow as increasingly large powers of k and a(t), the higher-order terms obtained within this method may remain of finite amplitude.

Secondly, we have described a 2PI effective action approach where all two-point functions are coupled. The numerical computation requires an iterative procedure but thanks to causality the equations can be solved by integrating forward over time with the help of a few iterations at each time-step (keeping an implicit scheme in the non-linear integro-differential equations). Then, we have found that the one-loop non-linear response again exhibits fast oscillations in the non-linear regime but its amplitude now shows a fast decay. We have explained how this damping is produced by the non-linearity of the evolution equation obeyed by the response function. This yields a power-law decay such as $R_{\rm L} J_1[2(a_1-a_2)k/k_*]/[(a_1-a_2)k/k_*]$. In the quasi-linear regime the non-linear correlation agrees again with the usual one-loop prediction but in the highly non-linear regime it no longer converges back to the linear power. It rather shows a moderate growth intermediate between the linear prediction and the result of numerical simulations, with an amplitude which depends on the redshift (whence on the slope of the linear power-spectrum). Therefore, this method appears to be superior to the direct steepest-descent approach which could not recover any non-linear damping (at one-loop order) but it requires more complex numerical computations.

Thirdly, we have described other approximation schemes built from these two approaches. Thus, we have found that using a simple explicit approximation for the non-linear response within either large-N scheme gives correct qualitative results but does not manage to reproduce accurately the quantitative results obtained from the full large-N approaches. Next, we have investigated the approximations obtained by separating the computations of the response and the correlation (but keeping a non-linear dynamics for the response). As expected we have found that this yields a non-linear response which closely follows the damping obtained by the full 2PI effective action approach and it gives good results (as compared with the full large-N schemes) for the correlation.

Finally, we have shown that these results also apply to more general cosmologies such as a $\Lambda$CDM universe, since our formalism extends to any expansion history of the Universe provided we use the correct time-redshift relation and linear growth factors.

Thus, we have shown in this article that large-N expansions provide an interesting systematic approach to the formation of large-scale structures in the expanding universe. They already show at one-loop order some damping in the non-linear regime for the response function and the amplitude of the two-point correlation remains well-behaved. Unfortunately, for practical purposes the accuracy provided by these large-N expansions at one-loop order is not significantly better than the usual perturbative result for weakly non-linear scales. Nevertheless, the improvement may become greater at higher orders since the expansion is likely to be better behaved, but this requires rather complex computations which are beyond the scope of this paper.

On the other hand, one of the goals of this article was to describe a different approach to the problem of non-linear gravitational clustering than the standard expansion schemes, namely a path-integral formalism. An important feature of this approach is that the statistical nature of the problem is already included in the definition of the system, that is the action $S[\psi,\lambda]$ describes both the equations of motion (here in the hydrodynamical limit) and the average over Gaussian initial conditions. This also means that one directly works with correlation functions (such as the power spectrum). By contrast, in usual expansion schemes (or in N-body numerical simulations) one first works with the density field associated with a specific initial condition, that is one tries to express the non-linear density field in terms of the linear density field (e.g. as a truncated power series) and next performs the average over the initial condition. In principles, this formalism can present several advantages. First, the statistical quantities such as the power-spectrum are precisely the objects of interest for practical purposes. Secondly, they exhibit symmetries (such as homogeneity and isotropy) which are not satisfied by peculiar realizations of the initial conditions. Moreover, they are smooth functions (such as power laws) rather than singular distributions. This could facilitate the numerical computations and the building of various approximations.

Of course, one can recover the standard perturbative results by expanding over the cubic part of the action, as recalled in this paper and discussed in more details in Valageas (2004). However, the advantage of this formalism is that it may serve as a starting point to other approximation schemes by applying the methods of statistical mechanics and field theory. Thus, we have described in details in this paper two large-N expansion schemes. In addition to the interesting properties of these two methods discussed above, we can note that it is useful to have several rigorous expansion schemes which only differ by higher-order terms than the truncation order. Indeed, since one can expect that they should be correct up to the scale where they start to depart from one another, this should allow one to estimate their range of validity without computing explicitly higher-order terms or performing N-body simulations. In fact, in the weakly non-linear regime where all expansion schemes agree this procedure should provide a better accuracy than N-body simulations which may still be inaccurate by up to $10\%$. Farther into the non-linear regime, it seems that the usual one-loop result and the 2PI effective action scheme give better results than the direct steepest-descent method. On the other hand, one can hope that other methods than these two expansion schemes could be applied to the action $S[\psi,\lambda]$.

In regard to recent works we can note that the Schwinger-Dyson equations obtained in Crocce & Scoccimarro (2006a) can also be obtained from the approach described in this paper, provided we use the equations of motion in their integral form (that is we first integrate the time-derivative), as also described in Valageas (2001). Whereas Crocce & Scoccimarro (2006a) used a diagrammatic technique and the resummation of all diagrams gives back the usual Schwinger-Dyson equations the path-integral formalism directly gives these equations by expanding the action (or its Legendre transform) about a saddle-point. On the other hand, the diagrammatic approach also provides a direct expression of the response function $R(\eta_1,\eta_2)$ (in the limit $\eta_2\rightarrow-\infty$) as a series of diagrams (which could be recovered here by writing the self-energy ${\Sigma }$ as a series of diagrams). Thus, Crocce & Scoccimarro (2006a) managed to resum these diagrams in the high-k limit and to obtain the asymptotic behavior of $R(\eta_1,\eta_2\rightarrow-\infty)$. It is not clear yet how this result could be directly obtained from the path-integral formalism without expanding over diagrams.

Among the points which deserve further studies, it would be interesting to apply these approaches to higher-order correlations beyond the two-point functions investigated here. On the other hand, the analysis presented in this article is based on the hydrodynamical equations of motion which break down at small scales beyond shell-crossing. It is possible to apply these large-N expansions to the Vlasov equation (Valageas 2004) but the numerical computations would be significantly more complex (since one needs to add velocities to space coordinates). Nevertheless, the results obtained within the hydrodynamical framework (damping and well-behaved non-linear asymptotics) can be expected to remain valid in the collisionless case and studies such as the present one may serve as a first step towards an application to the Vlasov dynamics. Such large-N expansions could also be applied to simpler effective dynamics which attempt to go beyond shell-crossing (based for instance on a Schroedinger equation, Widrow & Kaiser 1993). These works would be of great practical interest as several cosmological probes (such as weak lensing surveys and measures of the baryon acoustic oscillations) which aim to measure the recent expansion history of the Universe in order to constrain the dark energy equation of state (as well as other cosmological parameters) are very sensitive to the weakly non-linear scales where non-linear effects can no longer be neglected.

References

• Bernardeau, F. 1994, A&A, 291, 697 [NASA ADS] (In the text)
• Bernardeau, F., van Waerbeke, L., & Mellier, Y. 1997, A&A, 322, 1 [NASA ADS] (In the text)
• Bernardeau, F., Colombi, S., Gaztanaga, E., & Scoccimarro, R. 2002, Phys. Rept., 367, 1 [NASA ADS] [CrossRef] (In the text)
• Crocce, M., & Scoccimarro, R. 2006a, Phys. Rev. D, 73, 063519 [NASA ADS] [CrossRef] (In the text)
• Crocce, M., & Scoccimarro, R. 2006b, Phys. Rev. D, 73, 063520 [NASA ADS] [CrossRef] (In the text)
• Fry, J. N. 1984, ApJ, 279, 499 [NASA ADS] [CrossRef] (In the text)
• Goroff, M. H., Grinstein, B., Rey, S.-J., & Wise, M. B. 1986, ApJ, 311, 6 [NASA ADS] [CrossRef] (In the text)
• Jain, B., & Bertschinger, E. 1994, ApJ, 431, 495 [NASA ADS] [CrossRef] (In the text)
• Jain, B., & Bertschinger, E. 1996, ApJ, 456, 43 [NASA ADS] [CrossRef] (In the text)
• Koehler, R. S., Schuecker, P., & Gebhardt, K. 2007 [arXiv:astro-ph/0609354] (In the text)
• Liddle, A. R., & Lyth, D. H. 1993, Phys. Rep., 231, 1 [NASA ADS] [CrossRef] (In the text)
• Makino, N., Sasaki, M., & Suto, Y. 1992, Phys. Rev. D, 46, 585 [NASA ADS] [CrossRef] (In the text)
• McDonald, P. 2006, Phys. Rev. D, 75, 043514 [NASA ADS] [CrossRef] (In the text)
• Martin, P. C., Siggia, E. D., & Rose, H. A. 1973, Phys. Rev. A, 8, 423 [NASA ADS] [CrossRef] (In the text)
• Munshi, D., Valageas, P., Van Waerbeke, L., & Heavens, A. 2007 [arXiv:astro-ph/0612667] (In the text)
• Peebles, P. J. E. 1980, The large scale structure of the universe (Princeton University Press) (In the text)
• Peebles, P. J. E. 1982, ApJ, 263, L1 [NASA ADS] [CrossRef] (In the text)
• Phythian, R. 1977, J. Phys. A, 10, 777 [NASA ADS] [CrossRef] (In the text)
• Scoccimarro, R., & Frieman, J. 1996a, ApJS, 105, 37 [NASA ADS] [CrossRef] (In the text)
• Scoccimarro, R., & Frieman, J. 1996b, ApJ, 473, 620 [NASA ADS] [CrossRef] (In the text)
• Seljak, U. 2000, MNRAS, 318, 203 [NASA ADS] (In the text)
• Seljak, U., & Zaldarriaga, M. 1996, ApJ, 469, 437 [NASA ADS] [CrossRef] (In the text)
• Seo, H.-J., & Eisenstein, D. J. 2007 [arXiv:astro-ph/0701079] (In the text)
• Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311 [NASA ADS] (In the text)
• Valageas, P. 2001, A&A, 379, 8 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Valageas, P. 2002, A&A, 382, 477 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Valageas, P. 2004, A&A, 421, 23 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Vishniac, E. T. 1983, MNRAS, 203, 345 [NASA ADS] (In the text)
• Widrow, L. M., & Kaiser, N. 1993, ApJ, 416, L71 [NASA ADS] [CrossRef] (In the text)
• Zeldovich, Y. B. 1970, A&A, 5, 84 [NASA ADS] (In the text)
• Zinn-Justin, J. 1989, Quantum field theory and critical phenomena (Oxford: Clarendon Press) (In the text)