Kinetic Field Theory for Cosmic Structure Formation

We apply kinetic field theory to non-linear cosmic structure formation. Kinetic field theory decomposes the cosmic density field into particles and follows their trajectories through phase space. We assume that initial particle momenta are drawn from a Gaussian random field. We place particular emphasis on the late-time, asymptotic behaviour on small spatial scales of low-order statistical measures for the distribution of particles in configuration and velocity space. Our main result is that the power spectra for density and velocity fluctuations in ensembles of particles freely streaming along Zel'dovich trajectories asymptotically fall off with wave number $k$ like $k^{-3}$ for $k\to\infty$, irrespective of the cosmological model and the type of dark matter assumed, with the exponent set only by the number of spatial dimensions. This conclusion remains valid for density-fluctuation power spectra if particle interactions are taken into account in a mean-field approximation. We also show that the bispectrum of freely-streaming particles falls off asymptotically like $k^{-11/2}$ under the same general conditions.


Introduction
Our observable universe is permeated by structures on all scales. The Earth is part of the Solar System, located in one spiral arm of the Milky Way galaxy, which is a member of the Local Group of galaxies, which is part of the Virgo Supercluster. Galaxies identified in the Two-Micron All-Sky Survey mark large-scale filamentary structures, galaxy clusters and voids in the nearby universe [1]. These are structures characterizing the cosmic matter distribution at present, almost 14 billion years after the Big Bang. Temperature fluctuations in the cosmic microwave background on the other hand represent density fluctuations in the early universe, about 400000 years after the Big Bang. At that time, density fluctuations had an amplitude of ≈ 10 −5 relative to the mean [2]. By now, the matter density in the central regions of galaxy clusters exceeds the mean cosmic density by factors of ≈ 5 . . . 10 [3,4]. Within 14 billion years, the amplitude of density fluctuations has grown by a factor of 10 6 on such scales. The cosmological standard model requires dark matter, i.e. a dominant form of matter incapable of interacting electromagnetically, to explain this amount of growth [5]. With known forms of matter, the density-fluctuation amplitude could not have grown by more than a factor of ≈ 10 3 in the same time (see eg. [6]).
Based on the assumption of dominating dark matter, of cold dark matter in particular, more and more refined and extended numerical simulations have revealed over decades that the morphology and low-order statistical measures of the observed cosmic matter distribution can be well reproduced in detail [7][8][9]. Impressively highly resolved simulations have shown furthermore that the radial profiles of the darkmatter density in gravitationally bound objects from a wide range of masses have a self-similar, universal form, and measurements in galaxy clusters have confirmed this profile shape [8][9][10][11][12][13].
Cosmic structures are most commonly quantified by their power spectrum, which is the variance of the Fourier modes of their matter-density distribution as a function of the wave number k. As long as the relative amplitude of density fluctuations is smaller than unity, their evolution can be well described by the linearized system of the continuity, Euler, and Poisson equations on the expanding cosmic space-time (see e.g. [14] for a review). In this linear theory, the Fourier modes of the density field evolve independently and at a rate independent of wave number. The linearly evolved density-fluctuation power spectrum thus has the same shape as it had initially, but an amplitude expected to be approximately 10 6 times larger than right after the cosmic microwave background was released (see eg. [6]).
Reconstructions of the linearly evolved power spectrum from many different cosmological measurements confirm that its shape cannot be observationally distinguished from the simplest shape it could be expected to have, growing almost linearly with k at large scales (small k), reaching a broad maximum near k ≈ 10 −2 h Mpc −1 and turning over to a decrease approaching an asymptotic fall-off proportional to k −3 on small scales (large k) [5,[15][16][17], [18]. Numerical simulations show that this power spectrum is characteristically deformed by non-linear evolution on small scales: at wave numbers 1 h Mpc −1 , its amplitude is enhanced by way more than an order of magnitude, and it seems to approach an asymptotic fall-off proportional to k −3 on small scales [19,20]. Measurements of the non-linearly evolved power spectrum via the weak gravitational lensing effect confirm these numerical results on the scales accessible to observation [21,22]. Structures with power spectra approaching an asymptotic, logarithmic slope of −3 are distinguished because their power, defined as the number of density fluctuations times their variance in Fourier space, becomes scale-independent towards small scales: every decade in scale then adds an equal amount of fluctuation power.
While numerical simulations reproduce the statistical properties of the cosmic matter distribution very well and thus strongly support the cosmological standard model as well as the hypothesis of cold dark matter, they cannot identify any fundamental reasons for the self-similar density profiles of gravitationally-bound cosmic structures or the scale-independence of fluctuation power on small scales. Understanding the fundamental origin of this kind of universality of cosmic structures is necessary to decide whether it is caused by any specialities of the cosmological model, the assumed properties of dark matter, the functional shape of the gravitational law, or due to any other reason. Testing wide classes of theoretical possibilities with sufficiently detailed numerical simulations seems forbiddingly costly. The inevitable shot noise in numerically simulated matter distributions on small scales caused by the necessarily finite number of simulation particles adds another motivation to search for rigorous statements on the statistics of cosmic structures on small scales (see e.g. [12] for a recent review on cosmic N-body simulations).
Analytic approaches to cosmic structure formation exist. They fall into the two main classes of Eulerian [23][24][25][26][27][28][29][30][31] and Lagrangian perturbation theory [32][33][34][35][36][37][38][39][40], and enter into the non-linear regime of cosmic of structure formation by perturbative or effective methods [41][42][43][44][45][46][47][48]. These different approaches suffer from one essential problem: they describe the evolution of cosmic structures in terms of dynamical equations for the cosmic density and velocity fields, assumed to be smooth and differentiable. In this sense, dark matter is modelled as a fluid, based on the ideal or viscous hydrodynamical equations. Once convergent streams of dark-matter paricles cross, however, the velocity field is no longer uniquely valued, and the fluid description becomes inadequate. This is the origin of the notorious shell-or stream-crossing problem.
For this main reason, we choose a different approach here, based on kinetic field theory [49][50][51][52][53][54][55]. This theory is kinetic in the sense of describing the statistical properties of a large number of microsopic particles. These particles are classical, subject to Hamiltonian dynamics, and as an ensemble need not be in any kind of equilibrium. The field that the theory is acting on is the bundle of particle trajectories. The theory differs from conventional approaches to kinetic theory in that it does not assume a smooth phase-space density function of the particles to exist. Its central mathematical object is thus not a dynamical equation for a phase-space density, but a generating functional characterizing the initial statistical properties of the particle ensemble and the time evolution of the particle trajectories. Statistical properties of the ensemble are extracted from the generating functional by functional derivation. Since phasespace trajectories subject to Hamiltonian dynamics cannot cross, the theory avoids the shell-crossing problem by construction. By design, it is formally identical to a statistical quantum-field theory.
We shall focus in this paper on rigorous statements that can be derived on smallscale cosmic structures within the framework of kinetic field theory. Other aspects of the theory have been worked out elsewhere, most noticeable on its relation to conventional kinetic theory or cosmological perturbation theory [53], macroscopic reformulations including resummation [56], applications to mixtures of dark matter and gas [57,58], and others. We develop kinetic field theory in Sect. 2 where we describe in detail how the theory can be adapted to the expanding cosmological background space-time. In Sect. 3, we first characterize the statistical properties of the initial state of the particle ensemble and then describe how low-order statistical measures for the evolved particle distribution in configuration and velocity space can be derived. The asymptotic, small-scale behaviour of the density-fluctuation power spectrum, the density-fluctuation bispectrum, and the velocity power spectrum is derived for freely-streaming particles in Sect. 4. In this section, we also show how particle interactions can be included in a mean-field approximation, and demonstrate that the asymptotic behaviour of the density-fluctuation power spectrum is unchanged by such interactions. In Sect. 5, we summarize our results and present our conclusions.

Kinetic field theory in the cosmological context
Kinetic field theory studies the evolution of classical particle ensembles in and out of equilibrium. We are studying canonical ensembles of particles here whose motion is described by Hamiltonian mechanics. Their phase-space trajectories (q i (t), p i (t)) =: x i (t), with the index i = 1, . . . , N enumerating the particles, are completely determined once their initial values (q (i) , p (i) ) are given at some point in time. Phase-space trajectories cannot cross: if they did, they would have one point in common at a certain time, which would force them to have the same past and to continue identically since the solutions of the Hamiltonian equations are unique. This is one of the major advantages of kinetic field theory compared to other approaches: it avoids by construction any problems with matter streams crossing in configuration space.

Notation
We shall introduce two essential pieces of notation here: bundles of particle trajectories and low-order statistical measures for density fluctuations.

Bundles of particle trajectories
The phase-space trajectories for all N particles of the ensemble together are the fundamental field that kinetic field theory is concerned with. For notational convenience, we denote this field as where theê i are the Cartesian unit vectors of R N . Let now H(x, t) be the Hamiltonian on the single-particle phase space. With the symplectic matrix the Hamiltonian equation for a single trajectory can be written aṡ Introducing the matrix M = M ⊗ 1 N (4) for the entire particle ensemble, further the N-particle phase-space gradient 5) and the N-particle Hamiltonian H(x 1 , . . . , x N , t), the equations of motion for the entire trajectory bundle can be compactly written aṡ For brevity and convenience, we introduce the short-hand notations and fix the Fourier convention bỹ We assume three spatial dimensions unless explicitly stated otherwise.

