A&A 379, 8-20 (2001)
DOI: 10.1051/0004-6361:20011309

Dynamics of gravitational clustering

I. Building perturbative expansions

P. Valageas

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

Received 28 May 2001 / Accepted 13 September 2001

Abstract
We develop a systematic method to obtain the solution of the collisionless Boltzmann equation which describes the growth of large-scale structures as a perturbative series over the initial density perturbations. We give an explicit calculation of the second-order terms which are shown to agree with the results obtained from the hydrodynamical description of the system. Then, we explain that this identity extends to all orders of perturbation theory and that the perturbative series actually diverge for hierarchical scenarios. However, since the collisionless Boltzmann equation provides the exact description of the dynamics (including the non-linear regime) these results may serve as a basis for a study of the non-linear regime. In particular, we derive a non-perturbative quadratic integral equation which explicitly relates the actual non-linear distribution function to the initial conditions (more precisely, to the linear growing mode). This allows us to write an explicit path-integral expression for the probability distribution of the exact non-linear density field.

Key words: cosmology: theory - large-scale structure of universe

1 Introduction

In usual cosmological models large-scale structures in the universe have formed through the growth of small initial density perturbations by gravitational instability, see Peebles (1980). Moreover, in most cases of cosmological interest the power increases at small scales as in the CDM model (Peebles 1982). This leads to a hierarchical scenario of structure formation where smaller scales become non-linear first. Once they enter the non-linear regime high-density fluctuations form massive objects which will give birth to galaxies or clusters. Thus, it is important to understand the dynamics of these density perturbations in order to describe the large-scale structures we observe in the present universe.

The distribution of matter is usually described as a pressure-less fluid through the standard equations of hydrodynamics (continuity and Euler equation) coupled to the Poisson equation for the gravitational potential. Then, one often looks for a solution of the equations of motion under the form of a perturbative expansion (e.g., Fry 1984; Goroff et al. 1986; Bernardeau 1994). In fact, one cannot hope to describe the non-linear regime within such a framework since this hydrodynamical description becomes inexact as soon as shell-crossing appears. Then, there is no longer only one velocity at each point. Note that for hierarchical scenarios with no small-scale cutoff (i.e. the amplitude of density fluctuations diverges at small scales $R \rightarrow 0$) shell-crossing occurs as soon as t>0. Moreover, in the non-linear regime shell-crossing should play a key role (e.g., through virialization processes).

Therefore, in order to obtain a complete description of the dynamics one needs to go beyond the standard approximation of a pressure-less fluid. Then, one might add a pressure term to the hydrodynamical equations in order to describe the effects of multi-streaming, which provide some velocity dispersion (e.g., Buchert & Dominguez 1998; Buchert et al. 1999). However, this is not very satisfactory because collisionless and gaseous systems can exhibit very different behaviours starting from the same initial conditions. For instance, some physical processes are specific to collisionless systems (e.g., Landau damping, phase-mixing) and stability conditions for similar collisionless and gaseous systems are not always the same, see Binney & Tremaine (1987). Hence it is important to investigate the dynamics of density fluctuations within the framework of the collisionless Boltzmann equation, which provides an exact description of the system both in the linear and non-linear regimes.

Since this is a rather difficult task we may first try to devise a perturbative approach. This may then serve as a basis for further developments to tackle the non-linear regime. Thus, in this article we build a perturbative approach to the collisionless Boltzmann equation. One could also start from the BBGKY hierarchy, which can be derived from the Boltzmann equation. Such a study was performed in Bharadwaj (1996) where it was shown that the lowest non-linear order matches the results of the usual hydrodynamical approach. We recover this result here and we explain that indeed a perturbative method cannot handle shell-crossing. However, in the building of perturbative theory we also derive in this article some non-perturbative results which may prove to be useful for a study of the non-linear regime.

This paper is organized as follows. In Sect. 2 we write the equations of motion which describe the dynamics of the collisionless fluid. Then, in Sect. 3 we develop a systematic procedure to obtain the fluctuations of the distribution function as a perturbative series over the initial perturbations (or more precisely over the linear growing mode) in the case of a critical density universe. In doing so we also derive a non-perturbative quadratic integral equation which directly relates the actual non-linear distribution function to this linear mode. Next, we give in Sect. 4 the explicit calculations of the second-order terms. We explain in Sect. 5 why we recover the results of the hydrodynamical approach and why the perturbative series must diverge for hierarchical scenarios (with no small-scale cutoff). Then, we show in Sect. 6 how to extend this method to arbitrary cosmologies. Finally, we derive in Sect. 7 an explicit path-integral expression which describes the statistical properties of the exact non-linear density field.

   
2 Equations of motion

The gravitational dynamics of a fluid of collisionless particles of mass m is given by the collisionless Boltzmann equation (Peebles 1980), coupled with the Poisson equation:

