A microscopic, non-equilibrium, statistical field theory for cosmic structure formation

Building upon the recent pioneering work by Mazenko and by Das and Mazenko, we develop a microscopic, non-equilibrium, statistical field theory for initially correlated canonical ensembles of classical microscopic particles obeying Hamiltonian dynamics. Our primary target is cosmic structure formation, where initial Gaussian correlations in phase space are believed to be set by inflation. We give an exact expression for the generating functional of this theory and work out suitable approximations. We specify the initial correlations by a power spectrum and derive general expressions for the correlators of the density and the response field. We derive simple closed expressions for the lowest-order contributions to the nonlinear cosmological power spectrum, valid for arbitrary wave numbers. We further calculate the bispectrum expected in this theory within these approximations and the power spectrum of cosmic density fluctuations to first order in the gravitational interaction, using a recent improvement of the Zel’dovich approximation. We show that, with a modification motivated by the adhesion approximation, the nonlinear growth of the density power spectrum found in numerical simulations of cosmic structure evolution is reproduced well to redshift zero and for arbitrary wave numbers even within first-order perturbation theory. Our results present the first fully analytic calculation of the nonlinear power spectrum of cosmic structures.


Motivation and overview
In a sequence of pioneering papers, Mazenko and Das and Mazenko [1][2][3][4] have recently shown how the nonequilibrium kinetic theory of classical particles can be mapped to the path-integral approach familiar from statistical quantum field theory, in the spirit of [5,6]. Besides the unifying formal analogy, this approach has a multitude of advantages for the systematic development of perturbation theory and the calculation of correlators. Another substantial advantage is that the theory begins with the microscopic degrees of freedom of the individual particles, which usually follow structurally simple equations of motion such as the Hamiltonian equations. Collective fields are introduced as operators extracting the desired information when needed from the microscopic degrees of freedom in the generating functional of the theory.
This paper aims at building upon this approach to find a new access to the theory of cosmological structure formation. Despite heroic efforts and ingenious new developments, it has been notoriously difficult to calculate into the nonlinear regime of second-or higher-order statistics of cosmic structures, such as the power spectrum of the cosmic density field (see  for an inevitably eclectic list, and [29] for an extensive, impressive and complete review).
One of the cardinal difficulties in conventional approaches to cosmological structure formation, both Lagrangian and Eulerian, is that the perturbed quantities are the density and the velocity fields. With dissipationless dark matter dominating cosmic structures, however, particle trajectories can cross and form multiple streams, at which point of the evolution the description by unique and smooth density and velocity fields breaks down. This difference between fluid and particle dynamics seems to be the decisive reason for the difficulty in higher-order, standard cosmological perturbation theory.
Numerical simulations of cosmic structure formation do not encounter this problem because they follow the trajectories of large numbers of individual tracer particles. The density or other collective information is calculated when needed from the actual positions of these tracer particles.
A statistical, non-equilibrium field theory for classical particles can be seen as the analytic analog to a numerical simulation based on particles. As in simulations, following particle trajectories has the decisive advantages that the equations of motion are as simple as possible and that crossing trajectories pose no difficulty at all for the analytic treatment.
With the foundation of the theory worked out in [2,4], the next step towards cosmological structure formation is the definition of a suitable initial particle ensemble in phase space. Working this out is the first purpose of this paper. For rendering our discussion more accessible and self-contained, we shall begin in section 2 by summarising the non-equibrium field theory for classical particles. In section 3, we include operators for collective fields and for the particle interactions. We ask the expert readers for patience, but we believe that this paper is more useful if it contains an outline of the theoretical foundations even though they have already been developed and described in detail elsewhere. In section 4, we construct the free generating functional for particle ensembles initially correlated in phase space. Section 5 discusses first-order perturbation theory in canonical particle ensembles. In section 6, we derive specific expressions for low-order correlators of the density and response fields, which we then specialise in section 7 to derive density power spectra for cosmic structures at first order in the gravitational interaction. This section presents the first fully analytic calculation of the nonlinear power spectrum of cosmic density fluctuations. Section 8 summarises the paper and presents our conclusions.
Even though cosmic structure formation is our main motivation, we believe that both the approach and our central results, the generating functional for correlated classical particle ensembles and approximations to it, may be useful for other areas of statistical physics. We thus intend to lay out the formalism as generally as possible, with little or no reference to cosmology until section 7. To streamline the notation, we use the abbreviations where d is the number of spatial dimensions.