Correlation functions and power spectra
The cosmic matter density ρ is often written as a mean densityρ times a fluctuation, where δ(q) is called the relative density contrast at position q. The probability for finding one matter particle in a small volume dV around position q 1 and another one within dV around position q 2 is then given by where ξ δ (q 1 , q 2 ) is the correlation function of the density fluctuations δ between positions q 1 and q 2 . Thus, the correlation function quantifies the conditional probability for finding a particle at q 2 given another particle at q 1 . Due to statistical homogeneity, ξ δ may only depend on the separation vector q 2 − q 1 between the two points, and due to statistical istropy, it may only depend on the absolute value q = |q 2 − q 1 |, The Fourier transform of the correlation function is the power spectrum, Since the Fourier transform of the density is the two-point function of the Fourier-transformed density is related by to the two-point function of the Fourier-transformed density contrast. The latter is related to the power spectrum by we can write (14) as We shall return to this equation later in (110). On the other hand, the connected part of the two-point function of the Fourier-transformed density, i.e. the two-point density cumulant in Fourier space, is by combining (14) with (13) and (15).

Generating functional for a classical particle ensemble
A bundle of phase-space trajectories x(t) beginning at the initial phase-space points x (i) follows the classical Hamiltonian flow φ cl (x (i) , t). The probability to find the particles at the phase-space points x(t) at time t is where P(x (i) ) is the probability for the initial particle positions x (i) to be occupied, and P(x(t)|x (i) ) is the transition probability for the particle ensemble from there to x(t).
Since the particles are classical and follow deterministic trajectories, this transition probability is a functional Dirac delta distribution, Integrating P(x(t)) over the trajectory bundle x(t) and introducing a generator field J(t), we arrive at the generating functional where the parentheses in the exponent denote a suitably defined scalar product between J and x including a time integral, Evaluating the path integral in (20) over all possible trajectory bundles x(t), the delta distribution selects the classical solution of the Hamiltonian equations (6). The generating functional for the entire particle ensemble thus becomes with the integral measure dΓ on the initial phase space, The generating functional Z[J] in (23) is already the central mathematical object of kinetic field theory. By functional derivatives with respect to the generator field J, then setting J = 0, statistical information on the particle ensemble at any time t can be extracted from Z[J] [49,52]. Aiming at statistical information on the particle ensemble, we are not solving dynamical equations for density and velocity fields assumed to be sufficiently smooth. Instead, the time evolution of the particle trajectories implies the time evolution of the generating functional, from which statistical information on the particle distribution in phase-space can be extracted at any required time. Studying particle trajectories instead of smooth fields avoids the notorious shell-crossing problem, which arises in other approaches when a velocity field assumed to be smooth ceases to be uniquely valued after matter streams cross. It should be noted that the generating functional Z[J] contains the exact information on the time evolution of the entire particle ensemble if we insert the exact particle trajectoriesx(t) into (23). Approximation schemes for the trajectories will be introduced as we go along to arrive at tractable expressions. We will now turn to describing particle trajectories with suitably chosen Green's functions.

Particle trajectories
Beginning with the Lagrange function of point particles of equal mass in the expanding cosmic background, we shall derive Green's functions here solving the equations of motion and thus characterizing the particle trajectories through phase space. We shall pay particular attention to choosing appropriate reference trajectories and deriving the force acting between two particles relative to these trajectories (see eg. [52,59]).

Lagrange function, transformation to comoving coordinates
The trajectories of particles with massm on the expanding cosmic space-time are determined by their Lagrange functioñ where the gravitational potential Φ satisfies the Poisson equation with the matter density ρ and the cosmological constant Λ [60]. In terms of the comoving spatial coordinate q, the physical spatial coordinate r is r = a q. During the matter-dominated epoch, the scale factor a obeys Friedmann's equation containing the mean matter densityρ, which is a function of time only. Transforming to comoving coordinates, subjecting the Lagrangian to the gauge transformationL →L − d f dt with f =m 2 aȧ q 2 (28) and dropping the particle massm leads to the Lagrangiañ with the gravitational potential φ now satisfying the comoving Poisson equation containing the density contrast We now introduce the convenient time coordinate instead of the cosmic time t, and H i is the Hubble constant at some initial time t i . The growth factor D + (t) for cosmic density fluctuations is the growing solution of the linearized growth equation (40) below.
We set t i early in the matter-dominated era, e.g. right after the cosmic microwave background has decoupled. We normalize the growth factor to unity initially such that τ = 0 at t = t i .
The inverse Hubble constant at the initial time, H −1 i , sets an appropriate time scale. We express the mean densityρ by the critical density and the density parameter Ω i at the initial time,ρ where the scale factor is also normalized to unity at the initial time t i . Early in the matter-dominated phase, we can safely replace the density parameter by unity, Ω i = 1. If we finally drop the common factor H 2 i from both terms of the Lagrange functioñ L, we find the equivalent Lagrange function The dot now and in the following represents the derivative with respect to the time τ, the potential ϕ satisfies the Poisson equation and m is the effective, dimension-less, but time-dependent particle mass m = a 2 dτ dt .
Note that the time-dependence of the effective particle mass m is a mathematically exact and physically equivalent consequence of introducing comoving coordinates. In physical coordinates, particles are diluted by cosmic expansion such that the gravitational acceleration between them decreases over time. In comoving coordinates, their mean separation remains constant, and the decreasing gravitational acceleration is mapped to their increasing mass. We will need the time derivative of the mass later, i.e. the derivative of m with respect to the growth factor,ṁ where the prime indicates the derivative with respect to the scale factor a. We insert (36) into (37) and use dτ dt  (see eg. [54]). The growth factor itself satisfies the linear growth equation, which simplifies the right-hand side of (39). Early in the matter-dominated era, Ω i ≈ 1 to excellent approximation. Combining this with (35)-(40), we can bringṁ into the simple formṁ with A ϕ specified in (35). From now on, we shall write t instead of τ, expressing all times by the linear growth factor as defined in (32).

