Kinetic Field Theory: Effects of momentum correlations on the cosmic density-fluctuation power spectrum

In earlier work, we have developed a Kinetic Field Theory (KFT) for cosmological structure formation and showed that the non-linear density-fluctuation power spectrum known from numerical simulations can be reproduced quite well even if particle interactions are taken into account to first order only. Besides approximating gravitational interactions, we had to truncate the initial correlation hierarchy of particle momenta at the second order. Here, we substantially simplify KFT. We show that its central object, the free generating functional, can be factorized, taking the full hierarchy of momentum correlations into account. The factors appearing in the generating functional, which we identify as non-linearly evolved density-fluctuation power spectra, have a universal form and can thus be tabulated for fast access in perturbation schemes. In this paper, we focus on a complete evaluation of the free generating functional of KFT, not including particle interactions yet. This implies that the non-linearly evolved power spectra contain a damping term which reflects that structures are being wiped out at late times by free streaming. Once particle interactions will be taken into account, they will compensate this damping. If we suppress this damping in a way suggested by the fluctuation-dissipation relations of KFT, our results show that the complete hierarchy of initial momentum correlations is responsible for a large part of the characteristic non-linear deformation and the mode transport in the density-fluctuation power spectrum. Without any adjustable parameters, KFT accurately reproduces the scale at which non-linear evolution sets in. Finally, we further develop perturbation theory based on the factorization of the generating functional and propose a diagrammatic scheme for the perturbation terms.


Introduction
Based on the pioneering work [1][2][3][4] and in the spirit of [5], we have developed a new theory of structure formation in ensembles of classical particles assumed to be subject to Hamiltonian dynamics and initially correlated in phase space [6]. Structurally, the theory resembles a non-equilibrium quantum field theory. Its central object is a free generating functional describing how the initial phase-space distribution of the particles is transported forward in time. The symplectic structure of the Hamiltonian equations and the deterministic trajectories of the classical particles allow substantial simplifications compared to a quantum field theory. Particle interactions are taken into account by applying to the free generating functional an interaction operator which can be expanded into a power series reflecting increasing orders of the interaction. Cumulants of collective fields, such as the macroscopic mass density, can be read off the generating functional by repeated suitable functional derivatives.
Compared to other approaches [see 7-14, for a a review and a small selection of recent papers], KFT has several conceptual and methodical advantages. Most importantly, since KFT is based on the Hamiltonian flow in phase space, the problem of shell crossing does not occur. Instead, the flow in phase space is diffeomorphic and, due to the symplectic structure of the Hamiltonian equations, even volume-conserving. Difficulties with multiple streams, as they plague standard perturbation theories based on either the Boltzmann equation or the hydrodynamic equations which assume the existence of uniquely valued velocity fields, are thus avoided by construction in KFT because phase-space trajectories do not cross. Related difficulties caused in Lagrangian perturbation theories [see 15-19, for examples] by functional determinants developing singularities in convergent flows are absent from KFT because the functional determinant of the phase-space flow is constant at unity. Since KFT neither assumes the existence of smooth density or velocity fields nor their uniqueness and avoids taking moments over momentum space, it contains in principle the complete hierarchy of moments and of the particle correlations in configuration and momentum space. As we have shown in [20], the order of the perturbative approach to particle interactions in KFT controls the order of correlations taken into account in a related BBGKY hierarchy. Since no smooth and uniquely-valued velocity field is assumed to exist, the particle motions in phase space also trace the formation of vorticity on small scales and late times. Finally, the linearity of Hamilton's equations guarantees the existence of a Green's function. Since splitting the Hamiltonian into parts interpreted as unperturbed and perturbed contributions is to a large degree arbitrary, the Green's function can be chosen such that the interaction Hamiltonian becomes small. The latter is one of the main reasons for first-order perturbation theory to be highly successful in KFT, as shown in [6].
We note that the fundamental mathematical framework of KFT is adopted from statistical non-equilibrium field theory and thus routinely used in wide areas of theoretical physics [see 21, for a recent textbook]. Even though it is not common in cosmology yet, we emphasize that it is not the formalism of KFT that is new, but the application of this formalism to cosmological structure formation, using the particles' phase-space trajectories as fundamental fields. The possibly unfamiliar, but otherwise well established mathematical approach is outweighed by far by the substantial conceptual and methodical advantages of KFT listed above.
Applying this theory to dark-matter structures in cosmology, we showed that even at first order in the interactions between particles on Zel'dovich-type trajectories, the non-linear evolution of the cosmic density-fluctuation power spectrum known from numerical simulations is reproduced very well to redshift zero and to arbitrary wave numbers [6]. In addition to the linearization of the interaction operator, we further approximated the initial hierarchy of momentum correlations to second order, i.e. we truncated this originally exponential hierarchy after the quadratic momentum-correlation terms. This approach caused our formalism to become quite cumbersome.
Here, we show how the free generating functional of our theory can be fully factorized, taking the complete hierarchy of the initial momentum correlations into account. Besides being more accurate, this entails several major advantages: the formalism simplifies considerably, the development of perturbation theory becomes much more tractable, and the factors have a universal form that can be evaluated and tabulated for fast access in automated evaluations of perturbation terms.
In a first cosmological application of the free generating functional of KFT, we show that including the complete hierarchy of momentum correlations accurately reproduces the scale below which the characteristic non-linear deformations of the density-fluctuation power spectrum set in, that and under which circumstances it leads to mode transport from large to small scales, and that a good fraction of the amplitude of the fully non-linearly evolved power spectrum can already be recovered even in the free theory.
In Sect. 2, we briefly summarize the theory, focussing on its free generating functional, specify the initial momentum correlations, and show how the free generating functional can be factorized. First consequences of the theory for cosmological structure formation are described in Sect. 3. A systematic perturbative approach to the theory is developed in Sect. 4 together with a diagrammatic representation of the perturbation terms. We summarize our results in Sect. 5. Appendix A contains further detail on the correlation function of the initial velocity potential. In Appendix B, we present the details of the calculation leading to the factorization of the generating functional.