Summary of main concepts and results
This paper is quite technical. To provide a compact overview, we summarise here the concepts, the approximations made and the main results. Essentially, the theory laid out here on the foundations of [1][2][3][4][5][6] begins with an initial phase-space distribution of classical particles following Hamiltonian dynamics. Like in thermodynamics, the statistical properties of this ensemble are characterised by a free generating functional (or partition function) Z 0 given in (31). This generating functional assigns a probability ( ) P q p , for each initial phase-space position ( ) q p , to be occupied by a particle of the ensemble. The phase-space points are then evolved forward in time. A phase factor containing the retarded Green's function of the free Hamiltonian ensures that particles move along their free, classical trajectories.
This free generating functional Z 0 is then extended in two ways. First, in complete analogy to quantum field theory, the particle interactions are written in form of a multiplicative, exponential operator acting on Z 0 . Second, since the full microscopic information contained in the particle ensemble is rarely required, macroscopic or collective fields are introduced as superpositions of microscopic fields. The minimum set of collective fields consists of the number density ρ of the ensemble particles and the so-called response-field B, with the latter describing how the evolution of the particle ensemble responds to changes in the particle coordinates by means of an interaction potential. This leads to the generating functional Z in (58) which contains all interactions and the collective fields required.
So far, the theory is independent of the specific particle ensemble to be studied. With an eye on cosmological structure formation, the initial phase-space probability distribution ( ) P q p , is then constructed to incorporate the appropriate auto-correlations of particle positions and momenta, and the cross-correlations between particle positions and momenta required by continuity. This results in the probability distribution (65), which contains a correlation operator shown in (A.43).
Up to this point, the theory is exact. As (65) shows, the probability distribution contains the momentum auto-correlation matrix of the particles in the argument of an exponential. The dependence of the autocorrelation on the particle positions needs to be integrated out, which is not generally possible analytically. Therefore, we shall expand this exponential up to the second order in the momentum auto-correlations, which is justified in cosmology because the amplitude of these correlations is low.
Similarly, the particle interactions are expressed by an exponential interaction operator. Like in quantum field theory, the series expansion of this interaction operator leads to the Feynman graphs of the theory. We shall expand the interaction operator to first order only, leaving the (tedious, but ultimately inevitable) higher-order perturbation theory to later work.
Thus, apart from the general foundations of the theory, we apply two types of approximation in this paper, viz. the Taylor expansions in the momentum auto-correlations to second order, and in the interaction operator to first order. In addition, but without invoking further approximations, we describe the free particle dynamics in terms of an improved version of the Zel'dovich propagation [30]. Since the improved, free Zel'dovich trajectories already incorporate part of the gravitational interaction, it is interesting to see the nonlinear growth of the density power spectrum possible even before explicitly including the particle interactions. This result, which is of zeroth order in the interaction operator, is found in (182). The contributions to the nonlinear powerspectrum evolution to first order in the interaction are summarised in (190).
2. Non-equilibrium statistical theory for classical fields 2.1. Transition probability for classical fields a be a classical field with n components,   a n 1 , at time t and position  q , living in d space-time dimensions. Further, let the dynamics of this field be described by some equation of motion, written here symbolically as needs to satisfy (2) everywhere in the space-time domain considered. Among all classical field configurations satisfying (2), one particular configuration is singled out by specifying suitable initial conditions, , defined at some initial instant of time t i which we choose to be zero without loss of generality, = t 0 i . The initial field configuration is mapped on a later field configuration by the classical flow ( ) F t cl , Since a classical field evolves deterministically, a field configuration ( ) can be reached at  t 0 beginning with an initial field configuration . In analogy to quantum field theory, we aim to find the probability for the transition of an initial field configuration The meaning of expression (4) is straightforward: a path integration over all possible field configurations j a beginning with ( ) j a i is being carried out, but the functional delta distribution allows only that particular path to contribute which satisfies the equation of motion.
We now introduce a conjugate fieldĉ a to express the delta distribution by a functional Fourier transform where the integration within the exponential function proceeds in general over all d space-time coordinates that the fields j a andĉ a depend on, and a summation over a is implied. Note thatĉ a plays the role of the 'hatted' fieldŷ introduced by [5]. Even though we shall later connote operators with hats, we also add a hat here to emphasise the relation to [5]. The transition probability (4) then reads Identifying the integral in the exponential with the action S and its integrand with the Lagrange density , we define The functional derivative of S with respect to the conjugate fieldĉ a , set to zero, reproduces the equation of motion (2) a a 0 a 2.2. Generating functional for a classical theory A generating functional is now readily constructed from the transition probability P fi . Since the path beginning with a fixed initial field configuration is deterministic in a classical field theory, the only possible random element in such a theory is the configuration of the initial states. The configuration space to be summed or integrated over in the construction of the generating functional is thus the space of initial field configurations. Therefore, we integrate over all possible configurations of initial states, weighted by an initial probability distribution [ ] ( ) j P a 0 i . We shall abbreviate the path integral over the initial field configurations by Later, when we shall specify classical microscopic degrees of freedom for the fields, the initial states will be defined by a point set rather than by a set of functions. Then, the probability distribution P 0 for the initial conditions will be a function rather than a functional, and the path integrations over the initial states will turn into ordinary integrations. Finally, we introduce auxiliary source fields J a for j a and K a forĉ a into the Lagrangian and thus arrive at the generating functional   if Z is normalised, [ ] = Z 0, 0 1. Since we have obtained Z by integration over a functional delta distribution, normalisation is ensured. Likewise, as in quantum field theory, the functional = W Z ln is the generating functional for the connected correlators, i.e.the cumulants.