Hamiltonian Green's function
The Lagrange function L from (34) implies the canonical momentum p =˙ q/m with the dimension-less, time-dependent mass m. Due to our normalizing both the scale factor a and the growth factor D + to unity at the initial time, the effective particle mass m is also unity initially and increases from there. The Hamilton function splits into a free, kinetic part H 0 = p 2 /(2m) and an interacting part H I = mϕ (see eg. [49,52,54,59]). The free part implies the equation of motioṅ which is solved by the exponential SinceĀ(t, t ) is nilpotent,Ā 2 = 0, the solution (44) simplifies to Varying the constant x 0 leads to the solution for the phase-space trajectory or for the trajectory in configuration space, where is what we call the Hamiltonian propagator [59].

Motion relative to Zel'dovich trajectories
It is convenient in cosmology to replace this Hamiltonian propagator g H (t, t ) by the time difference t−t , which we have chosen to be the difference between linear growth factors (see eg. [59]). To achieve this, we re-write the trajectory (47) as where the dot denotes the time derivative with respect to the first time argument of the Hamiltonian propagator, and define a function A p (t) implicitly by  Differentiating (50) twice with respect to t gives A p =ṁ = mA ϕ D + ; cf. (41) and the definition of g H in (48).
We replace the potential ϕ by the shifted potential where ψ is an initial velocity potential defined to satisfy By assuming a scalar potential for the initial particle velocities, we neglect any initial vortical flows, which does however not mean that the flow remains non-vortical in the course of the further evolution. Since initial velocities and density fluctuations have to satisfy the continuity equation, the velocity potential must be related to the initial density fluctuations δ (i) by the Poisson equation Considering (50) and replacing ϕ by φ, we can bring the trajectory (49) into the form with the potential φ satisfying the Poisson equation where δ (lin) = D + δ (i) is the linearly evolving density contrast. This is an important result in our context. The potential φ is sourced by the difference between the actual density contrast δ and its linearly evolving representative δ (lin) , i.e. φ is the potential created by the non-linear part of the density contrast only. Initially, therefore, φ = 0 and the particles follow the inertial trajectories representing the Zel'dovich approximation [61]. As the density contrast develops a non-linear contribution at later times and on small scales, particles are deflected from the Zel'dovich trajectories. This is the essential reason for introducing the Zel'dovich trajectories (56) as reference trajectories here: the force with respect to these trajectories initially vanishes, and it only builds up as density fluctuations become non-linear. We can finally replace the Hamiltonian propagator g H under the integral in (54) by implicitly defining an effective force f (t) acting relative to the fiducial, Zel'dovich trajectories (56) (see eg. [54]). Then, Differentiating this equation twice with respect to t results in the effective force With the Green's function the full phase-space trajectory can be written as

Effective interaction potential
The Poisson equation (55) for the potential is interesting in its own right. In Fourier space, it readsφ At the same time, φ is the convolution δn * v of the fluctuation δn of the mean particlenumber densityn and the one-particle potential v, thus φ =nδṽ (62) according to the convolution theorem. The effective one-particle potential is then determined byṽδ Multiplying this equation once withδ and once withδ (lin) , taking the ensemble average and eliminating δδ (lin) between the resulting two equations gives At large wave numbers, P (lin) δ P δ , and the one-particle potential approaches the expected Newtonian form ∝ k −2 . At small wave numbers, P (lin) δ = P δ and the potential drops to zero. Relative to the Zel'dovich trajectories (56), the potential thus acquires a large-scale cutoff. In good approximation, its shape resembles the Yukawa form with a time-dependent cutoff scale k 0 delineating the regimes of linear and non-linear structure formation [54].

Remarks concluding this section
We finish this subsection by five remarks which appear important here also in view of wrong statements frequently repeated about the kinetic field theory of cosmic structure formation. 1. For completeness, we repeat that the time-dependent, effective particle mass m is due to the exact transformation from physical to comoving coordinates. It is not at all an approximation. 2. The trajectories (60) are not approximate either, but exact. By a sequence of transformations, we have introduced inertial trajectories with respect to a time coordinate t given by the linear growth factor as introduced in (32). This time coordinate is non-uniform in cosmic time. These inertial, or Zel'dovich, trajectories are chosen to incorporate part of the gravitational interaction between the particles, as quantified by the Poisson equation (56) and the definition (58) of the effective force. Despite their simple form, the Zel'dovich trajectories are subject to gravity by linear density fluctuations because the initial particle momenta are correlated with the initial density fluctuations via the velocity potential ψ: where particles are overdense, their flow converges. 3. The trajectories (60) show that the kinetic field theory of cosmic structure formation goes beyond the Zel'dovich approximation inasmuch as the effective force is taken into account. We shall show later how the gravitational interaction relative to Zel'dovich trajectories can be incorporated in a mean-field approach. 4. Introducing reference trajectories of the Zel'dovich form and an effective force relative to them is nothing mysterious or unusual, and not at all any limitation of the kinetic field theory of structure formation. Newton's axioms introduce reference trajectories defined by inertial motion in coordinate time, and forces as reasons for deviations from inertial motion. The equivalence principle introduces the trajectories of freely-falling particles as a reference, in consequence of which gravitational force is transformed away altogether, and gravitational interaction is reduced to the gravitational tidal field. This illustrates that reference trajectories, and thus Green's functions, can be chosen at will if forces acting relative to them are suitably adapted. The discussion of the effective gravitational potential above shows that the interaction relative to the Zel'dovich trajectories has not an infinite range any more. 5. Deviations from the Zel'dovich trajectories are exactly quantified by the time integral over the force term in (61). Since this force is sourced by non-linear density fluctuations only, it is initially small even at small scales, and remains so at large scales. The integral in (61) thus precisely determines a quantity suitable for a perturbative approach, if perturbation theory is what is wanted. Other approaches, such as the mean-field approximation described below, are possible, less tedious, and more efficient.

Statistics of the particle ensemble
Having set up the generating functional of kinetic field theory including the particle trajectories, we shall proceed in this section by characterizing the initial state of the particle ensemble in the cosmological situation, and by describing generally how statistical information on the evolved particle ensemble can be extracted from the generating functional.

The initial state
Since this section is about the initial particle configuration in phase space only, we drop here the superscript (i) on the initial density contrast δ (i) and the initial phasespace positions The set of the initial phase-space coordinates for all particles is what we call the initial state of the ensemble. Two assumptions are crucial for this initial state: initial velocities are non-vortical, and the initial velocity potential is a Gaussian random field [62].

Probability distribution of initial phase-space positions
As discussed before, we assume that an initial velocity field exists which is the gradient ∇ψ of a velocity potential ψ. By continuity, the initial density contrast δ must then be the negative Laplacian of this potential, δ = −∇ 2 ψ; see (52) and (53). This initial density field is now sampled by N point particles, placed by a Markov process such that the probability P(q|δ) of finding particles at positions q = q i ⊗ê i in the density field characterized by the density contrast δ is (see eg. [49]). Given the joint distribution P(δ, p) of the density-contrast values δ = δ(q i ) ⊗ê i at the particle positions and the momenta p, the probability P(q, p) of the initial phase-space positions for the particle ensemble is Assuming that ψ is a Gaussian random field, the joint probability P(δ, p) is a multi-variate Gaussian. Its characteristic function is then given by characterized by the covariance matrix The joint probability P(q, p) is then With this, we return to (67), where we introduce for later convenience a source field R associated with δ, This allows us to elevate the conditional probability P(q|δ) from (67) to an operator acting on the source field R, which we can pull in front of the integral in (72). This, together with (71), turns (67) into Performing the integral over δ results in the Dirac delta distribution δ D (r − R) such that the ensuing integration over r replaces r by R, with the differential operatorD As we shall discuss in detail later, this differential operator can in cosmological applications at late times safely be approximated bŷ such that is the correlation matrix for the entire set of particle momenta, and C p i p j is the correlation matrix of the momenta of particles i and j [49].