Generating functional, cumulants, and momentum correlations
We give a brief overview of our non-equilibrium kinetic theory for correlated classical particle ensembles here as we have developed it in [6]. For further detail, we refer the reader to that paper.

Generating functional and cumulants
The central object of the theory is the free generating functional Z 0 [J, K] with generator fields J and K coupled to the phase-space coordinates x = (q, p) and the one-particle response fields, respectively. The bold-faced symbols denote tensors bundling contributions from all N particles in the ensemble. Let x j = ( q j , p j ) be the phase-space coordinates of particle j, and e j an N-dimensional unit vector with components ( e j ) k = δ jk , then We define a scalar product between two such tensors by with Einstein's summation convention implied.
This generating functional is the integral over the N-particle phase space at the initial time t 0 = 0. In terms of the Green's function G of the free equations of motion, the particle trajectoriesx in phase space arē The phase-space measure in (3) is with a suitable probability distribution P to be adapted to the initial conditions at hand. For N particles initially correlated in phase space, the probability distribution P(q (i) , p (i) ) will be given in (21) below. More explicitly, the N-particle Green's function G is the tensor product of the matrix-valued, one-particle Green's function and the unit matrix I N in N dimensions. Particle interactions are included by applying an interaction operator to the free generating functional, producing the full generating functional Here, v is the interaction potential, andρ andB are the density-and response-field operators, respectively. The many-particle response field B describes how a particle ensemble responds to a change in the phase-space coordinates of one of its particles. The arguments abbreviate 1 := (t 1 , k 1 ) and −1 := (t 1 , − k 1 ). In a Fourier-space representation, these operators are sums over one-particle operatorŝ and the integral in (9) is taken over 1 := (t 1 , k 1 ). The response-field operatorB thus contains a density operatorρ. Correlators of order n in the density, say, are obtained by applying n density operators to As usual in statistical field theory, each of the generator fields J and K is set to zero once all functional derivatives with respect to J or K have been applied. An approach to a perturbative evaluation of (11) begins with expanding the exponential interaction operator exp(iŜ I ) into a power series, introducing two density operatorsρ and one response-field operatorb per power ofŜ I . Thus, for an n-th order correlator with m-th order particle interaction, we need to evaluate one-particle expressions of the form with r = n + 2m, and with the particle indices j s = 1 . . . N. The density operators thus replace the generator field J by the shift tensor We define the position and momentum components of the shift tensor L out by the projections and abbreviate L q j := L q j (0) , L p j := L p j (0) (15) with t 0 = 0. The symbol I 3 denotes the unit matrix in three dimensions. The subsequent application of a single, one-particle response-field operatorb j c (c) to Z 0 [L, K], taken at K = 0, simply returns a response-field pre-factor b j c (c), given by in terms of L p jc as defined in (14). Inserting (14) results in which has two important consequences for our later considerations. First, the positionmomentum component g qp (t s , t c ) appears here, evaluated at the two times t s and t c . Due to causality, the particle ensemble can only respond to causes prior to the response, expressed by g qp (t s , t c ) = 0 for t c ≥ t s . Thus, only response-field factors with t c < t s do not vanish. Second, the Kronecker symbol δ j c j s identifies two particle indices. We shall return to these two properties of the response-field factors later. Our next goal will now be to evaluate the free generating functional Z 0 [L, 0] after application of an arbitrary number r of density operators. Constructing the phase-space probability distribution P(q, p) in [6], we have assumed that a statistically homogeneous and isotropic, Gaussian random velocity potential ψ exists such that the momentum p at an arbitrary position is its gradient, Continuity then demands that the density contrast δ is its negative Laplacian, Then, the density-fluctuation power spectrum P δ (k) specifies the initial phase-space probability distribution P(q, p) completely. As we have shown in [6], it is given by where C pp is the momentum-correlation matrix to be defined in (28) and discussed in the following Section. The correlation operator C( p) appearing here and defined in [6] can safely be approximated by unity, for correlators evaluated at sufficiently late times if the q-p-component of the Green's function is unbounded. In the cosmological application we are aiming at here, sufficiently late means that the cosmological scale factor a needs to be much larger than the scale factor a i corresponding to the time when the phase-space distribution of the particles is initially set, a a i . Adopting a i ≈ 10 −3 according to the release of the cosmic microwave background, a > 0.01 seems safe for approximation (22) to hold. With the probability distribution (21), the integrations over the momenta in (3) can be carried out straightforwardly. Then, after applying an arbitrary number r of density operators and setting the generator fields to zero, the free generating functional of our microscopic, nonequilibrium, statistical field theory for canonical ensembles of N classical particles enclosed by the volume V can be written in the form valid at all sufficiently late times. We note that this expression for Z 0 [L, 0] needs to be summed over all different particle configurations, expressed by the indices j 1 , . . . , j r appearing in (12). The integral in (23) is carried out over all particle positions q, and L q and L p are the position and momentum shift tensors resulting from applying the density operators, with components defined in (14). If only density operators are applied, the shift tensors L q and L p will be sums over as many terms as density operators have been applied, with each term representing the contribution of a particle to the density. Since the possible later application of response-field operators will identify two particles per operation, the number of particles involved will be lowered by one for each response-field operator applied. For each particle, one pair of shift vectors ( L q j , L p j ) will remain. Thus, if r density operators and m response-field operators have been applied, the number of different particles involved will be l = r − m = n + m.

