Kinetic Field Theory: Higher-Order Perturbation Theory

We give a detailed exposition of the formalism of Kinetic Field Theory (KFT) with emphasis on the perturbative determination of observables. KFT is a statistical non-equilibrium classical field theory based on the path integral formulation of classical mechanics, employing the powerful techniques developed in the context of quantum field theory to describe classical systems. Unlike previous work on KFT, we perform the integration over the probability distribution of initial conditions in the very last step. This significantly improves the clarity of the perturbative treatment and allows for physical interpretation of intermediate results. We give an introduction to the general framework, but focus on the application to interacting $N$-body systems. Specializing the results to cosmic structure formation, we reproduce the linear growth of the cosmic density fluctuation power spectrum on all scales from microscopic, Newtonian particle dynamics alone.


Introduction
In a 1988 paper [1], E. Gozzi outlined how the successful path integral formalism of quantum field theory can be used to describe classical systems as well. In the following years, Gozzi further developed and studied this approach in collaboration with M. Reuter and W.D. Thacker [2,3] and others [4,5]. Relatively recently, G.F. Mazenko and S.P. Das [6,7] developed a perturbation scheme and applied the formalism to phase transitions in glass. Based on this, the group of M. Bartelmann began to apply the path integral approach of classical mechanics to the study of cosmic structure formation [8]. Since then, the approach has been subject to significant further development [9,10,11,12,13] and is now known as Kinetic Field Theory (KFT) [14].
Conceptually, KFT mimics the approach of numerical simulations to describe the time evolution of an interacting N-particle system. Via use of path integrals over phase space trajectories, one keeps track of the positions and momenta of all particles over time. Mathematically, this information is stored in the generating functional Z of the theory. Collective quantities like, e.g., the density field can be extracted via functional derivatives of Z. Explicit expressions for the generating functional are obtained either via a perturbative approach [8] or through the use of suitable approximations [9].
There is one significant difference between the formalism of KFT and numerical N-body simulations.
Namely, the latter always have to start out with some explicit initial conditions, i.e., initial positions and momenta for all particles. A priori this is also true for KFT. However, in an analytic setting we can abstractly integrate over a probability distribution of initial conditions to obtain expectation values. This would be prohibitively expensive for numerical N-body simulations which therefore have to resort to other methods of extracting macroscopic quantities. The integration over initial conditions in KFT yields an ensemble average which in many applications is the observable we are interested in.
In this article we revisit to the initial perturbative approach to KFT in [8]. We give a self-contained and pedagogical introduction to the framework of KFT in a general setting and then specialize it to interacting N-body systems. Unlike the perturbative treatment in [8], we perform the integration over the initial conditions not on the level of the generating functional, but on the level of observables. This allows us to interpret intermediate results describing observables of a system with explicit initial conditions. Upon integrating them over a suitable probability distribution, we obtain statistically averaged expectation values, i.e., macroscopic observables describing the collective system. While ultimately mathematically equivalent, this change in procedure makes the perturbative treatment significantly more transparent and simpler.
Our ultimate aim is to supply a simple and consistent procedure for obtaining observables in cosmic structure formation, in particular the non-linear density fluctuation power spectrum. The initial probability distri-bution for this particular application is relatively complicated. For the sake of self-containedness and clarity, we repeat the discussion from appendix A of [8]. In particular, we explain how the integration over initial conditions fits together with the perturbative treatment. Having done so, we will be able to obtain analytic parameter-free perturbative expressions for the non-linear density fluctuation power spectrum. These expressions are expanded in both the initial correlations and the interactions. More precisely, on the one hand we expand the density and momentum correlations imprinted in the initial conditions of the N-particle system and on the other hand we treat the particle interactions perturbatively by expanding the forces relative to the free particle trajectories.
Surprisingly, we find that the expansion in the initial correlations is significantly more difficult than the expansion in the interactions. Indeed, for higher-order correlations our expressions acquire integrals over wave vectors which are challenging to solve even numerically. Conversely, each interaction order merely supplies a simple time integral. To conclude this article, we calculate the density fluctuation power spectrum in first order (i.e. linear) correlations and up to very high interaction order. Combining the contributions for linear correlations and arbitrarily high interaction orders precisely reproduces the linear growth of the power spectrum at all scales. This reveals that any non-linearities in the power spectrum are solely due to higher-order correlations in the initial conditions and thus can be studied completely independently of the linear evolution. Results for higher-order correlations will be presented in a future article [15].

KFT Framework
Let us consider a generic classical physical system. Denote its space of possible states by X, such that any point x ∈ X uniquely characterizes a state of this system. 1 In practice, we usually take X to be phase space consisting of all (generalized) coordinates and associated momenta. Any evolution of such a system can be described as a map ϕ : R → X which associates to any point in time t ∈ R a state of the system ϕ(t) ∈ X.
However, the vast majority of such maps describe unphysical evolution. We are only interested in the maps which satisfy the equations of motion which we symbolically write as E(ϕ(t)) = 0. In almost all cases, the equations of motion are differential equations which require that the change of the state is determined only by the state of the system at the same instance in time. For classical systems, given a state x ∈ X at initial time t i , there is exactly one solution of the equations of motionφ(t; x) which satisfiesφ(t i ; x) = x.
In a more general setting, we might not precisely know the initial state x ∈ X, but instead our knowledge of the system at initial time t i may be encoded in form of a probability distribution P(x) on the space X. In case P(y) = δ D (y − x) is given by a Dirac δ-function, we recover the deterministic setting above. However, 1 All bold faced symbols take values in X, i.e., they are generally some high-dimensional vectors or tensors. 4 if the probability distribution is more complicated, the state of the system is stochastic for all times t. This means that for any time t there is a unique probability distribution P(x; t) describing the system. In fact, it is simple to give an explicit expression for this probability distribution by using the solution to the equations of motion, where on the right-hand side we have the initial probability distribution. Obviously it is P(x; t i ) = P(x).
The probability distribution P(x; t) allows us to extract expressions for observables at any point in time. As an observable we regard in this context any function O : X → R, i.e., a map which assigns to a state of a system some number (this can easily be generalized to tensor-valued observables). For such an observable, we define its expectation value at time t as In the deterministic case P(y) = δ D (y − x), we simply have Here we defined the notationÕ(t; x) for the observable evaluated along the solution of the equations of motion for initial value x at time t.
From the discussion above one may arrive at a seemingly simple recipe to calculate observables: determine the solutionφ and use it to express any quantity of interest. However, as is well known, even for seemingly simple classical systems one can prove that there is no elementary expression forφ(t; x). In many cases one can determineφ numerically instead, but this can become computationally costly for systems with a highdimensional space of states X. A prime example for this is an interacting N-particle system. The formalism of KFT provides an alternative approach for calculating observables by analytical methods which we want to introduce below.

Path Integral Formulation of Classical Mechanics
The formalism of KFT is based on the path integral formulation of classical mechanics which to the best of our knowledge was first proposed by E. Gozzi [1]. Readers familiar with the path integral approach to Quantum Field Theory (QFT) will notice many similarities between KFT and QFT. In fact, these similarities are one of the main theoretical motivations to use KFT, as one can use the powerful tools developed in the framework of QFT to describe the evolution of a classical physical system.
Let us consider a classical physical system which at initial time t i is in a state x ∈ X. Denote the solution of the equations of motion E(ϕ) = 0 subject to these initial conditions byφ(t; x). Further, let y ∈ X be some possible state of the system at a later time t > t i . Then the transition probability T f i (y, t; x, t i ), i.e., the probability of observing y at t > t i given initial conditions x at t = t i , is given by The interpretation of this quantity is rather simple: While the path integral adds contributions from all possible paths in state space X from x to y, the δ D -functional (a generalization of the Dirac δ-function) picks out a single one of these. This is reminiscent of taking the classical limit in the path integral formulation of quantum mechanics, where the classical path supplies the dominant contribution as → 0.
Clearly, the transition probability T f i (y, t; x, t i ) is one if y =φ(t; x) and zero otherwise. This is not surprising given that the evolution of classical systems with explicitly specified initial conditions is deterministic.
Using the fact that the path selected by the δ D -functional is the solution of the equations of motion, we can replace condition ϕ(t) =φ(t; x) encoded in its argument by the condition E(ϕ) = 0 which avoids explicit reference to the solutionφ. Generally, such a change in the argument of a δ D -functional changes the normalization according to This is a generalization of the identity for the usual Dirac δ-function which is valid if the function f has exactly one root x =x. However, as proven in equation (3.51) of [2], the determinant in equation (5)  While it is possible to keep the determinant in equation (5), it would unnecessarily complicate our treatment.
Therefore, we shall restrict ourselves to the case of Hamiltonian dynamics, i.e., the space of states X actually is phase space and the equations of motion E(ϕ) = 0 are of Hamiltonian form 6 Here, J is the usual symplectic matrix and H is the Hamiltonian. Using equation (5), the transition probability can be expressed as In the second line, we used a standard functional Fourier transform to write the δ D -functional as an exponential function. The object χ is a purely auxiliary abstract mathematical function and does not have an obvious physical interpretation. In analogy to the path integral approach to QFT, we may regard the exponent as the action of the system under consideration.
Further building on the analogy to QFT, we now define the so-called generating functional which is the central mathematical object of KFT. Integrating over the final state y of the system (or, equivalently, setting y φ(x, t)) to get rid of the endpoint of the path integral over ϕ and introducing a source field J for ϕ, it is The generating functional contains complete information about the physical system under consideration.
Any observable depending on the phase space position at a single time t can be derived from it by means of functional derivation with respect to J(t). 2 Indeed, for a given observable O : We used the notationÕ introduced in equation (3) to denote the observable as evaluated along the solution of the equations of motionφ(t; x) for initial value x -which evidently is exactly what we obtain via this differentiation. In fact, this can also be seen by performing the integration over χ directly in the generating functional which yields We remark that we implicitly assume that the observable O is an analytic function of the state x ∈ X, i.e., it is possible to write O as convergent power series in x (to make equation (11) mathematically rigorous).
As a simple example for an observable, take the phase space position O(y) = y for y ∈ X. It has values in phase space X itself (as opposed to being real-valued), but this is no issue at all. Using the prescription (11), we obtainÕ Thus, the observable O precisely extracts the phase space trajectory which is indeed the phase space position as evaluated along the solution of the equations of motion.