$$\left\{ \begin{array}{l} {\displaystyle \frac{\partial f}{\p... ...int f({\bf x},{\bf p},t) \; {\rm d}{\bf p}} \end{array}\right.$$ (1)

where $f({\bf x},{\bf p},t)$ is the distribution function (phase-space density), $\phi$ is the gravitational potential, $\rho({\bf x},t)$ is the comoving density, $\overline{\rho}$ is the mean comoving density of the universe (hence it is independent of time), a(t) is the scale factor and the momentum ${\bf p}$ is related to the comoving coordinate ${\bf x}$ of a particle by ${\bf p}= m a^2 \dot{{\bf x}}$. Thus, Eq. (1) describes the growth of gravitational structures in an expanding universe which is homogeneous on large scales.

Note that in order to derive the collisionless Boltzmann equation one needs to neglect the two-particle correlations. Then, the gravitational interaction is described by a smooth gravitational potential which neglects the effects of the discrete character of the distribution of matter (e.g., Peebles 1980). This becomes exact in the continuous limit $m \rightarrow 0$ which we perform below in Eq. (2). However, one must also add some small-scale cutoff to the power-spectrum of the density fluctuations in order to guarantee small-scale smoothness.

The definition itself of the problem we investigate here also involves an implicit "regularization'' of the gravitational interaction. Indeed, within an infinite uniform universe the gravitational force is not well defined since the pull due to the matter located within a small solid angle ${\rm d}\Omega$ diverges as we integrate over all shells up to infinity. This problem is cured by first integrating the relevant integrals over angles, and next over radial distances (e.g., Peebles 1980). This can also be seen as the introduction of a large-scale cutoff L for the gravitational interaction which preserves rotational symmetry. More precisely, this means that in Fourier space we integrate on the wavenumber ${\bf k}$ over $\vert{\bf k}\vert>k_{\rm c}$ and in the final results we take the limit $k_{\rm c} \rightarrow 0$. In fact, except in Eq. (76) in Sect. 7 this feature will not show up in the equations we encounter in this article so that we can directly take $k_{\rm c}=0$.

Next, we can absorb the mass m into the distribution function f and the momentum ${\bf p}$: $f \rightarrow f/m^4$, ${\bf p}\rightarrow m {\bf p}$, and we obtain:

$$\left\{ \begin{array}{l} {\displaystyle \frac{\partial f}{\p... ...int f({\bf x},{\bf p},t) \; {\rm d}{\bf p}} \end{array}\right.$$ (2)

which is also valid in the continuous limit $m \rightarrow 0$. Here $f({\bf x},{\bf p},t) {\rm d}^3x {\rm d}^3p$ is the mass enclosed in the phase-space element ${\rm d}^3x {\rm d}^3p$ and the momentum ${\bf p}$ verifies:

$${\bf p}= a^2 \dot{{\bf x}} .$$ (3)

In order to single out the deviations of the distribution function from the case of a perfectly homogeneous universe we define the perturbation $\delta f({\bf x},{\bf p},t)$ by:

$$f({\bf x},{\bf p},t) \equiv \overline{\rho}\; \delta_{\rm D}({\bf p}) + \overline{\rho}\; \delta f({\bf x},{\bf p},t)$$ (4)

where $\delta_{\rm D}$ is Dirac's function (here in 3 dimensions). Then, the density contrast $\delta({\bf x},t)$ is simply given by:

$$\delta({\bf x},t) = \int \delta f({\bf x},{\bf p},t) \; {\rm d}{\bf p}.$$ (5)

From the definition (4) and the equation of motion (2) we obtain the evolution equation for $\delta f$:

$$\left\{ \begin{array}{l} {\displaystyle \frac{\partial\delta... ...}{a} \overline{\rho}\; \delta({\bf x},t) } .\end{array}\right.$$ (6)

Note the term $\delta_{\rm D}'({\bf p})$ which is sometimes forgotten in the literature. Next, we define the spatial Fourier transform of the field $\delta f$ by:

$$\left\{ \begin{array}{l} {\displaystyle \delta f({\bf x},{\bf... ....{\bf x}} \; \delta f({\bf x},{\bf p},t) }. \end{array}\right.$$ (7)

To simplify the notations, we note the distribution function by $\delta f$, both in real space ${\bf x}$ and in Fourier space ${\bf k}$. Its argument removes any ambiguity. We also define the Fourier transforms $\delta({\bf k})$ and $\phi({\bf k})$ of the density contrast and of the gravitational potential as in Eq. (7). Then, Poisson's equation reads:

$$-k^2 \phi({\bf k}) = \frac{4\pi {\cal G}}{a} \overline{\rho}\; \delta({\bf k}) .$$ (8)

After substitution of Eq. (8) into the equation of motion (6) we obtain:

 

$$\frac{\partial\delta f}{\partial t} + i \frac{{\bf k}.{\bf p}}{a^... ... \cdot \frac{\partial\delta f}{\partial{\bf p}}({\bf k}-{\bf k}',{\bf p}) = 0 .$$ (9)

This non-linear and non-local equation describes the dynamics of the collisionless fluid.

   
3 Building a perturbative expansion

   
3.1 General framework

Here we consider the case of a critical density universe $\Omega_{\rm m}=1$. Then, the average comoving density is given by:

$$\overline{\rho}= \frac{a^3}{6 \pi {\cal G}t^2} \cdot$$ (10)

We normalize the time scale such that the Hubble time today t₀ (i.e. at z=0) satisfies t₀=1, while the scale factor is a=1. Then, we define a new time coordinate $\tau$ by:

$$\tau \equiv \frac{1}{3}\ln (t/t_0) = \frac{1}{3}\ln t , \hspace{0.2cm} t={\rm e}^{3\tau} , \hspace{0.2cm} a = {\rm e}^{2\tau}$$ (11)

and Eq. (9) reads:

 

$${ \frac{\partial\delta f}{\partial\tau} + 3 i ({\bf k}.{\bf p}) {... ... \cdot \frac{\partial\delta f}{\partial{\bf p}}({\bf k}-{\bf k}',{\bf p}) = 0 .$$ (12)

This equation is quadratic over the distribution function $\delta f$. In order to make apparent the structure of this equation of motion it is convenient to separate the linear and quadratic terms. Thus, we can also write Eq. (12) as:

$$({\cal O}. \delta f) ({\bf r}) = g \int {\rm d}^7r_1 {\rm d}^... ... r}_1,{\bf r}_2) \; \delta f({\bf r}_1) \; \delta f({\bf r}_2)$$ (13)

where we note ${\bf r}$ the 7-dimensional coordinate ${\bf r}=({\bf k},{\bf p},\tau)$ and we introduced the linear operator ${\cal O}$ and the kernel K defined by:

 

$${ ({\cal O}. \delta f) ({\bf r}) \equiv \frac{\partial\delta f}{\... ...artial{\bf p}}({\bf p}) \int {\rm d}{\bf p}' \; \delta f({\bf k},{\bf p}',\tau)$$ (14)

and:

 

$$K({\bf r};{\bf r}_1,{\bf r}_2) \equiv 2 i \; {\rm e}^{\tau} \; \d... ..._1^2} \cdot \frac{\partial\delta_{\rm D}}{\partial{\bf p}}({\bf p}_2-{\bf p}) .$$ (15)

In the r.h.s. in Eq. (13) the integration over ${\bf r}_1$ and ${\bf r}_2$ runs over all space. In particular, the time coordinates $\tau_1$ and $\tau_2$ are integrated over the whole range $]-\infty,+\infty[$. Note however that the Dirac factor $\delta_{\rm D}(\tau_1-\tau) \delta_{\rm D}(\tau_2-\tau)$ which appears in the kernel K implies that the r.h.s. only depends on the distribution function $\delta f$ at the time $\tau$, as in Eq. (12). This also ensures that causality is not violated. In the relation (13) we also introduced a "coupling constant'' g which is equal to unity: g=1. The advantage of generalizing the quadratic Eq. (13) to arbitrary g will appear later on.

   
3.2 Perturbative expansion

The dynamics of the distribution function $\delta f$ is fully determined by the non-linear Eq. (13) supplemented by initial conditions. At early times $\tau \rightarrow -\infty$ the universe becomes homogeneous and the velocity is given by the Hubble flow (i.e. ${\bf p}\rightarrow 0$) so that the perturbations vanish: $\delta f\rightarrow 0$. Thus, at early times (or at large scales for CDM-like initial conditions) the equations of motion can be linearized and the distribution function is equal to the linear solution $\eta({\bf r})$. As time goes on, the perturbation grows and it eventually enters the non-linear regime. In the mildly non-linear regime, one can seek a solution of Eq. (13) in the form of a perturbative expansion over the linear mode  $\eta({\bf r})$. Thus, we write the expansion:

$$\delta f({\bf r}) = \sum_{n=1}^{\infty} \delta f^{(n)}({\bf r})$$ (16)

with:

$$\delta f^{(1)}({\bf r}) \equiv \eta({\bf r}) \hspace{0.2cm} \mbox{and} \hspace{0.2cm} \delta f^{(n)}({\bf r}) \sim \eta^n .$$ (17)

That is, the term of order n is of the order of $\eta$ to the power n. Substituting the expansion (16) into the equation of motion (13) we obtain:

$${\cal O}. \eta = 0$$ (18)

and for $n \geq 2$:

$${\cal O}.\delta f^{(n)} = g \int {\rm d}{\bf r}_1 {\rm d}{\bf... ...^{n-1} \delta f^{(m)}({\bf r}_1) \delta f^{(n-m)}({\bf r}_2) .$$ (19)

Thus, the linear mode $\eta({\bf r})$ must belong to the kernel of the operator ${\cal O}$. Since we require $\eta \neq 0$ this implies that ${\cal O}$ is not invertible. We shall check below that indeed the kernel of the operator ${\cal O}$ defined by Eq. (14) is not reduced to $\{0\}$. Then, the higher-order terms $\delta f^{(n)}$ are given by the recursion (19).

However, in order to obtain an explicit expression for $\delta f^{(n)}$ we must "invert'' the operator ${\cal O}$ in Eq. (19). More precisely, we need to build a kernel ${\tilde K}({\bf r};{\bf r}_1,{\bf r}_2)$ which satisfies:

$$({\cal O}.{\tilde K})({\bf r};{\bf r}_1,{\bf r}_2) = K({\bf r};{\bf r}_1,{\bf r}_2) .$$ (20)

Indeed, if such a kernel ${\tilde K}$ exists, we can solve the recursion (19) by choosing for the term of order n the value:

$$\delta f^{(n)} = g \int {\rm d}{\bf r}_1 {\rm d}{\bf r}_2 \; ... ...^{n-1} \delta f^{(m)}({\bf r}_1) \delta f^{(n-m)}({\bf r}_2) .$$ (21)

This relation fully defines the expansion (16): it provides an explicit procedure to compute any order term $\delta f^{(n)}$, in a recursive fashion. Note that it is not obvious a priori that a solution ${\tilde K}$ to Eq. (20) exists. On the other hand, if a solution exists it is not unique since we can add to any peculiar solution ${\tilde K}$ any factorized kernel of the form $\eta({\bf r}) H({\bf r}_1,{\bf r}_2)$ where H is arbitrary, as shown by Eq. (18). This amounts to add a factor $\eta$ to all terms $\delta f^{(n)}$. We shall come back to this point below.

Thus, we now need to construct a solution ${\tilde K}$ to Eq. (20). To do so, we use the following trick. First, we define a new operator ${\cal O}_{-1}$ by:

$${\cal O}_{-1}({\bf r}) \equiv \int_0^{\infty} {\rm d}\sigma \; {\rm e}^{-\sigma {\cal O}({\bf r})} .$$ (22)

Since the operator ${\cal O}_{-1}$ is only an intermediate tool Eq. (22) should be interpreted in a formal sense. In particular, one may introduce a large cutoff $\Lambda$ in the intermediate steps of the calculation for the integration over $\sigma$ in the definition (22) and take the limit $\Lambda \rightarrow \infty$ in the final results. The reason behind the definition (22) is that this operator ${\cal O}_{-1}$ is formally the "inverse'' of the operator ${\cal O}$. Indeed, expanding the exponential e $^{-\sigma {\cal O}}$ we can write:

$$\frac{\partial}{\partial\sigma} \left( {\rm e}^{-\sigma {\cal... ...\; {\cal O}^q . f = - {\cal O}. {\rm e}^{-\sigma {\cal O}} . f$$ (23)

where $f({\bf r})$ is an arbitrary function of ${\bf r}$. Using Eq. (23) we obtain:

 

$${{\cal O}. {\cal O}_{-1}. f = \int {\rm d}\sigma \, {\cal O}. {\r... ...c{\partial}{\partial\sigma} \left( {\rm e}^{-\sigma {\cal O}} . f \right) = f }$$ (24)

if all integrals converge. Here, we must point out that Eq. (24) is not a rigorous proof. In particular, using the same arguments we could also obtain ${\cal O}_{-1}. {\cal O}. f=f$. However, we know that ${\cal O}_{-1}. {\cal O}$ cannot be the identity operator since the kernel of ${\cal O}$ is not reduced to $\{0\}$, as shown by Eq. (18) which implies ${\cal O}_{-1}. {\cal O}. \eta=0$. Nevertheless, we can still work with the operator ${\cal O}_{-1}$. In particular, we define a kernel ${\tilde K}({\bf r};{\bf r}_1,{\bf r}_2)$ by:

$${\tilde K}({\bf r};{\bf r}_1,{\bf r}_2) \equiv {\cal O}_{-1}({\bf r}) . K({\bf r};{\bf r}_1,{\bf r}_2) .$$ (25)

If the operator ${\cal O}$ were invertible and all integrals were convergent the kernel ${\tilde K}$ would be the unique solution of Eq. (20). Here, the operator ${\cal O}$ is not invertible hence it is not obvious a priori that the kernel ${\tilde K}$ defined by Eq. (25) is a solution of Eq. (20). However, if this procedure provides a finite kernel ${\tilde K}$ we could expect that it obeys Eq. (20). We compute in Appendices A and B the kernel ${\tilde K}$ defined by Eq. (25). This yields the equivalent expressions (B.1) and (B.3). As discussed in Appendix B we also explicitly check that the kernel ${\tilde K}$ written in Eq. (B.3) is a solution to Eq. (20). Then, we obtain the explicit perturbative expansion of the distribution function $\delta f$ through Eq. (21). As noticed above, at each step n we could add to $\delta f^{(n)}$ an element of the kernel of ${\cal O}$ like $\eta$. However, since we require that $\delta f^{(n)}$ be of the order of $\eta$ to the power n (so that a perturbative expansion makes sense) this additional term is zero and the high-order terms $\delta f^{(n)}$ are indeed given by Eq. (21) for $n \geq 2$.

Here we note that instead of the operator ${\cal O}_{-1}$ defined by Eq. (22) we could have tried to use the operator ${\tilde {\cal O}}_{-1}$ defined by:

$${\tilde {\cal O}}_{-1}({\bf r}) \equiv - \int_0^{\infty} {\rm d}\sigma \; {\rm e}^{\sigma {\cal O}({\bf r})} .$$ (26)

Indeed, we can easily see that this operator also satisfies ${\cal O}. {\tilde {\cal O}}_{-1}= 1$ in a formal sense, following the steps (23) and (24). However, only one of the two choices ${\cal O}_{-1}$ and ${\tilde {\cal O}}_{-1}$ is valid and leads to finite results while the other one involves divergent quantities. In fact, it is easy to see a priori that the correct choice is ${\cal O}_{-1}$ from causality requirements. Indeed, from Eq. (14) we note that ${\cal O}= \partial/ \partial\tau + ..$. Then, since we have:

$${\rm e}^{ - \sigma \; \partial/ \partial\tau} . f(\tau) \equi... ... \frac{\partial^q}{\partial\tau^q} \; f(\tau) = f(\tau-\sigma)$$ (27)

we can see that $\delta f^{(n)}(\tau)$ defined by Eq. (21) using ${\cal O}_{-1}$ involves the linear mode at earlier times $\eta(\tau-\sigma)$ (with $\sigma \geq 0$), as required by causality. On the contrary, using ${\tilde {\cal O}}_{-1}$ would imply that $\delta f^{(n)}(\tau)$ depend on $\eta(\tau+\sigma)$, that is on the linear mode at later times, which violates causality. As a consequence, the correct choice is the operator ${\cal O}_{-1}$. However, the actual justification of Eq. (25) is the explicit check (discussed in Appendix B) that the kernel ${\tilde K}$ derived in Appendices A and B from the definition (25) is a solution to Eq. (20).

From the recursion (21) we can also express $\delta f^{(n)}$ explicitly in terms of the linear mode $\eta$. Indeed, we can easily check that we can write Eq. (21) as:

$$\delta f^{(n)}({\bf r}) = g^{n-1} \int {\rm d}{\bf r}_1 ... {... ...{\bf r}_1, ..., {\bf r}_n) \eta({\bf r}_1) ... \eta({\bf r}_n)$$ (28)

for $n \geq 1$. Here we defined the kernels F_n by:

$$F_1({\bf r};{\bf r}_1) \equiv \delta_{\rm D}({\bf r}-{\bf r}_1)$$ (29)

and for $n \geq 2$:

 

$${ F_n({\bf r};{\bf r}_1, ..., {\bf r}_n) \equiv \int {\rm d}{\bf ... ...1';{\bf r}_1, ...,{\bf r}_m) F_{n-m}({\bf r}_2';{\bf r}_{m+1}, ..., {\bf r}_n).$$ (30)

The expression (28) clearly shows that the expansion (16) is a perturbative expansion over powers of the linear mode $\eta$. It is also apparent from Eq. (28) that $\eta$ should only consist of the linear growing mode. Indeed, the decaying mode would lead to divergences in Eq. (28) as $\tau \rightarrow -\infty$. This is not surprising since at early times the decaying mode becomes increasingly large so that a perturbative approach fails. Of course, in principle it is possible to include decaying modes if we set up the initial conditions at a finite time $\tau_{\rm i}$ (i.e. $t_{\rm i}>0$). Note that in this case we cannot use the kernel ${\tilde K}$ derived from Eq. (25). This can be seen from the fact that if we only include the growing mode but set up the initial conditions at $t_{\rm i}>0$ we should recover the same results for $t>t_{\rm i}$. However, this cannot be achieved from Eq. (21) using the same kernel ${\tilde K}$ since the integration range over $\tau_1$ and $\tau_2$ is changed to $\tau_1>\tau_{\rm i}$ and $\tau_2>\tau_{\rm i}$.

Fortunately, for practical purposes the decaying modes should play no role: $t_{\rm i} \ll t_0$ so that the decaying modes become negligible before one enters the non-linear regime. Hence it is convenient to restrict $\eta({\bf r})$ to the linear growing mode and to take $t_{\rm i}=0$.

   
3.3 Non-perturbative integral equation

Here, we note that the expression (28) also shows that the expansion (16) can be formally reinterpreted as a perturbative expansion over powers of the "coupling constant'' g. Thus, the recursive solution defined by Eqs. (17) and (21) is also the perturbative solution of the integral equation:

$$\delta f({\bf r}) = \eta({\bf r}) + g \int {\rm d}{\bf r}_1 {... ...\bf r}_1,{\bf r}_2) \; \delta f({\bf r}_1) \delta f({\bf r}_2)$$ (31)

where the perturbative parameter is g. Moreover, if we apply to both sides of Eq. (31) the operator ${\cal O}$, using Eq. (20), we obtain Eq. (13). This means that the solution $\delta f[\eta]$ of Eq. (31) is also a solution of the non-perturbative Eq. (13). Hence Eq. (31) is actually equivalent to Eq. (13) supplemented by the initial condition $\delta f\rightarrow \eta$ at early times $\tau \rightarrow -\infty$. Indeed, Eq. (13) admits only one solution (i.e., the distribution function $\delta f$ is completely specified by its value at some initial time $t_{\rm i}$ since Eq. (13) involves a first-order time derivative). Thus, the problem of solving the collisionless Boltzmann equation with this initial condition comes down to look for the solution $\delta f$ of Eq. (31) for a given linear mode $\eta$. We stress here that Eq. (31) is non-perturbative, since it is equivalent to Eq. (13). Note that in our case we shall eventually put g=1. In particular, the relation (31) does not assume that the distribution function $\delta f$ can be written as a perturbative series over the linear mode $\eta$, that is Eq. (31) remains valid even if the perturbative series is only asymptotic (i.e. divergent).

The Eq. (31) is similar to the integral form of usual stochastic differential equations (e.g., Langevin equation) where $\eta$ would be called a noise term. Moreover, this formulation is very convenient since it reduces the three equations we started from (Boltzmann equation, Poisson equation, initial conditions) to the only one Eq. (31). Thus, the initial conditions $\eta({\bf r})$ are encoded within the "evolution equation'' itself, which fully determines the properties of the distribution function $\delta f$ once the characteristics of $\eta({\bf r})$ are specified. Here we note that a quadratic equation of the form (31) was derived in Scoccimarro (2000) from the usual hydrodynamical description. However, the latter formulation breaks down when shell-crossing occurs while Eq. (31) remains valid in the non-linear regime.

   
3.4 Diagrams

The perturbative expansion (21) can be conveniently written as a sum of tree-diagrams with a three-leg vertex. To do so, we first define the symmetrized kernel ${\tilde K}_{\rm s}$ by:

$${\tilde K}_{\rm s}({\bf r};{\bf r}_1,{\bf r}_2) \equiv \frac{... ...{\bf r}_2) + {\tilde K}({\bf r};{\bf r}_2,{\bf r}_1) \right] .$$ (32)

Thus, the relations (21), (30) and (31) remain valid after we replace the kernel ${\tilde K}$ by ${\tilde K}_{\rm s}$. Then, one can write the term $\delta f^{(n)}$ of order n as the sum over all trees built from (n-1) three-leg vertices with n external points ${\bf r}_1, ...,{\bf r}_n$ over which we integrate. For instance, we show in Fig. 1 the two diagrams which give the fourth-order term $\delta f^{(4)}$. To construct such diagrams one simply applies a recursive procedure of (n-1) steps which builds the tree upward. To go from step k to step k+1 one merely chooses an end-point (i.e. with no links upward) and connects to this point two "sons'' through a three-leg vertex. At the first step one starts from the "root'' ${\bf r}$ while the final end-points which are left are labeled ${\bf r}_1, ...,{\bf r}_n$ (the order is irrelevant). Then, each vertex contributes a weight $g {\tilde K}_{\rm s}$ while to each external point is affected a weight $\eta({\bf r}_i)$. Finally, to obtain $\delta f^{(n)}({\bf r})$ one simply integrates over the points ${\bf r}_1, ...,{\bf r}_n$ and takes the sum over all diagrams. Note that each diagram of topology $\alpha$ needs to be multiplied by a multiplicity factor  $a_{n,\alpha}$.

	Figure 1: The two diagrams which give the fourth-order term $\delta f^{(4)}({\bf r})$. The filled circles correspond to the three-leg vertices, which contribute a weight $g {\tilde K}_{\rm s}$ for each one, and the circles represent the external points ${\bf r}_1, ...,{\bf r}_n$. The numbers (1) and (4) give the multiplicity $a_{4,\alpha }$ of each diagram.
Open with DEXTER

It is instructive to consider the simple case where $\delta f$ and $\eta$ are mere numbers, rather than functions. Then we can absorb ${\tilde K}_{\rm s}$ into g and the integral Eq. (31) becomes:

$$\delta f= \eta + g \; \delta f^2 .$$ (33)

Let us note c_n the sum over all multiplicity factors associated with the diagrams of order n, that is: $c_n = \sum_{\alpha} a_{n,\alpha}$ (in particular, c₄=5). Then, the perturbative expansion (16) reads:

$$\delta f= \sum_{n=1}^{\infty} g^{n-1} \; c_n \; \eta^n$$ (34)

with $c_1 \equiv 1$. However, in this simple case we can directly solve the non-linear Eq. (33) which yields (the choice among the two roots of the quadratic Eq. (33) is determined by the constraint $\delta f=\eta$ in the limit $g \rightarrow 0$):

$$\delta f= \frac{1-\sqrt{1 - 4 g \eta}}{2 g} \cdot$$ (35)

Then, expanding the square root which appears in (35) we obtain the coefficients c_n:

$$c_1=1 \hspace{0.2cm} \mbox{and for} \hspace{0.2cm} n \geq 2: \; c_n= \frac{2(2n-3)!}{n! (n-2)!} \cdot$$ (36)

Note that these coefficients c_n are simply the Catalan numbers. From the expression (35) we can see that the perturbative expansion (34) diverges for $\vert g \eta\vert>1/4$. Note that this occurs rather early as the second-order term is only one fourth of the first-order term. This means in particular that we cannot use a local approximation for Eq. (31), which would consist of taking the mean of Eq. (31) over a small phase-space volume V and approximating the integral in the r.h.s. by its value when ${\bf r}_1$ and ${\bf r}_2$ are restricted to V. Of course, this result does not imply that the exact perturbative series diverges for large $\eta$ (i.e. in the non-linear regime), however it suggests that such a behaviour would be quite likely.

   
4 Explicit calculations up to the second-order

We have described in the previous section the general structure of the evolution equation which governs the dynamics of the distribution function $\delta f({\bf r})$. In particular, we have shown how we can obtain a perturbative expansion of $\delta f({\bf r})$ over powers of the linear growing mode $\eta({\bf r})$. In this section, we apply these results to the specific case of gravitational clustering in an expanding universe (with $\Omega_{\rm m}=1$).

   
4.1 Linear growing mode

First, we have to obtain the linear growing mode $\eta({\bf r})$ (as explained in Sect. 3.2 we do not need the decaying modes). This function $\eta({\bf r})$ must belong to the kernel of the operator ${\cal O}$, see Eq. (18), and grow with time. To get $\eta({\bf r})$ one can simply consider the linear growing mode of the hydrodynamical approximation. Thus, in this approximation all particles at a given point ${\bf x}$ have the same peculiar velocity ${\bf v}$ defined by:

$${\bf v}= a \dot{{\bf x}} , \hspace{0.2cm} \mbox{hence} \hspace{0.2cm} {\bf p}=a {\bf v}.$$ (37)

In the linear regime the growing modes of the density and velocity fields are related by (Peebles 1980):

$${\bf v}_{\rm L}({\bf k},\tau) = i \frac{a}{D_+} \dot{D}_+ \fr... ...\; \dot{a} \frac{{\bf k}}{k^2} \; \delta_{\rm L}({\bf k},\tau)$$ (38)

where the growth factor D₊(t) verifies $D_+(t) \propto a(t)$ since we consider a critical density universe. The dot over "$\dot{D_+}$'' or "$\dot{a}$'' stands for the time-derivative " ${\rm d}D_+/{\rm d}t$'' or " ${\rm d}a/{\rm d}t$''. The subscript L refers to the "linear'' regime. Thus, using Eq. (37) we have in the linear regime:

$${\bf p}_{\rm L}({\bf k},\tau) = i \; \dot{a} a \; \frac{{\bf k}}{k^2} \; \delta_{\rm L}({\bf k},\tau) .$$ (39)

Moreover, since $\delta_{\rm L}({\bf k},\tau) \propto a$ we can write:

$$\delta_{\rm L}({\bf k},\tau) = {\rm e}^{2\tau} \delta_{\rm L0... ... e}^{3\tau} \delta_{\rm L0}({\bf k}) \frac{{\bf k}}{k^2} \cdot$$ (40)

Thus the linear growing mode which is solution of the hydrodynamical equations is:

$$\delta f_{\rm hydro}^{(1)}({\bf x},{\bf p},\tau)=(1+\delta_{\... ...}- {\bf p}_{\rm L}({\bf x},\tau)] - \delta_{\rm D}({\bf p}).\,$$ (41)

The Dirac distribution $\delta_{\rm D}[ {\bf p}- {\bf p}_{\rm L}({\bf x},\tau)]$ corresponds to the hydrodynamical description. However, one can check that this function $\delta f_{\rm hydro}^{(1)}$ is not a solution of Eq. (18). To obtain the correct linear mode $\eta({\bf r})$ which belongs to the kernel of ${\cal O}$ one simply needs to keep only the linear term in Eq. (41). Thus, we write:

$$\delta f^{(1)}({\bf x},{\bf p},\tau) = \delta_{\rm L}({\bf x}... ...cdot \frac{\partial\delta_{\rm D}}{\partial{\bf p}}({\bf p}) .$$ (42)

In Fourier space we obtain for $\eta({\bf r}) \equiv \delta f^{(1)}({\bf r})$:

$$\eta({\bf r}) = {\rm e}^{2\tau} \delta_{\rm L0}({\bf k}) \del... ...cdot \frac{\partial\delta_{\rm D}}{\partial{\bf p}}({\bf p}) .$$ (43)

Then, one can easily check that ${\cal O}.\eta=0$ with $\eta$ given by Eq. (43) and ${\cal O}$ by Eq. (14).

Here, we can note that up to first-order our system can be described as a pressure-less fluid. Indeed, we have not included any velocity dispersion in the initial conditions. As is well-known, caustics will form as the density field evolves and multi-streaming regions appear. The perturbative approach should break down at this stage since it cannot go beyond these singularities. Thus, we shall check in the following sections that our perturbative treatment recovers the predictions of the standard perturbative approach applied to the equations of motion derived within the hydrodynamic approximation (see Peebles 1980 or Goroff et al. 1986 for a presentation of this method). This clearly implies that the perturbative method cannot handle the multi-streaming regime. However, we stress that the equations of motion (1), (12) and (31) remain valid in the non-linear regime after multi-streaming appears.

In order to handle the multi-streaming regime in a perturbative way, one may try to add a small velocity dispersion to the initial conditions. Indeed, this will smooth out the caustics encountered for a pressure-less fluid. Such an approach is investigated in Buchert & Dominguez (1998) and Buchert et al. (1999) for instance. However, it is clear that this is not sufficient to go beyond shell-crossing. Indeed, as non-linear structures form, shocks will appear around virialized objects. This can be seen from the analytic solution of gaseous spherical collapse in Bertschinger (1985). Then, a perturbative approach cannot give any information on the inner region beyond the shock. In particular, one cannot rely on a simple equation of state of the form $p=p(\rho)$ for the pressure. Moreover, it is well-known that collisionless and gaseous systems can exhibit qualitatively different behaviours (e.g., existence of Landau damping and phase-mixing for the former, see Binney & Tremaine 1987). This is why it is important to devise methods which are based on the collisionless Boltzmann Eq. (1). This article is a first step in this direction.

   
4.2 Second-order term. Comparison with the hydrodynamical approach

Once the linear growing mode $\eta({\bf r})$ is specified we can derive the higher-order terms from the recursion (21). Thus, using the expression (43) and Eq. (B.1) obtained in Appendix B for the kernel ${\tilde K}$ we get after a simple calculation:

 

$${ \delta f^{(2)}({\bf k},{\bf p},\tau) = \int {\rm d}{\bf k}_1 \; \delta_{\rm L0}({\bf k}_1) \delta_{\rm L0}({\bf k}-{\bf k}_1) }$$

$$\times \biggl \lbrace \left(1+\frac{2}{3} \frac{{\bf k}.({\bf k}-... ...\rm e}^{4 \tau} \frac{{\bf k}.{\bf k}_1}{k_1^2} \delta_{\rm D}({\bf p}) \right.$$

$$\left. - \frac{2i}{5} {\rm e}^{5 \tau} \frac{{\bf k}_1}{k_1^2} \c... ... k}}{k^2} \cdot \frac{\partial\delta_{\rm D}}{\partial{\bf p}}({\bf p}) \right]$$

$$- \frac{2}{9} {\rm e}^{6 \tau} \frac{{\bf k}_1}{k_1^2} \cdot \fra... ...{\partial\delta_{\rm D}}{\partial{\bf p}}({\bf p}) \right) \biggl \rbrace \cdot$$ (44)

To derive Eq. (44) we used identities of the form:

$$F({\bf k}.{\bf p}) \; \delta_{\rm D}({\bf p}) = F(0) \; \delta_{\rm D}({\bf p})$$ (45)

and:

 

$$F({\bf k}_1.{\bf p}) \frac{{\bf k}}{k^2} \cdot \frac{\partial\del... ...f p}}({\bf p}) - F'(0) \frac{{\bf k}.{\bf k}_1}{k^2} \; \delta_{\rm D}({\bf p})$$ (46)

which are valid for any function F. This allows us to eliminate the functions $\psi$ which appear in ${\tilde K}$ in Eq. (B.1) since we can express all factors in terms of $\psi(s,0)=1/s$, see Eq. (B.2).

It is interesting to compare this result for $\delta f^{(2)}$ with the results obtained within the usual hydrodynamical approach. In this latter case, the system is fully defined by the density contrast and the mean fluid velocity at each point ${\bf x}$ (there is no velocity dispersion). Thus, we need to derive the density and the mean collective velocity from Eq. (44).

   
4.2.1 Density contrast

From Eq. (5) and Eq. (44) we obtain the second-order term for the density contrast:

$${ \delta^{(2)}({\bf k},\tau) = {\rm e}^{4 \tau} \int {\rm d}{\bf ... ...\frac{{\bf k}.({\bf k}-{\bf k}_1)}{\vert{\bf k}-{\bf k}_1\vert^2} \right) \cdot$$ (47)

Going back to real space, we get after a change of variables:

$${ \delta^{(2)}({\bf x},\tau) = {\rm e}^{4 \tau} \int {\rm d}{\bf ... ... \frac{2}{7} \left( \frac{{\bf k}_1.{\bf k}_2}{k_1 k_2} \right)^2 \right] \cdot$$ (48)

Thus, we recover the result obtained by the hydrodynamical approach (Peebles 1980).

We can note that this result agrees with the calculations performed in Bharadwaj (1996) from the BBGKY hierarchy. Indeed, in that paper it was shown that the lowest order non-linear correction to the two-point correlation function derived from the BBGKY hierarchy (which remains valid in the non-linear regime after shell-crossing) matches the prediction of the standard hydrodynamical approach. Note that the BBGKY hierarchy can be derived from the collisionless Boltzmann Eq. (1), see Peebles (1980). Therefore, both works explicitly show that multi-streaming effects cannot be handled through a perturbative method.

   
4.2.2 Peculiar velocity

From the distribution function $\delta f$ we can also derive the mean velocity of the fluid. Indeed, the hydrodynamics equations are obtained from moments of the Boltzmann equation. In particular, the hydrodynamic peculiar velocity is defined by:

$$\rho {\bf v}\equiv \int {\rm d}{\bf p}\; f({\bf x},{\bf p},\tau) \; \frac{{\bf p}}{a} \cdot$$ (49)

Since $\rho=\overline{\rho}(1+\delta)$ we obtain using Eq. (4):

$$\left[ 1 + \delta({\bf x},\tau) \right] {\bf v}({\bf x},\tau)... ...int {\rm d}{\bf p}\; \delta f({\bf x},{\bf p},\tau) \; {\bf p}$$ (50)

which yields for the second-order term:

$${\bf v}^{(2)} = {\rm e}^{-2\tau} \int {\rm d}{\bf p}\; \delta f^{(2)} \; {\bf p}\; - \; \delta^{(1)} {\bf v}^{(1)} .$$ (51)

The first-order term is simply ${\bf v}^{(1)}={\bf v}_{\rm L}={\rm e}^{-2\tau}{\bf p}_{\rm L}$ which is given by Eq. (40). Then, we obtain the second-order term from Eq. (44). This yields:

 

$${ {\bf v}^{(2)}({\bf k},\tau) = i \; {\rm e}^{3 \tau} \int {\rm d... ...k}-{\bf k}_1)}{\vert{\bf k}-{\bf k}_1\vert^2} \frac{{\bf k}}{k^2} \right] \cdot$$ (52)

Within the hydrodynamical approach, when one expands the density and velocity fields over the linear growing mode (i.e. there are no decaying modes) the velocity field is nonrotational (e.g., Peebles 1980). Then it is fully defined by its divergence and one usually introduces:

$$\theta({\bf x},\tau) \equiv \frac{\nabla . {\bf v}}{\dot a} ,... ...c{3}{2} \; {\rm e}^{\tau} \; i {\bf k}.{\bf v}({\bf k},\tau) .$$ (53)

In real space, using Eq. (52) we recover the result of the hydrodynamical approach (Peebles 1980):

 

$${ \theta^{(2)}({\bf x},\tau) = - {\rm e}^{4 \tau} \int {\rm d}{\b... ... \frac{4}{7} \left( \frac{{\bf k}_1.{\bf k}_2}{k_1 k_2} \right)^2 \right] \cdot$$ (54)

In a similar fashion, we can consider the normalized vorticity ${\vec \omega}\equiv (\nabla \times {\bf v})/{\dot a}$. Using Eq. (52) we can easily check that we obtain ${\vec \omega}^{(2)}=0$ (and of course the linear mode as defined in Eq. (40) is nonrotational). Thus, we see that up to the second-order we recover the results of the hydrodynamical approach.

   
5 Divergence of perturbative series

In the previous section we noticed that the perturbative treatment of the collisionless Boltzmann equation yields the same results at second-order as the hydrodynamical approach when the linear mode is given by Eq. (43). In fact, it is easy to see that this equality extends to all orders of perturbation theory. Indeed, let us consider an initial condition of the form (43) which is very smooth. That is, $\delta_{\rm L0}({\bf k})$ shows no power at small scales (i.e. large wavenumbers k) and the flow is nonrotational, defined by the linear growing mode of hydrodynamics. Then, there is no shell-crossing until a finite time $\tau_{\rm c}$ (i.e. $t_{\rm c}>0$). This means that the hydrodynamical description is actually exact at early times $t<t_{\rm c}$: at each point there is only one velocity. Then, until the non-zero time $t_{\rm c}$ the solution of hydrodynamics (i.e. continuity and Euler equations) is also a solution of the collisionless Boltzmann equation from which they derive. This implies that the density contrast and velocity fields obtained by the perturbative theory in both frameworks are identical. Moreover, the form of the perturbative expansion does not depend on the properties of $\delta_{\rm L0}({\bf k})$, hence this identity holds for any initial conditions $\delta_{\rm L0}({\bf k})$. As a consequence, if the linear mode is of the form of Eq. (43) the perturbative results obtained from the collisionless Boltzmann equation and from the hydrodynamical description are identical. Here, we note that this identity was explicitly derived at all orders within the framework of the Zel'dovich approximation in Bharadwaj (1996). Therefore, in order to study the non-linear regime one needs to devise non-perturbative methods applied to the collisionless Boltzmann Eq. (1).

On the other hand, in the case of CDM-like power-spectra (or power-law power-spectra) there is some power at all scales. In particular, the density fluctuations increase with no limit at small-scales (if we do not include a cutoff at small scales). Then, as soon as t>0 there is some shell-crossing. This implies that the hydrodynamical solution is no longer a solution of the collisionless Boltzmann equation (actually, the former is not well defined). Then, the perturbative expansion must diverge (otherwise the hydrodynamical solution and the Boltzmann solution would be identical since they would be equal to the same perturbative solution). Hence, we conclude that for hierarchical initial conditions the perturbative expansion is only asymptotic. If we include a cutoff at small scales for the initial power-spectrum the perturbative series may converge until the finite time $t_{\rm c}$ when shell-crossing occurs for the first time (though this is not guaranteed a priori). However, it must diverge at later times which implies that the perturbative series diverges at all times of practical interest.

In fact, the system we investigate here can actually be even more pathological than these results suggest. First, as we show in a companion paper (Paper II) where we develop a non-perturbative method to study the quasi-linear regime, the generating function of the cumulants of the density contrast obtained at leading order has a radius of convergence which is zero independently of shell-crossing for $P(k) \propto k^n$ with n<0. This means that the high-density tail of the probability distribution function of the density contrast cannot be evaluated by perturbative means in this case. Secondly (and more importantly), the terms encountered in perturbative expansions actually diverge beyond some finite order for power-law power-spectra (e.g., Scoccimarro & Frieman 1996, Paper V). This implies that one can only use perturbation theory up to a finite order at best.

   
6 Other cosmologies

In this section, we show how we can extend the perturbative calculations developed in Sect. 3 for a critical density universe to arbitrary cosmological parameters. In the general case, the scale factor a(t) is no longer a power-law hence it is not useful to introduce the time coordinate $\tau$ as in Eq. (11). Thus, from Eq. (9) and keeping the coordinate t we define the operator ${\cal O}$ by:

 

$${ ({\cal O}. \delta f) ({\bf k},{\bf p},t) \equiv \frac{\partial\... ...partial{\bf p}}({\bf p}) \int {\rm d}{\bf p}' \; \delta f({\bf k},{\bf p}',t) .$$ (55)

Next, as in Sect. 4.1 we obtain the growing mode $\eta({\bf k},{\bf p},t)$ by keeping only the linear terms of the linear hydrodynamical growing mode which can still be written as in Eq. (41). This yields:

$$\eta = D_+ \delta_{\rm L0}({\bf k}) \delta_{\rm D}({\bf p}) -... ... \cdot \frac{\partial\delta_{\rm D}}{\partial{\bf p}}({\bf p})$$ (56)

where D₊(t) is the usual linear growing mode, which is no longer proportional to a(t). Then, one can easily check that ${\cal O}.\eta=0$ with $\eta$ given by Eq. (56) and ${\cal O}$ by Eq. (55), using the fact that D₊(t) obeys (Peebles 1980):

$$\ddot{D}_+ + 2 \frac{\dot a}{a} \dot{D}_+ = \frac{4\pi{\cal G}\overline{\rho}}{a^3} \; D_+ \cdot$$ (57)

Then, in order to get the second-order term we simply keep the structure of the result (44) obtained for the critical density universe but we replace all time-dependent factors by unknown functions. Thus, we write:

 

$${ \delta f^{(2)}({\bf k},{\bf p},t) = \int {\rm d}{\bf k}_1 \; \delta_{\rm L0}({\bf k}_1) \delta_{\rm L0}({\bf k}-{\bf k}_1) }$$

$$\times \biggl \lbrace g_1 \frac{{\bf k}.{\bf k}_1}{k_1^2} \delta_... ... k}.({\bf k}-{\bf k}_1)}{\vert{\bf k}-{\bf k}_1\vert^2} \delta_{\rm D}({\bf p})$$

$$- i g_3 \frac{{\bf k}_1}{k_1^2} \cdot \frac{\partial\delta_{\rm D... ...{\bf k}_1}{k_1^2} \cdot \frac{\partial\delta_{\rm D}}{\partial{\bf p}}({\bf p})$$

$$- i g_5 \frac{{\bf k}.{\bf k}_1}{k_1^2} \frac{{\bf k}}{k^2} \cdot... ...rac{{\bf k}}{k^2} \cdot \frac{\partial\delta_{\rm D}}{\partial{\bf p}}({\bf p})$$

$$- g_7 \frac{{\bf k}_1}{k_1^2} \cdot \frac{\partial}{\partial{\bf ... ... \frac{\partial\delta_{\rm D}}{\partial{\bf p}}({\bf p}) \right) \biggl \rbrace$$ (58)

where g₁(t),...,g₇(t) are functions of time. Next, we substitute Eq. (58) into Eq. (9), using Eq. (56) for the first-order term of $\delta f$, and keeping only the second-order terms we obtain a system of differential equations for the unknown functions g₁(t),...,g₇(t). This yields:

$$a^2 {\dot g}_1 = g_3+g_5, \;\; a^2 {\dot g}_2 = g_4+g_6, \;\; a {\dot g}_3 = 4\pi{\cal G}\overline{\rho}D_+^2$$ (59)

$$a^2 {\dot g}_4 = 2 g_7, \;\; a {\dot g}_5 = 4\pi{\cal G}\overline{\rho}g_1$$ (60)

$$a {\dot g}_6 = 4\pi{\cal G}\overline{\rho}g_2, \;\; {\dot g}_7 = 4\pi{\cal G}\overline{\rho}D_+^2$$ (61)

which determines the functions g₁(t),...,g₇(t). The boundary condition at t=0 is given by the constraint that for $t \rightarrow 0$ we must recover the result (44) obtained for the critical density universe. We shall not go further here since the discussion of Sect. 5 remains valid for arbitrary cosmologies hence we already know that we shall eventually recover the results obtained through the hydrodynamical approach, described in Bouchet et al. (1992).

Thus, in order to derive the perturbative expansion of the distribution function $\delta f$ in the general case one first computes the terms $\delta f^{(n)}$ up to the required order for a critical density universe, using the method developed in Sect. 3. This provides an explicit expression for $\delta f^{(n)}$. Then, one looks for a solution of Eq. (9) in the case of interest (which differs from the critical density universe) by keeping the same structure for $\delta f^{(n)}$ but replacing all time-dependent factors by unknown functions of time. Then, these functions are determined by substitution into Eq. (9). The need for this two-step procedure comes from the fact that in the general case the calculation of the operator ${\cal O}_{-1}$ is not straightforward. More precisely, in order to derive ${\cal O}_{-1}$ one can still follow the procedure outlined in Appendix A up to Eq. (A.12). However, the recursion obeyed by the functions $R_n(\sigma)$ is no longer a simple Laplace convolution so that it is not obvious to find a simple analytic solution.

   
7 Functional formulation

Finally, we briefly describe how one can use the integral Eq. (31) to obtain an explicit expression for the statistical properties of the distribution function. In this section, we shall work in real space ${\bf x}$, and we note ${\vec \omega}$ the 7-dimensional coordinate ${\vec \omega}=({\bf x},{\bf p},\tau)$. Thus, the fields $\eta({\vec \omega})$ and $\delta f({\vec \omega})$ are real. Taking the Fourier transform of Eq. (31) we see that $\delta f({\vec \omega})$ and $\eta({\vec \omega})$ are related by the same quadratic equation where the kernel ${\tilde K}({\bf r};{\bf r}_1,{\bf r}_2)$ is replaced by its real-space Fourier transform ${\tilde K}({\vec \omega};\bom_1,\bom_2)$. In practice we are not really interested in the exact solution $\delta f$ associated with a given field $\eta$. Indeed, since the initial condition $\eta$ is a random field we are only interested in the statistical properties of the distribution function $\delta f$. These can be derived from the functional Z[j] of the test field $j({\vec \omega})$ defined by:

$$Z[j] \equiv \bigg\langle{\rm e}^{\int {\rm d}{\vec \omega}\; j.\delta f} \bigg\rangle$$ (62)

where $\langle... \rangle$ expresses the average over the initial conditions. Expanding the exponential we get:

 

$$Z[j] = 1 + \sum_{n=1}^{\infty} \frac{1}{n!} \int {\rm d}\bom_1 ..... ...j(\bom_1) ... j(\bom_n) \; \langle\delta f(\bom_1) ... \delta f(\bom_n) \rangle$$ (63)

which clearly shows that Z[j] yields all moments $\langle\delta f(\bom_1) ... \delta f(\bom_n) \rangle$ of the distribution function $\delta f$. If we assume the initial conditions to be Gaussian we can write the average (62) as the path-integral:

$$Z[j] = {\cal N}_0\int [{\rm d}\eta({\vec \omega})] \; {\rm e}^{j.\delta f- \frac{1}{2} \eta . \Delta_{\eta}^{-1}. \eta}$$ (64)

where we introduced the normalization constant:

$${\cal N}_0\equiv \left( \mbox{Det}\Delta_{\eta}^{-1}\right)^{1/2} .$$ (65)

Here $\delta f[\eta]$ must be understood as the solution of Eq. (31) (in real ${\vec \omega}$-space) associated with a given $\eta$. We also used the short-hand notations: $j.\delta f\equiv \int {\rm d}{\vec \omega}\; j({\vec \omega}) . \delta f({\vec \omega})$ and $\eta . \Delta_{\eta}^{-1}. \eta \equiv \int {\rm d}\bom_1 {\rm d}\bom_2 \; \eta(\bom_1) . \Delta_{\eta}^{-1}(\bom_1,\bom_2) . \eta(\bom_2)$. The kernel $\Delta_{\eta}^{-1}$ is the inverse of the kernel $\Delta_{\eta}$ defined by:

$$\Delta_{\eta}(\bom_1,\bom_2) \equiv \langle\eta(\bom_1) \eta(\bom_2) \rangle$$ (66)

which fully determines the statistics of the random field $\eta({\vec \omega})$ since the latter is assumed to be Gaussian. Here we must note that the kernel $\Delta_{\eta}$ is not invertible. This is related to the fact that the linear mode $\eta({\vec \omega})$ is restricted to the form (43), which also ensures that ${\cal O}.\eta=0$. Thus, in Eq. (64) it is understood that the kernels $\Delta_{\eta}$ and $\Delta_{\eta}^{-1}$ are regularized through a small modification described by an infinitesimal parameter $\epsilon$, which makes $\Delta_{\eta}$ and $\Delta_{\eta}^{-1}$ invertible. For instance, one can add to the kernel $\Delta_{\eta}$ a term of the form $\epsilon \times \delta_{\rm D}(\bom_1-\bom_2)$. Then, the final results are obtained by taking the limit $\epsilon \rightarrow 0$ in the calculations. Using Eq. (31) we can make the change of variable $\eta \rightarrow \delta f$ in the expression (64). This yields:

 

$$Z[j] = {\cal N}_0\int [{\rm d}\delta f({\vec \omega})] \; \vert \... ...{\tilde K}\delta f^2) . \Delta_{\eta}^{-1}. (\delta f- g {\tilde K}\delta f^2)}$$ (67)

where we used the notation:

$$({\tilde K}\delta f^2)({\vec \omega}) \equiv \int {\rm d}\bom... ...vec \omega};\bom_1,\bom_2) \delta f(\bom_1) \delta f(\bom_2) .$$ (68)

In Eq. (67) we introduced the matrix M defined by:

 

$${ M(\bom_1,\bom_2) \equiv \frac{\delta \eta(\bom_1)}{\delta ( \delta f(\bom_2) )} } = \delta_{\rm D}(\bom_1 - \bom_2) - 2 g \int {\rm d}{\vec \omega}' \; {\tilde K}_{\rm s}(\bom_1;\bom_2,{\vec \omega}') \delta f({\vec \omega}') .$$ (69)

Here we used Eq. (31) (in real ${\vec \omega}$-space) and ${\tilde K}_{\rm s}$ is the symmetrized part of ${\tilde K}$, as in Eq. (32). Next, we need to compute the Jacobian $\vert\mbox{Det}M\vert$. To do so, we use the relation:

$$\mbox{Det}M = {\rm e}^{{\rm Tr} \; \ln M}$$ (70)

and from Eq. (69) we have:

$$\mbox{Tr}\; \ln M = - \sum_{n=1}^{\infty} \frac{(2 g)^n}{n} \; \mbox{Tr}\; A^n$$ (71)

where we introduced the matrix A:

$$A(\bom_1,\bom_2) \equiv \int {\rm d}{\vec \omega}' \; {\tilde... ...s}(\bom_1;\bom_2,{\vec \omega}') \; \delta f({\vec \omega}') .$$ (72)

Hence a simple recursion yields for $n \geq 2$:

 

$${ A^n({\vec \omega},{\vec \omega}') = \int {\rm d}\bom_1 ...{\rm ... ...m_1;\bom_2,\bom_2') ... {\tilde K}_{\rm s}(\bom_{n-1};{\vec \omega}',\bom_n') .$$ (73)

Next, we note from Eq. (B.1) that the kernel ${\tilde K}_{\rm s}({\vec \omega};\bom_1,\bom_2)$ involves a factor $\theta(\tau-\tau_1) \delta_{\rm D}(\tau_2-\tau_1)$ (which also ensures that causality is satisfied). Then, we see that the integrand in Eq. (73) involves a factor $\theta(\tau-\tau_1) ... \theta(\tau_{n-1}-\tau')$. Thus, the Heaviside factors imply that the time-coordinates must obey: $\tau>\tau_1>...>\tau_{n-1}>\tau'$ in order that the integrand be non-zero. This also implies that $A^n({\vec \omega},{\vec \omega})=0$ for $n \geq 2$. On the other hand, the trace $\mbox{Tr}A^n$ is defined by:

$$\mbox{Tr}\; A^n \equiv \int {\rm d}{\vec \omega}\; A^n({\vec \omega},{\vec \omega}) .$$ (74)

Hence we see that $\mbox{Tr}A^n=0$ for $n \geq 2$. Note that this result only relies on causality requirements (which translate into the Heaviside factors) and it also applies to usual stochastic differential equations (e.g., Langevin equation, see Zinn-Justin 1989). Thus, we are left with $\mbox{Det}M = \exp(-2 g \mbox{Tr}A)$, with:

$${ 2 \mbox{Tr}A = 2 \int {\rm d}{\vec \omega}{\rm d}{\vec \omega}'... ...rm s}({\vec \omega};{\vec \omega},{\vec \omega}') \; \delta f({\vec \omega}') }$$

$$= \int {\rm d}{\vec \omega}{\rm d}{\vec \omega}' \; \left[ {\tild... ...vec \omega};{\vec \omega}',{\vec \omega}) \right] \; \delta f({\vec \omega}') .$$ (75)

Then, from Eq. (B.1) we can see that:

$$\int {\rm d}{\vec \omega}\; {\tilde K}({\vec \omega};{\vec \o... ...; {\tilde K}({\vec \omega};{\vec \omega}',{\vec \omega}) = 0 .$$ (76)

Indeed, the second term in Eq. (B.1) vanishes since $\tau=\tau_1$ and we can push the integration path over s₁ to the right: Re $(s_1) \rightarrow \infty$. The first term vanishes after we integrate over s₁ and next over ${\bf p}$ or over the angles of ${\bf k}'$ (using the fact that ${\tilde K}({\vec \omega};{\vec \omega}',{\vec \omega}) \propto \delta_{\rm D}({\bf k}')$). Note that this latter integration involves the "regularization'' of the gravitational interaction discussed in Sect. 2 below Eq. (1). Indeed, as discussed in Sect. 2 apparently undetermined integrals are given a meaning by the prescription that one must first perform the integration over angles (or introduce a large-scale cutoff $k_{\rm c} \rightarrow 0$). Thus, we get $\mbox{Det}M=1$ and we can write Eq. (67) as:

 

$$Z[j] = {\cal N}_0\int [{\rm d}\delta f({\vec \omega})] {\rm e}^{j... ...}^{-1}. \delta f+ \frac{g}{3!} K_3 \delta f^3 - \frac{g^2}{4!} K_4 \delta f^4}.$$ (77)

Here we used the short-hand notations:

 

$$K_3 \; \delta f^3 \equiv \sum_{\small {\rm perm.}} \int {\rm d}{\... ...Delta_{\eta}^{-1}(\bom_1,{\vec \omega}) {\tilde K}({\vec \omega};\bom_2,\bom_3)$$ (78)

and:

 

$${ K_4 \delta f^4 \equiv \frac{1}{2} \sum_{\small {\rm perm.}} \in... ...ta}^{-1}({\vec \omega},{\vec \omega}') {\tilde K}({\vec \omega}';\bom_3,\bom_4)$$ (79)

where the sums run over the permutations of $(\bom_1,\bom_2,\bom_3)$ and $(\bom_1,\bom_2,\bom_3,\bom_4)$. We symmetrized the kernels K₃ and K₄ in Eqs. (78) and (79) because it simplifies perturbative expansions. Thus, the path-integral (77) yields an explicit expression for the functional Z[j], hence for all statistical properties of the distribution function $\delta f({\vec \omega})$. Indeed, all terms in this expression are known. The kernel ${\tilde K}$ which yields K₃ and K₄ is given in Eq. (B.1) while the kernel $\Delta_{\eta}^{-1}$ is obtained from $\Delta_{\eta}$ (after a suitable regularization). The latter is given by Eqs. (66) and (43) and it can be expressed in terms of the power-spectrum P(k) of the linear density field $\delta_{\rm L0}({\bf k})$.

The correlation functions can be obtained from Eq. (77) as a perturbative series over the coupling constant g. Indeed, we noticed in Sect. 3.3 that this is equivalent to the usual expansion over the linear growing mode $\eta({\vec \omega})$. To do so we simply develop the exponential term in Eq. (77) over g K₃ and g² K₄ as a series over g. This corresponds exactly to the Feynman diagrams in standard Quantum Field Theories for a " $\phi^3+\phi^4$'' theory (e.g., Zinn-Justin 1989) (note that we defined the kernels K₃ and K₄ in such a way that they are symmetric). One must merely pay attention to the fact that the kernels K₃ and K₄ are not pure numbers. Of course, this is due to the long-range character of gravity which leads to a non-local theory (while usual Quantum Field Theories are local). Note that we do not need the explicit expression of the inverse kernel $\Delta_{\eta}^{-1}$. For instance, in order to compute the skewness some terms which depend on $\Delta_{\eta}^{-1}$ appear but they actually cancel one another.

In practice, the expression (77) is not very convenient for perturbative calculations and it is easier to use the direct method described in Sect. 3, or simply use the standard hydrodynamical approach. However, we stress that Eq. (77) presents the advantage to be an explicit non-perturbative expression for the functional Z[j]. Thus, one may hope to be able to extract some useful information from this relation for the non-linear regime. However, this is beyond the scope of this article.

   
8 Conclusion

Thus, in this article we have developed a systematic method to obtain the solution of the collisionless Boltzmann equation, which describes large-scale structure formation in the universe, as a perturbative series. Moreover, we have explained that these perturbative results are identical to those obtained from the hydrodynamical description of the fluid, if we start with the same initial conditions. Thus, multi-streaming cannot be handled through perturbative methods. Besides, for hierarchical scenarios the perturbative expansions diverge. The interest of this approach is that, contrary to the hydrodynamical model, the collisionless Boltzmann equation provides a rigorous description of the dynamics, even in the non-linear regime. Hence in order to tackle the non-linear regime where shell-crossing plays a key role one should use these equations of motions. Then, the results presented in this paper may serve as a basis for a study of the non-linear regime. In particular, we have obtained a non-perturbative integral Eq. (31) which explicitly relates the non-linear distribution function to the initial conditions. Besides, we have derived an explicit path-integral expression (77) which yields the correlation functions of the actual non-linear density field at all orders. Although this result provides a formal description of the statistics of the density field one still needs to check whether it is possible to extract some practical information from this expression. Thus, we shall describe in companion papers an alternative non-perturbative method, based on the path-integral (64), which allows one to derive rigorous results in the quasi-linear regime for the probability distribution of the density contrast for Gaussian initial conditions (Paper II) as well as for some non-Gaussian scenarios (Paper III). Such a formalism can also be used in the non-linear regime (Paper IV).

   
Appendix A: Calculation of the exponential ${\rm e}^{-\sigma {\cal O}}$

In order to obtain the kernel ${\tilde K}$ defined in Eq. (25) we first need to compute the exponential e $^{-\sigma {\cal O}}$. Thus, for a given function $F({\bf r})=F({\bf k},{\bf p},\tau)$ we look for the function $B({\bf r},\sigma)$ defined by:

$$B({\bf r},\sigma) \equiv {\rm e}^{-\sigma {\cal O}({\bf r})} . F$$ (A.1)

To derive $B({\bf r},\sigma)$ we note that this function satisfies the initial-value problem:

$$\frac{\partial B}{\partial\sigma} = - {\cal O}. B \hspace{0.3cm} \mbox{and} \hspace{0.3cm} B({\bf r},0) = F({\bf r}) .$$ (A.2)

Moreover, we are only interested in the restriction of $B({\bf r},\sigma)$ to $\sigma \geq 0$. Using Eq. (14) we obtain the first-order integro-differential equation:

 

$${ \frac{\partial B}{\partial\sigma} + \frac{\partial B}{\partial\... ...tial{\bf p}}({\bf p}) \int {\rm d}{\bf p}' \; B({\bf k},{\bf p}',\tau,\sigma) .$$ (A.3)

This equation is formally solved by the method of characteristics which yields the integral form:

 

$${ B({\bf k},{\bf p},\tau+\sigma,\sigma) = B({\bf k},{\bf p},\tau,0) {\rm e}^{3i ({\bf k}.{\bf p}) {\rm e}^{-\tau} [ {\rm e}^{-\sigma}-1 ] } } - 2 i \frac{{\bf k}}{k^2} \cdot \frac{\partial\delta_{\rm D}}{\pa... ...{3i ({\bf k}.{\bf p}) {\rm e}^{-\tau} [ {\rm e}^{-\sigma}-{\rm e}^{-\sigma'} ]} \times \; {\rm e}^{\tau+\sigma'} \int {\rm d}{\bf p}' \; B({\bf k},{\bf p}',\tau+\sigma',\sigma') .$$ (A.4)

The form of Eq. (A.4) shows that we can consider in the next intermediate steps ${\bf k}$ and $\tau$ to be fixed constants and we are led to introduce the functions:

$$G({\bf p},\sigma) \equiv B({\bf k},{\bf p},\tau+\sigma,\sigma)$$ (A.5)

and:

$$U({\bf p},\sigma) \equiv B({\bf k},{\bf p},\tau,0) \; {\rm e}... ... ({\bf k}.{\bf p}) {\rm e}^{-\tau} [ {\rm e}^{-\sigma}-1 ] } ,$$ (A.6)

as well as the kernel:

$${\cal K}({\bf p},\sigma\vert{\bf p}',\sigma') \equiv \theta(\sigm... ...{3i ({\bf k}.{\bf p}) {\rm e}^{-\tau} [ {\rm e}^{-\sigma}-{\rm e}^{-\sigma'} ]}$$ (A.7)

where $\theta(\sigma-\sigma')$ is Heaviside's function. Then, Eq. (A.4) reads:

$$G({\bf p},\sigma)=U({\bf p},\sigma) + \nu \int {\cal K}({\bf ... ...}',\sigma') G({\bf p}',\sigma') {\rm d}{\bf p}' {\rm d}\sigma'$$ (A.8)

where we introduced the constant:

$$\nu \equiv - 2 i \; {\rm e}^{\tau}$$ (A.9)

and we integrate over $\sigma'$ from 0 to $+\infty$: $\int {\rm d}\sigma' \equiv \int_0^{\infty} {\rm d}\sigma'$. The linear integral Eq. (A.8) is a standard Fredholm equation of the second kind and it can be solved as a perturbative series over $\nu$:

 

$${ G({\bf p},\sigma) = U({\bf p},\sigma) } + \sum_{n=1}^{\infty} \... ...\sigma\vert{\bf p}',\sigma') U({\bf p}',\sigma') {\rm d}{\bf p}' {\rm d}\sigma'$$ (A.10)

where the kernels ${\cal K}_n$ are defined by the recursion:

$${\cal K}_n \equiv \int {\rm d}{\bf p}'' {\rm d}\sigma'' {\cal... ...gma'') {\cal K}_{n-1}({\bf p}'',\sigma''\vert{\bf p}',\sigma')$$ (A.11)

with ${\cal K}_1({\bf p},\sigma\vert{\bf p}',\sigma') \equiv {\cal K}({\bf p},\sigma\vert{\bf p}',\sigma')$. Thus, in order to obtain $G({\bf p},\sigma)$ we simply need to sum up the series in the r.h.s. of Eq. (A.10). One can check by substitution into Eq. (A.11) that these kernels ${\cal K}_n$ are given for $n \geq 2$ by:

 

$${\cal K}_n = \frac{{\bf k}}{k^2} \cdot \frac{\partial\delta_{\rm ... ... e}^{-\tau} [ {\rm e}^{-\sigma}-{\rm e}^{-\sigma_{n-1}} ]} \; R_n(\sigma_{n-1})$$ (A.12)

where we introduced:

 

$${R_n(\sigma_{n-1}) \equiv \int_0^{\infty} {\rm d}\sigma_1 ... {\r... ...\sigma_{n-2}} \right] ... \left[ {\rm e}^{-\sigma_1}-{\rm e}^{-\sigma'} \right]$$ (A.13)

for $n \geq 3$ and $R_2(\sigma_1) \equiv \theta(\sigma_1-\sigma') \; \left[ {\rm e}^{\sigma'}-{\rm e}^{\sigma_1} \right]$. One can see from Eq. (A.13) that the functions $R_n(\sigma)$ actually satisfy the recursion:

$$R_n(\sigma) = \int_0^{\sigma} {\rm d}l \; \left( 1 - {\rm e}^{\sigma-l} \right) \; R_{n-1}(l)$$ (A.14)

which has the form of a simple Laplace convolution. This recursion is easily solved by taking the Laplace transform of Eq. (A.14) which yields, after going back to $\sigma$-space through an inverse Laplace transform:

$$R_n(\sigma) = {\rm e}^{\sigma'} \int_{-i\infty}^{+i\infty}\fr... ...igma-\sigma')} \; \left[ \frac{-1}{s (s-1)} \right]^{n-1}\cdot$$ (A.15)

Here and in the following calculations the integration path over the Laplace variables like s will always be to the right of the singularities of the integrand. Thus, in Eq. (A.15) the integration path for s obeys: $\mbox{Re}(s)>1$. Then, Eq. (A.15) provides an explicit expression for the kernels ${\cal K}_n$, see Eq. (A.12). Next, we can now compute the sum over n in Eq. (A.10), which eventually provides an explicit expression for $B({\bf r},\sigma)$ using Eq. (A.5). Thus, a simple calculation gives:

 

$${ B({\bf k},{\bf p},\tau,\sigma) = {\rm e}^{3i ({\bf k}.{\bf p}) {\rm e}^{-\tau} [ 1-{\rm e}^{\sigma} ]} F({\bf k},{\bf p},\tau-\sigma) }$$

$$- 2 i \frac{{\bf k}}{k^2} \cdot \frac{\partial\delta_{\rm D}}{\pa... ...t \frac{{\rm d}s}{2\pi i} {\rm e}^{s(\sigma_1-\sigma')} \theta(\sigma-\sigma_1)$$

$$\times \; {\rm e}^{3i ({\bf k}.{\bf p}) {\rm e}^{-\tau} [ 1-{\rm ... ...bf k}.{\bf p}') {\rm e}^{-\tau} [ {\rm e}^{\sigma-\sigma'}-{\rm e}^{\sigma} ] }$$

$$\times \; {\rm e}^{\tau+\sigma'-\sigma} \; \frac{s(s-1)}{(s-3)(s+2)} \; F({\bf k},{\bf p}',\tau-\sigma) .$$ (A.16)

It is convenient to replace the integrations over $\sigma'$ and $\sigma_1$ in Eq. (A.16) by integrals over the conjugate Laplace coordinates s' and s₁. To do so, we simply compute the integrals over $\sigma'$ and $\sigma_1$ in Eq. (A.16) in terms of the function $\psi(s,x)$ which we define by:

$$\psi(s,x) \equiv \int_0^{\infty} {\rm d}\sigma \; {\rm e}^{- s \sigma + x [1-{\rm e}^{\sigma}]} .$$ (A.17)

In the following the argument x of the function $\psi(s,x)$ will be imaginary. Then, we have the asymptotic behaviour for large positive s:

$$\mbox{Re}(x)=0 , \; \mbox{Re}(s) \rightarrow +\infty: \;\; \psi(s,x) \rightarrow 0$$ (A.18)

since for imaginary x and $\mbox{Re}(s)>0$ Eq. (A.17) implies: $\vert\psi(s,x)\vert \leq 1/\mbox{Re}(s)$. The integral representation (A.17) only applies to $\mbox{Re}(s)>0$ however the function $\psi$ can be continued to $\mbox{Re}(s) \leq 0$. Indeed, from Eq. (A.17) one can check that $\psi(s,x)$ can be expressed in terms of the incomplete Gamma function $\Gamma$ or the confluent hypergeometric function $\psi_{\rm K}$ (here we add a subscript "K'' to distinguish Kummer's function $\psi_{\rm K}$ from $\psi$ defined in Eq. (A.17)) by (e.g., Bateman 1953):

$$\psi(s,x) = {\rm e}^x x^s \Gamma[-s,x] = \psi_{\rm K}(1,1-s;x) .$$ (A.19)

This implies that $\psi(s,x)$ is an entire function of s. In particular, one obtains from Eq. (A.19) the asymptotic behaviour for large negative s at fixed x (see Bateman 1953, vol. I, Sect. 6.13.2 (6)):

$$\mbox{Re}(s) \rightarrow -\infty: \;\; \psi(s,x) \sim x^s {\rm e}^{s-(s+1/2) \ln(-s)}$$ (A.20)

which shows that $\psi(s,x)$ diverges for $\mbox{Re}(s) \rightarrow -\infty$. On the other hand, from Eq. (A.17) we see that the function $\psi(s,x)$ also satisfies the relation (inverse Laplace transform):

$${\rm e}^{- s_1 \sigma + x [1-{\rm e}^{\sigma}]} = \int_{-i\in... ...rac{{\rm d}s}{2\pi i} \; {\rm e}^{s \sigma} \; \psi(s+s_1,x) .$$ (A.21)

Then, we can write Eq. (A.16) as:

 

$${ B({\bf k},{\bf p},\tau,\sigma) = \int \frac{{\rm d}s_1}{2\pi i}... ...bf k}.{\bf p}) {\rm e}^{-\tau}) \biggl \lbrace F({\bf k},{\bf p},\tau-\sigma) }$$

$$- 2i \; {\rm e}^{\tau} \frac{{\bf k}}{k^2} \cdot \frac{\partial\d... ... d}s_3}{(2\pi i)^2} {\rm e}^{(s_2+s_3)\sigma} \frac{s_1(s_1-1)}{(s_1-3)(s_1+2)}$$

$$\times \; \frac{1}{s_1+s_2-1} \int {\rm d}{\bf p}' \; \psi(s_2,-3i ({\bf k}.{\bf p}') {\rm e}^{-\tau})$$

$$\times \; \psi(s_3+1,3i ({\bf k}.{\bf p}') {\rm e}^{-\tau}) \; F({\bf k},{\bf p}',\tau-\sigma) \biggl \rbrace \cdot$$ (A.22)

Using the definition (22) this result provides the operator ${\cal O}_{-1}$ through the formal integral:

$$({\cal O}_{-1}. F)({\bf r}) = \int_0^{\infty} {\rm d}\sigma \; B({\bf r},\sigma) .$$ (A.23)

   
Appendix B: Calculation of the kernel ${\tilde K}({\bf r},{\bf r}_1,{\bf r}_2)$

From Eq. (A.22) and the expression (A.23) for the operator ${\cal O}_{-1}$ we can derive the explicit form of the kernel ${\tilde K}({\bf r};{\bf r}_1,{\bf r}_2)$ defined in Eq. (25), using Eq. (15). First, we write the factor $\delta_{\rm D}(\tau_1-\tau) \delta_{\rm D}(\tau_2-\tau)$ which appears in Eq. (15) as $\delta_{\rm D}(\tau_1-\tau) \delta_{\rm D}(\tau_2-\tau_1)$. Then, using Eq. (A.22) we see that $B({\bf r},\sigma) = {\rm e}^{-\sigma {\cal O}} . K$ exhibits a factor $\delta_{\rm D}(\tau_1-\tau+\sigma)$. Thus, the integration on $\sigma$ over the range $[0,\infty[$ is now straightforward and it yields a Heaviside factor $\theta(\tau-\tau_1)$. Next, a simple calculation yields for ${\tilde K}({\bf r};{\bf r}_1,{\bf r}_2)$:

 

$${ {\tilde K}= 2i \; {\rm e}^{\tau_1} \; \theta(\tau-\tau_1) \; \delta_{\rm D}(\tau_2-\tau_1) \; \delta_{\rm D}({\bf k}_1+{\bf k}_2-{\bf k}) }$$

$$\times \int \frac{{\rm d}s_1}{2\pi i} {\rm e}^{s_1(\tau-\tau_1)} ... ...{k_1^2} \cdot \frac{\partial\delta_{\rm D}}{\partial{\bf p}}({\bf p}_2-{\bf p})$$

$$- \; 6 \; \frac{{\bf k}_1 . {\bf k}}{k_1^2} \; \frac{{\bf k}}{k^2... ...nt \frac{{\rm d}s_2 {\rm d}s_3}{(2\pi i)^2} \; {\rm e}^{(s_2+s_3)(\tau-\tau_1)}$$

$$\times \; \frac{s_1(s_1-1)}{(s_1-3)(s_1+2)(s_1+s_2)(s_1+s_2-1)}$$

$$\times \; \psi(s_2,-3i ({\bf k}.{\bf p}_2) {\rm e}^{-\tau}) \psi(s_3,3i ({\bf k}.{\bf p}_2) {\rm e}^{-\tau}) \biggl \rbrace$$ (B.1)

where the integration paths for the variables s₁ and s₂ obey: $\mbox{Re}(s_1)>3$ and $\mbox{Re}(s_2)>0$. To get Eq. (B.1) we used the following properties of $\psi(s,x)$ which are obvious from Eq. (A.17):

$$\psi(s,0) = \frac{1}{s} , \hspace{0.2cm} \frac{\partial\psi}{\partial x}(s,x) = \psi(s,x) - \psi(s-1,x) .$$ (B.2)

Note the factor $\theta(\tau-\tau_1)$ in Eq. (B.1) which explicitly shows that the expansion (21) satisfies causality requirements. This agrees with the discussion in Sect. 3.2.

Finally, as discussed in Sect. 3.2 we need to check that the kernel ${\tilde K}({\bf r};{\bf r}_1,{\bf r}_2)$ obtained in Eq. (B.1) satisfies Eq. (20). To do so, it is convenient to first integrate over s₁ in Eq. (B.1). The integration of the first term is straightforward, using the inverse Laplace transform (A.21). In order to integrate the second term over s₁ we use the identity (46). This allows us to express this factor in terms of $\psi(s_1,0)=1/s_1$ and $\psi(s_1-1,0)=1/(s_1-1)$, so that the integration over s₁ becomes quite simple. In particular, the contribution from the poles s₁=-s₂ and s₁=1-s₂ vanishes when we take the integration path over s₂ to be $\mbox{Re}(s_2)>3$. Then, we can write the kernel ${\tilde K}$ as:

 

$${ {\tilde K}= 2i \; {\rm e}^{\tau_1} \; \delta_{\rm D}(\tau_2-\tau_1) \; \delta_{\rm D}({\bf k}_1+{\bf k}_2-{\bf k}) }$$

$$\times \biggl \lbrace \theta(\tau-\tau_1) \; \frac{{\bf k}_1}{k_1... ...bf p}) \; {\rm e}^{3i ({\bf k}.{\bf p}) [ {\rm e}^{-\tau} -{\rm e}^{-\tau_1} ]}$$

$$- \; \frac{6}{5} \; \frac{{\bf k}_1 . {\bf k}}{k_1^2} \int_{-i\in... ...ty}\frac{{\rm d}s_2 {\rm d}s_3}{(2\pi i)^2} \; {\rm e}^{(s_2+s_3)(\tau-\tau_1)}$$

$$\times \; \psi\left(s_2,-3i ({\bf k}.{\bf p}_2) {\rm e}^{-\tau}\right) \; \psi\left(s_3,3i ({\bf k}.{\bf p}_2) {\rm e}^{-\tau}\right)$$

$$\times \biggl [ \frac{{\bf k}}{k^2} \cdot \frac{\partial\delta_{\... ...(s_2+2)(s_2+3)} + \frac{3 \; {\rm e}^{-2(\tau-\tau_1)}}{(s_2-2)(s_2-3)} \right)$$

$$+ 3 i {\rm e}^{-\tau} \delta_{\rm D}({\bf p}) \left( \frac{{\rm e... ...frac{{\rm e}^{-2(\tau-\tau_1)}}{(s_2-2)(s_2-3)} \right) \biggl ] \biggl \rbrace$$ (B.3)

where the integration path obeys $\mbox{Re}(s_2)>3$. Note that we removed the Heaviside factor $\theta(\tau-\tau_1)$ from the second term in Eq. (B.3) because the integrals over s₂ and s₃ vanish for $(\tau-\tau_1)<0$ (since for $(\tau-\tau_1)<0$ we can push the integration path to $\mbox{Re}(s) \rightarrow +\infty$). Next, we can apply the operator ${\cal O}$ defined in Eq. (14) onto the kernel  ${\tilde K}$ written in Eq. (B.3). After a tedious calculation we finally obtain Eq. (20). The only technical trick in this calculation is to use Eq. (B.2) as well as the property:

$$s \; \psi(s,x) = 1 - x \; \psi(s-1,x) .$$ (B.4)

This identity is obtained by multiplying Eq. (A.17) by s and integrating by parts. This explicit check of Eq. (20) fully justifies the procedure developed in Sect. 3.

References

 
• Bateman, H. 1953, Higher transcendental functions (McGraw-Hill Book Company, Inc.) In the text  
• Bernardeau, F. 1994, A&A, 291, 697 In the text NASA ADS  
• Bertschinger, E. 1985, ApJS, 58, 39 In the text NASA ADS  
• Bharadwaj, S. 1996, ApJ, 460, 28 In the text NASA ADS  
• Bharadwaj, S. 1996, ApJ, 472, 1 In the text NASA ADS  
• Binney, J., & Tremaine, S. 1987, Galactic dynamics (Princeton University Press) In the text  
• Bouchet, F. R., Juszkiewicz, R., Colombi, S., & Pellat, R. 1992, ApJ, 394, L5 In the text NASA ADS  
• Buchert, T., & Dominguez, A. 1998, A&A, 335, 395 In the text NASA ADS  
• Buchert, T., Dominguez, A., & Perez-Mercader, J. 1999, A&A, 349, 343 In the text NASA ADS  
• Fry, J. N. 1984, ApJ, 279, 499 In the text NASA ADS  
• Goroff, M. H., Grinstein, B., Rey, S.-J., & Wise, M. B. 1986, ApJ, 311, 6 In the text NASA ADS  
• 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 In the text NASA ADS  
• Scoccimarro, R., & Frieman, J. A. 1996, ApJ, 473, 620 In the text NASA ADS  
• Scoccimarro, R. 2000, A new angle on gravitational clustering, Proc. of the 15th Florida Workshop in Nonlinear Astronomy and Physics [astro-ph/0008277] In the text  
• Valageas, P. 2001, A&A, submitted, Paper II [astro-ph/0107126] In the text  
• Valageas, P. 2001, A&A, submitted, Paper III [astro-ph/0107196] In the text  
• Valageas, P. 2001, A&A, submitted, Paper IV [astro-ph/0107333] In the text  
• Valageas, P. 2001, in preparation, Paper V In the text  
• Zinn-Justin, J. 1989, Quantum Field Theory and Critical Phenomena (Clarendon Press, Oxford) In the text