Momentum-correlation matrix
The momentum-correlation matrix is where the matrix E jk singles out the particles j and k from the ensemble of N particles, The amplitude σ 2 1 is defined as the first moment of the power spectrum P ψ of the velocity potential ψ. Generally, the moments σ 2 n are defined by and the velocity-potential power spectrum is related to the density-fluctuation power spectrum P δ by due to (20). By definition, the momentum-correlation matrix is given by where q jk is the separation vector between particles j and k. Expression (28) is equivalent to where ξ ψ (q) is the correlation function of the velocity potential, taken at distance q and normalized to σ 2 1 . Due to isotropy, ξ ψ can only depend on q, but not on the direction of q. Accordingly, the tensor of second derivatives ∇ ⊗ ∇ needs to be expressed in terms of derivatives with respect to the particle separation q. Letq be the unit vector in the direction of q, and by its means define the projectors π :=q ⊗q ,π ⊥ := I 3 −π (30) parallel and perpendicular toq. Then, and with the prime denoting the derivative with respect to the argument. The correlation function ξ ψ of the velocity potential and its derivatives ξ ψ and ξ ψ are worked out in Appendix A together with accurate fit formulae convenient for fast numerical evaluations. The quadratic form remaining in (23) splits into two terms by inserting (24), Replacing the sum of squares by a squared sum, we can write instead with the damping terms We shall see below that Q 0 will vanish identically in important cases and that Q D has an intuitive and important effect on the time evolution of the density-fluctuation power spectrum. We shall refer to the Q 0 and Q D as dispersion and diffusion terms, respectively.

Factorization of the generating functional
We now turn to factorizing the generating functional in the form (23), which is a lengthy procedure detailed in Appendix B. The essential ideas are that, in a statistically homogeneous field, only relative particle coordinates q j − q i must matter, and that all these coordinate differences must be statistically indistinguishable. The final result of the calculations presented in Appendix B is that the free generating functional (23) for a shift tensor L with contributions from l different particles can be completely factorized, The index pairs (a, b) and ( j, k) are defined in (B.5) and (B.8), respectively. The function P jk appearing in each of the factors in (37), is a non-linearly (and non-trivially) time-evolved density-fluctuation power spectrum, as we shall discuss in detail in the next Section. The expression ∆ jk abbreviates The wave vectors k jk are defined in (B.15), with the indices ( j, k) given in (B.5) and the indices (a, b) defined in (B.8). We shall call the wave vectors k j1 external because they contain the shift vectors L q j , and the remaining wave vectors k jk with k ≥ 2 internal because they can entirely be integrated out. The quantities λ ,⊥ jk are defined by according to (B.30), with the projectors π jk and π ⊥ jk being defined with respect to k jk . Letk be the unit vector in the direction of k jk , then The factorization (37) of the free generating functional and the expression (38) for the timeevolved density-fluctuation power spectrum are the first main results of our paper. The power spectrum P jk is particularly relevant for cosmological structure formation.