Free Generating Functional
If we are able to find an explicit analytical expression for the generating functional, i.e., express Z[J; x] as an elementary functional of J, we would know everything about the physical system under consideration.
Indeed, we could simply obtain the phase space trajectory of the system by functional derivation. Hence, it is expected that in any non-trivial case the path integrals in the generating functional cannot be evaluated analytically without explicitly referring to the solution of the equations of motion. As such, the formalism as developed so far is not all that useful for obtaining explicit results for observables of a classical physical system.
To remedy this, we make the crucial assumption that the equations of motion can be split into a free part and an interaction part, i.e., where we demand that the free part E 0 (ϕ) is given by a linear differential operator acting on ϕ. A simple example would be E 0 (ϕ) = ∂ t ϕ which shows that this decomposition is always possible for Hamiltonian equations of motion. Ultimately, we will perform a perturbative expansion of the interaction part E I (ϕ), such that it is advantageous to include as many terms as possible into the free part of the equations of motion.
Any linear differential operator admits a Green's function (also known as propagator) which can be used to write down a solution to the associated differential equation. In the specific case of E 0 (ϕ) = 0, it is 8 where we denote the free solution byφ and the Green's function by G(t, t ). Note that the free solution usually describing straight paths has been overset by a bar and should not be confused with the phase space trajectoryφ generally describing curved paths and accordingly featuring a tilde.
Using the decomposition of the equation of motion, we can simplify the generating functional significantly.
Suppressing the time dependence for compactness of notation, we find We introduced another auxiliary function K which acts as a source field for χ in order to move the interaction term in front of the path integrals. In the last line we isolated the so-called free generating functional The second term in the last line is to be understood as a matrix product and precisely cancels the additional factor K in the equation E 0 (ϕ) + K = 0. We remark that in performing the integration over ϕ we assumed that this equation has symplectic structure such that the associated functional determinant is unity. This is generically the case for Hamiltonian equations of motion.
Using the explicit form obtained for the free generating functional, we can immediately perform the functional derivatives with respect to K in equation (25). 3 This yields Evidently, there are no path integrals left in these expressions and there is no explicit reference to the solution of the equations of motion. Instead, only the free trajectory appears which can be easily calculated using the Green's function. The free generating functional takes exactly the same form as the generating functional (16) withφ replaced byφ.

Interactions
The expression in equations (29) and (30) for the generating functional strongly suggests a perturbative treatment of the interactions by means of a Taylor expansion. We pursue this strategy momentarily, but as a preparatory step we need to investigate the term E I δ i δJ(t ) in slightly more detail. As usual, the derivative in the argument is to be understood by the virtue of an analytic power series for the function E I . As such, we have the identity This comes as no surprise as this is exactly the purpose for which we introduced the functional derivative in equation (24) -it extracts the trajectory from the free generating functional. This identity can be generalized to the case we will be facing in the perturbative treatment. There, in addition to the generating functional, the derivative in E I can also act on single factors of J. We can use where A is a placeholder for terms of the form G(t , t ) E I δ i δJ(t ) and · · · denotes further terms on which the derivatives may act. The appearance of two terms on the right-hand side is due to the product rule. 3 Note that the term J (t ) G(t , t ) obtained from the differentiation can be moved past the functional derivative in E I δ i δJ(t ) . Indeed, G(t , t ) = 0 if t ≤ t such that the following procedure is possible: We write the exponential as a power series and in each term we sort the products by their time variable t (this is accomplished by splitting the integral into time-ordered regions and considering each of them separately). In each product, we move all the functional derivatives with respect to K to the right and let them act on the free generating functional. The resulting terms then can be moved back to the original positions of the corresponding functional derivatives, allowing the power series to be converted back into an exponential. 10 The second term warrants some further remarks as it may seem counter-intuitive that we obtain a derivative acting on E I . However, this can once again be understood by virtue of power series. Indeed, for a monomial where the factor n in the second term accounts for the choice of which of the derivatives acts on the term ax.
Observe that resulting expression has the form of the derivative of the monomial. The generalization to a power series is straightforward as the same procedure can be applied term-by-term.
It is possible to proceed in full generality. In higher-orders we obtain higher derivatives of the interaction term E I and it is possible to devise a diagrammatic representation of terms order by order to keep track of how the derivatives (which mathematically are tensors) are contracted with the remaining terms. We are happy to supply such prescription upon request, but refrain from detailing it here. Instead we momentarily specify our physical system to be an interacting N-body system in an attempt to sacrifice generality for clarity. Nonetheless, below we give the explicit expressions for the generating functional in full generality up to second order in a Taylor expansion of the exponential in equation (29): Note that we denote with Z 1 and Z 2 , respectively, the additional contributions from first and second order.
The full second order generating functional is given by the sum Z 0 + Z 1 + Z 2 . Observables may be obtained via functional derivatives with respect to the source field J. As an example, the contributions to the phase space trajectoryφ(t; x) up to second order is obtained by acting with δ i δJ(t) (evaluated at J = 0) on the expressions above. We obtaiñ Note that the term containing two factors of J at second order vanishes for this particular observable because the leftover factor J is set to zero (for the same reason, the functional derivative cannot act on Z 0 at any positive order). The first order phase space trajectory, given byφ 0 (t; x) +φ 1 (t; x), coincides with the so-called Born approximation which adds the interactions as evaluated along the free trajectory and is a common ad hoc approximation used, e.g., in scattering problems. 4 Here it is derived as the first order perturbative approximation in the KFT framework.

Choice of Free Motion
As a last topic for the general treatment, we want to briefly discuss the splitting of the equations into a free and an interaction part. It should be stressed that this is a choice which is to a certain degree arbitrary. Indeed, one can always take the trivial splitting E 0 = 0 and E I = E regarding the entirety of the equations of motion as interactions. This way, the Green's function is the identity and the free trajectory is constantφ(t; x) = x. However, this is a poor choice because the interaction part of the equations of motion is only treated perturbatively, while the free part is solved exactly. Hence, it is advantageous to include as many terms of the equations of motion as possible into the free part E 0 . The ability to do so is of course limited by the requirement of the existence of a Green's function for the free motion. Nonetheless, minimizing the interaction part E I generally is expected to improve the accuracy and convergence of the perturbative treatment.
This situation is familiar from many other areas of physics. Fundamentally, what we consider as free (or background) evolution and what we regard as interactions (or perturbations) on top of it, is a choice. As an example, given a density field ρ(q) we may regard it as the sum of a background componentρ(q) and a perturbation δρ(q), i.e., ρ(q) =ρ(q) + δρ(q). However, there is a large freedom in performing this split and depending on the situation different choices might be advantageous. In the example, the most natural split is obtained by introducing a condition like 1 V dq ρ(q) =ρ. But if we, e.g., describe the density field of gas in a galaxy, one might want to consider the dark matter halo as part of the background density.
Likewise in our situation, there is a natural choice for the free motion. Newton' first law asserts that the reference motion is uniform along straight lines. His second law can then be used to determine the acceleration relative to this reference motion. Commonly we would claim that this acceleration is due to a force caused by interactions. However, in some situations it may be advantageous to include part of this acceleration into the free motion. For example, on an expanding background, we usually prefer to use comoving coordinates. In doing so, the physical trajectories of the particles are no longer straight lines.
While this makes the free motion a bit more complicated, it has the advantage that the forces between particles are no longer sourced by the density field, but only by the density excess over the mean density (cf. Appendix A).
There is an even more sophisticated choice. It is possible to use the Zel'dovich approximation for the free motion and model the interactions relative to it [16,17]. In this case it is less clear what the interaction potential should be and, in fact, it is not even clear whether a simple two-body interaction can reproduce the forces between the particles. Recent work suggests, however, that this is indeed possible and demonstrates that a Yukawa-like interaction potential describes the force relative to Zel'dovich trajectories well [18].
Here, we limit ourselves to the more conventional case of using uniform motion along straight paths in comoving coordinates as our free motion.
We emphasise that the choice of free motion does not have any effect on observables. The trajectoryφ and all other observables are exactly the same for any split of the equations of motion E = E 0 + E I . The free trajectoryφ and the interaction part E I can be altered significantly, though, which can have significant effects on computational complexity and the validity of approximations. For example, the perturbative treatment of interactions presented in this work may have very different convergence properties and speed depending on the choice of free motion. It might be an interesting task for future works to formalize and study the symmetry (E 0 , E I ) → (E 0 + ξ, E I − ξ) and find optimal choices for ξ.

Interacting N-body Systems
The discussion above was completely general, but also rather abstract. The only assumption we made was that the equations of motion are of Hamiltonian form and split into a free and an interaction part. We needed their symplectic structure to omit a functional determinant which otherwise would have been a bit of a nuisance. As such, we have been studying a Hamiltonian system governed by an equation E 0 (ϕ)+ E I (ϕ) = 0 evolving from initial conditions x ∈ X along a phase space trajectoryφ(t; x).
In this section we specialize our abstract treatment and notation to interacting N-body systems. We stress that it is possible, albeit notationally laborious, to keep the entire perturbative expansion abstract and general. In particular, it is possible to define Feynman rules for the general setting which yield the Feynman rules derived below as special cases. We refrain from listing them here or even in a dedicated appendix as this would require to introduce an excessive amount of non-trivial notation or, alternatively, work with equations featuring a zoo of indices.