Momentum correlations
Due to the definition p = ∇ψ of the initial momentum field in terms of the velocity potential ψ, we have where ξ ψ is the auto-correlation function of the initial velocity potential taken at the separation q = |q i − q j | of the particles i and j. Thus, Notice that C p i p j = C p j p i because a sign change in q can be cancelled by an irrelevant sign change in k in (82). Carrying out the integral over the directions of the wave vector k, we obtain with the correlation functions and the projector π =q ⊗q (85) parallel to the line connecting the two points which are being correlated [50]. The vectorq is the unit vector in q direction,q = q/|q|. We have used in (84) that P (i) δ = k 4 P ψ due to the Poisson equation (53) between the velocity potential ψ and the initial density contrast δ (i) . The functions a 1,2 (q) have the limits where σ 2 1 is one of the moments of the initial density-fluctuation power spectrum; see also the asymptotic expansions (156) below. Of course, we implicitly need to assume here that these moments exist up to the order n needed in later expressions.

Statistics of the evolved particle distribution
Statistical information can be extracted from the time-evolved generating functional by applying suitable operators (see eg. [49,52]). We shall first introduce density operators here and discuss density correlation functions before we proceed to power spectra for density and velocity fluctuations. −a 1 (q)/(σ 1 2 /3) The initial momentum correlation functions a 1,2 (q) as defined in (86) normalized to σ 2 1 3 (87) are shown as a function of particle separation.

Density operators
The number density of our particle ensemble is a sum over delta distributions, We Fourier transform this expression, and elevate it to an operator by replacing the particle positions q i (t) by the respective functional derivative with respect to J q i (t), We thus obtain the density operator which is a sum of one-particle density operatorsρ i (k, t).

Density correlation functions
Since the operatorsρ i (k) contain a functional derivative in the exponential, they create a translation. Thus, any sequence of n one-particle density operators applied to the generating functional gives, setting J = 0 at the end, where the short-hand notationρ( j) =ρ(k j , t j ) was introduced. The shift tensor is Applying finally n density operators to the generating functional results in where the first equality in the second line holds because the particles of the ensemble are indistinguishable, so each of the N(N −1) · · · (N −n+1) particle tuples ( j 1 , . . . , j n ) must contribute the same statistical result. The final approximate equality holds if N n, as will usually be most safely the case. If the thermodynamic limit N → ∞ could be taken or would be inappropriate, shot-noise terms could occur [49]. In cosmology, however, we can assume that any reasonably large volume will be filled with an extremely large number of particles. We shall therefore replace (94) by an equality from here on.

Low-order statistics of the free particle distribution
After this general discussion, we shall now specify the results obtained so far to loworder statistical measures for the free particle distribution, i.e. for the distribution of the particle ensemble flowing along force-free trajectories.

Free power spectrum
For power spectra, we have with n = 2 in (94) For synchronous power spectra, t 2 = t 1 , and the tensor L in (95) is such that the scalar product (L,x(t)) turns out to be If we neglect the contribution to the particle trajectories (54) due to the particle interactions and insert the Zel'dovich trajectories (56) into (97), we obtain from (23) what we call the free generating functional, where (q i , p i ) are now meant to be initial particle positions and momenta. Due to spatial homogeneity, we can refer the positions of all particles to the position of any particular particle, say particle 1, replacing q i → q i −q 1 , and integrate over q 1 . Since the momentum-correlation matrix depends only on particle separations, but not on absolute particle positions, this results in a delta distribution, leaving the free generating functional in the form where the Gaussian integral measure now does not contain dq 1 any more, the momentum integration in (99) can be carried out directly, leading to (103) KFT for cosmic structure formation Combining (102) with (79), the quadratic form in the exponential of (103) is with the definition a (q 2 , µ) = a 1 (q 2 ) + µ 2 a 2 (q 2 ) (105) containing the cosine µ =k ·q 2 of the angle between k 1 and q 2 (see eg. [50,52]). Notice the limit of the momentum-correlation function a (q, µ).
Since the integrand in (103) depends only on the (relative) position q 2 , we can integrate over all initial particle positions q j with j > 2, which results in a factor V N−2 . We can further pull part of (107) in front of the integral, arriving at where q now abbreviates q 2 for simplicity, and Q D is Returning to (105), the free two-point density correlator is now As shown in (16), the free power spectrum P(k) is defined in terms of this expression, except for the preceding delta distribution and the prefactorρ 2 = (N/V) 2 , Combining (110), (109) and (107), we thus find It is evident from the argument of the first exponential in the integrand of (111) that multiplying the amplitude of the initial power spectrum, and thus the correlation function a , by a certain factor is equivalent to leaving this amplitude unchanged but multiplying the time coordinate by the root of the same factor. This will later be reflected by our results containing the time coordinate only in combination with σ 2 , suggesting to introduce τ 2 = tσ 2 as the time coordinate relevant for structure formation.
The exponential prefactor exp(−Q D ) seems to indicate that the spectrum will be damped exponentially on small scales, but we will show later that it will be compensated exactly on small scales. If the argument of the exponential in the integrand of (111) is small enough, we can approximate where we have used in the last step that our time coordinate t is the linear growth factor. Since a is bounded, this approximation holds for sufficiently small k 2 t 2 , i.e. for any given time on sufficently large scales, as expected [50].

Free bispectrum
For the bispectrum, we proceed in an exactly analogous way (see eg. [49]). The shift tensor L contains three wave vectors (k 1 , k 2 , k 3 ) now, thus (97) is replaced by Consequently, the free generating functional, evaluated at L, is extended to and (116) Similar to (110), the free bispectrum B is defined by for non-degenerate configurations of the three wave vectors k 1,2,3 , i.e. if none of them vanishes. The free bispectrum depends on only two wave vectors because the delta distribution in (114) ensures that k 1 = −(k 2 + k 3 ). Taking account of (94) with n = 3, then comparing (117) to (114), we find the expression understanding again that the wave vector k 1 in Q (3) C and Q (3) D is fixed by the condition In the large-scale limit, when the norm of all wave vectors involved is getting small, k 1 , k 2 , k 3 → 0, the expression resulting from (114) can be brought into the familiar form