Cosmological consequences
For illustrating the cosmological consequences of our results (37) and (38), we will now reduce (37) to the simplest possible case. For l = 2, the free generating functional (37) contains information on the evolution of the density-fluctuation power spectrum, but neglecting any particle interactions. Thus, this corresponds to the free evolution of the density-fluctuation power spectrum. The free generating functional then shrinks to The single remaining wave vector is k 21 = L q 2 , and the Dirac delta function in (42) ensures that L q 2 = − L q 1 . Then, according to (14), and therefore, by the definitions (40) This brings (38) into the form where we have dropped all indices because only the single index pair 21 now remains to be considered. According to (B.41), this expression turns into in the limit of early times or small wave numbers, which reflects the linear growth of the power spectrum. This emphasizes once more that P(k) is a density-fluctuation power spectrum, with its time evolution modified by the onset of non-linear evolution.
Since the two momentum shift vectors L p 1 and L p 2 are equal in length and opposite in sign, the dispersion term Q 0 from (36) vanishes, but the diffusion term Q D does not, We show in Fig. 1 the power spectrum P(k) from (45) times the diffusion factor exp(Q D /2) after different periods of time expressed by the propagator g qp (τ, 0). This product expresses the free evolution of the power spectrum and will be denoted by P (0) δ .  Figure 1. The solid curves show the the freely evolved power spectrum P (0) δ , i.e. the non-linearly evolved spectrum P times the diffusion factor exp Q D , for five different times, expressed by the propagator g qp (τ, 0). The dashed lines show the initial density-fluctuation power spectrum, linearly evolved to the same times. At large scales, all spectra P (0) δ coincide with the linearly evolved spectrum. On smaller scales, diffusion sets in and suppresses the structure, beginning with small-scale structures at early times and progressing towards larger scales as time proceeds. Figure 1 illustrates two aspects of the evolution. First, the non-linearly evolved spectrum P coincides with the linearly evolving density-fluctuation spectrum for all times at sufficiently small wave numbers, i.e. for sufficiently large structures. Since the free generating functional does not contain interactions between individual particles, structures are damped by free streaming, provided they are small enough. Thus, diffusion proceeds from small to larger scales as time progresses. Starting from (38) and including the diffusion term, we can write the expression for the damped power spectrum briefly as with y := g 2 qp k 2 (a λ + a ⊥ λ ⊥ ) and y → y 0 = Q D /2 for q → 0. Our analysis allows us to largely remove this damping effect from the evolution of the power spectrum. We do so by introducing a parameter 0 ≤ α ≤ 1 meant to gradually switch on the momentum correlation in a thought experiment. We can then write e y 0 P = q e y 0 (e y − 1) e i k· q = e y 0 We have arranged terms such that one contribution to the damping term builds up as α increases and the correlation strengthens, while the other contribution is strongest for α = 0 and decreases as the correlation grows. We ignore this latter term, thus suppressing the contribution to damping present without particle correlations, and replace (49) by If the remaining damping contribution was absent,P = P. By construction, this operation removes that part of the power suppression due to uncorrelated particle momenta, which dissipate the newly formed structures by free streaming. The resulting power spectraP are shown in Fig. 2 for the same times as the power spectra P (0) δ in Fig. 1. The undamped power spectra shown in Fig. 2 illustrate how structure builds up beginning at small and proceding towards larger scales as time progresses. At the latest time used in Figs. 1 and 2, approximately corresponding to the present time in the standard ΛCDM cosmology, non-linear structure formation has reached wave numbers of k ∼ 0.2 − 0.3 h Mpc −1 , in excellent agreement with the transition to non-linearity revealed by numerical simulations.
We can of course not expect our simple procedure to return the correct amplitude of the power spectrum because we have not included particle interactions yet and thus have not gone beyond the lowest-order approximation of the generating functional, expressed by setting l = 2. Nonetheless, the shape of the undamped, non-linearly evolved density-fluctuation power spectrum at late cosmic times, which is obtained in a completely parameter-free way in our theory, reflects the shape of the fully non-linear power spectrum found in numerical simulations very well. We should like to emphasize that the procedure according to (50) applied to the non-linearly evolved power spectrum finds a deeper explanation in the fluctuation-dissipation theorems following from our kinetic theory of structure formation. We shall discuss this in detail in a separate paper now in preparation, dedicated to fluctuation-dissipation relations in KFT. The response to fluctuations in the particle density is mediated by particle interactions. This supports the intuition that the excess diffusion seen in the curves of Fig. 1 will be compensated by particle interactions, which are not included yet in our analysis of the free generating functional. Based on the factorization of the free generating functional demonstrated in this paper, and using the diagrammatic approach to the factorization to be proposed in the following Section, we shall develop a systematic perturbative approach to non-linear structure formation including particle interactions in another dedicated paper.
The fact that the undamped, non-linearly evolved power spectraP shown in Fig. 2 already resemble the shape of the non-linearly developed power spectra of fully numerical simulations shows that this shape is not set by the particle-particle interactions, but rather by their only property we have included here, which is the correlation of their initial momenta. ] g 2 qp = 5 · 10 1 g 2 qp = 5 · 10 2 g 2 qp = 5 · 10 3 g 2 qp = 5 · 10 4 g 2 qp = 5 · 10 5 Figure 2. Non-linearly evolved power spectra are shown here for the same times as in Fig. 1, but now with the late-time diffusion suppressed as specified in (50). In this representation undamped by free streaming, the non-linearly evolving power spectra indicate how structure formation proceeds from small to large scales.
The deformation of the power spectrum beginning at k 0.2 − 0.3 characteristic for the late-time non-linear evolution of cosmic structures is therefore determined by the statistical initial conditions of the particle ensemble, specifically by the initial correlation properties of the particle momenta, presumably imprinted by the Gaussian random density-fluctuation field generated by cosmological inflation.
We notice also that the curves shown in Fig. 2 do not diverge for large wave numbers k. Rather, they approach an asymptotic slope close to that of the linearly evolved densityfluctuation power spectrum, which tends to k −3 for cold dark matter. This indicates that there is no inherent limit to applying our theory to arbitrarily large k, which was to be expected because our phase-space approach does by construction not suffer from any problems with multiple streams. Of course, the amplitude of the density-fluctuation power spectrum will be affected by particle interactions. However, the factorization (37) of the free generating functional together with the asymptotic behaviour of the curves shown in Fig. 2 both indicate that non-linearly evolved density-fluctuation power spectra evaluated at an arbitrary order of particle interactions will be a convolution of curves with an asymptotic slope near k −3 , multiplied by the subsequent application of damping factors. It is plausible that the result will retain the same asymptotic slope, but we will have to demonstrate that. However, the regular asymptotic behaviour of our non-linearly evolved power spectraP(k) for large k suggests that KFT can in this sense be extended to arbitrarily large wave numbers. Figure 3 shows the derivative ofP from (50) with respect to g 2 qp , divided by the initial density-fluctuation power spectrum P δ (k). As the limit (46) shows, for large scales or early times. Since g qp = D + , i.e. the linear growth factor of density perturbations, any deviations of this expression from unity indicates structure growth different from the growth expected in linear theory.  . This figure shows the derivative of the non-linearly growing power spectrumP with respect to the squared propagator g 2 qp = D 2 + , divided by the initial power spectrum P δ (k). If the power spectrum grew linearly, the curves would be flat at unity. The figure shows that, for late times or small-scale structures, the power spectrum grows far beyond the linear growth. The decrease of the curves at late times reflects the remaining damping, which will be compensated by gravitational interaction.
The curves in Fig. 3 illustrate how the growth factor of cosmic density fluctuations becomes scale dependent. They show that, for small-scale structures or at late times, the power spectrum grows substantially more rapidly than expected from linear theory. The decrease of the curves for small-scale structure at late times indicates the remaining damping, which will be counteracted by gravity once particle interactions will be included.
So far, we have evaluated the free generating functional for the simplest possible case, i.e. l = 2, where λ jk = −1 and λ ⊥ jk = 0. This happens if the two momentum shift vectors L p j and L p k appearing in P jk are equal and opposite to each other, and if they are aligned with the wave vector k jk . If the momentum shift vectors are however not aligned with the wave vector k jk , the parallel projection λ jk > −1. This expresses a configuration in which the momenta of two particles are misaligned with the separation vector between the two particles. Power spectrā P jk for such cases are shown in Fig. 4.  Figure 4. Power spectraP jk for late times and for configurations with λ jk ≥ −1, indicating that the momentum shift vectors L p j and L p k are misaligned with the separation vector connecting the two particles. In this case, the power at large scales is reduced compared to the linearly evolved power spectrum, and enhanced at small scales. This shows that such a misalignment causes mode transport from large to small scales.
The figure shows that such a misalignment between particle momenta and their separation vector leads to a reduction of power relative to the linearly evolved power spectrum at large scales, and to an enhancement on small scales. This causes mode transport from large to small scales, adding to the characteristic deformation of the power spectrum by non-linear evolution.
The non-linear growth of the power spectrum and the mode transport contained inP jk and illustrated by various examples in Figs. 1 through 4 thus leads to a characteristic deformation of the power spectrum compared to its linearly evolved shape. As our derivation shows, this deformation does not reflect the gravitational interaction between the particles, but is a consequence exclusively of the initial correlations of the particle momenta. Particles with momenta aligned with their separation vector lead to an enhancement of power and thus to structure growth on small and intermediate scales. The more these particle momenta are misaligned with their separation vector, the more power is transported from large to small scales. Even before any gravitational interaction between the particles is taken into account, the full hierarchy of initial momentum correlations gives rise to a characteristic deformation of the power spectrum. Initial momentum correlations of the dark-matter particles thus contribute substantially not only to the amplitude of cosmic structures, but also to characteristic re-distribution of power from larger to smaller scales.