Generating Functional
Let us consider a system of N interacting particles in three-dimensional Euclidean space. The state of this physical system can be described by a point in phase space x ∈ X R 6N . Let us take one such point as our initial conditions and definẽ to be the phase space trajectory of the system. Here,φ j (t; x) is a six-dimensional vector describing the position and momentum of the j th particle at time t. As we have to be able to access the position and momentum parts of this and other objects frequently below, let us introduce the following notation: For any point x in single-particle phase space X R 6 we write x (q) for its position and x (p) for its momentum, Similarly, the 3N-dimensional vector containing all the positions is given by Using this notation and particularizing to an interacting N-body system, the equations of motion for the j th particle are given by The potential V is given by featuring the two-particle interaction potential v. We allow for an explicit time-dependence of this potential for later convenience. It is useful to perform a Fourier transform resulting in where we use the symbol v for the Fourier transform of the two-particle interaction potential, too. We abuse notation in this way throughout and distinguish functions and their Fourier transforms only via their arguments. For Newtonian gravitational interactions in Newtonian gauge we have v(k ) ∝ − 1 k 2 and the splitting of the equations of motion is given by A Green's function for the free equations of motion is given by in Newtonian gauge on a static background. We keep the propagators g qp (t, t ) and g pp (t, t ) as well as the interaction potential v(k , t ) general below to allow for gauge changes in the final results. Moreover, this leaves sufficient freedom to transfer our results to an expanding background necessary for the cosmological case (cf. Appendix A).
Using the explicit forms for the free and interaction parts of the equations of motion, we can simplify the expressions in equations (29) and (30) for the generating functional. Indeed, we obtain Note that we have set the momentum part of the source field J (p) to zero as we are exclusively interested in observables featuring positions of particles. We remark for later reference that this avoids any appearance of the function g pp (t, t ) here and below.

Density Correlation Functions
Aside from the trajectories themselves, the density and its r-point correlation functions are the most interesting observables for an interacting N-body system. In fact, the collective information encoded in the density correlation functions tends to be much more useful in practice compared to the microscopic information encoded in the specific positions of the individual particles. In our applications we integrate the observables over a probability distribution of initial conditions making the specific realization of the system irrelevant.
Instead, we are interested in statistical quantities like the density fluctuation power spectrum. In this section we explain how density r-point functions can be obtained from the generating functional.
The particle number density at position q, given the positions of particles ϕ (q) (t) at time t, is Performing a Fourier transformation which introduces a wave vector k conjugate to q, the density takes the We remark that these equations are valid both for the comoving and the Eulerian densities (cf. Appendix A) because the normalization of the Dirac δ-function, resp. the wave vectors, compensate tacitly the factors of the scale factor. Thus, the equations written down here and below are valid both for a static and an expanding background as long as we regard the wave vector as comoving in the expanding case.
One of the advantages of working in Fourier space is that it is very easy to construct the density correlation functions. Indeed, the (Fourier space) density r-point correlation function simply is G ρ · · · ρ r times (k 1 , . . . , k r ; ϕ (q) (t)) = r s=1 ρ(k s ; ϕ(t)) .
Written out, this yields where the sum in front of the exponential runs over all r-tuples of particles indices 1 ≤ j s ≤ N.
The observables defined above are all given in terms of the positions of the ensemble of particles ϕ (q) (t). In order to calculate these observables for the specific system under consideration, we ought to replace these positions by the solution of the equations of motionφ (q) (t; x) for intial values x. However, we do not attempt to find this solution explicitly, but instead use the identity in equation (11), i.e., 16 Explicitly, we replace the positions ϕ (q) (t) by a functional derivative with respect to J (q) (t) in the expressions for the density and the r-point correlation function above. This yields the operator with which we can act on the generating functional Z[J; x] to obtain the r-point correlation function. The density operatorρ(k, t) is obtained for the special choice r = 1.

KFT Perturbation Scheme
In this section we present a perturbation scheme for obtaining the density correlation functions for an interacting N-body system. We stress again that the entire procedure can be done for any Hamiltonian system allowing for a separation of free motion and interactions. For clarity, concreteness and simplicity of notation we focus on interacting N-body systems here which are also our main application.
Above we already hinted at the possibility of performing a Taylor expansion of the exponential containing the interaction part of the equations of motion in the generating functional (29). In doing so, we obtain a polynomial expression containing in each term factors of the source field J and functional derivatives with respect to it. To the very right-hand side of the expression is the free generating functional which also contains J. Evidently, in each term a certain amount of combinatorics is required to keep track of all the ways the derivatives could act on the various factors of J.
These combinatorics can be easily done algorithmically via successive product rules, but this is neither enlightening, nor efficient. Instead we take inspiration from the path integral formulation of QFT which arrives at very similar expressions for observables. There, the method of diagrammatically representing mathematical expressions has been applied with great success. While our expressions are somewhat unusual from a QFT point of view -most of our functional derivatives come inside Fourier phase factors -it is nonetheless possible to obtain a diagrammatic representation in our case, too. Below we derive our Feynman rules and use them to find perturbative expressions for the observables we are interested in.

Taylor Expansion
Let us consider the generating functional for interacting N-body systems given in equation (47). Replacing the exponential by its Taylor expansion we have the expression with the free generating functional as given in equation (48). It is worthwhile to point out that up to this point the entire treatment has been exact. Perturbatively, the n th order contribution, denoted below by Z n , is simply given by the n th term in this sum. The full expression for the generating functional to n th order is the sum of the first n terms in the sum. Evidently, this truncation constitutes an approximation which may have a limited validity. Taylor's theorem comes with estimates of the truncation error, but it is not trivial to apply them to our case of a derivative operator present in the expansion. We intend to investigate this in future works.
Just like one would do in QFT, we intend to derive Feynman rules to diagrammatically represent the perturbative expressions for observables. The class of observables we want to focus on are density correlation functions as introduced in section 3.2. The usefulness of the Fourier transform becomes evident now, as the phase factors coming from the interaction are similar to the ones coming from the observables. This allows for a unified and hence simplified treatment. We remark that similar Feynman rules can be obtained for, e.g., the trajectories of particles yielding a general recipe for the expressions in equations (37) to (39).
Let us write down the expression for the n th order density r-point correlation function. It is obtained by acting with the operatorĜ ρ···ρ given in equation (55) on the n th term of the sum in equation (56). We obtain 18 Despite the seeming complexity, this expression has actually quite a pleasant structure. The key mathematical identity to use is which can be seen as a special case of the discussion around equation (33). Again the dots · · · denote further terms on which the derivative may act. Given that the phase factors are not modified, they ultimately act in this form on the free generating functional Z 0 . Hence, in the final expression we can replace all the functional derivatives by the free trajectoryφ (q) . The only difficulty is to keep track of which of the j is acted upon by which of the phase factors. Note that this action happens exactly once as any remaining factors of J (q) j would cause the entire term to vanish upon setting the source field to zero.

Feynman Rules
We want to find a diagrammatic representation of equation (57). The first step in achieving this is to understand the term in the square bracket which is taken to the n th power. It contains both factors of the source field J (q) as well as derivatives with respect to it. Hence, the ordering of the individual factors is crucial. Back in equation (25) the ordering was irrelevant because the two derivatives commute. But as soon as we performed the functional derivative with respect to K in equation (29) we had to be careful with this issue. In fact, in footnote 3, we explained in detail how it is even possible to perform these derivatives -by time ordering all terms and exploiting the fact that the Green's function is causal. This allowed us to move the factors of J back to the left past some of the functional derivatives.
We can postpone the step in which we move the factors of J back past the functional derivatives for the moment. Then all these factors -here in the form of J (q) j -are towards the very right of the expression just next to the free generating functional from which they were obtained via differentiation. All Fourier phase factors containing the functional derivatives are left of it. Hence, a priori the situation is such that any phase factor could act on any J We choose this as the starting point for formulating our Feynman rules. Let us use the phase factors as the vertices of our diagrams: We call the first of them density vertex and the second interaction vertex. Note that we replaced the functional derivatives with the free trajectory in accordance with our earlier comments. The remaining parts of the expression in equation (57) are encoded in the edges of the diagrams. These are where the uncoloured vertices mean that the corresponding vertex factors are not yet included.
We claim the following: 5 The expression (57)  j is acted upon exactly once -yet another way of seeing that we obtain collections of trees with density vertices as roots. As a final remark, the lower bound for the time integration is due to our initial conditions x being set at t = t i . 6 For concreteness, we list the diagrams and expressions for the density, i.e., the case r = 1, up to first order below. In zeroth order we have zero interaction vertices and therefore simply obtain the free densitỹ As this corresponds to the interactions being completely absent, it is no surprise that we evaluate the (Fourier transform) of the density function along the free trajectory. The first order correction is given by the diagram In addition, let us list the diagrams which appear in second and third order. The corresponding expressions follow from the Feynman rules and become rather lengthy quite quickly.
We explain in detail why these are the correct sets of diagrams momentarily. Note that there is a certain amount of freedom in defining equivalence of diagrams which in turn induces some symmetry factors. The convention we choose here avoids symmetry factors altogether and allows for easy numerical evaluation of the resulting expressions. There is a minimal change in the Feynman rules to achieve this. If one were to use the rules as written down above, symmetry factors 1 2 and 1 for the diagrams ofρ 2 and factors 1 6 , 1 3 , 1 3 , 1 3 , 1 2 and 1 forρ 3 would have to be put into the expressions.