Free velocity power spectrum
The generating functional of kinetic field theory contains the complete statistical information on the phase-space trajectories of the particle ensemble and thus also allows calculation velocity power spectra [52]. The momentum p j of a particle j needs to be localized at the spatial position q j by the expression  Fig. (9)).
which we elevate to the velocity operator The density operator locates the particle j and the derivative with respect to the source field J p j extracts this particle's momentum.
Since the derivative with respect to J q j contained in the density operator and the derivative with respect to J p j commute, we can apply multiple operatorsΠ 1 . . .Π n by shifting all involved density operatorsρ 1 . . .ρ n to the right to apply them first to the generating functional, which they will translate by an amount L, If we specialize to synchronous spectra, all operators are applied at the same time t j = t. The free generating functional Z 0 [J + L] can then be brought into a form similar to (98), with the components of the shift tensor L, generalizing the definition of L p in (102). Expression (124) shows that we can extract information on particle momenta p j by taking derivatives with respect to L p j rather than J p j . Therefore, for a velocity power spectrum with n = 2, analogous to (109). Taking the free generating functional from (103) and defining the velocity power spectrum P Π (k) in analogy to (110) results in which is a second-rank tensor because it correlates all momentum components at one position with all momentum components at another position. Here, D j = ∂/∂L p j abbreviates the derivative with respect to L p j and Q is the quadratic form The derivatives of Q with respect to L p 1,2 are while the second derivative is Thus, the velocity power spectrum reads with the matrix 4 Asymptotic small-scale behaviour of power spectra We shall now turn to deriving rigorous statements on the small-scale behaviour of the main quantities derived above, in particular the free density-fluctuation power spectrum, the free density-fluctuation bispectrum, and the free velocity power spectrum. We shall then generalize some of our results by including the complete set of initial correlations and by including interactions in the mean-field approximation. at redshift z = 0 is shown together with P (ini) δ k 2 (light gray line), which matches the trace at large scales, and the free power spectrum, multiplied by σ 2 1 (dark gray line), matching the trace at small scales. The initial power spectrum is for WIMP dark matter with a Gaussian cut-off at k s = 10 6 Mpc −1 .

Free density-fluctuation power spectrum
We shall focus on density-fluctuation power spectra first, deriving their asymptotic small-scale behaviour in two different ways [55].

Leading-order asymptotics from Morse's lemma
We are interested in general statements on the power spectrum on small scales. For this reason, we first study the asymptotic behaviour of the free power spectrum P(k) as derived in (111).
We begin with Laplace's method for d-dimensional integrals of the form over functions f, g ∈ C ∞ on a domain Ω ⊂ R d . The integral is supposed to converge absolutely for sufficiently large λ ∈ R, the function f is assumed to have a minimum at x 0 ∈ Ω and only there, and the Hessian A of f in x 0 is supposed to be positive definite. Then, J(λ) has the asymptotic expansion for λ → ∞ (see eg. [63,64]). For specifying the coefficients c n , we introduce some elements of notation. Under the given conditions, Morse's lemma ensures that neighbourhoods U, V of y = 0 and x 0 and a diffeomorphism h : U → V exist such that With the Jacobian determinant det H of h, we define the function G : U → R by We further introduce the multi-index α = (α 1 , . . . , α d ) and agree on the notation where µ in the last equation is a d-dimensional vector. Moreover, we define the symbol and the derivative operator Then, the coefficients c n are given by Applying the result (134) to integrals of the Laplace-Fourier type, and resumming the coefficients, we find for λ → ∞. Next, we specialize this statement to integrals of the form with s ≥ 2. We have shown in [55] that although the kernel does not meet the aforementioned conditions, the theorem can still be applied. We thus obtain to leading asymptotic order for |k| → ∞.
Comparing the integrals in (143) and (111), we now set s = 2, d = 3 and f (q) = t 2 a (q), use the limit (106) of a (q, µ) and the Hessian Since the inverse of the A is and its determinant is we immediately find with Recently, this result has also been derived in the framework of Lagrangian perturbation theory [65]. This is the leading-order asymptotic term, from which important conclusions can be drawn. Before we get to them, we derive the full asymptotic series for the free power spectrum in a different manner from Erdélyi's theorem.

Asymptotic series from Erdélyi's theorem
We begin with the integral (111), introduce spherical polar coordinates withk as the polar axis, and introduce a finite upper limit q max > 0 for the integration over q, for k → ∞. This is asymptotically correct since the asymptotic behaviour of P for large wave numbers is determined by the behaviour of the exponent next to its critical point, so the contribution to the integral from q max to ∞ can be ignored [55,66]. Then, we can use Erdélyi's theorem (see eg. [67][68][69]) to derive the asymptotic series for P [55].
This theorem states that one-dimensional integrals of the form over functions f and g admitting the asymptotic series have the asymptotic expansion with the coefficients provided the function f has a global minimum at a, the function f can be term-wise differentiated, f and g are continuous in a neighbourhood of a, except possibly at a itself, and I(λ) converges absolutely for sufficiently large λ.
We need to set f (q) = a (q) here with a (q) from (105) and a 1,2 (q) from (84). Using the series expansions for the spherical Bessel functions, we find first with the moments σ 2 n from (87). These imply the asymptotic series Setting g(q) = q 2 exp(ikqµ), we further find Comparing (157) and (159) to (152), we can read off α = 2 and β = 3. Inserting the coefficients a 2m and b m into (154), setting ν = (n + 3)/2, using (153) and integrating the resulting expressions over µ, we can derive the complete asymptotic series The two lowest-order terms are The leading-order term in (160) reproduces the result (148), as it should. Generally, the functions P (m) (t) are proportional to the moments σ 2 n of the initial power spectrum, and also depend on lower order moments. Explicit expressions for the coefficient functions P (m) (t) are given in [55].

Conclusions from the asymptotics of the free power spectrum
Our conclusions on the asymptotic behaviour of the free power spectrum rest upon rather general assumptions. We have assumed that the initial momenta of the particles in our ensemble are drawn from a Gaussian random velocity-potential field and correlated in such a way as to satisfy the continuity equation between the initial density and velocity fields. For deriving the asymptotic series of the free power spectrum, we did not have to specify the shape of the initial power spectrum, but only had to assume that its moments σ 2 2 for the leading-order and σ 2 3 for the next-to-leading order terms exist.
Moreover, none of our results obtained so far depends on anything specific for cosmology. We have transformed the Hamiltonian equations of motion to the expanding spatial background, which resulted in a time-dependent particle mass and a specific form of the Poisson equation. Having focussed on the free power spectrum, however, we could describe the particle trajectories purely kinematically, without invoking any specific dynamics. All we have assumed in this regard is that a time coordinate t exists in terms of which particle trajectories take on the inertial form (56). The formation of the asymptotic k −3 tail of the free power spectrum is thus an effect of collective free streaming of classical particles with phase-space positions drawn from an initially Gaussian random field, and the exponent −3 is set solely by the number of spatial dimensions.
Most importantly, the exponential damping factor exp(−Q D ) appearing in (107) is exactly cancelled in the asymptotic terms. Freely-streaming particles with correlated initial momenta thus do not lead to exponential damping of small-scale structures.
The amplitude P (0) (t) of the leading-order asymptotic term starts at zero for t = 0, reaches a maximum value of when Σ max = 3/2 and then decreases again. The increase is due to the fact that freelystreaming particles create structures where their flow is locally convergent, while the decrease is due to the fact that they fly past each other and erase these structures again after they pass the point of convergence. The values of Σ max and P (0) max have an absolute meaning, irrespective of the cosmological background and the shape of the initial density-fluctuation power spectrum. The definition of Σ in (149) shows that τ 2 = tσ 2 is the relevant time coordinate for structure formation by collective streaming. The lower the moment σ 2 of the initial power spectrum is, the more time t it will take the asymptotic k −3 tail to reach its maximum amplitude, with This time scale t max is expected to set an important scale for structure formation. The ratio of the asymptotic terms of the next-to-leading and the leading order is When this value drops below unity above a certain wave number k 0 , the free power spectrum attains its asymptotic behaviour ∝ k −3 for scales smaller than k 0 .