Expansion of the generating functional
The main goal of this Section is to develop a systematic approach to evaluating (37). It begins by realising that the product over the index pair ( j, k) in (37) can be expanded into a sum ordered by an increasing power of power-spectrum factors. Since there are factors in this product, we can write formally If the power-spectrum factors P are small enough, this sum can be truncated at low powers f . According to (37), all internal wave vectors k ab = k ab are then to be integrated over. For evaluating the result of this sequence of integrals applied to the sum terms in (53), it is important to see how many of these integrals can trivially be carried out due to the ∆ factors appearing in these sum terms. For f = 0, for example, the remaining integrals set all of the internal wave vectors to zero, leaving k j1 = L q j and a>b k ab j>k This is a pure shot-noise term. Identifying and counting the power-spectrum factors with external or internal wave numbers thus allows to simplify the remaining expressions substantially. Let us illustrate this simplification with one more example. For f = 1, the single powerspectrum factor can depend either on an external or an internal wave number. If it is external, the ∆ factors set all internal wave numbers to zero as for f = 0 before, leaving again the external wave numbers k j1 = L q j . All of these wave vectors k j1 except one are also set to zero by the remaining ∆ factors. Without loss of generality, we can label the one non-zero wave vector by the index j = 2 because the particles are indistinguishable. Then, the power-spectrum factor P( L q 2 ) appears besides the remaining ∆ factors. If, however, the single power-spectrum factor depends on an internal wave number, all except one internal wave numbers are integrated out. Again without loss of generality, we can label the only remaining internal wave number k 32 = k 32 . The external wave vectors are then and k j1 = L q j for j > 3. All external wave vectors appear as arguments of ∆ factors in this case. The remaining integration over k 32 sets k 32 = L q 3 by means of the factor ∆ 31 and leaves k 21 = L q 2 + L q 3 . There are l−1 external and F −l+1 internal wave vectors. Thus, there are l−1 possibilities for choosing an external and F − l + 1 possibilities for choosing an internal wave number, allowing us to write a>b k ab where the overall Dirac delta distribution in (37) allows us to replace ∆ 21 ( L q 2 + L q 3 ) by ∆ 21 ( L q 1 ).
From these examples for f = 0 and f = 1, a straightforward scheme emerges for evaluating the sum terms in (53) and the subsequent integrations over all internal wave vectors: (i) Expand each term into a sum ordered by the number of power-spectrum factors depending on external wave numbers, where the subscripts 'ext' and 'int' indicate that respective factors depend on external or internal wave vectors.
(ii) Use the ∆ factors depending on internal wave vectors to set (F − l + 1 − f + e) of the internal wave vectors to zero. Label the remaining ( f − e) internal wave vectors beginning with k 32 and determine the external wave vectors k j1 according to (B.15).
(iii) Use the ∆ factors depending on external wave vectors and the integrals over the remaining internal wave vectors to eliminate as many of the internal wave numbers as possible from the external wave vectors k j1 and from the arguments of the power-spectrum factors depending on internal wave numbers.
This procedure lends itself to be evaluated by a symbolic computer code returning the terms appearing in the sum (58), given l and f . We are now in the process of developing such a code, aiming at possibly completely automatizing the evaluation of perturbation terms.