Symmetry Factors and Shot-Noise
As pointed out above, any discussion of symmetry factors needs to start out with defining what we mean by equivalence of diagrams. We shall use here the simple notion that two diagrams are equivalent if and only if they only differ by the labels of the vertices. In particular, we distinguish between topologically equal graphs with different horizontal ordering of the vertices (an example are diagrams two to four in (67)).
Equivalent diagrams only need to be listed once in the sum over all diagrams. Naively, the cardinality of each such equivalence class appears to n!, where n is the number of vertices. This seems to conveniently cancel the factor 1 n! coming from the Taylor expansion. However, one has to be slightly more careful.
As an example, consider the first diagram in equation (66). This graph actually only appears once in the Taylor expansion, because there is only one way of connecting all interaction vertices to the density vertex.
As a consequence, the cardinality of the corresponding equivalence class is less than n! in such cases.
However, this can precisely be compensated for by a minor change in the Feynman rules. Namely, we demand that all vertices are time-ordered, i.e., the time coordinates should satisfy t i ≤ t 1 ≤ t 2 ≤ · · · ≤ t n ≤ t for diagrams of order n. This can easily be accomplished by adjusting the integral boundaries in the rule for the interaction vertex. The combinatorial factor from the equivalence of diagrams is exactly the same as the one from the possible time-orderings.
In summary, we obtain observables through the sum of all different time-ordered diagrams of the respective order. No symmetry factors are necessary. Note that the time ordering for vertices connected via edges is redundantly encoded by the propagator g qp (t , t). Hence, we can omit the Heaviside function in g qp (t , t) in the following without altering the result which can be advantageous when numerically integrating the ex- There is a second issue which suggests another slight change to our Feynman rules. Namely, the discrete nature of the system yields shot noise terms which do not automatically vanish once we integrate over a probability distribution of initial conditions. Instead, they would require to explicitly perform the thermodynamic limit N → ∞. This is discussed in detail in [8, section 5.2]. It is convenient to remove these ultimately irrelevant terms right away by demanding that density correlation functions are taken only for different particles. Hence, instead of equation (53), we should use where the sum runs over all combinations of different indices. Note that the terms removed are precisely the shot noise terms found in [8].

Particle Trajectories and Interpretation
It is possible to obtain a very similar set of Feynman rules for the particle trajectories, i.e., the observable . Indeed, compared to the case of the density, the only thing to change is to replace the 22 density vertex by the result of acting with the operator on a single factor J (q) j . It can easily be verified that the Feynman rules for the interaction vertex and the edges between them remain unchanged. The discussion of symmetry factors and time ordering remains unchanged, too. The density vertex is gone and instead all diagrams are trees ending in an outgoing edge There is a special case in zeroth order, where we only have an outgoing edge. There (and only there) the operator O j can act on the free generating functional Z 0 directly and hence yields the free trajectoryφ (q) We thus obtain the contributions to the trajectory of the j th particle up to third order as Again we only gave the diagrams for third order as the expressions become rather lengthy. The expressions up to second order are more transparent once we perform the integrations over k a . Indeed, the first order trajectory obtained from summing the contributions of zeroth and first order, respectively, is just the Born approximation. This is in complete agreement to our considerations in the general case and this expression is a special case of equation (38). For the second order contribution we can perform the integrations, too, and arrive at is the Hessian matrix of the two particle interaction potential v. Again, we reproduce one of our general expressions, namely equation (39), specialized to an interacting N-body system.
For the interpretation of these expressions and in particular the associated diagrams, consider first the outgoing edge (69). Here, we want to extract the trajectory of the j th particle which is accomplished by propagating the particle j a from the interaction vertex a at time t a to the final time t. The fact that we use δ j j a for this, suggests to take the particle with index j a as the outgoing particle from the interaction vertex. Then, the form of the expression for the interaction vertex (60) becomes transparent: We consider the interaction of particle j a with all other particles l a at possible times t a . In order to calculate this interaction, we consider the force as evaluated along the free trajectories of the particles. Consequently, in first order we obtain the Born approximation.
Matters become more interesting once we have two interaction vertices connected together by an edge (61). Now, the outgoing particle j a from the left-hand interaction vertex a is identified with one of the particles taking part in the right-hand interaction vertex b. More precisely, the factor δ j a j b − δ j a l b yields two terms.
The first corresponds to the case where the outgoing particle j b is acted upon by the other particles twice at times t a and t b . The second term takes into account the possibility that one of the particles l b taking place in the interaction of j b at time t b had been acted upon already at time t a by all other particles l a . In this case right-hand interaction vertex b determines the forces acting on particle j b at time t b by adding the contributions of all particles l b as if all but one of them would be moving on a free trajectory, while this one particle j a is moving on a Born trajectory. 7 While this is quite an instructive interpretation, it is slightly oversimplified. Indeed, all diagrams only give corrections to the free trajectory. Hence, interaction vertices give corrections to free particle trajectoriessuch that, e.g., in the case of the interaction vertex b what we actually calculate is the difference between 7 Continuing the comparison with the iterative Born approximation from footnote 4: In the second order iterative Born approximation we would have that all other particles move on (first order) Born trajectories, not just the particle j a . Again we can see that KFT is adding in these contributions slower -order by order more particles on Born trajectories are considered for the calculation of the force acting on the outgoing particle. Diagrammatically, the second order iterative Born approximation corrects the first order trajectory by the sum over all diagrams of the form the force as given by all particles l b moving freely and the force given by one of the particles moving on a Born trajectory. Only upon adding in the diagrams of lower orders we are in the situation described above.
The interpretation of diagrams involving density vertices is similar. The key difference is that the density vertex collects the density contributions of all incoming particles similarly to an interaction vertex. Depending on the precise form of the diagram, there might be one or multiple such incoming particles and they might have different orders (in the sense that the subdiagrams have a different number of interaction vertices). As such, the n th order density calculated in KFT perturbation theory is not the density of a system of particles moving on n th order KFT trajectories as one would naively expect. Instead, contributions are added more gradually.

Density Correlation Functions and Power Spectra
The sections above contain all necessary ingredients to derive expressions for density r-point correlation functions. However, since these are the key observables in the cosmological applications we are aiming at, we want to discuss them explicitly here. As stated above, the diagrams to consider for the r-point correlation function in n th order KFT perturbation theory are the ones with n internal vertices and r density vertices. In accordance with our discussion in section 4.3 we avoid symmetry factors and shot-noise terms by enforcing a time ordering t a ≤ t a+1 ∀ a and remove tuples with common indices from the sums in the density vertices.
Then, the contributions to the density 2-point correlation functions up to second order are given by

25
Let us write down the mathematical expression encoded by one of the diagrams to provide a further example of the usage of our Feynman rules. It is For our applications we are specifically interested in the density fluctuation power spectrumP δ . Given a probability distribution of initial conditions P(x) which is spatially homogeneous, the power spectrum is whereρ is the mean density. As evident from this equation, the power spectrum fully describes the density two-point correlation function. If the distribution P(x) is also statistically isotropic, the power spec-trumP δ (k 1 , t) only depends on the norm k 1 . We discuss the process of integrating observables over initial conditions in the following section.

Integration over Initial Conditions
The KFT formalism and the associated perturbation theory as presented in the previous sections can be used to study classical Hamiltonian systems with explicit initial conditions x ∈ X. However, in applications we rarely have access to, nor are interested in the specific initial conditions of the system. Rather, we often prefer to deduce statistical properties of a system. KFT allows for this generalization in a natural way.
Indeed, instead of having our generating function depend explicitly on initial phase space coordinates x, we can integrate it over a probability distribution P(x) of them. We define the averaged generating functional It has been common practice to use this object as the starting point in KFT [8,14,18]. If one performs the integration in equation (84), the observables obtained via functional derivatives with respect to J are averaged over initial conditions, i.e., Without doubt it is convenient to have a single object Z[J] describing the averaged system and being able to deduce observables from it. However, it is quite difficult to perform the integration in equation (84) We opt for the latter option here which in our opinion is advantageous for the perturbative treatment. In this section we apply this strategy to increasingly complicated physical systems.