Generating functional for the non-interacting theory
Suppose now that the equation of motion can be brought into the form where E 0 represents the free motion while E I is due to any interaction. We can then split the action and the Lagrangian into a free partˆ[˙( x a a a 0 0 0 0 and an interacting partˆ( 1 5 x a a I I I I We shall proceed with the free part first, ignoring for now any interactions between the fields j a themselves or between the fields j a and any external field. Restricting the action to its free part S 0 , we obtain the generating functional for the free theory from (10), e x pi   exp i  .  17   a  a  a  a  x  a  a  a  a  a a   a  a  a  a  x  a a   0  i  0   i  D  0 For any given initial field configuration ( ) j a i , the delta distribution in (17) singles out the solution¯( ) j x a of the free equation of motion, augmented by the inhomogeneous source term K a . Let ab be the propagator (Green's function) of the free equation of motion, then this solution is Absorbing a constant functional determinant into the normalisation of the generating functional, we can replace the delta distribution in (17) by and write the free generating functional as where j a was replaced by the free solutionj a from (18) by integrating over the delta distribution.

Microscopic and collective fields
3.1. Introductory remarks So far, the formalism chosen for classical fields is independent of the specific equations of motion and of the general properties of the fields. For the following discussion, the distinction between macroscopic and microscopic fields or degrees of freedom will be important. Instead of a macroscopic field such as an electromagnetic field, we can also use the formalism for describing the kinematics of point particles under the influence of Hamiltonian dynamics in three spatial dimensions. Then, delta distributions at the phase-space coordinates replace the fields j a . The equations of motion of the phasespace points  x j are Hamilton's equations where  is the Hamiltonian,  is the symplectic matrix 3 and the derivative ¶ j acts upon all six phase-space coordinates x j of the jth particle. The matrix  d is the unit matrix in d dimensions. With microscopic degrees of freedom, the action S in (7) simplifies to a time integral, as in classical mechanics. Similarly, the Green's function will then also depend on time only.

Data structure
For notational as well as conceptual simplicity, we follow [31] and organise the positions { }  q j and the momenta { }  p j of N microscopic particles by means of the tensor product into the phase-space coordinate tensors where summation over repeated indices is implied, and  e j is the N-dimensional column vector whose only nonvanishing entry is 1 at component j. Recall that the tensor product has the convenient properties Tr Tr . 24 We further introduce the scalar product where the sum over the repeated indices is again implied. Bundling the phase-space points accordingly Like the phase-space coordinates, we bundle the source fields J and K as We then need to introduce an analogous tensor product for the propagators where G is a 6×6 dimensional matrix describing the free propagation of an individual phase-space point.
With this notation, we can write the free solution (18) as for all particles together, and the free generating functional assumes the form Note that the integral over the initial phase-space configurations is now an ordinary rather than a path integral.

Collective fields
If the fields j a represent microscopic degrees of freedom, such as the phase-space coordinates of point particles, it will be appropriate to introduce collective fields in addition, i.e.fields representing collective properties of the particle ensemble. Perhaps the most obvious example of such a collective field is the density ( ) here assumed to be composed of N point-particle contributions.
The potential ( )  V t q , experienced by any particle at time t and at position  q is the sum over all pointparticle potentials v, which we can re-write in terms of an integral over the density (32), According to Hamilton's equations, and given the interaction potential ( )  V t q , , the interaction contribution to the equations of motion of the particle ensemble is where the last step was taken by partial integration to remove the gradient from the potential for later convenience. With this result and with the conjugate fieldˆ≔ˆ c c Ä e j j , we can thus write the interaction contribution to the Lagrange density  I from (15) as The term in brackets defines the response field in terms of the conjugate 'hatted' fieldĉ. Introducing this into (37), we can write the interaction part of the Lagrange density as Expressing the response field B, the potential v and the density ρ by their Fourier transforms, we can re-write the interaction Lagrangian as where we have assumed that the potential v in (39) is translation invariant and thus depends on the difference   q y only. We now combine the two collective fields ( ) and write the interaction part of the action in the compact form where the conventional abbreviations = t k d1 d d 1 3 where the delta distribution ( where the product is understood as an implicit sum over the collective-field indices and integral over the spacetime coordinates

Operator expressions for the collective fields
The collective fields Φ typically contain the field variables x orĉ. These are obtained from the free functional [ ] J K Z , 0 by functional derivatives with respect to the sources J or K , For introducing the values of x andĉ into Φ, we replace Φ by an operatorF acting on the free functional , with all occurrences of x andĉ replaced by functional derivatives according to (46). Then, the free generating functional including the collective fields is expressed by The minimum set of collective fields that we require are the density ρ and the response field B. We shall now construct their operator expressions. The density ρ is assumed to be composed of delta contributions, see (32). In Fourier space, the one-particle contribution of particle j to the density at the space-time position ( ) where ( )  q t j 1 is the position of particle number j in configuration space at time t 1 . In this expression for the density, we replace the particle position  q j by a functional derivative with respect to ( )  J 1 q j , obtaining the oneparticle density operatorˆ( The action of the density operator (49) on the free generating functional (31) becomes clear by expanding the exponential into a series. This immediately leads tô Thus, the application of the density operatorˆ( ) F r 1 j amounts to a shift of the source field J in the free generating functional by the tensor ( ) L 1 j . According to (38), the one-particle contribution ( ) B 1 j of particle j to the collective field B(1) is determined by the gradient of the density Taken into Fourier space, the one-particle response-field operator thus turns intô SinceF r j involves functional derivatives with respect to  J q j while the response-field operator takes functional derivatives with respect to  K p j , the relevant functional derivatives commute. We can thus reorder the operators and apply all required response-field operators after all density operators.

Operator expression for the interaction part of the action
We have seen in (42) and (43) that the interaction between the particles can be included by adding the expression to the free action. By means of the operator expressions (49) and (53) for the two collective fields ρ and B, we can write this action contribution in the operator formˆ( where the operatorsF r andF B will now be responsible for acquiring the respective collective-field values from the free generating functional . Therefore, the interaction part of the action iŝ as noted in (45) before in more general form.

Generating functional for correlated initial conditions
Having developed the general formalism for a set of Hamiltonian point particles, the final step to be taken towards defining the generating functional (58) completely concerns the initial phase-space distribution. In this section, we consider points in a domain of phase space at an initial time only. For convenience, we shall drop the superscript (i) on any quantity here, understanding that all quantities are to be taken at the initial time throughout this section.

Initial phase-space distribution
So far, the construction of the non-equilibrium field theory for classical particles has been completely general, with the one exception that we required the microscopic degrees of freedom to follow the Hamiltonian equations of motion. To specify the generating functional completely, we now have to define the initial phasespace measure that is, we have to construct the probability distribution ( ) q p P , for initial particle positions in phase space. Having cosmological structure formation in mind, we need the particles to be spatially correlated such that their number density is a homogeneous and isotropic Gaussian random field. By continuity, spatial correlations imply correlations also in momentum space as well as cross-correlations between spatial and momentum coordinates. Our main goal here is thus to derive the probability distribution for the initial phase-space coordinates under these requirements.
A central (and, as we shall see, a sufficient) quantity characterising all required correlations is the power spectrum of density fluctuations. Calling the number density of particles ρ and its meanr, the density contrast is and its power spectrum ( ) d P k is defined as where the density contrast ( )  d k as a function of the wave vector  k implicitly denotes the Fourier transform. The delta distribution ensures translation invariance and thus the statistical homogeneity of the density contrast. If it is statistically isotropic as well, the power spectrum depends on the wave number k only and not on the direction of the wave vector  k . The cosmological motivation aside for now, we are thus aiming to derive the probability distribution for correlated phase-space points drawn from a statistically homogeneous and isotropic Gaussian random field characterised by the power spectrum of the density fluctuations. Such initial conditions may be interesting far beyond cosmology.
For clarity of the discussion in the main part of this paper, we shall develop the probability distribution ( ) q p P , in the appendix. The central variables in the derivation of this probability distribution will be the values d j of the density contrast and  p j of the momentum at the positions of all particles j. We organise these variables at all N particle positions into a data tensor by means of the tensor product with the vectors  e j defined in (23). A major intermediate result will be the covariance matrix¯≔ of this data tensor, which contains the density-contrast and momentum auto-correlations d d á ñ j k and   á Ä ñ p p j k , respectively, and the density-momentum cross-correlation  d á ñ p j k . The entries of the covariance matrix are detailed in the appendix. As derived there, the initial phase-space probability distribution is with q and p as defined in (23). The correlation operator ( )  p appearing here is given in (A.43). It contains the correlation matrices dd C and d C p introduced above and defined in (A.31). If the correlations dd C and d C p are weak, as we can expect them to be early in time, the probability distribution (65) can be approximated by

General expressions for density correlators
We are aiming at calculating m-point correlators of collective fields, such as the density and response fields. Equation (53) shows that the response-field operator contains the density operator as a factor. Thus, for an mpoint correlator, m density operators will have to be applied to the free generating functional first. Since no further derivatives with respect to J will be required afterwards, the source field J can then be set to zero. The operator for the density contributions by N particles is the sum over the one-particle density operatorŝ As illustrated in (51) for a single one-particle density operator, the result of applying m one-particle density operators to the free generating functional iŝ Following (53), a single one-particle response-field operatorˆ( ) b l j l applied subsequently then leads tô According to (68) and (69), an m-point density correlator is found by summing over all particle indices This shows that all we have to evaluate for m-point correlators of the density and response fields is the free generating functional taken at = where the free phase-space trajectory¯( appears because the term containing the source K disappears. With (70), the phase in the exponential in (73) is for the free generating functional evaluated at = J Land = K 0. The shift tensors L q p , have non-vanishing components only for the particles specified by the indices j s set by the one-particle density operators applied to the free generating functional. For any shift tensor specified by a complete set of The integral over the initial phase-space configuration remaining in the free generating functional [ ] L Z , 0 0 still has to be carried out.

Integration over the initial phase-space distribution
Since the first and the third term in parentheses in the second line do not depend on the momenta p, they can easily be integrated over p using The second term in parentheses can be integrated after pulling a factor p down by applying the derivative (67), this leads to Combining the results (79), (80) and (82), we find where we have used (A.31) to split up the quadratic formQ into¯= , and . 84 For later convenience, we introduce one-particle shift vectors  L q j and  L p j by the projection In terms of  L p j , we can write the quadratic forms in (84) as with the quadratic form Q D from (83) or (87) appearing in the free generating functional (83) requires a separate consideration. As the derivation of [ ] J K Z , 0 shows, it originates from the initial one-point momentum variance and thus arises from the free streaming of the particles with the initial root mean square velocity quantified by s 1 . In absence of momentum correlations, it would lead to a Maxwellian or thermal velocity distribution of the particles.
In the cold-dark matter cosmogony, free streaming is suppressed by the long-ranged gravitational interaction between the massive particles. In the free generating functional, gravitational interaction is not included by definition. Later, we shall introduce gravitational interaction between the particles in a perturbative approach. As long as we neglect gravitational two-particle interaction, the damping expressed by ( ) -Q exp 2 D will be unrealistic for cold dark matter because it is counteracted by the gravitational interaction.
The effect of the damping term depends on its relation to the quadratic form Q, also defined in (84), which contains the initial momentum correlations between different particles. As we shall see later, these initial momentum correlations will be mainly responsible for the growth of structures. Realising that Q will commonly be a small number, we shall approximate i.e. we shall expand in powers of the initial momentum correlations. For appropriately suppressing the damping term relative to the growth of structures, we shall approximate the damping term ( ) -Q exp 2 D consistently at one order less than the term (88). This implies that damping will only be included at loop order, but not at tree order. While this may appear arbitrary here, we shall show in a follow-up paper how the damping term is counteracted when the complete hierarchy of momentum auto-correlations is taken into account. Since this calculation is quite involved, we postpone it here.
Thus, when we derive results from the free generating functional restricted to linear initial momentum correlations, we shall ignore the damping term completely, approximating At the next higher order of the initial momentum correlations, we shall include the damping term at linear order, approximating where the second approximation is advantageous because it remains positive definite.

One-point functional and normalisation
If we consider deriving free one-point 'correlators' of a collective field, e.g. of the density, from the generating functional, we can ignore all correlation terms because they appear only if two or more points are involved. Then, Q=0, and the generating functional (83) shrinks to If a single point is involved, L q will have a single non-vanishing component, which we can without loss of generality label with the index j=1. Then , . 92 because the delta distribution further ensures that  = L 0 q 1 and thus also Q D =0. The factor N in (93) takes into account that there are N possibilities to select a particle from the ensemble. Since the remaining delta distribution is the Fourier transform of unity we see that (93) simply reproduces the mean particle density, as it should.

Low-order approximations
We can now Taylor-expand the factore Q 2 in (83) in powers of Q, for example to first or second order. This gives the two contributions . We shall now evaluate these expressions in detail.
The first integral in L q i , q has a particular meaning. It simply returns a delta distribution or products thereof, depending on L q . Since it does not incorporate any correlations, it can only represent shot-noise terms or powers of the mean particle densityr. For a general discussion of shot-noise terms, see section 5 below. Besides, all other terms involve double sums over = ¼ Note that the ( ) Z jk 1 are not necessarily symmetric in ( j, k) because of the d C p j k correlation between densities and momenta. In contrast, the ( ) Z jklm 2 are symmetric under the permutations (

Linear momentum correlations
With the explicit expressions (A.21)-(A.23) for the components of d d C j k , d C p j k and C p p j k , and using that the power spectra for the density and the velocity potential are related by (A.6) we immediately obtain from (99) the result where the abbreviations were defined. The prime on  ¢ jk indicates that  q j and  q k are excluded from q here.

Quadratic momentum correlations
For taking quadratic initial momentum-momentum correlations into account, we need to evaluate different terms contributing to ( ) Z jklm 2 in (100). In view of the symmetries of ( ) Z jklm 2 , we find it convenient to distinguish terms with two equal index pairs, ( ) Z jkjk 2 , terms with one double index, ( ) Z jkkl 2 , and terms with four different indices, ( ) Z jklm 2 .
For two equal index pairs, we find and for four different indices Given those expressions, the free generating functional is approximated by (95) and (96). Further progress can be made once the shift tensor L is specified, for example when the density correlators are to be calculated; see section 6 below. with the interaction part of the action given by the operator

First
defined with slightly more explicit notation in (57). As we have noted before, this expression for the interaction part of the action contains the two assumptions that the potential is assumed to be translation invariant and acts instantaneously.
Since the functional derivatives with respect to H in (110) act only on the collective-field operator ·F e H i , the effect of the interaction operator can be brought into the form The density and response-field operators,F r andF B , in the interaction part S I of the action now act directly on the free generating functional and produce correlators introduced in (72) before. To lowest non-trivial order, the interaction operator isˆ( The corrections to the one-and two-point density correlators in first non-trivial order are then Note that we mark with primes the internal vertices of the interaction, which are integrated over in the interaction operator.
For calculating the first-order approximation of the nonlinear density evolution and the nonlinear power spectrum, we thus have to work out the three-and four-point correlators from the free generating functional.
In view of our later cosmological application, we anticipate that the potential satisfies a Poisson equation of the form with a function g v (t) to be specified, where δ is again the number-density contrast of the particles. Since we are here aiming at the potential caused by a single particle, we can place this particle without loss of generality into the origin of a coordinate system and write its contribution to the density contrast as 1 D withr being the mean particle number density. Fourier transforming (116) then gives The Fourier-transformed unity1 can be neglected later because the zero mode of the potential will not contribute to any correlators. We can thus insert for the Fourier-transformed, one-particle potential v. Notice in particular that this potential scales inversely with the mean particle densityr. This is because, for a fixed mean mass per volume, the particle mass has to decrease in inverse proportion to the particle number N if that number is increased.

Shot noise and the relevance of terms
In our microscopic approach, shot-noise terms appear because the density field is composed of discrete particles. To identify these terms and to clarify their relevance, consider a statistically homogeneous density field In terms of the density contrast δ, the power spectrum of a continuous density field is On the other hand, calculating the variance of (121) results in Obviously, only the second term in (124) corresponds to the result (123) for the continuous density field, while the first arises only because the density field is composed of discrete particles. Thus, the first term in (124) is a shot-noise term which arises from summing over pairs of identical particles, as the calculation shows.
More generally, for m-point correlators of density fields composed of discrete particles, an analogous calculation shows that terms proportional to all powers ofr occur,r s , with   s m 1 . Only the term proportional tor m is not a shot-noise term. It is the only term arising from summing over combinations of particles which are all different. Terms proportional to powers ofr s with < s m are all shot-noise terms in the sense that they arise because of the discrete nature of the density field. In the thermodynamic limit  ¥ N , the shot-noise terms can be neglected relative to the dominant term proportional tor m .
In the case of gravitational interaction between the microscopic particles, the interaction potential scales with the particle mass. Resolving the density field into an increasing number of particles while keeping the mass density constant, the particle mass must be decreased proportional to -N 1 . This repeats the argument made following (119): the Poisson equation then implies that the gravitational interaction potential must scale inversely with the mean number density of particles, i.e. liker -1 .
According to (112), the interaction operator from the interaction part S I of the action increases the order of the density ρ and the response field B in the free correlators by one each and multiplies with a potential. As (71) shows, the response field identifies two particles, as expressed by the Kronecker symbol d j j m s there. Comparing this with our earlier result on the origin of shot-noise terms, we see that the identification of particles by the response field only selects shot-noise terms from the free density correlators because the only non-shot noise term in the free density correlators arises from combinations of different particles, for which d = 0 j j m s . Specifically, for an m-point density correlator in nth order perturbation theory, free correlators of order up to + m n 2 need to be calculated which are of ( ) + m n th order in the density and nth order in the response field. In these free correlators, terms proportional to all powers ofr up tor + m n 2 will occur. Their subsequent multiplication by v n will reduce the power ofr by n tor + m n . Each response field will identify particles pairwise and will thus further reduce the power of the leading term tor m , as expected for an m-point density correlator.
This shows that only such terms in the free correlators of order + m n 2 need to be considered which are proportional tor + m n . Terms proportional to lower powers ofr will vanish in the limit  N 1, while terms proportional to higher powers ofr disappear because of the identification of particles by the response fields. If we wish to calculate density correlators taking momentum correlations into account to first or second order, we need to evaluate the expressions ( ) Z jk 1 from (102)  . Then , 126 If we set = t t taking into account that the remaining delta distribution ensures  As discussed below (100), the terms ( ) Z jklm 2 are symmetric in the first and second index pairs and under exchanges of the two index pairs, and the indices in the first and in the second index pairs must be different. Under these requirements, the term (132) appears 4 times in the sum over particle indices: terms with the index combinations ( ) 1212 , ( ) 2112 , ( ) 1221 and ( ) 2121 are all equivalent, and others do not appear. Furthermore, we have to multiply with the number ( ) -» N N N 1 2 of ways for selecting a pair from the N particles. Thus, by (98), we arrive at the contribution of quadratic momentum correlations to the two-point density correlator. The damping term is set to unity here as discussed in section 4.4 before.

One-and two-point response-field correlators
The effect of a single, one-particle response-field operator on the free generating functional was shown in (53). That expression, valid for a single response-field operator applied to the free generating functional, is easily generalised. Suppose we apply m operators in total, of which n are density and m−n are response-field operators. Since each response-field operator contains a density operator to be executed first, we will have to apply m density operators in total. The result will have to be multiplied by m−n response-field factors. Thus, we haveˆ( Applying a single response-field operator to the generating functional, we obtain the mean response field. In the general approach outlined above, we set m=1 and n=0. Then, from (135), we have qp vanishes for = ¢ t t , as it usually will. Then, the mean response field vanishes identically.
For m=2, we have the density-response correlator after summing over all   j N 1 1 . Changing the order of B and ρ in (137) only changes the ordering of the times t 1 and t 2 , thus leading to and multiply the results with the number of ways to choose particle pairs and particle triples from an ensemble of N particles.
We begin with the correlators derived from the generating functional [ ] ( ) L Z , 0 0 1 from (96), which contains momentum correlations to linear order only. For m=3, the one-particle response-field factor in (135) reduces to the single term , 0 qp qp 1 2 2 2 . Moreover, we have replaced  ¢ k 2 by  -¢ k 1 , expressing the translation invariance of the potential v. Since the Kronecker symbol in the response-field factor identifies the particles j 1 and ¢ j 2 , only two particle indices are free, which we set without loss of generality to For the two-particle term (102), we can label the two particles by ( ) ( ) = j k , 1,2 and thus write . 145 q q 1 1 1 1 2 We can stop here: the delta distribution in the two-particle term in (102) shrinks to to first order in the interaction and to linear order in the momentum correlations: to this order, the interaction does not change the mean density. For the two-particle term (106) contributing to the quadratic momentum correlation, we can also set ( ) ( ) = j k , 1,2 and arrive at the same conclusion: the delta distribution ensures  = k 0 1 and thus sets the response to zero. The three-particle term (107) cannot contribute because  = L 0  in (96), we first notice that again neither  L q 1 nor  L q 2 must vanish because otherwise individual wave vectors would be set to zero, causing the power spectrum to disappear. Therefore, at least one each of the particle indices ( ) j j j , , 1 2 3 must be set to 1 and 2, while the third particle index available for the threepoint correlator remains free.
If we set this third index to 1 or 2 as well, the multiplicity of the resulting term is ∝N 2 , which is lower by a factor of N than the multiplicity ∝N 3 required for the three-point correlator. This term is thus negligible. Only terms with the third index set to >2 will remain. Adopting ( Moreover, for a synchronous three-point correlator, = = t t t 1 2 3 . Thus where the latter step follows because one of the remaining delta distributions ensures The index combination ( ) ( ) = j j j , , 2, 1, 3 1 2 3 adds the same expression. Taking the remaining cyclic index permutations into account leads to the contribution to the three-point correlator from linear momentum correlations. The terms of second order in the momentum correlation can be read off (106) and (107). The two-particle . 152 taken from (133). Finally, the three-particle term ( ) Z jkkl 2 from (107) gives to the three-point correlator. Since there are no interactions included at this level, we set the damping factor to unity.

Four-point correlator
from linear momentum correlations Turning to the effect of first-order interactions on the density power spectrum, we need to work out the fourpoint correlator . The response-field factor is again. Other terms do not appear here because 0 , qp qp 1 2 2 2 . We shall further consider synchronous correlations only and thus set = t t 1 2 . Of the two terms remaining in (156), we now focus on the first, in which the Kronecker symbol ensures that = ¢ j j 1 2 . The second term will then be obtained from the result by interchanging the indices j 1 and j 2 or, equivalently, the wave vectors without loss of generality. The shift vectors are then The three particles need to be placed on three different positions to achieve the largest possible multiplicity. We choose three positions labelled by ( ) ( ) = j k l , , 2, 3, 1 , obtain the shift vectors (157) and A k  2  2  159   jk  23   1  3  D 2  1  3  D 1  1  2  2  2 from (102). The second delta distribution arises from the factor  ¢ jk in (105). Since it ensures   ¢ = k k 1 1 , it allows us to write (159) as 2 due to the delta distributions. For permutations of ( ) j k l , , with ¹ l 1, the factor  ¢ jk results in a delta distribution setting one individual wave vector to zero, which causes the result to vanish. The only other permutation leading to a non-vanishing result is thus ( ) ( ) = j k l , , 3, 2, 1 , for which After collecting results, the summation over particle indices multiplies the result by ( ) 3 , and the relevant three-particle contribution to the four-point density correlator turns out to be . It is quite straightforward to see that the contribution for = ¢ j j 2 2 is identical, multiplying the correlator by two. Thus, the four-point correlator required for the firstorder perturbation theory according to (115) is According to (115), this implies the contribution  Since the response-field prefactor in (156) identifies particle pairs, only three particle indices are free, which immediately implies that no four-particle terms can contribute. The two-and three-particle terms from (106) and (107) are thus the only ones to consider. Again, we label the particles by ( ) ( ) ¢ = j j j , , 1, 2, 3 1 2 1 without loss of generality. Regarding the three-particle term ( ) Z jkkl 2 , the position-index combination ( ) ( ) = j k l , , 1, 2, 3 leads to , 169 . 170 Finally, for the two-particle term in (106) to contribute, the factor  ¢ jk returns a delta distribution for an individual wave number except for the particle-index combinations ( ) ( ) = j k l , , 2, 3, 1 or ( ) 3, 2, 1 . For these follows again by multiplying with the response-field factor (156), taking into account that both terms lead to same result. Thus Inserting this into (115), we find 7. Cosmological power spectra 7.1. Power-spectra contributions from the free generating functional All results obtained so far for the generating functional, for the initially correlated phase-space distribution and the low-order density-and response-field correlators are generally valid for systems of classical particles. The free phase-space trajectories of these particles are characterised by a known retarded Green's function and they interact with a two-particle potential v.
In this section, we shall specialise these results to classical point particles in cosmology. The essential difference to common classical particle systems is that space is expanding with time. The physical distance  r between any two particles thus grows with time in proportion to a scale factor a(t). The spatial coordinates  q are taken to be comoving coordinates, defined by ( )   = r a t q. The Green's function for particles moving freely in such an expanding space has been derived in [30]. Specifically, the free propagator has been shown to be 1 was introduced as a time coordinate more convenient than the cosmological time t or the cosmological scale factor a. The function ( ) + D a is the linear growth factor, describing the increase in densityfluctuation amplitudes as long as they remain linear. The growth factor is assumed to be normalised to unity at the initial time such that t = 0 initially. Moreover, ( ) t g is defined by into the free linear power-spectrum term (130) cannot reproduce the result well-known from ordinary cosmological perturbation theory that the matter power spectrum evolves linearly as ( ) We can, however, achieve this behaviour by mimicking the Zel'dovich approximation, which implies extrapolating the first-order solution of Lagrangian perturbation theory (LPT) beyond the linear regime. LPT describes the motion of fluid elements in terms of a displacement field ( ) ( ) i 0 which maps the initial Lagrangian coordinate ( )  q 0 of any fluid element to its final position  q at a later time t. (See also appendix for the notation.) Applying this map to the evolution of fluid mass elements leads to the continuity equation , assuming that the initial density field is nearly uniform, ( )  d 1 i . Using the Jacobian determinant of this mapping and linearising, one finds the first-order relation , , . 178 Together with the equation of motion for  Y [32] and assuming an irrotational flow, this equation is solved by Thus, in this approximation, particles simply move on straight trajectories with the time dependence given by the linear growth factor. The Zel'dovich approximation now lies in extrapolating these trajectories to the present day. In our approach, this entails replacing the Hamiltonian propagator (175) with a Zel'dovich propagator qp Z Then, the time-evolution factor in the free two-point density cumulant (130) derived from linear initial correlations turns into Consequently, the free linear power-spectrum contribution scales as , as expected from Eulerian standard perturbation theory (SPT). This time evolution is due to the fact that the equation of motion for has the same form as that for the linear density contrast in SPT. Since this equation of motion contains the gravitational potential, the Zel'dovich trajectories already include part of the gravitational interaction between particles. This interaction and the actual deviations from inertial motion it causes are hidden in the time dependence of straight Zel'dovich trajectories.
To first and second order in the initial correlations, with the Hamiltonian propagator replaced by the Zel'dovich propagator and suitably approximated damping terms, we thus find the following contributions to the power spectrum in our free theory where we have specified for clarity that the power spectra on the right-hand sides are the power spectra ( ) d P i characterising the initial particle distribution.
As one would expect, these expressions are also found in linear LPT by going to quadratic order in the initial correlations of ( )  Y 1 (see the ( ) C ij 11 terms in equation (35) of [32]). However, one should be aware of the fact that the LPT formalism only includes pure initial momentum correlations due to the assumption of a uniform initial density in (177). In our approach, quadratic power-spectrum contributions coming from density autocorrelations and density-momentum cross-correlations are also present, as they should be. As mentioned before, we dropped them here since they scale with lower powers of the propagator g qp . The assumption of a uniform initial density field is also responsible for the slightly different time dependence of LPT when compared with our approach.
Strictly speaking, the choice of LPT in (179) is inconsistent with this assumption and the boundary condition ( ) if the linear growth is normalised as ( ) = + D 0 1at the initial time t=0. However, it is necessary in LPT to achieve the same growth of the linear power spectrum with time as in Eulerian SPT, since LPT lacks the density auto-correlation and the density-momentum cross-correlations which together lead to the correct time evolution factor in (130) Figure 2 shows these four terms A, B, C and D together with the linearly evolved density power spectrum. Dashed curves indicate negative contributions. Clearly, term A is positive throughout and expresses how structures grow on small scales by gravitational collapse. Term C is negative. At large scales, its amplitude is larger than that of term A, while it falls below at small scales. This reflects two important aspects of cosmological structure formation. Gravitational contraction removes power from large scales and transports it to smaller scales. However, the reduction of power on small scales by term C is exaggerated here because the improved Zel'dovich propagators still overshoot and are only partly compensated by the first-order interaction. The power on small scales is thus suppressed too strongly. The terms B and D are much lower in amplitude, except on the smallest scales. There, however, the negative contribution by term B almost exactly cancels the positive contribution by term D. While term D contributes to structure growth on small scales, term B adds power on intermediate scales, but removes power on small scales, albeit at a low level.
The combined effect of terms A-D is shown in figure 3, together with the sum of terms A-D and the linearly evolved power spectrum. For comparison, the approximation by the Coyote cosmic emulator [34] of the nonlinearly evolved density power spectrum obtained in fully numerical simulations is also shown. Clearly, at first order perturbation theory, our analytic result falls below the numerical result. The reduction of power on large scales and the increase on small scales is clearly shown. Yet, the overshooting due to the improved second-order perturbation theory may already yield fully satisfactory results with the Newtonian interaction potential (189). This will be worked out in a follow-up study.

Conclusions
Mazenko and Das, Mazenko have recently described how the evolution of classical point-particle ensembles can be described as a non-equilibrium statistical field theory [2,4]. The theory uses the phase-space coordinates of the N particles in the ensemble as elementary microscopic fields. The central object of the theory is a generating functional. This is a path integral over phase-space trajectories, weighed by a probability distribution for the initial particle positions, in which each particle is represented by a phase factor containing the free particle trajectory.  To reduce the dynamic range, we normalise the linearly evolved power spectrum, the power spectrum evolved with the cosmic emulator and our first-order analytic solutions with and without viscosity to the linear spectrum according to [31], which has no BAOs. The spectra are are evolved to s = 0.8 For expressing the microscopic properties of the N-particle ensemble in a compact way, we have introduced a notation bundling these properties in tensor-valued structures which appear straightforward to calculate with (see also [31]).
Collective fields, such as the density and the response field, are described by operators extracting correlators of the collective fields from the generating functional by functional derivatives. Likewise, interactions between particles are described by an interaction operator containing the collective density and response fields. This theory closely resembles statistical quantum field theory. Since the form for the equation of motion of the microscopic degrees of freedom is assumed to be very general, the theory should be widely applicable to ensembles of classical point particles.
Using the theory of cosmic structure formation as a motivation, we have derived in section 4 a probability distribution for the initial phase-space positions of the particles which accounts for auto-correlations of the spatial positions and the momenta, as well as for cross-correlations between spatial positions and momenta.
We argue that, under very general conditions, continuity requires the momenta to be correlated if the positions are. If the particles are supposed to sample an initial Gaussian density field, a single initial power spectrum, e.g. for the initial density field, suffices to specify the initial probability distribution. We derive its exact form, which contains a correlation operator containing the complete correlation hierarchy of the ensemble. We then approximate this correlation operator to low order in the correlations and give explicit expressions for the generating functional containing momentum correlations to first and second order.
The main results of this derivation in section 4 are the exact equation (65) for the initial probability distribution in phase space, the contributions (102) to the generating functional with linear momentum correlations, and the terms (102), (107) and (108) for the generating functional with quadratic momentum correlations. From these results, low-order density-and response-field correlators for correlated ensembles of classical particles can be readily determined, as shown in section 6.
Based on this general formalism of the theory and on a free generating functional for initially correlated, canonical particle ensembles, we have derived two-and three-point correlators of the cosmic density field without interaction, and the two-point density correlator of the cosmic density field with interaction included in first-order perturbation theory. Our results of section 7 can be summarised as follows: • If momentum correlations are taken into account to linear order, and if the free particle propagator is taken from the Zel'dovich approximation, the density power spectrum (130) reproduces the linear growth well known from SPT.
• Evolving quadratic momentum correlations with the free Zel'dovich propagator leads to a first contribution to the nonlinear evolution of the power spectrum, for which the simple, closed expression (133) can be given. This contribution is a convolution of the initial power spectrum with itself, multiplied by a mode-coupling kernel.
• Deriving the bispectrum, we obtain the connected term (183), containing the kernel (184). In this form, it resembles the bispectrum result from Eulerian perturbation theory, but with a small difference in two coefficients. The reasons for this difference are subtle and will be explained in future work.
• Proceeding to first-order perturbation theory, using the improved Zel'dovich propagator and the appropriate interaction potential derived in [30], we showed that the shape of the first-order nonlinear terms reproduce the shape of the nonlinear density power spectrum known from numerical simulations rather well, while the amplitude at large wave numbers turns out to be substantially too low. This reflects the fact that, in the improved Zel'dovich approximation, the re-expansion of cosmic structures is still not fully suppressed by the first-order interaction.
• While this calls for higher-order perturbation theory, an effective remedy is provided even at first order by an adapted version of the adhesion approximation [35][36][37]. If we modify the interaction potential accordingly, the shape as well as the amplitude of the nonlinear corrections to the power spectrum reproduce the numerical results very well.
Our calculation extends to redshift zero and to arbitrary wave numbers. Apart from the viscosity introduced to strengthen gravity in our first-order calculation, our theory has no free parameters once the power spectrum of the initial phase-space particle distribution is fixed and normalised. We expect that, once we can proceed to second-order perturbations, the theory will have no free parameters. The form of the nonlinear terms and the inevitable damping factor suggest that the expected asymptotic behaviour of the nonlinear power spectrum for large wave numbers will be retained in higher-order calculations.
The main difference to conventional, Eulerian or LPT of cosmic-structure evolution is that we do not require, solve or perturb a dynamical equation for the cosmic density. Rather, we study the statistical evolution of a particle ensemble in phase space, weakly perturbing their trajectories, and read out any collective information, such as the density, from the evolved phase-space distribution when needed. Since even small perturbations of trajectories can lead to large increases in density, our approach is able to extend into the regime of highly nonlinear density perturbations even at low perturbative orders. It also appears crucial to keep the complete phase-space information of the particles because this allows us to use the Hamiltonian equations of motion with their simple structure and their equally simple Green's function.
While our main motivation for this work has been the extension of this theory to cosmic structure formation from dark-matter particles, the results given here may be useful for a wider class of problems involving the nonequilibrium statistics of correlated ensembles of classical particles. implying that the initial density contrast needs to satisfy the Poisson equation The power spectra ( ) y P k for the velocity potential and ( ) d P k for the density contrast must thus be related by At the same time, the initial particle velocity iṡ˙( Finally choosing the time coordinate such that˙( ) = = b t 0 1and setting the particle mass to unity, we can identify the gradient of the velocity potential with the initial momentum The initial density contrast and the initial momentum are thus related by the velocity potential. For the remainder of appendix, we shall drop the superscript (i) for the initial positions and momenta, understanding that all positions and momenta are to be taken at the initial time for now.
A.2. Probability for particle positions in phase space We aim at the probability for finding a point particle j at position  q j with momentum  p j . The point particles need to sample the density ρ, hence the probability for finding a particle at position  q j , given the density ( )  r r = q j j , must be proportional to r j , ( | ) ( )  r r = -P q N , A . 9 j j j 1 with the normalisation -N 1 set by the requirement ò r = V N d j .
Expressing the initial density ρ by the density contrast δ and introducing the velocity potential ψ, we can substitute¯( ) r r y = - 1 2 for the density. Note that ρ must be chosen to be a number density here since we intend to draw particles from it.
In a similar manner, we need to account for the conditional probability for a particle at position  q j to have momentum  p j . By (A.8), we have understanding that ( )  y y  =  q j j . For ease of notation, we now abbreviate the negative Laplacian and the gradient of the potential ψ at position  q j by The probability for finding a particle at position  q j with momentum  p j can thus be related by