Summing over particles
As often in perturbation theories, diagrams are useful to keep an overview of the terms involved.
We shall now develop suitable diagrams for the perturbative approach to our non-equilibrium field theory for correlated classical particle ensembles. We return to the free generating functional Z 0 [L, 0], evaluated at a specific shift tensor L, with the generator field K set to zero. The position and momentum components of L are given in (14) or, after inserting (13) there, by In these expressions, the index 1 ≤ j ≤ l identifies the l particles whose positions are being correlated, while the indices j s assign particles to wave vectors with field labels s. The Kronecker symbols appear in (59) because only such phase-space positions contribute to the density which are occupied by particles. Given a specific set of particle indices { j 1 , . . . , j r }, we write to denote the particle indices explicitly. As indicated above, the response-field factors given by (18) have two crucial properties affecting the selection of terms that can or cannot contribute to the perturbation series: the Kronecker symbol appearing there identifies the two particles with indices j s and j c and thus assigns the same particle to the wave vectors labelled by c and s. Furthermore, since the propagator g qp (t s , t c ) vanishes if t s ≤ t c , only such terms can contribute for which t s > t c . Applying r one-particle density and m one-particle response-field operators to the free generating functional thus leaves us with the expression Several aspects are important to note at this point. First, each response-field factor with index c identifies two particles with the indices j c and j s c . Second, no terms can contribute to any response-field factor which belong to the same times because g qp (t, t) = 0. Therefore, only such particles may be identified which are assigned to wave vectors at different times. Fourth, since g qp (t 1 , t 2 ) = 0 for t 2 ≥ t 1 , any response-field factor implies a time ordering in the sense Beginning with a set of particle indices { j 1 , . . . , j r }, we thus proceed as follows to evaluate the terms in (61): We first identify as many particle pairs as response fields occur, i.e. we identify m pairs of particle indices, taking care that particle pairs must not identify positions with equal times and that a time-ordering is involved in all response-field factors. The remaining r − m particles form a reduced set { j} of particle indices. All possible reduced sets { j} must finally be summed over.
Let us illustrate this procedure with the simple example r = 4 and m = 1, corresponding to the terms contributing to a two-point density correlator calculated at first-order in the particle interaction. Since the interaction is instantaneous, t 3 = t 4 , and if the correlator is chosen to be simultaneous, t 1 = t 2 .
In this case, (61) simplifies to the two possible terms Other terms would identify particles at the same time and thus return zero.

Diagrams
Combining (37) and (53) can profitably be represented by diagrams which greatly help constructing and ordering the terms appearing in the generating functional. The essential point of the diagrammatic representation which we are going to construct now is to systematically construct all wave vectors k jk according to (B.15) which enter into the factors (P jk + ∆ jk ) appearing in the generating functional. The diagrams are constructed according to the following rules: (i) Mark the free generating functional Z 0 [J, K] by a circle. According to (10) and (12), each one-particle density operatorρ j s applied to Z 0 corresponds to a functional derivative with respect to a component J q js of the generator field J. According to (12) and (13), each of these operations adds a wave vector k s to the shift tensor L that we need to find for each term in a perturbation series. Thus, for each of the r = n + 2m density operators applied for a perturbation term appearing in an n-th order correlator at m-th order in the particle interaction, we attach a wave vector to the circle symbolizing Z 0 , pointing outward.
(ii) The times when these density operators act are represented by filled dots on the circumference of the Z 0 symbol. Each k vector thus begins at a filled dot. Since each interaction is instantaneous, all internal k vectors need to be pairwise attached to the same time. If the correlator to be calculated is simultaneous, the external k vectors are also attached to the same time.
Thus, the internal wave vectors representing interactions are pairwise attached to the same points in time. Times are ordered counter-clockwise, with the latest times appearing on top. The external wave vectors appearing in the final correlator are attached to the same time if the correlator is simultaneous.
(iii) Between the internal wave vectors corresponding to the density and the response-field operator of an interaction term (9), a further factor appears representing the interaction potential v. If this potential depends only on the separation between the two particles interacting, it requires the two k vectors of the density and the response-field operator it connects to be equal in magnitude and opposite in sign.
Internal wave vectors are marked with primes. To include the interaction potential into the diagram, we connect any two internal wave vectors with a circled v which enter into a single interaction term. For translation-invariant potentials v, the two internal wave vectors connected by a potential are equal and opposite.
(iv) According to (18), each response field identifies two particles at different times, i.e. it assigns the same particle to two different wave vectors at two different times. Adopting the convention in (9), we assign response fields to the negative internal wave vectors appearing in the interaction terms.
Thus, response fields are represented by dashed circle segments between two different wave vectors attached to different times. Each response field must begin at a negative internal wave vector and must be connected to exactly one other wave vector, internal or external.
(v) Equivalent diagrams can appear multiple times. For example, in the diagram shown in Fig. 5 representing a contribution to the density power spectrum at first order in the particle interactions, the response-field arrow can end on k 1 or k 2 . In a homogeneous random field, these two wave vectors are equivalent, and the diagram with the response-field arrow attached to the vector k 2 corresponds to an identical perturbation term. Thus, diagrams are assigned a multiplicity corresponding to the number of equivalent configurations they express. Figure 5 shows the single diagram according to these rules representing the first-order interaction contribution to the second-order density power spectrum. Figure 5. Diagram representing the contribution of the first-order interaction term to a twopoint density correlator. The diagram has a multiplicity of 2 because the wave vectors k 1 and k 2 are indistinguishable in a homogeneous random field.
Diagrams for higher-order interactions or correlators of higher order are now easily constructed. To give an example, we show in Fig. 6 the four non-equivalent diagrams contributing to a second-order density correlator at second order in the particle interactions. Figure 6. Diagrams representing the terms contributing to a two-point correlator at second order in the interaction. Both diagrams in the top row have a multiplicity of 2 because k 1 and k 2 can be interchanged in a homogeneous random field. The multiplicity of the diagrams in the bottom row is 4 because the time order of the two interactions can be exchanged.
Position and momentum shift vectors according to (59) can now be read off these diagrams as follows: Randomly assign particle indices to the k vectors extending from Z 0 , thereby assigning the same particle index to such k vectors connected by a dashed line representing a response field. These will be r − m indices in total, for which the numbers 1 . . . r − m can be chosen without loss of generality. Any permutation of these indices will result in an equivalent set of shift vectors.
For example, we can use the diagram in Fig. 5 to assign the particle indices (1, 2, 3) clock-wise to the k vectors extending from Z 0 , beginning with k 1 . This results in the assignment of particles (1, 2, 3, 1) to the wave vectors ( k 1 , k 2 , k 1 , − k 1 ) because k 1 and − k 1 are connected by a response field. With (59), this implies the position shift vectors and the momentum shift vectors with t denoting the time of the interaction and t > t the time where the correlator is to be evaluated. Any permutation of the particle indices would merely permute the labels on these shift vectors. These diagrams allow a quick construction of all terms contributing to the generating functional Z for given orders of correlators and of particle interactions. From these diagrams, the shift vectors L q and L p as well as the response-field factors can be read off. They can then be inserted into the complete factorization of the generating functional to evaluate the perturbation terms. The procedures involved can now be implemented in a symbolic computer code.