Free density-fluctuation bispectrum
We now turn to the asymptotic behaviour of the free bispectrum, given in (118).
For finding an expression valid on small scales, an approach based on Laplace-type integrals is possible. For the lowest order term of an asymptotic series, we require an approximate expansion for the exponent Q (3) C around its relevant critical point. While C has a critical point at zero, its Hessian with respect to q is non-invertible at q = 0 [70]. The previously applied method based on Morse's lemma thus fails, but it can be suitably adapted using the splitting lemma of functions, which extends the validity of the Morse lemma to cases of Hessians with positive corank (see eg. [71]). The splitting lemma states that, if f is a polynomial of order ≥ 2 with a critical point at x 0 and a Hessian at x 0 with rank k, then with either g = 0 or g a polynomial of order ≥ 3 which is uniquely determined up to an equivalence transformation, provided the critical point is isolated. Under this assumption, the proof of the splitting lemma can be adapted to the specific form of Since we are studying a statistically isotropic situation, we can without loss of generality orient the coordinate system such that k 2 points into e z direction and k 3 falls into the x-z plane. Then, the only non-vanishing components of the wave vectors k 2,3 are k 2z , k 3x and k 3z , in terms of which the asymptotic behaviour of the bispectrum is given by for k 1 , k 2 , k 3 → ∞ with the coefficients c 0 = 9000 4 √ 6.3 π 5/2 Γ(5/4) , as shown in [70]. In the triangle formed by k 1,2,3 , let α and β be the angles between k 2 and k 3 and between k 1 and k 3 , respectively. Then, by the sine theorem, the norm of k 3 is and its components are k 3x = −Ak sin α and k 3z = −Ak cos α. We can then bring the function c 4 (k) into the form allowing us to write the asymptotic expression for the bispectrum as Values of the function f 1/4 (α, β) are tabulated for some special cases in Tab. 1. Table 1 Values of the function f 1/4 (α, β) for some special configurations of the triangle k 1,2,3 . case

Free velocity power spectrum
The asymptotic expansions (156) of the functions a 1,2 (q) show that thus terms proportional to k 2 show up only at O(q 4 ) in the terms in parentheses in (131). The dependence on the wave number k of the leading-order asymptotic term of the velocity power spectrum P Π (k) for k → ∞ thus remains unchanged compared to that of the density-fluctuation power spectrum P(k). Only its amplitude changes  because we need to replace the coefficient b 0 in the asymptotic expansion (159), and thus also the coefficient function P (0) in (160), as This takes us to the leading-order asymptotic expression for k → ∞, with The time dependences of the leading-order asymptotic terms in the densityfluctuation and the velocity power spectra are thus the same, their maximum amplitudes are reached at Σ max = 3/2, but the maximum amplitude of the leading-order asymptotic term in the velocity power spectrum is compare (163). The velocity correlation on small scales is negative and has the same amplitude as the density-fluctuation power spectrum multiplied by the initial velocity dispersion σ 2 1 . This explains the origin of the k −3 asymptotics, i.e. structures on small scales: initially convergent particle streams cross which leads to caustics.

Density-fluctuation power spectrum for interacting particles
So far, we have neglected particle interactions altogether. This does however not mean that no gravity was included because the Zel'dovich inertial trajectories are subject to the large-scale, linear part of the gravitational interaction, as specified in Sect. 2. We shall now proceed to include particle interactions in the mean-field approximation [54].

Forces in the mean-field approximation
So far, we have neglected particle interactions beyond those that are already contained in the reference trajectories modelled with the Zel'dovich approximation. The density-fluctuation two-point function written down in (95) is still exact, however. The shift tensor L, multiplied with the actual particle trajectories, is according to (97), taking into account that k 1 + k 2 = 0 due to statistical homogeneity. The spatial trajectoriesq j (t) can be split into the free partq (0) j (t) described by the Zel'dovich approximation, and a partȳ(t) describing the deviations from the Zel'dovich reference trajectories caused by the particle interactions, q j (t) = q j + tp j We are quoting the result (54) here which contains the Hamiltonian propagator g H given in (48), and the gradient of the potential φ which satisfies the Poisson equation (55). In (98), we continued withq (0) 1,2 , neglectingȳ 1,2 . We shall now includeȳ 1,2 in a mean-field approximation. For completing the generating functional Z[L], we require the scalar product We now approximate the potential gradients ∇ j φ by a mean-field average, which we construct in the following way. Let V(q) be an interaction potential linearly superposed by contributions v(q) due to individual particles at positions q i , then because the sum of delta distributions is the (number) density ρ of the particles. The potential gradient at the position q j of another particle is where ρ j (q) is the contribution of particle j to the particle number density. Thus, the potential gradient acting on particle j is where v is the interaction potential between a pair of individual particles separated by q − y. The potential gradient contributed by particles at a distance q − y from particle j is In the cosmological situation, we imagine a test particle moving through a collection of other particles, exposed to their collective gravitational field. The number density of these particles can be considered to be arbitrarily high. The force exerted by any individual particle on the test particle is then an arbitrarily small contribution to the total force. Fluctuations of this force caused by individual particles can then be ignored. This situations is typically well described by a mean-field approximation.
Averaging the potential gradient ∇φ leads to containing the two-point function of the (number) density field. By definition of the two-point density auto-correlation function ξ, wheren is the mean number density. The division by N takes into account that we cross-correlate the one-particle density contribution ρ j with the density ρ. Inserting (187) into (186) leaves because the contribution by uncorrelated particles averages to zero in a statistically homogeneous and isotropic random field. Now, both factors in the integrand of (188) depend on the fixed separation q − y only, but not on the point q any more. This is a consequence of statistical homogeneity. The integral over q thus only results in a volume factor which, together with the prefactor N −1 , cancels one of then factors. Thus, the mean potential gradient is By the Fourier convolution theorem, the Fourier transform of the average potential gradient is the convolution of the Fourier transforms of the individual factors. The Fourier transform between particles is given by (64), approximated by (65). Thus, the Fourier transform of its gradient is The Fourier transform of the two-point density auto-correlation function ξ is the density-fluctuation power spectrum P δ (k). Since we select the motion of particles along Zel'dovich inertial trajectories (56) as a reference, the most appropriate power spectrum to insert here would be the free power spectrum P(k). We shall further simplify this choice below and keep the symbol P δ (k) for now. We thus find the Fourier transform of the mean potential gradient to be We thus replace the expressions k 1 · (∇ 1 φ − ∇ 2 φ) appearing in the scalar product (180) by the averaged expression where we have used that ∇φ (−k 1 ) = − ∇φ (k 1 ). For carrying out the convolution, we introduce the cosine µ of the angle between k 1 and k and define κ = |k |/|k 1 | as well as κ 0 = |k 0 |/|k 1 |. This enables us to write where σ 2 J is the moment of the density-fluctuation power spectrum with the filter function

Density-fluctuation power spectrum in the mean-field approximation
We can now write the mean-field averaged scalar product between L andȳ as Since this expression does not depend on the initial phase-space coordinates of the particle ensemble any more, it can be pulled in front of the integral over dΓ in the generating functional. This leads to containing the mean-field interaction term Invoking (109) and (110) once more, the non-linear, density-fluctuation power spectrum, including particle interactions in the mean-field approximation, can be written as P (nl) δ (k) = e S I (k) P(k) .
(199) This is a convenient and simple expression whose merits need to be assessed by comparing it with, e.g., density-fluctuation power spectra derived from numerical simulations. In [54], we have further simplified it by replacing the free, non-linear power spectrum P by the linearly evolved power spectrum P (lin) δ , and the power spectrum P δ (k) to be inserted into the moment σ 2 J given in (194) by a suitably damped for the mean-field averaged, non-linear power density-fluctuation power spectrum, which reproduces the results from numerical simulations remarkably well. This result may give sufficient credit to the mean-field approximation for the particle interactions.
In our present context, we are aiming at a different conclusion, however. As (198) shows, the scale dependence of the mean-field interaction term is determined by the moment σ 2 J , which we now repeat in a different way, The filter function J(k/k 1 , k 0 /k 1 ) flattens for k < k 1 at a level decreasing with increasing k 0 /k 1 , and falls off ∝ k −2 for k > k 1 . For k 0 /k 1 → 0, the filter function approaches J → 2 for k k 1 . In the small-scale limit, k 1 → ∞, the filter function is well approximated by J ≈ 2 almost everywhere in the integration range in (201). Then, the integration covers the entire power spectrum, and the result becomes scaleindependent. While the power spectrum in (201) is growing linearly, comparing with (87) shows that σ 2 J (k 1 ) → 2t 2 σ 2 2 (202) for k 1 → ∞. The mean-field interaction term thus becomes independent of scale in the small-scale limit, which implies that the asymptotic behaviour of the non-linear power spectrum P (nl) δ (k) for k → ∞ will be the same as that of P, while its amplitude will be enhanced compared to (148) or the identical expression in (161) by the exponentiated mean-field interaction term.