Toy Model: Gaussian Distribution
In describing physical systems with KFT there are two independent aspects of complexity: Firstly, if we have a complicated Hamiltonian with non-trivial interactions, obtaining perturbative expressions for observables may be laborious. Secondly, if we have complicated initial conditions, the final integration in equation (87) may become rather challenging. While the first aspect generally is less of a problem in itself (the Feynman rules only require derivatives), the more complicated the pertrubative expressions become, the more difficult is the integration over initial conditions. Hence, even if we have a relatively simple probability distribution P(x), it might be infeasible to perform this integration. The second aspect is more direct: If P(x) is a difficult distribution, averaging any non-trivial function over it requires some effort.
In order to get used to the procedure and to introduce some calculational techniques, we consider a series of toy models. The physical system we want to investigate is an ensemble of particles of mass m moving in one dimension and interacting via a two-particle interaction potential v(k) without time-dependence. We can regard this as a one-dimensional interacting N-body system on a static background and use the onedimensional analogues of the expressions derived above. For the initial values at time t i = 0 let us use the probability distribution i.e., a Gaussian distribution both in position and momentum. An elementary calculation shows that the expectation value for the zeroth order KFT trajectory for any particle j reproduces the center of mass trajectory as it should. None of higher-order perturbative contributions to the trajectory alter that result. We show this explicitly for the diagram coming from first order KFT perturbation theory, taken from equation (71): In the last line we already performed the integration over all components of x except x j and x l 1 . This is one key simplification of performing the integration of perturbative observables compared to integrating the generating functional directly -the latter necessarily contains all particle indices which prevents this step.
Inserting the expression of the free trajectory forφ (q) and a similar factor with indices l 1 . This clearly yields Fourier transforms of the position and momentum distributions. In the present case, these Fourier transforms are again Gaussians multiplied with some conveniently cancelling phase factors and therefore their product is an even function. Together with the tacit assumption that the interaction potential v is symmetric, the integration over k 1 yields zero.
Above we have encountered a generic feature of the averaging process for perturbative observables. Because the free trajectory and hence the initial values only appear inside Fourier phase factors, the integration over x always yields Fourier transforms of the initial probability distribution P(x). In probability theory, such Fourier transforms of a probability distribution P(x) of a random variable x is known by the name of characteristic function and usually denoted by Φ x . It can be defined in terms of expectation values via where we point out the unusual sign convention for the Fourier transform. We use the letter f for the conjugate variable for x instead of the in this context usually used letter t to avoid confusion with time coordinates.
In terms of characteristic functions, the expression in equation (91) can be compactly written as where we abbreviated f k 1 , g qp (t 1 , 0) k 1 . Given that in our toy model the distributions for position and momentum are uncorrelated, the characteristic functions factorize into characteristic functions for position and momentum, respectively. In our applications this is generally not possible and we need to work with the characteristic functions of x or even x.
As a more quantitative example, we consider the same Gaussian initial distributions, but determine the expectation value of the density. In zeroth order perturbation theory, the expectation value of the density is Let us immediately perform a Fourier transform to obtain the density as a function of position q: Using the definition of the inital probability distribution from equation (88) we obtain i.e., a Gaussian with mean given by the center of mass position φ (q) 0 j (t) and variance σ 2 + g qp (t, 0) ξ 2 which is monotonically increasing (for standard choices for the propagator g qp ). While this is the expected result, the free density could have been obtained directly without using KFT perturbation theory.
For an actual test of the predictions of KFT perturbation theory, let us consider the first order contribution to the density given by 29 The expression encoded in the diagram can again be simplified and integrated mostly analytically. For simplicity and transparency, we assume a static background, for which the propagator takes the simple (46). Then, the resulting evolution of the density for zeroth and first order is shown in figure 1. Clearly, attractive interaction encoded by the first order correction counteracts the diffusion. Going to higher-orders, this behaviour is enhanced and prolonged. 8 initial spatial distribution at t = 0 has mean value µ = 2 and standard deviation σ = 2 (leftmost Gaussian). The momentum distribution has mean value λ = 2 and standard deviation ξ = 1, such that in every unit of time the density distribution moves to the right by two units and broadens. This precisely is the behaviour in zeroth order perturbation theory (blue). The first order correction counteracts the diffusion due to the encoded attractive interaction (red). [19] Note that not all physical properties of observables are preserved at finite interaction orders. In the specific example, there can be regions where the density is negative in the first order approximation. This is to be expected because the first order correction needs to be negative in some regions in order not to change the normalization. By increasing the interaction strength, the negative terms can then dominate over the positive zeroth order contribution. However, higher-order terms correct for this and restore positivity of the density. 8 Since the interaction is evaluated along free trajectories, at no finite order a bound system can be modelled, i.e., diffusion always eventually wins in this example. With respect to our intended application to cosmic structure formation, this means that at any finite order of perturbation theory we are unable to accurately describe objects below a certain size. Specifically, we expect that for any finite order the resulting non-linear density fluctuation power spectrum falls below the observed one above a certain wave number.

Toy Model: Shell-Crossing
In the previous toy model, the distributions for position and momentum were completely independent. Here we study an example where the momenta and positions of the particles are correlated in the sense that the momentum assigned to particles depends on their position. This is a preparation for the upcoming section on correlated initial conditions which generalizes this. However, the toy model calculations are also of interest in their own right because they show that KFT is not affected by the notorious shell-crossing problem one encounters in cosmic structure formation. Common hydrodynamic approaches break down when streams of dark matter particles cross due to the assumption of the existence of a (single-valued) velocity field. KFT considers the particle dynamics in phase space where this conceptual limitation is circumvented.
Let us consider an interacting N-body system in one dimension subject to the probability distribution of initial values x at time t i = 0. Evidently, the particle positions initially form a Gaussian distribution, but their momenta are assigned such that a particle at initial position x j . This linear relationship implies that there is a time t c , where simultaneously the free trajectory of each particlē This is a Gaussian distribution with variance 1 − g qp (t, 0) σ which transitions through a δ D -distribution as t → t c . Note that the density evolution for t > t c is completely unaffected by the singularity at t c in the KFT formalism. As can be seen in figure 2 (again evaluated for static background), the singularity in the density is a result of the projection of the phase space distribution of the system onto the spatial axis and therefore does not constitute a problem when working in phase space. is a result of the spatial projection rather than an inherent singularity of the system. [19]

Correlated Initial Conditions
As our next step, we want to consider a case of initial conditions where all particle positions and momenta are correlated with each other. More precisely, we sample a (statistically homogeneous and isotropic) density field and assign velocities according to an associated velocity field. This is no longer a mere toy model, but rather is the way in which standard cosmological N-body simulations set up their initial conditions [20].
Indeed, standard cosmological perturbation theory is applicable in the early stages of cosmic structure formation in good accuracy. 9 Stopping the evolution prior to shell-crossing we obtain a density and a velocity field which we can then discretize and follow the particle dynamics as structure formation becomes nonlinear. This way, we use the hydrodynamic picture in its range of validity and switch to a description via N-body dynamics as this becomes necessary. We follow Appendix A of [8] here, but adapt it somewhat.
Let us attempt to formalize the ideas of the previous section. We want that the initial positions of our particle ensemble sample a (number) density field ρ(q) at initial time t i . This can be encoded in the equation The meaning of this equation can be read off: Given the density field, the probability of finding a particle j at the position x (q) j should be given by the density at that precise point. The prefactor 1 N accounts for the fact that in a given volume V we deposit a total particle number N. Having assigned a position to a particle, we want to equip it with a momentum. Again, we can write this down in a rather transparent equation Here, v(q) is the velocity field at the initial time t i . It can clearly be seen from this equation, that individually for each particle its momentum is correlated with its position.
So far, we have treated the density and velocity fields completely independently. However, in our case of cosmic structure formation, we can regard them as the result of a hydrodynamic description valid for times t < t i . Then, the equations of standard cosmic perturbation theory yield a relationship between the two. As a first step, we assert that the velocity field is irrotational. Even if the primordial perturbations have a vector component in the usual scalar-vector-tensor decomposition (which for standard inflationary scenarios they have not), it decays during cosmic expansion. Therefore there exists a velocity potential ψ such that the velocity field is its gradient, v(q) = ∇ψ(q).
Assuming that the density and velocity fields indeed arise via a hydrodynamic description, the are related via the equations of standard cosmic perturbation theory. The relevant equation is which is the final result of the discussion in Appendix B. Here, δ(q) is the density contrast,ρ the mean density and κ a constant which depends on the choice of time variable. Using this equation as well as the definition of ψ as the velocity potential, we can replace the density and velocity fields in equations (102) and (103). This way, we promote the velocity potential ψ to the fundamental object controlling the statistics of our initial particle distribution. Indeed, the initial positions and momenta of our particles are subject to the conditional probability distributions Formally, in order to combine these conditional probabilities into a probability distribution P(x), we have to average them over P(ψ). Written as an equation, it is where the integration runs over all possible smooth field configurations ψ.
On the surface, the integration in equation (106) seems rather difficult. However, our task is simplified substantially by the fact that the velocity potential ψ is a Gaussian random field if we choose our initial time t i sufficiently early. A good choice for t i is the release of the Cosmic Microwave Background (CMB) at redshift z ≈ 1090. This is well in the linear regime of cosmic structure formation and CMB measurements give us very precise statistics for the density contrast δ(q). In particular, we know that the density contrast at time t i is very well described by a Gaussian random field [21]. Hence, the velocity potential ψ(q) likewise is a Gaussian random field and therefore its statistics are solely described by its power spectrum P ψ (k) = κ −2 k −4 P δ (k).
In addition to this, the conditional probabilities only depend on the derivatives of ψ at the particle positions x (q) j . Since ψ is a Gaussian random field, whenever we pick a set of positions follow a multivariate Gaussian distribution. More precisely, It is a standard result that the statistics of the derivatives of a Gaussian random field are related to the statistics of the field itself in a simple manner. Indeed, it is Here we abbreviated the velocity field ∇ψ ∇ψ x and contains initial density-density, density-momentum and momentum-momentum correlations. Using the joint probability P ∇ψ, ∇ 2 ψ , we can rewrite equation (106) into Equation (112) is the initial phase space probability distribution we use for the cosmic structure formation.
The initial positions and momenta of the particles are correlated via density and momentum correlations.
Instead of the probability distribution P(x) we work with its characteristic function Φ x ( f ) below. Therefore, let us derive a compact expression for it: In the second equality we inserted equations (112), (109) and (105) and immediately used the Dirac δfunction to cancel the integral over ∇ψ. We also inserted a factor equal to one which allows to perform the Gaussian integral by introducing an auxiliary Fourier conjugate ξ for m ∇ 2 ψ. Note that in practice the form of f can be read off from the exponential factors in our expressions for observables. A simple example is the expectation value of the zeroth order density After exchanging the sum with the integration, each term with index m is precisely given by the character- while all other components are zero. Inserting this into equation (116), we obtain The Gaussian exponential factor remains difficult even in this simple example.