Summary and conclusions
In [6], we have developed a kinetic non-equilibrium field theory for cosmic structure formation. The central object of this theory is a free generating functional which describes how an initially correlated ensemble of classical particles propagates in time under Hamiltonian dynamics. Particle interactions are included by an exponential interaction operator whose series expansion suggests a natural perturbative approach. The initial correlation of the particle momenta was shown in [6] to be described by a Gaussian in which the momentum-correlation matrix enters as a quadratic form. In [6], we expanded this Gaussian to second order in the quadratic form and showed that, to first order in the particle interactions, the non-linear cosmic density-fluctuation power spectrum known from numerical simulations could be well reproduced. To a large part due to the expansion in the order of the momentum correlations, the notation as well as practical calculations became quite cumbersome.
In this paper, we have substantially simplified the theory by factorising the free generating functional into factors of universal shape, taking the complete hierarchy of momentum correlations into account. Our first main results are thus the expressions (37) for the free generating functional and (38) for the non-linearly evolving density-fluctuation power spectra P. The factorization suggests the expansion (53) of the free generating functional in powers of P, which contain the full hierarchy of momentum correlations. This expansion in terms of power-spectrum factors together with the expansion of the interaction operator suggest a diagrammatic representation of perturbation terms, which we have developed here.
Our second main result is that the initial momentum correlations of the dark-matter particles, prior to any gravitational interaction, lead to a characteristic deformation of the density-fluctuation power spectrum compared to its initial, linearly evolved shape. For particle momenta aligned with the separation vector between the particles, power is enhanced on moderate and small scales, while any misalignment between these two vectors leads to a substantial transport of power from large to smaller scales.
Our third main result is the development of a diagrammatic approach to the perturbative terms developed from the complete factorization of the generating functional. Taken together, the factorization of the free generating functional and the diagrammatic approach to perturbation theory open a way to have perturbation terms automatically calculated and evaluated by a combination of a symbolic and a numerical computer code, which we are now planning to develop.
Appendix A.1. Relation to the density-fluctuation power spectrum The potential correlation function ξ ψ (x) introduced in (29) is where P ψ the power spectrum of the velocity potential. With and using the recursion relations for the spherical Bessel functions and their derivatives, we find where the density-fluctuation power spectrum P δ (k) = k 4 P ψ (k) was introduced. The asymptotic behaviour for small arguments of the spherical Bessel functions of the first kind ensure that lim q→0 ξ ψ (q) The functions ξ ψ (q)/q and ξ ψ (q) are shown in Fig. A1. For the ΛCDM density-fluctuation power spectrum, it turns out that these two functions a 1 and a 2 allow simple fits by rational functions, both of which fall like q −2 asymptotically for q q 1 or q q 3 . For a ΛCDM power spectrum according to [22] with cosmological parameters Ω m0 = 0.3 and Ω Λ0 = 0.7, we find the best-fitting parameters listed in Tab. A1. Figure A2 shows the functions a 1,2 (q) together with the fits (A.7) specified by the parameters in Tab. A1.

Appendix B. Factorization of the generating functional
In this Section of the Appendix, we factorize the free generating functional (23). We have shown in [6] that, due to the scaling of the interaction potential with the inverse mean particle density and the identification of particles with each other by the response field, the calculation Table A1. Parameters of the rational functions (A.7) fitting the functions a 1 and a 2 defined in (A.6) for a ΛCDM power spectrum according to [22]  -a 1 (q) a 2 (q) fit to -a 1 (q) fit to a 2 (q) Figure A2. The functions a 1,2 (q), normalized to σ 2 1 /3, are shown here together with the fit functions (A.7) specified by the parameters given in Tab. A1. of a non-shot noise contribution to an n-point correlator at m-th order in the particle interactions requires the contribution of order r = n + 2m to the free density correlator due to l = n + m particles. Thus, there are l entries in L q,p and uncorrelated particles can simply be integrated out from (23). Generally, we denote the number of entries in L q,p by l and enumerate them consecutively from 1 . . . l, which is possible without loss of generality because the enumeration of the particles is unimportant.
Appendix B.1. Introducing relative particle coordinates Returning to (23), we introduce the coordinate differences with respect to the arbitrarily chosen particle 1, Then, the scalar product of the positions q with the spatial shift tensor L q can be written as Since the correlation matrix C pp depends on coordinate differences only and not on absolute positions, the integration over q 1 can be carried out. This results in 3) The momentum-correlation matrix depends on the absolute values of all pair-wise coordinate differences not just on the coordinate differences q j1 with respect to particle number 1. Since there are l particles to consider in total, the number of coordinate differences that C pp depends on is Of these coordinate differences, (l − 1) are taken into account by the (l − 1) difference vectors q j1 . We now extend the integration in (B.3) to all these coordinate differences by introducing the remaining difference vectors additional, dependent coordinates related by (B.7), which must be ensured by including appropriate delta distributions δ D q ab − q a1 + q b1 (B.10) into the integrand, with a and b from (B.8). Thus, the free generating functional turns into Replacing all of these delta distributions by their Fourier expansions, δ D q ab − q a1 + q b1 = k ab e i k ab ·( q ab − q a1 + q b1) (B.12) with auxiliary wave vectors k ab conjugate to the coordinate differences q ab , we arrive at the expression for the free generating functional. Reordering terms, we can rewrite the phase factor as l j=2 L q j · q j1 + a>b k ab · q ab − q a1 + q b1 and a, b from (B.8).
With these results and the definition of Q 0 and Q D in (35), the generating functional (B.13) factorizes completely in the integrations over the coordinate differences q jk , since the correlation matrix C p j p k depends on the distances | q jk | between the particles only. We can thus decompose the generating functional into independent factors for all particle pairs, which are then to be convolved in Fourier space by integrating over all auxiliary wave vectors k ab .