Lifting limitations
Up to this point, we have made several simplifying assumptions. In particular, we have assumed that the moments σ 2 n of the initial density-fluctuation power spectrum exist up to the order necessary, and we have neglected density-density and densitymomentum correlations in the initial state. We shall now show how these assumptions can be lifted, and that they do not change our conclusions about the asymptotic behaviour of the density-fluctuation and velocity power spectra.

Asymptotic behaviour for strictly Cold Dark Matter
In deriving the leading-order asymptotic behaviour of the free power spectrum, we assumed that the second moment σ 2 of the initial velocity potential power spectrum exists as defined in (87). This is a save assumption also for cold dark matter since measurements of the spectral index n s return values smaller than unity. For this reason, the initial cold dark matter power spectrum, which is assumed to be accurately described by linear theory, has a tail at large wave numbers falling off like k n s −4 log 2 k. The asymptotic series that we derived above from Erdélyi's theorem, however, requires that the moments σ 2 n of arbitrary high order should exist. This is indeed only the case for power spectra that are exponentially cut off above some possibly very large wave number k s , i.e. if structures smaller than a typical length scale of k −1 s do not exist at the initial time. Should this not be the case, and it is unknown whether arbitrarily small darkmatter structures can be formed initially, we need stronger methods than previously applied to study the asymptotic behaviour of the free power spectrum. Such methods are provided by analytic continuations of inverse Mellin transforms. Moreover, it turns out that the asymptotics of a power spectrum for strictly cold dark matter correctly describes the intermediate regime of dark matter with large k s , which is not accessible from our previously derived asymptotic series due to its divergent nature. In this section, we begin with an assumed wide class for the asymptotic behaviour of the initial density-fluctuation power spectrum and use the Mellin-transform technique to derive the asymptotic behaviour of the momentum correlation function a (q) for q → 0. We then insert the resulting expression into the free, non-linear power spectrum P from (111) and obtain its asymptotic behaviour in the small-scale limit.
We thus begin by assuming that the initial density-fluctuation power spectrum admits an asymptotic expansion of the form for k → ∞ with real coefficients c mn . This represents a very wide class of asymptotic behaviour (see eg. [72,73]). Next, we need the Mellin transform of a function f , which is defined by where the integral typically converges on a strip of the complex plane with a < Re z < b. Equipped with the Mellin transform, we study integral transforms of functions f with a kernel h, called H-transforms and defined by Under sufficently general conditions, allowing I to be not absolutely convergent, integral transforms of this type can be expressed by Mellin transforms of the functions involved, where c is a real number falling into the common strip of analyticity of both Mellin transforms. Defining the function and assuming G can be analytically continued to the right half-plane, the integral I(q) can be expanded asymptotically as for q → 0 + , where G(z) has no pole for Re R such the remaining integral is o(q R ) (see eg. [63]). In order to derive the asymptotic behaviour of a (q) on small scales, we apply this method to the integrals (84) and obtain the results for q → 0, where M (n− j) denotes the (n − j)-th derivative of the Mellin transform with respect to the function argument. The detailed derivation can be found in [SK: 2022Konrad]. The moments σ 2 n of the initial velocity potential can be expressed by the Mellin transform as which also allows us to identify the analytic continuation of these moments for nonconverging integrals. Since a (q) = a 1 (q) + µ 2 a 1 (q), we conclude from (210) that we can write as q → 0, where the function ξ(µ 2 , log q) is implicitly given by (210).
We can now insert this asymptotic, small-scale expression for a (q) into the free power spectrum (111) and use a suitable Taylor expansion of the exponential to arrive at the asymptotic expansion with time-dependent coefficient functions P nm (τ 2 ) to be specified, turning into constants at late times, τ 2 1 [74].

Including all initial correlations in the density-fluctuation power spectrum
We have so far simplified the probability distribution P(q, p) from (75) by approximating the differential operatorD by (77), leading to the Gaussian (78) in the momenta p. We shall now demonstrate that this approximation is very well justified at late times [75]. With (75), the free generating functional (23) reads For calculating a power spectrum, we apply two density operators directly and integrate over q 1 as before. This results in = N 2 (2π) 3 δ D (k 1 + k 2 )D dqˆ1dp s Φ(R, s)e i(s+L p )· p+ik 1 ·q 2 (215) with the shift tensor L p from (102). The integral over p results in a delta distribution δ D (s + L p ) which allows carrying out the s integration immediately. This leads to ρ(1)ρ(2) = N 2 (2π) 3 δ D (k 1 + k 2 )D dqˆ1Φ(R, L p )e ik 1 ·q 2 .
In view of our further calculation, it is important to note that L p contains entries for two particles only, conveniently labelled as particles 1 and 2. If we setD = V −N1 and R = 0 here and inserted the characteristic function Φ from (69), we reassuringly returned to (103). The characteristic function, evaluated at (R, L p ), reads The density-density auto-correlation and density-momentum cross-correlation functions are determined by the initial density-fluctuation power spectrum, The dimensionless free power spectrum k 3 P for WIMP dark matter with small scale smoothing wave number k s = 10 6 Mpc −1 (purple line) at redshift z = 0 is shown together with the various asymptotics that we derived. At large scales, i.e. up to wave numbers k ≈ 0.02 Mpc −1 , k 3 P is well described by the damped linearly evolved power spectrum (gray line). At small scales, above k ≈ 2 · 10 4 Mpc −1 , k 3 P becomes constant, accurately described by the k −3 asymptotics (green line). Below k 1.2 · 10 4 , the asymptotics of the free power spectrum of cold dark matter (golden line) aligns with the initially smoothed spectrum, where we chose M = 1000. Note, that this is remarkable, as this line contains only information of the unsmoothed spectrum up to q 3−ns order of the initial momentum correlation functions a 1 and a 2 . The crossing point of the k −5 asymptotics (blue line) with the k −3 asymptotics marks the transition from k −3 asymptotics to the CDM asymptotics. To summarize, strictly CDM Zel'dovich power spectra approximate the intermediate regime of WIMP Zeldovich power spectra.
Since the differential operator acts on R only, we can continue writing (216) in the form The differential operatorD, given by (66) and (76), iŝ Since M is quadratic in R, at most second-order derivatives of M with respect to R can appear. The differential operatorD applied to M will thus contain the two types of term a 2 (q) q 2 asymptotics −q 3−n s asymptotics sum q 2 and q 3−n s asymptotics Fig. 23 The functions a 1 + σ 2 1 3 (left, purple line) and a 2 (right, purple line) with their q 2 (green lines) and q 3−ns (blue dashed lines) asymptotics are shown. The sum of the two asymptotic orders (golden lines) shows that we indeed need the q 3−ns order to accurately describe those functions on relevant scales.
only. If j 1, 2 in the first type of term, the correlation matrix C δ j p depends only on relative coordinates that appear nowhere else in the integrand. Integrating −iC δ j p L p over the spatial coordinates then results in zero. Likewise, if j, k 1, 2 in the second type of term, the subsequent integration over C δ j δ k vanishes. Therefore, only those contributions toDM remain which contain ∂ R 1 or ∂ R 2 . We can thus replaceDe M in (219) bŷ Furthermore, since C δ j p j = 0, and likewise for C δ 2 p L p . Therefore, we havê with The expression for C δp in (218) shows that C δ j p k = −C δ k p j because the exchange of indices implies q → −q. Moreover, we know from (102) that L p 2 = −L p 1 = k 1 t. Thus, we can somewhat simplify F(k, q, t) to read C δδ 2 |C δ p | |C δ p | 2 Fig. 24 The terms that enter the function F as defined in (226). The initial density-density correlation term (C δδ , purple line) enters linearly. The linear density-momentum correlation term (2C δp , green line) enters linearly in kt, while the quadratic density-momentum correlation term (C d p 2 , blue line) enters quadratically in kt.
Integrating finally over the free particle positions q 3 . . . q N and comparing (219) to (110), we can identify the expression for the free power spectrum including not only momentum auto-correlations, but all correlations between density and momentum fluctuations [75]. The expressions (218) for C δδ and C δp and the series expansions (155) for the spherical Bessel functions lead to the asymptotic expansions for q → 0. Therefore, the factor (1 + F) has the asymptotic expansion The modified asymptotic amplitude P (0) , where initial density-density, density-momentum and momentum-momentum correlations are considered (see (227)), are shown as a function of scale factor a in colored lines for three different initial small-scale smoothing wave numbers k s . The comparison to the asymptotic amplitudes where only initial momentum-momentum correlations C pp are considered shows that the other correlations have a barely noticeable impact on the small scale structure evolution.
Searching for the effect of the factor (1 + F) on the leading-order asymptotic behaviour of the free density-fluctuation power spectrum, we thus need to multiply the asymptotic expansion for the function g in (159) by each term in the asymptotic expression (229) and change m to m + 1 and m + 2, respectively. The dependence on the wave number of the leading order terms remain unchanged. The leading-order asymptotic behaviour of the free power spectrum is thus given by P(k) ∼ 1 + 4 9 σ 2 2 − 10 3t + 25 9t 2 with the coefficient function P (0) from (161), and the k −3 behaviour remains [75].