Correlation Hierarchy
We have seen in the previous section that already the zeroth order density expectation value yields a complicated expression (119). If we go to higher-order perturbation theory or consider higher-point density correlation functions, the expressions become even more intimidating. However, structurally they are quite similar. Therefore, in this section, we study in detail how the expression in equation (119) can be simplified. Afterwards, we show how to deal with the general case of r-point density correlation functions in higher-order perturbation theory. Ultimately, we are interested in the density fluctuation power spectrum.
We start our analysis by simplifying the derivatives with respect to ξ equation (119). Using the decomposition of the correlation matrix C in equation (110) into density-density, density-momentum and momentummomentum correlations, we have Here and in the following we refrain from explicitly writing the dependence of C and its components on x (q) .
The products of 1 + i κ m ∂ ∂ξ j applied to the exponential can be expressed as Since we have a quadratic expression in the exponential in equation (121) and ξ is set to zero after taking the derivative, derivatives act either alone We used C δδ jl = C δδ l j to cancel the factor of 1 2 (we implicitly assume that P ψ (−k) = P ψ (k) which is justified in a statistically isotropic setting). In the second term, we wrote the products on the right-hand side in brackets to remind ourselves that they are actually scalar products since C δp jk are 1×3-matrices for all particle indices j, k. By combining the terms from equations (122) and (123) appropriately, one can construct the expressions for higher derivatives with respect to ξ in a systematic manner. As an example, for the third order derivative in equation (121) we obtain In order to treat the momentum-momentum correlations C pp on the same footing as C δδ and C δ p , we can perform a Taylor expansion of the exponential containing them in equation (120). This way we obtain a product of two infinite sums in which we order terms in powers of the initial power spectrum or, equivalently, the number of C •• for • ∈ {δ, p}. This is in contrast to the treatment in [8], where the different kinds of correlations are treated differently. 10 In our case, the correlation hierarchy is build up from the three types of terms At m th order in the correlation expansion, we take a product of m of these terms. Each such product acquires a prefactor 1 r! , where r is the number of momentum-momentum correlation factors, coming from the Taylor expansion of the exponential. Moreover, there is a prefactor 1 s! , where s is the number of indices not paired with a factor f (p) · . Finally, we sum over all appearing indices (of which there are 2m), demanding that the set formed by both indices of factors D δδ jl and first indices of factors D δp jl has no duplicates (the cardinality of this set is s). This is due to the way the derivatives appear in equation (121)  A priori, all indices in the correlation terms constructed above are summed over the total number of particles N. Once again, if we were to integrate the generating functional Z over the initial probability distribution P(x), this would be the best we could do. However, given that we are considering observables, we can reduce the number of terms substantially. Let us consider the expectation value ρ 0 (k, t) which has been our example throughout. There we know that f (p) has exactly one non-vanishing component f where m is the index being summed over to obtain the contributions of each particle to the density. Thus, we immediately see that all the sums in our correlation expansions over indices appearing in D pp jl or as second index in D δp jl collapse into a single term. Generally, for any density r-point correlation functions, these sums 10 In particular, in [8] the authors treat the auto-and cross-correlations between the particles differently. This leads to a damping j which subsequently needs to be treated with special care. We do not see a strong reason for treating particle auto-correlations differently and therefore do not obtain a separate damping term. Instead, the contribution Q D is kept as part of the momentum-momentum correlation matrix C pp and automatically appears at the appropriate orders of our correlation hierarchy. receive contributions only for indices which already appear in the expression prior to integration over the initial conditions.
Perhaps surprisingly, the same statement is true for the remaining sums. Note that generally components of f (q) are non-zero if and only if the corresponding component of f ( p) is non-zero. This is because they originate from Fourier factors of the form (126) Thus, the general form for density r-point correlation functions in equation (116)  once. This, however, implies that these variables can not appear at all. To see this, recall that the correlation factors have the form for certain functions α •• (k ) and • ∈ {δ, p}. If the variable x (q) l does not appear elsewhere in the expectation value, then the integral over it results in which vanishes for usual choices for the initial density fluctuation power spectrum P δ (k).
Summarizing this discussion, we cannot have any indices in the correlation expansion which are not present in the expression for the density r-point correlation function prior to integration over the initial conditions. This is a substantial simplification, because it means that we can perform the integration over almost all components of x (q) in equation (116) immediately. Specifically for the zeroth order density expectation value in equation (119) we find that only the zeroth order correlation gives a contribution:

Correlations vs. Interactions
We have seen above that for the expectation value of a general density r-point correlation function all the integrals over initial phase space coordinates except for the position variables x (q) l which explicitly appear 38 in the expression for the diagrams are trivial. In the following, we consider the remaining integrals in more detail. This will lead us to the surprising observation that initial correlations between particles are much more difficult to treat computationally than the interactions during the evolution of the system.
Once again, we use the density function as an example in order to make the discussion more concrete. This time, we consider the first order correctioñ and focus on the contribution coming from the density-momentum correlation D δp l 1 j . Recall that the free trajectory is given byφ j , such that once we combine the exponential factors, we obtainρ where the only non-vanishing components of f are the position and momentum parts with indices j and l 1 .
Hence, the expectation value of the first order correction to the density can be written compactly as where we used the characteristic function Φ x ( f ). We remark that this procedure is generally possible for all density r-point correlation functions for any perturbation order.
As the next step, we rewrite the characteristic function as in equation (116) and expand in the correlations.
Leveraging the discussion in the previous subsections, we know that we obtain an infinite sum of contributions with certain correlation factors (125). Crucially, the possible particle indices for the correlation factors are severely restricted, as they must be picked from the set {l 1 , j}. Therefore, the possible contributions are proportional to D δδ l 1 j , D  11 Note that while C δp ·· itself is symmetric in the indices, only the second index is contracted with a factor f (p) · , as can be seen in (125). Thus the contribution for exchanged indices is generally different. at least one D pp ·· factor. The restriction that the set formed by the indices of D δδ ·· and the first indices of D δp ·· must not have duplicates, removes two contributions in first and roughly halves the number of contributions in second correlation order.
Let us return to our example, where we specifically wanted to analyze the contributions coming from D δp l 1 j which is given by Note that the position variables only appear in Fourier phase factors and thus the integrals over them can be performed analytically. Using the definition each of the two integrals yields a Dirac δ-function. These encode the linear system of equations Adding these two equations, we immediately see that the condition k = 0 is implied. However, the entire contribution (134) is proportional to k and thus vanishes. In fact, it can be shown that for any perturbation and correlation orders the expectation value of the density function does not receive corrections. 12 Consequently, the expectation value of the density ρ (k, t) = ρ 0 (k, t) = (2π) 3 δ D (k)ρ is the (Fourier transform of the) mean density as it should be for statistically homogeneous and isotropic initial conditions.
The appearance of such a system of linear equations is a generic feature. In fact, its form can be read off quite easily from the diagram and the indices of the correlations. Let us consider the general case of a density r-point correlation function G (r) ρ (k 1 , . . . , k r ) at n th order perturbation theory and at m th order in the expansion in the initial correlations. We temporarily ignore any initial particle auto-correlations, i.e., terms C •• j j with repeated indices j. For those we do not have any exponential factor such that they do not contribute to the system of linear equations and the corresponding integrals over wave vectors can be evaluated independently.
Generally, the number of equations is given by r + n, as this is the number of particle indices appearing in the expression for the diagrams. Note that while any interaction vertex has two indices l a and j a , the latter is always identified with another index by the Kronecker deltas in the edges. The number of variables in the equations is equal to r + n + m. Specifically, these are the arguments k 1 , . . . , k r of the correlation function, the wave vectors integrated over in the interactions k 1 , . . . , k n and those coming from the initial correlations k 1 , . . . , k m . Because our diagrams are acyclic, the system of linear equations always has full rank. Indeed, the (r + n) columns corresponding to (k 1 , . . . , k r , k 1 , . . . , k n ) are linearly independent. 13 Consequently, we have r + n independent Dirac δ-functions which can be used to cancel r + n integrals over wave numbers k j and k l . Of these integrals there are n + m such that naively the remaining number of integrals is m − r. In practice, however, it is possible that after solving the linear system of equations some of the Dirac δ-functions only contain the arguments (k 1 , . . . , k r ). This is the case, e.g., for the zeroth order density expectation value ρ 0 (k, t) as we saw above. In these cases an equal number of additional integrals remain.
Thus, denoting the number of remaining Dirac δ-functions by p, the actual number of remaining integrals is m − r + p. As an explicit example, the density fluctuation power spectrum P δ has r = 2 and p = 1 (corresponding to a factor δ D (k 1 + k 2 ), where k 1 and k 2 are the wave vectors coming from the density vertices), and hence there are m − 1 integrals remaining. In particular, there is no contribution coming from the zeroth order term in the correlation expansion.
Even in the case of a maximum of r remaining Dirac δ-functions, there is only m integrals over wave vectors remaining. In particular, it is always the case that all integrals coming from the interactions can be cancelled. Conversely, in general the integrals over wave vectors coming from the correlations remain and need to be solved numerically. In practice, this makes the expansion in the initial correlations significantly more difficult than the expansion in the interactions. Indeed, generically we expect one additional simple time integral for any interaction order, but one additional difficult multidimensional wave vector integral for any correlation order.

Reproducing Linear Cosmic Structure Growth
We have seen above that going to high orders in the expansion in the correlations is difficult in practice.
Generically, in the resulting expressions there are multidimensional integrals over wave vectors remaining which are numerically challenging to solve. It is beyond the scope of this paper to discuss these computational issues in detail. Instead, we restrict ourselves to the first order in the correlation expansion in the following. We will see momentarily that this way we reproduce exactly the linear growth of perturbations as predicted by standard hydrodynamic cosmic perturbation theory. The non-linear growth of structures is exclusively encoded in the higher-order initial correlations and will be studied in detail in a follow-up paper [15]. Below we study the density two-point correlation function from which we extract the density fluctuation power spectrum -an important observable in the context of cosmic structure formation. In an analogous fashion, one can derive higher correlation functions and obtain, e.g., the density fluctuation biand trispectra.
Before we can interpret the results we obtain for the power spectrum, we need to lay the groundwork by specifying the form of the propagator g qp (t, t ) and the two-particle interaction potential v(k, t) in the context of cosmic structure formation. In particular, this requires us to transfer our construction of KFT on an expanding background. Similarly to numerical simulations of cosmic structure formation, Einstein's theory of General Relativity is only incorporated on the background level, while the particle interaction is given by Newtonian gravity, ignoring baryonic, radiative and relativistic effects at t > t i . This is a well-established approximate treatment on sub-horizon scales and can be justified both theoretically and numerically (cf. e.g. [22]).