Appendix B.2. Examples
For two-point density correlations, l = 2 + m, hence the number of coordinate pairs is are dependent. For first-order perturbation theory of a two-point spectrum, we have l = 3, and we need to introduce one auxiliary wave vector k 32 . According to (B.15), we then have k 21 = L q 2 + k 32 , k 31 = L q 3 − k 32 , k 32 = k 32 , (B.20) and L q 1 = −( L q 2 + L q 3 ) because of the δ distribution in (B.3). For second-order perturbation theory of a two-point spectrum, l = 4, and three auxiliary wave vectors k 32 , k 42 and k 43 need to be introduced. Then, k 21 = L q 2 + k 32 + k 42 , and L q 1 = −( L q 2 + L q 3 + L q 4 ) .

Appendix B.3. Evaluation of the generic factors
We now need to evaluate the integrals I jk defined in (B.17), which are all of the type with an independent wave vector k 21 representing any of the vectors k jk defined in (B.15). The momentum-correlation matrix C p 2 p 1 defined in (32) depends on the particle separation q only. Repeating (32), we can write C p 2 p 1 as C p 2 p 1 = −π ξ ψ (q) −π ⊥ ξ ψ (q) q (B.23) with the projectorsπ andπ ⊥ defined in (30). We now expand the projectorsπ ,⊥ with respect toq into the respective projectors π ,⊥ with respect to k 21 . Letk be the unit vector in the direction of k 21 , then π 21 =k ⊗k , π ⊥ 21 = I 3 − π 21 . (B.24) For doing so, we expand the Hessian of the potential-correlation function into the projectors π 21 and π ⊥ 21 , D 2 ξ ψ (q) =π ξ ψ (q) +π ⊥ ξ ψ (q) q = a π 21 + a ⊥ π ⊥ 21 , (B.25) multiply this equation by π 21 and π ⊥ 21 and take the trace of the resulting two equations to find a = ξ ψ (q) trπ π 21 + ξ ψ (q) q trπ ⊥ π 21 , 2a ⊥ = ξ ψ (q) trπ π ⊥ 21 + ξ ψ (q) q trπ ⊥ π ⊥ 21 . (B.26) Introducing the cosine µ :=q ·k of the angle between q and k 21 , the traces on the right-hand sides are trπ π 21 = µ 2 , trπ ⊥ π 21 = 1 − µ 2 = trπ π ⊥ 21 , trπ ⊥ π ⊥ 21 = 1 + µ 2 . (B.27) We thus have a = µ 2 ξ ψ (q) + (1 − µ 2 ) ξ ψ (q) q , q , (B.28) and the quadratic form in (B.22) turns into L p 2 C p 2 p 1 L p 1 = − L p 2 π 21 L p 1 a − L p 2 π ⊥ 21 L p 1 a ⊥ . (B.29) It is now convenient to quantify the remaining projections by λ ,⊥ 21 defined by L p 2 π 21 L p 1 = g 2 qp (τ, 0) k 2 21 λ 21 , L p 2 π 21 L p 1 = g 2 qp (τ, 0) k 2 21 λ ⊥ 21 . (B.30) With these definitions, we arrive at the form qp (τ,0) k 2 21 a λ 12 +a ⊥ λ ⊥ 21 +i k 21 · q (B.31) for the integral to be solved. This integral has an intuitive physical meaning. To clarify it, we first split off a delta distribution, Due to the large dynamic range of the argument of the exponential in (B.33) and the fast oscillations of the Fourier phase, the integration required to evaluate P 21 is numerically difficult in some regions of parameter space. We use Levin collocation [23][24][25] for a fast and reliable integration scheme. The meaning of P 21 becomes best visible in the limit of early times. Then, the argument of the exponential in (B.33) is small, the exponential can be approximated by a first-order Taylor expansion, and we are left with This is the density-fluctuation power spectrum, linearly evolved with g 2 qp (τ, 0), times −λ 21 . As we shall see in the main text, λ 21 = −1 in the most straightforward applications. This shows that the function P 21 generalises the linearly evolved density-fluctuation power spectrum and thus that P 21 is a non-linearly time-evolving density-fluctuation power spectrum. Its Fourier transform is the corresponding generalization of the spatial density correlation function.