Including all initial correlations in the velocity power spectrum
As shown before, taking the initial density-density and density-momentum correlations into account in addition to the momentum auto-correlations inserts a factor (1 + F) into the integrand of the free power spectrum, with the function F given by (225). Therefore, the velocity power spectrum taking the complete set of initial correlations into account is obtained by inserting (1 + F) into (127), P Π (k) = −D 1 ⊗ D 2 dq 1 + F(k, q, t) e −Q e ik·q = dq F e −Q e ik·q (231) with the matrix As before, we need the leading asymptotic order in q of the matrix F , which is implying that the velocity power spectrum including all initial correlations must have the leading-order asymptotic behaviour with P (0) Π from (176).

Conclusions
We have reviewed kinetic field theory for classical particle ensembles, putting particular emphasis on the evolution of cosmic structures in collision-less dark matter. Kinetic field theory dissolves the cosmic density field into particles subject to Hamiltonian dynamics and studies the evolution of an initial phase-space probability distribution under the Hamiltonian flow. Compared to other analytic approaches to cosmic structure formation, the essential advantage of kinetic field theory is that trajectories in phase space do not cross. The notorious shell-crossing problem occurring inevitably in methods building upon uniquely valued density and velocity fields in configuration space is thus avoided by construction. The central mathematical object of kinetic field theory is a generating functional encapsulating the statistically defined initial state of the particle ensemble together with the dynamics of phase-space trajectories. Kinetic field theory does not set up and solve a dynamical equation for smooth density or velocity fields or a phasespace distribution function. Rather, it evolves this generating functional in time and allows extracting statistical information on the evolved particle ensemble by functional derivatives. Formally, it resembles a statistical quantum field theory, however some aspects of it are considerably simpler.
After the introduction in Sect. 1, we set up in Sect. 2 the generating functional of kinetic field theory in an expanding cosmic space-time and chose a suitable Green's function to solve the Hamiltonian equations of motion. We introduced the growth rate of linear density fluctuations as a time coordinate and argued that inertial motion with respect to this time, the so-called Zel'dovich approximation, is particularly appropriate for describing cosmic structure formation. We emphasized that the gravitational potential between particles on Zel'dovich inertial trajectories is sourced only by nonlinear fluctuations of the cosmic matter density, which implies that the effective gravitational potential is short-ranged and is approximately of Yukawa form.
We defined density operators in Sect. 3, showed how they can be used to extract statistical information from the generating functional, and derived low-order statistical measures for distributions of particles freely streaming along Zel'dovich reference trajectories. We emphasize once more that this kind of reference motion, even though being referred to as free, does include gravitational interaction at early times, and long-range gravitational interaction at all times. In particular, we derived equations for the non-linear, free power spectrum of cosmic density fluctuations and their bispectrum as well as the free velocity power spectrum.
Aiming at rigorous statements on the small-scale behaviour of low-order statistical measures, Sect. 4 is the core of this paper. There, we used extensions of Laplace's method to derive the small-scale asymptotics of the free density-fluctuation power spectrum, the free bispectrum, and the free velocity power spectrum. Our central results obtained there are that 1. the free density-fluctuation power spectrum and the free velocity power spectrum asymptotically fall off proportional to k −3 for wave number k → ∞, and that 2. the free bispectrum asymptotically falls off like k −11/2 .
These results assume only that the initial particle momenta are drawn from a Gaussian random field whose power spectrum has finite low-order moments. The nature of the dark matter and the cosmological model are irrelevant. The k −3 tail of the free power spectra and the k −11/2 tail of the bispectrum evolve by collective streaming of the matter particles, and the exponents are set by the number of spatial dimensions only.
These results persist if some simplifying limitations are given up. We have shown further that gravitational interaction between particles in the mean-field approximation and strictly cold dark matter with an initial power spectrum without small-scale cut-off do not affect the asymptotic k −3 dependence of the matter-fluctuation power spectrum. In addition, including initial density-density and density-momentum correlations together with momentum auto-correlations do not change the asymptotic behaviour of the density-fluctuation and the velocity power spectra either.
It would thus appear that the asymptotic behaviour of low-order statistical measures of the cosmic matter and velocity distribution is a consequence of the initial particle momenta being drawn from a Gaussian random field, without further assumptions entering. Late-time, non-linear density-fluctuation power spectra would therefore necessarily develop a k −3 tail on small scales, irrespective of the cosmological model and the nature of the dark matter.
Characteristic and universal time scales for structure formation can be derived from the leading-order asymptotic term reaching its maximum amplitude, and from the next-to-leading order term falling below the leading order term. Possible observational consequences and the potential significance of these time scales need to be worked out.
The k −3 behaviour of power spectra implies that the power, i.e. the product of the power spectra times the number of Fourier modes, flattens off and becomes constant on small scales. Each fixed scale interval will then contribute the same amount of power. This result should provide clues for the universal density profile of gravitationally bound structures. Of course, many questions remain to be addressed and answered, but the conclusions on the small-scale behaviour of non-linear cosmic structures obtained from kinetic field theory here appear promising.