Free Linear Power Spectrum
In our quest for the linear power spectrum, we start with the expression for the density two-point correlation function in zeroth order perturbation theory given by The aim of this subsection is to obtain the expectation value of this observable and to extract the density fluctuation power spectrum from it. In order to do so, we integrate the expression in equation (137)  The calculations of these contributions are analogous to the ones in the previous section. Firstly, we rewrite the integral over G ρρ 0 (k 1 , k 2 , t; x) in terms of the characteristic function ϕ x ( f ) and utilize its expansions in the correlations. Then we consider each contribution individually. As an example, for the correlation In the limit of large particle numbers N → ∞, the fraction in the last line approaches unity. Moreover, we can replace the initial power spectrum of the velocity potential P ψ (k) by the initial density fluctuation power spectrum P δ (k) = κ 2 k 4 P ψ (k).
The full list of non-zero contributions is given by We remark that the factor of 1 2 in the expression from contributions proportional to D δδ ·· is the prefactor of the second order term in the expansion in equation (121). Combining these four contributions, the expectation value for the density two-point correlation function is given by This is precisely the form of a general two-point correlation function of statistically homogeneous density field as given in equation (83). In particular, we can extract the density fluctuation power spectrum The subscript 0 in P δ 0 (k 1 , t) indicates that we are working in zeroth order perturbation theory for the interactions between the particles. The term O C 2 symbolizes that we are working up to first order in the expansions in the initial correlations between the particles. We observe that the power spectrum in equation (146) is proportional to the initial power spectrum. This should come as no surprise given that we argued earlier that the first order correlation contributions to the power spectrum do not come with any integrals over k.

Higher-Order Linear Power Spectrum
Our next step is to compute the first order correction to the two-point density correlation function. It is given by for explicit initial conditions x. The expression for the second diagram is the same as for the first diagram, but the wave vectors k 1 and k 2 as well as the indices m 1 and m 2 are swapped. We have already familiarized ourselves with the procedure of integrating such expressions over our probability distribution of initial values P(x). As a first step, we read off the non-zero components of f (p) which are given by for the first diagram. Then, we rewrite the entire expression for G ρρ 1 in terms of the characteristic function and expand to first order in the correlations.
We know already that the only contributions D •• jl are those where indices j and l already appear in the expression of G ρρ 1 . However, there is a very useful further condition. Namely, any index m corresponding to the left-most vertices of the diagram must appear as an index of some correlation factor D •• jl . This is because these left-most vertices always come with a single exponential phase factor containing the particle position x To provide an example for the calculation, let us consider the contribution of D δδ m 2 l 1 from the first diagram. It is given by (cf. footnote 14) Like for the zeroth order case, we take the limit of the particle number N going to infinity and replace the power spectrum for the initial velocity potential P ψ by the initial density fluctuation power spectrum P δ .
Then, the non-zero contributions are given by Combining the contributions listed above, we obtain the first order correction to the expectation value to the density fluctuation power spectrum

45
where the factor 2 in front accounts for the contribution of the second diagram. In order to identify the two contributions we assumed that the two-particle interaction potential v(k 1 , t 1 ) is symmetric in the wavevector, i.e., it is v(−k 1 , t 1 ) = v(k 1 , t 1 ). We observe that if the potential is of Newtonian form, i.e., v(k 1 , we again only obtain a rescaling of the initial density fluctuation power spectrum.
Going to higher-orders, there is a crucial simplification when we only consider the first order contributions in the correlation expansion. Let us consider the diagram which is part of the second order two-point density correlation function G ρρ 2 (k 1 , k 2 , t; x). We argued above that all left-most vertices must have their respective particle indices appear as indices of the correlation contributions D •• ·· . This is because these vertices -here both interaction vertices and the lower density vertex -feature a Fourier phase factor of a particle position (with indices l 1 , l 2 and m 2 , respectively) which does not appear elsewhere in the expression. However, in linear order of the correlation expansion, we have only two slots for indices, hence there always is one of these phase factors left which subsequently sets the corresponding wavevector to zero and leads to a vanishing result for the entire term. The same is even true if we consider the contributions proportional to δ j a j b in the Feynman rule for the propagator given in equation (61), because then a Fourier phase factor involving x (q) l b remains. Hence, there are yet more contributions which may be ignored when considering linear initial particle correlations only.
Since all diagrams for the density two point function are forests consisting of two trees with density vertices as their roots, these trees must never branch out to make sure that there are only two left-most vertices.
Hence, the n th order contribution to the density two-point correlation function is given by combine (after relabeling 1 ↔ 2 in, e.g., the second diagram) into a single expression in which both t 1 and t 2 are integrated from t i to t.
Let us examine the contributions to the power spectrum in zeroth and first interaction order given in equations (146) and (157) once more. In both cases, the combined contributions from D δδ ·· , D δp ·· and D pp ·· can be written as a product of two terms of which one depends on the time variables of the upper and the other on the time variables of the lower tree. 15 This remains valid at higher-orders and allows us to completely decouple the two trees in the diagrams contributing to the power spectrum. Indeed, we arrive at Note that the quantity f r (k, t) is defined iteratively. As an example, it is For the calculation of f r (k, t) on an expanding background, we use the expressions derived in Appendix A.
In this case, it is for r > 0 and k 0. 16 Aside from exposing that f r (k, a) is dimensionless and scale-invariant for the case of Newtonian gravitational interaction, this also shows that the particle mass m and the mean densityρ do not affect the result of this calculation. The numerical evaluation of this expression can be done iteratively, such that the computational complexity is linear in the perturbation order. Depending on the desired numerical accuracy, the calculation takes about two minutes per perturbation order on a laptop. 15 This is somewhat concealed in these cases by the fact that both density vertices are evaluated at the same time t. We can write these contributions from the different kinds of initial correlations as a product because it happens to hold that 2D δδ jl D pp jl = D δp jl D δp l j upon integrating over x (q) j and x (q) l (the integrals identify the wave vectors and cancel the integrals in the correlation factors C •• ·· ). 16 Note that on an expanding background, the factor m in equation (162)   Like the left panel, but the relative difference is plotted against the perturbation order at final scale factor a = 1. Percent-level accuracy (indicated by the grey dashed curve) is achieved by going to at least 12 th perturbation order. [19] In figure 3 we plot the growth of the amplitude of the power spectrum using equation (161) up to 15 th order in perturbation theory as a function of scale factor for a standard ΛCDM-model. The colorful curves correspond to the orders in perturbation theory, starting from zeroth order with the lowest blue curve. The orange curve above adds the first order correction, the green curve additionally the second order correction, and so forth. Visually, the curves converge against the dashed black line which is simply the square of the linear growth factor. Going to high enough orders in KFT perturbation theory for linear initial correlations, we fully recover the linear growth of the power spectrum. This is quantitatively verified by figure 4. Its left panel shows the residuals of the previous plot, i.e., the relative difference between the expected linear growth and the perturbative approximations. The zeroth order curve (leftmost blue) is accurate to percent-level only up to scale factor a ≈ 0.0010 which is almost immediately after the scale factor at which we set our initial conditions a(t i ) ≈ 0.0009. However, as we include more and more orders in perturbation theory, our approximation remains accurate for longer. If we go up to about 12 th order, we retain percent-level accuracy for the entire evolution. The right panel of figure 4 shows the relative mismatch at scale factor a = 1, i.e., at the end of evolution, plotted against the order of perturbation. It can clearly be seen that the perturbation series converges towards the correct value.

Conclusions and Outlook
Based on the path integral approach to classical mechanics, Kinetic Field Theory (KFT) provides a framework to obtain observables of a classical physical system subject to probabilistic initial conditions. In this work we have disentangled the deterministic evolution and the averaging over the initial conditions. In particular, we have written down the generating functional with explicit dependence on the initial conditions in equations (29) and (30). These general expressions -valid for any Hamiltonian system which allows for a split into free evolution and interactions -can be taken as a starting point for obtaining perturbative expressions for observables. Indeed, it is possible to generalize the Feynman rules we derived in section 4.2 to this general setting.
For the sake of concreteness, we presented the perturbative approach to KFT specialized to interacting Nbody systems. Conceptually, the perturbative approach expands in the force between particles and thus, effectively, in the deviation from the free trajectories. We have argued that the perturbative expansion of KFT adds the interactions in a relatively slow fashion. Therefore, in systems for which the free evolution yields a good approximation to an observable, we expect that the perturbation expansion of KFT provides stable convergence to this observable. We point out that our main application, cosmic structure formation, has precisely this property -free evolution already provides a reasonable estimate of the density fluctuation power spectrum.
The main result of our work is the systematic perturbation scheme presented in section 4. Using our Feynman rules we can write down the mathematical expressions for any density r-point correlation functions for arbitrary order in the perturbative expansion. They also have a clear physical interpretation: At interaction vertices the outgoing particle is perturbed by the entire system via the two-particle interaction potential. If there are any ingoing edges, some of the particles participating in this interaction might already have been perturbed by earlier interactions. The expansion parameter is the interaction potential such that the n th order contribution to an observable accounts for all combinations of n-fold perturbation of particles of the system.
Up to this point, the entire scheme was explicitly dependent on a specific initial phase space position for the system x. The true potential of KFT lies in the averaging process of observables over a probability distribution of the initial values P(x). Crucially, this averaging can be done in large part analytically. In particular, even in the limit of the particle number becoming infinite, all integrals over initial particle positions and momenta can be performed analytically for the observables of interest. This is due to the fact that the initial values x only appear inside Fourier phase factors in our Feynman rules which allows us to express expectation values for observables in terms of the characteristic function Φ x ( f ) of the probability distribution P(x). We explained this using a toy model in section 5.1.
The initial conditions for cosmic structure formation feature a Gaussian random field ψ which acts as the velocity potential and determines the density contrast δ via its Laplacian. This gives rise to density-density, density-momentum and momentum-momentum correlations C •• between different particles. The complicated nature of this probability distribution and its characteristic function (given in equation (116)) suggests an expansion in the initial correlations C •• . Indeed, without this expansion the momentum-momentum correlations appear in an exponential and thus the initial particle positions appear as a Fourier phase factor within this exponential.
The initial conditions for cosmic structure formation feature a Gaussian random field ψ which acts as the velocity potential and determines the density contrast δ via its Laplacian. This gives rise to density-density, density-momentum and momentum-momentum correlations C •• between different particles. The complicated nature of this probability distribution and its characteristic function (given in equation (116)) suggests an expansion in the initial correlations C •• . In doing so, the particle positions appear exclusively inside Fourier phase factors and thus can be integrated analytically. In fact, these integrations yield Dirac δfunctions which cancel the wave vector integrals from the expansion in the interactions.

50
We have given a general recipe to deduce the form of any r-point density correlation function in m th order in the expansions in initial correlations and in n th order in the expansion in the interactions. Combining the expressions for explicit initial conditions obtained from the Feynman rules with appropriate combinations of the correlation factors given in equation (125) we obtain the relevant terms. We have automated these process and are able to produce all terms of any order in our two expansions. The remaining difficulties are in, firstly, the rapidly growing number of diagrams and correlation combinations at higher orders and, secondly, the numerical evaluation of the remaining integrals in their mathematical expressions. In section 5.5 we provided a precise count of the remaining integrals in our expressions for expectation values of density r-point functions in cosmology.
From this count we observe that for any interaction order there remains one time-integral, while for any correlation order we have one leftover integral over a wave vector. The precise count of the wave vectors we need to integrate over is lower if we calculate higher-point correlation functions, but higher if there are leftover Dirac δ-functions. Some simple algebra reveals that if we consider linear correlations, then the terms contributing to the density fluctuation power spectrum do not feature any integrals over wave vectors.
In fact, we have seen that the power spectrum can be calculated via an iterative method in this case which makes the numerical evaluation rather efficient.
One key feature of the case of linear initial correlations is that the density fluctuation power spectrum receives an enhancement which is independent of wave number, i.e., we obtain linear growth. By computing the relevant contributions, we show that KFT fully reproduces the expected linear growth of the power spectrum. Setting up initial conditions at redshift z = 1090, expanding these up to linear order in the initial density-and momentum-correlations and solving the Hamiltonian dynamics perturbatively up to 12th order in the particle interactions reproduces the amplitude of the power spectrum to percent accuracy.
In figure 3 we show how the various perturbation orders approximate the expected linear growth P δ (k, a) ∝ D + (a) 2 up to increasingly low redshift. Given that the perturbative expansion is essentially in the number of interactions per particle, this means that surprisingly few gravitational encounters are sufficient to accurately reproduce linear structure growth.
The fact that linear initial correlations fully reproduce linear growth of the power spectrum also implies that the entire non-linear growth of the power spectrum is due to non-linear initial correlations. In particular, this implies that the momentum-autocorrelations appearing at higher correlation orders are exactly cancelled, most likely order-by-order in the expansion. We therefore caution against treating these autocorrelations separately as a damping term. The formalism we have presented can be used to study non-linear correlations, too, and we will present the results in a future paper [15]. We also intend to study the relationship of our perturbative approach to the Resummed-KFT (RKFT) formalism [13].

Appendix A. Particles on an Expanding Background
When incorporating an expanding background into KFT, the first decision to make is whether our coordinates x ∈ X are comoving or absolute. While both choices are viable, there are a couple of advantages of working with comoving coordinates. Firstly, in comoving coordinates free motion is along straight paths and, secondly, we immediately obtain the commonly used comoving wave vectors via Fourier transform.
Therefore, let x be comoving coordinates and ϕ(t) the comoving phase space trajectory, i.e., the physical coordinates of the j th particle are given by r j (t) = a(t) ϕ We remark that by using an uniformly expanding universe described by a scale factor a(t), we implicitly assumed homogeneity and isotropy of our N-particle system when averaged over large enough volumes.
Since we are working in a flat and thus infinite universe, this necessitates taking the limit of the particle number N going to infinity. Preferring to work with a finite dimensional phase space for the moment, we postpone taking this limit. We note, however, that working with finite N can also be viewed as considering a finite volume V of an infinite universe filled with statistically homogeneously and isotropically distributed particles. When we later send N to infinity for constant average densityρ, we implicitly take simultaneously the limit V → ∞.
Labelled by physical coordinates, the particles are subject to the Lagrangian L (r, ∂ t r, t) = m 2 ∂ t r(t) · ∂ t r(t) − m 2 j k j ν E r k (t) − r j (t), t , (A.1) where we introduced the Eulerian two-particle interaction potential ν E . We can rewrite the second term of L as m 2 j k j ν E r k (t) − r j (t), t = m V E (t; r) , where the Eulerian potential V E (t; r) is defined in analogy to equation (42). This expression should be thought of as the potential energy of the N-particle system. If we regard one of the particle positions r r l (t) as a free variable, we can write down the cosmological Poisson equation Here, G is the gravitational constant and Λ the cosmological constant. The (Eulerian) particle number density field ρ E (r, r(t)) is simply given by a sum of Dirac δ-functions, ρ E (r; r(t)) = j δ D r − r j (t) .

(A.4)
Proceeding similarly to [17] and section 7 of [23] we translate the above to comoving coordinates ϕ. Up to a total derivative, the above Lagrangian is equal to where the peculiar Lagrangian potential V L t; ϕ (q) V E t; a ϕ (q) + 1 2 a(t) ∂ 2 t a(t) ϕ (q) (t) · ϕ (q) (t) (A. 6) appears. Using equation (A.3) as well as the cosmic acceleration equation, one can show that the Lagrangian potential V L satisfies In this equation we introduced the comoving number density ρ q; ϕ (q) (t) a(t) 3 ρ E a(t) q; a(t) ϕ (q) (t) which henceforth, as well as in the main document, is used exclusively. Note that the mean comoving number densityρ is time-independent. Note, that the expression g qp (t, t ) = 1 m (t − t ) θ(t − t ) derived in the special case of a static background in equation (46) is recovered upon replacing a → t and m eff (a) → m.
The second object we are interested in is the two-particle interaction potential. The potential V satisfies the Poisson equation We can define the (comoving) two-particle interaction potential ν implicitly via a decomposition of the potential V as V a; ϕ (q) 1 2 j k j ν ϕ  where in the last step we introduced the matter density parameter Ω m (a). Interestingly, however, the Dirac δ-function originating from the mean density does have an effect for k = 0. In that case, we observe that using the identity (2π) 3  instead of the correct result ρ (q, a) =ρ.

Appendix B. Poisson Equation for the Initial Density Contrast
In order to set up cosmological initial conditions, we need a relationship between the density and the velocity potential. This can be achieved by using the continuity equation of hydrodynamics which is given by ∂ρ E (r, t) ∂t + ∇ r · (ρ E (r, t) v E (r, t)) = 0 . (B.1) Here, the subscript E for the number density field ρ E and the velocity field v E indicates that this equation is valid in the Eulerian description of hydrodynamics which relies on the physical coordinates r. We can go over to the Lagrangian description by splitting the velocity and density fields into a background and a perturbative component. Specifically, we write v E (r, t) = H(t) r + v(r, t) and (B.2) ρ E (r, t) =ρ E (t) (1 + δ(r, t)) , where H(t) = a(t) −1 ∂ t a(t) is the Hubble function. We assume that the Lagrangian peculiar velocity field v(r, t) and the density contrast δ(r, t) are small compared to the Hubble flow and unity, respectively.  δ(r, t))) + ∇ r · ((ρ E (t) (1 + δ(r, t))) (H(t) r + v(r, t))) = 0 (B.5) (t) + 3H(t) (1 + δ(r, t)) + dδ dt (r, t) + ∇ r · v(r, t) + O (δ(r, t) v(r, t)) = 0 . (B.6) In the last term in the second line we collected all terms quadratic in the perturbation. In accordance with the usual treatment in linear hydrodynamic perturbation theory -which is valid in the very early universe where we set up our initial conditions -we will neglect these terms in the following. Furthermore, the first term in the second line yields zero by means of the Friedmann equation. Indeed, it is dρ E dt (t) + 3H(t)ρ E (t) = 0 (B.7) for pressureless cold dark matter even in a universe filled by a mixture of non-interacting fluids. Finally, we introduce comoving coordinates by r(t) = a(t) q(t) and obtain the Lagrangian version of the continuity equation dδ dt (q, t) + 1 a(t) ∇ q · v(q, t) = 0 . (B.8) As is common in linear hydrodynamic perturbation theory, we choose a product ansatz for the density contrast, i.e., we assume that δ(q, t) = D + (t) δ(q). The object D + (t) is the so-called linear growth factor which we normalize via the condition D + (t i ) = 1. This implies that δ(q) = δ(q, t i ) is simply the initial density contrast. Moreover, we recall that we made the assumption that the velocity field is irrotational at initial time t i , i.e., we can write v(q, t i ) = ∇ r ψ(q) = a(t i ) −1 ∇ q ψ(q), where ψ is the velocity potential. Lastly, we change our time parameter from cosmic time t to scale factor a and obtain δ(q) = −κ ∇ 2 q ψ(q) with κ = a(t i ) 2 dD + dt (t i ) .

(B.9)
Here we introduced the logarithmic derivative f (a) = d ln D + d ln a (a) and used D + (a i ) = 1.
The derivation of the Poisson equation for the initial density contrast relied on a number of assumptions which we want to collect here. Firstly, we assumed the validity of linear hydrodynamic perturbation theory at (and only at) the initial time t i . This includes the assumption of the existence of a velocity field v.
Moreover, we assumed that this velocity field is irrotational at initial time t i such that there exists a velocity potential ψ(q). For the density contrast we assumed that we can at initial time factorize the spatial and temporal dependencies, making use of the linear growth factor D + (a). Quantitatively, we only need the inputs f (a i ) and H(a i ). We stress that we only need to know these quantities at initial scale factor a i .