Resummed Kinetic Field Theory: general formalism and linear structure growth from Newtonian particle dynamics

In earlier work, we have developed a nonequilibrium statistical field theory description of cosmic structure formation, dubbed Kinetic Field Theory (KFT), which is based on the Hamiltonian phase-space dynamics of classical particles and thus remains valid beyond shell-crossing. Here, we present an exact reformulation of the KFT framework that allows to resum an infinite subset of terms appearing in the original perturbative expansion of KFT. We develop the general formalism of this resummed KFT, including a diagrammatic language for the resummed perturbation theory, and compute the lowest-order results for the power spectra of the dark matter density contrast and momentum density. This allows us to derive analytically how the linear growth of the largest structures emerges from Newtonian particle dynamics alone, which, to our knowledge, is the first time this has been achieved.


Introduction
The statistical properties of cosmic large-scale structures can be used to study the nature of gravity and test our cosmological model. This, however, requires a precise and accurate theoretical understanding of the evolution of self-gravitating systems, which has proven to be an extremely challenging task.
While conventional Eulerian and Lagrangian analytic descriptions of structure formation provide precise predictions on linear and mildly nonlinear scales, it has been notoriously difficult for them to reach further into the nonlinear regime [1]. This is largely due to the fact that they assume matter to move along a unique velocity field. With collisionless dark matter dominating cosmic structures, though, particle trajectories can cross and form multiple streams, at which point the assumption of uniqueness breaks down. Numerical N -body simulations, on the other hand, are able to describe the nonlinear small-scale dynamics, as they follow the phase-space trajectories of individual tracer particles. But on the downside, they offer only little insight into the underlying principles of structure formation and are computationally expensive.
To overcome these problems, we have developed a new analytic description of cosmic structure formation, dubbed Kinetic Field Theory (KFT), which is built on the full Hamiltonian dynamics of classical particles in phase-space and can thus be seen as an analytic analogue to an N -body simulation. At its very core, KFT is based on the Martin-Siggia-Rose formalism for describing classical systems in terms of a path integral approach [2]. This was first applied to Hamiltonian dynamics by Gozzi et al. in [3,4] and further developed into a full field-theoretical description of kinetic theory by Mazenko and Das in [5,6]. Building on these pioneering works, we extended the KFT framework to describe the dynamics of an initially correlated ensemble of particles and applied it to cosmic structure formation in [7][8][9].
An unusual aspect of KFT, when compared to other nonequilibrium statistical field theories, is that the fundamental dynamical fields in KFT are of microscopic nature even though we are actually interested in macroscopic fields like the density contrast or the momentum density. This renders the application of many well-established and powerful field-theoretical tools, e. g. renormalization group techniques [10][11][12], rather difficult. To facilitate the use of such tools in the future, we develop here an exact reformulation of the KFT path integral in terms of macroscopic fields while preserving all information on the microscopic dynamics.
We show how this reformulation also gives rise to a new macroscopic perturbation theory that resums an infinite subset of contributions to the original microscopic perturbative expansion in orders of the interaction potential presented in [7]. To systematise the calculation of macroscopic-field cumulants within this resummed KFT (RKFT) framework, we introduce a Feynman diagram representation similar to the one used in [13].
Using RKFT enables us to treat dark matter particles in terms of their fundamental Newtonian dynamics rather than the modified version of Zel'dovich dynamics we have been using so far [14]. In previous works, we exploited that the inertial motion of particles on Zel'dovich-type trajectories already contains part of the gravitational interaction. Due to this, even a first-order calculation in the microscopic perturbation theory reproduces the nonlinear evolution of the density contrast power spectrum known from numerical simulations over a wide range of scales remarkably well [7]. For Newtonian dynamics this is not the case and any expansion to finite order in the Newtonian gravitational potential does not even reproduce the linear growth of structures correctly. Using RKFT allows to overcome this limitation. In particular, we show how the lowest order in the new macroscopic perturbation theory precisely recovers the linear large-scale growth by partially resumming the gravitational interactions between particles following Newtonian dynamics. To our knowledge, this is the first time this has been achieved.
In section 2, we briefly review the general framework of KFT by setting up the path integral for classical particles and explaining how cumulants of macroscopic fields can be obtained from it. We also summarise how these cumulants can be calculated using the microscopic perturbative expansion in orders of the interaction potential. Afterwards, we show in section 3 how the KFT path integral can be reformulated in terms of macroscopic fields, yielding the new RKFT framework. Some physical intuition for the resummation process at the level of the microscopic particle dynamics is given and the properties of the new macroscopic perturbation theory are discussed in this section. For this purpose, we introduce a diagrammatic language and derive the respective Feynman rules. In section 4, the RKFT framework is used to describe dark matter structure formation in a standard ΛCDM cosmology. We specifically compute the lowest-order results for the power spectra of the density contrast and momentum density, demonstrating how the linear growth of the largest structures emerges from Newtonian particle dynamics. Finally, we summarise our results and conclusions in section 5. Some of the more technical aspects of our derivations can be found in the appendices.
2 General framework of KFT To render this paper as self-contained as possible, we give a brief review of the general framework of KFT here, focussing on the key aspects necessary to understand the subsequent sections. For a more detailed description we refer to [6,7,9].

Path integral for classical particles
Consider a classical N -particle system described by the phase-space trajectories x j (t) = q j (t), p j (t) of its individual particles, j = 1, . . . , N . Here q j (t) and p j (t) are the position and momentum of the j-th particle, respectively. To condense the notation, we will combine these single-particle vectors into N -particle tensors x(t), q(t) and p(t), adopting the conventions introduced in [7] and [9], where e j is the N -dimensional Cartesian base vector with entries e j i = δ ij and t i denotes the initial time.
The initial phase-space configuration x (i) := x(t i ) of such an N -particle system is assumed to be characterised by a probability distribution P i x (i) and its dynamics at times t ≥ t i shall be governed by some equations of motion E[x(t)] = 0. A central idea of KFT is to encapsulate both of these in the partition function Here, the dynamics is incorporated by functionally integrating over all possible phase-space trajectories x(t) starting from an initial configuration x (i) , where a functional Dirac delta distribution ensures that only the classical N -particle trajectory x (cl) x (i) , t , which solves the equations of motion, contributes. The stochastic initial conditions are then taken into account by averaging over the initial probability distribution P i x (i) .
Since the full solution of the equations of motion is generally not known, we re-express the delta distribution in (2.2) in terms of the equations of motion themselves, .
In this work, we will only consider Hamiltonian dynamics for which the functional Jacobian determinant on the right-hand-side is just a constant that can be absorbed into the normalisation of the path integral measure [4,6]. 1 The partition function can then be brought into the more convenient form by expressing the delta distribution in terms of a functional Fourier integral with respect to an auxiliary field χ(t) with components χ j (t) = χ q j (t), χ p j (t) conjugate to x j (t) and defining the short-hand notation dΓ i := dx (i) P i x (i) for the initial phase-space average. The equations of motion are given by Hamilton's equations, where I d denotes the d-dimensional identity matrix and 0 d a d-dimensional matrix with all entries being zero. The Hamiltonian shall take the form H = H 0 + H I , with H 0 describing the free evolution and H I the interactions. In addition, we assume on the one hand that the free equations of motion are linear in x and can hence be solved in terms of a retarded Green's function, On the other hand, we restrict ourselves to systems without external forces and assume that H I is the total potential energy generated by a superposition of single-particle potentials v depending only on the spatial distances between particles and possibly on time, Let us further introduce the combined microscopic field ψ := (x, χ) and define the microscopic action Like the Hamiltonian, this action splits into a free and an interacting part, S ψ = S ψ,0 + S ψ,I , given by

Collective fields
Usually, one is not interested in the microscopic state of the system but rather in its macroscopic properties. For this purpose, it is useful to introduce some collective fields Φ[ψ] whose cumulants, i. e. connected correlation functions, will contain the desired macroscopic information. Following [9], we choose Φ = Φ f , Φ B with Φ f being the Klimontovich phasespace density and Φ B the phase-space response field which encodes how the particle momenta are changed by a given interaction potential. Working in the Fourier space conjugate to phase-space, with s = k, l denoting the Fourier vector conjugate to x = q, p , these two collective fields are given by (2.14) They allow us to express the interacting part of the action defined in (2.11) as where we used the conventional abbreviations (±r) := ± s r , t r as well as and defined the interaction matrix element The delta functions of l 1 and l 2 appearing in (2.17) are a consequence of the interaction potential being independent of the particle momenta. As we will frequently encounter integrals over 1-and 2-point functions similar to those in (2.15), introducing a short-hand notation for these will greatly improve the clarity of the calculations. For general 1-point functions A 1 (1), B 1 (1) and 2-point functions A 2 (1, 2), B 2 (1, 2) we thus define their dot products as Using this notation, (2.15) condenses to Φ f · σ f B · Φ B . Finally, we construct the generating functional of collective-field correlators by introducing a source field H = H f , H B for the combined collective field Φ into the partition function (2.12), Any macroscopic physical observable of the system can be derived from the pure phasespace density cumulants G f ···f . Of particular interest are their momentum moments. In real space, these can be computed by multiplying an n f -point phase-space density cumulant with an appropriate function of the momenta, F p ( p 1 , . . . , p n f ), before integrating over these arguments. In Fourier space, this translates to applying the derivative operator In the simple case F p = 1 this reduces a phase-space density cumulant to a cumulant of the spatial density field Φ ρ , (2.24)

Microscopic perturbation theory
For interacting particles, the collective-field cumulants can generally only be computed perturbatively. This assumes that the non-interacting theory, described by the free generating functional is exactly solvable in the sense that the free collective-field cumulants If this is the case, one possible way to proceed is to pull the interacting part of the action in front of the path integral by replacing its Φ-dependence with functional derivatives with respect to H, acting on the free generating functional Z Φ,0 [H], (2.27) Expanding the exponential then yields a perturbative series that allows to calculate the interacting collective-field cumulants up to a finite order in σ f B , and hence in the interaction potential v, only requiring the computation of a finite number of free cumulants. This approach has been discussed in detail in [7,8] and we will refer to it as the microscopic perturbation theory.
The reformulation of KFT that we will present in section 3 also requires the knowledge of the free collective-field cumulants. However, for that purpose it will prove more natural to work with free cumulants involving the dressed response field rather than the bare response field Φ B . The reason for this is that in the calculation of pure phase-space density cumulants G f ···f every occurrence of the field Φ B will be dressed with an interaction matrix element σ f B anyway. This is because Φ B appears in the generating functional (2.20) only via the interaction term once we set H B to zero. Let us thus define the dressed combined collective fieldΦ := Φ f , Φ F and its source fieldH := H f , H F to construct the generating functional of dressed free collective-field correlators, The respective dressed free cumulants can then either be obtained from the cumulantgenerating functional WΦ ,0 H := ln ZΦ ,0 H or via their relation to the bare cumulants, In the case of initially correlated particles, the explicit computation of the free cumulants, be it dressed or bare, is rather involved. A thorough treatment of this can be found in [9], where a diagrammatic approach inspired by the Mayer cluster expansion [15] was developed to compute them systematically. For the considerations in the present work, though, it is completely sufficient to be aware of the free cumulants' physical interpretation and a few of their properties derived in [9], which we will only summarise here: (1, . . . , n f ) is the connected n f -point correlation of the phase-space density observed at times t 1 , . . . , t n f that emerges from the ensemble average over the initially correlated, freely evolving particles.

The cumulant G
(0) f ···f F ···F (1, . . . , n f , 1 , . . . , n F ), on the other hand, is the n F -th order response of this n f -point correlation to perturbations of the non-interacting system at times t 1 , . . . , t n F . If n F = 1, this is a linear response. Microscopically, these perturbations correspond to particles being deflected from their free trajectories by the forces − ∇v acting between them and other particles. In Fourier space, this leads to the proportionality 2. Causality tells us that the phase-space density Φ f (u) evaluated at some time t u can only respond to perturbations experienced at earlier times t r ≤ t u . A more detailed analysis further shows that the spatial density Φ ρ k u , t u = Φ f (u) lu=0 can in fact only respond to perturbations experienced at strictly earlier times t r < t u . On the level of the free cumulants this manifests itself in the property . . , n f , 1 , . . . , n F ) = 0 if ∃ r ∈ {1 , . . . , n F } such that t r > t u or t r = t u and l u = 0 ∀ u ∈ {1, . . . , n f } . (2.35) 3. In statistically homogeneous systems, a free collective-field cumulant containing n f phase-space density fields is given by a sum of individual free -particle cumulants with 1 ≤ ≤ n f , describing the contribution from correlations between particles. They satisfy whereρ denotes the constant mean number density and the delta function arises due to the translational invariance in statistically homogeneous systems. Note that all -particle cumulants with < n f describe shot-noise contributions, arising due to the discrete nature of the density field. In the thermodynamic limit N → ∞ these become negligible relative to the dominant term proportional toρ n f .
In appendix A we exemplarily list the resulting general expressions for the free 1-and 2-point cumulants in the case of statistically homogeneous and isotropic systems with Gaussian initial correlations.

Resummed KFT
From the perspective of most quantum and statistical field theories, the KFT partition function (2.12) is rather unusual in the sense that the path integral is expressed in terms of the microscopic fields ψ even though we are actually interested in macroscopic fields like the phase-space density. To facilitate the application of already established perturbative as well as non-perturbative field-theoretical techniques, we would thus like to reformulate the KFT partition function as a path integral directly over the macroscopic fields we are interested in. This can be achieved by exploiting the exact solvability of the free theory and the fact that the interacting part of the action (2.29) depends on ψ only implicitly via the collective fields Φ f [ψ] and Φ F [ψ].

Macroscopic action
We begin by introducing a functional delta distribution to replace the explicitly ψ-dependent field Φ f [ψ] with a new formally ψ-independent field f , Note that this new field f effectively still carries all the information of Φ f [ψ] due to the delta distribution. Most importantly, f -and Φ f -correlation functions are precisely the same. But to emphasise the different origin we will deliberately call f the macroscopic and Φ f the collective phase-space density field in the following.
Similar to what we did in (2.4), we can now express the delta distribution as a functional Fourier transform with respect to a new macroscopic auxiliary field β conjugate to f , Pulling all ψ-independent parts to the front, we find that the remaining microscopic part of the path integral precisely takes the form of the free generating functional ZΦ ,0 , where we defined the combined macroscopic fields φ := (f, β) andφ = (β, f ). It isφ which now plays the role of the dressed collective source fieldH introduced in section 2.3. Finally, using WΦ ,0 = ln ZΦ ,0 , we arrive at our desired result, with the macroscopic action We want to emphasise that this reformulation is exact and hence (3.5) still contains the complete information on the microscopic dynamics, even though S φ does not depend on ψ any more. The microscopic information is now encoded in the free generating functional WΦ ,0 φ and thus, by means of a functional Taylor expansion, in the dressed free collective-field cumulants, (3.8)

Macroscopic perturbation theory
The path integral in the form of (3.5) allows us to set up a new perturbative approach to KFT following the standard procedure in quantum and statistical field theory, i. e. in terms of propagators and vertices. For this purpose, we first split up the action (3.6) into a propagator part S ∆ , collecting all terms quadratic in φ, and a vertex part S V , containing the remaining terms, introducing the inverse macroscopic propagator ∆ −1 and the macroscopic ( Furthermore, we define the macroscopic generating functional Z φ by introducing a source field M = M f , M β for the combined macroscopic field φ into the partition function, Then, the vertex part of the action can be pulled in front of the path integral by replacing its φ-dependence with functional derivatives with respect to M ,Ŝ V := S V δ iδM , acting on the remaining path integral, (3.14) In the last step we dropped a constant functional determinant, as it can be absorbed into the normalization of the path integral. Expanding the first exponential in (3.14) in orders of the vertices now gives rise to a new perturbative approach that we will refer to as the macroscopic perturbation theory. Within this approach, the interacting macroscopic-field cumulants are obtained via In particular, this implies that the lowest perturbative order of the 2-point phase-space density cumulant G f f is given by the f f -component of the macroscopic propagator, For a systematic computation of the higher-order contributions involving vertices, we first need to specify an explicit expansion scheme. We will return to this at the end of section 3.3 after investigating the general properties of the macroscopic perturbation theory in more detail first.
To find explicit expressions for ∆ −1 and V β···βf ···f , we insert (3.8) into (3.6) and identify terms with (3.10) and (3.11), respectively, Here, I denotes the identity 2-point function, The propagator ∆ is then obtained by a combined matrix and functional inversion of (3.18), defined via the following matrix integral equation, (3.21) The matrix part of this inversion can be performed immediately and yields where we defined the abbreviations for the remaining functional inverses. In appendix B we describe how these as well as ∆ f f can be computed explicitly for a given physical system. But even without specifying the system it is always possible to express (3.23) formally in terms of a Neumann series by expanding the functional inverse in orders of iG we can see that ∆ R and ∆ A contain terms of arbitrarily high order in σ Bf and hence in the interaction potential v. Inserting them into the f f -component of (3.22) yields This demonstrates that the lowest perturbative order within the macroscopic approach already captures some features which could only be described at infinitely high order within the microscopic perturbation theory. Consequently, we expect the macroscopic perturbation theory to be generally more powerful than the microscopic one.
We can make this statement more precise by recalling from (3.17) that ∆ f f is exactly the result we get for G f f by taking no vertices into account. Due to (3.11) and (3.19) this is equivalent to computing G f f while ignoring all contributions involving free n-point cumulants with n = 2. Thus, the transition from the micro-to the macroscopic formulation leads to a resummation of all contributions appearing in the microscopic perturbative series that only contain free 2-point cumulants. In a statistically homogeneous system this corresponds to a resummation of all contributions which do not lead to any mode-coupling beyond the one already introduced by the free evolution. Due to this property, we will refer to the macroscopic reformulation of the KFT framework as resummed KFT (RKFT).
To understand what this resummation corresponds to at the level of the microscopic particle dynamics, we need to take a closer look at the meaning of the free cumulants appearing in the propagator. From the properties 1 and 2 discussed in section 2.3 we know that iG (0) f F (1, 2) describes the linear response of the phase-space density Φ f (1) at time t 1 to particles being deflected from their free trajectories by some perturbing forces acting on them at an earlier time t 2 ≤ t 1 . The product iG 2), appearing under the integral in (3.25), accordingly describes the linear response of Φ f (1) to a situation in which the forces generated by those already deflected particles perturb the trajectories of so far still freely evolving particles at an intermediate time This process can be repeated for any possible number of iterations n, with the deflections occurring continuously at all possible intermediate times t1, . . . , tn with t 2 ≤ t1 ≤ · · · ≤ tn ≤ t 1 , as illustrated in figure 1. ∆ R and ∆ A capture this by summing over all possible numbers of linear response cumulants iG , following from (2.34). Consequently, we find (3.27) and we will hence call ∆ R and ∆ A the retarded and advanced (macroscopic) propagators, respectively. If we do not want to explicitly distinguish between the two, we will refer to them more generally as the causal (macroscopic) propagators. In f f · ∆ A the forces generated by the freely evolving correlated particles associated with G (0) f f cause the perturbation needed to initiate the iterative deflection process encoded in ∆ R and ∆ A . Accordingly, ∆ f f describes the connected 2-point phase-space density correlation emerging from the ensemble average over initially correlated particles that have undergone this process. We will thus call ∆ f f the statistical (macroscopic) propagatoradopting the nomenclature used in nonequilibrium quantum and statistical field theory.
We want to stress that even though ∆ R and ∆ A only contain linear response cumulants, ∆ f f will generally contain contributions that are nonlinear in the initial phase-space density correlation, as the free-streaming of particles in G (0) f f itself builds up those nonlinearities as long as the initial particle momenta are correlated [9].

Feynman diagrams
Before applying the macroscopic perturbation theory to cosmic structure formation in section 4, let us first examine some of its general properties. For this purpose, it is convenient to introduce diagrammatic representations for the propagators and vertices. We want to do this in such a way that the diagrams indicate the causal structure inherited from the free collective-field cumulants. Combining (2.34) and (3.19) yields V β···βf ···f (1, . . . , n β , 1 , . . . , n f ) = 0 if ∃ r ∈ {1 , . . . , n f } such that t r > t u or t r = t u and l u = 0 ∀ u ∈ {1, . . . , n β } , (3.28) meaning that every f -argument of a vertex must always be evaluated at an earlier time than at least one of its β-arguments. We visualise this by placing incoming arrowheads on the f -legs and outgoing arrowheads on the β-legs, indicating the direction of time. A general vertex is then represented as where the dots indicate a leg to be drawn for each argument of the vertex. If n β = 0, it follows from (2.35) and (3.19) that the vertex vanishes identically, V f ···f = 0.
The propagator (3.22) contains the building blocks G (0) f f , ∆ R and ∆ A . For the first one of these we can directly adopt the representation (3.29) by means of G The causal propagators ∆ R and ∆ A , however, require a new kind of diagram. In accordance with their respective retarded or advanced causal structure we represent them as a line with an arrowhead in the middle, pointing towards the later time, additionally introducing a conventional factor −i. Whenever we connect the diagrams (3.29), (3.30) and (3.31) in such a way that there are consecutive arrowheads on a line pointing into the same direction, we further join these arrowheads into one. This yields the following diagrammatic representations for the different components of the propagator, Every term in the macroscopic perturbative expansion (3.14) can now be represented by combining the diagrammatic expressions (3.29) and (3.32) for the vertices and propagators appearing in that term. Furthermore, it is a well-established fact in quantum and statistical field theory that the terms appearing in the cumulant-generating functional W φ [M ] := ln Z φ [M ] always correspond to connected diagrams. We will thus only consider those in the following.
In such a diagram every f -leg of each vertex will be attached to an f -end of a propagator and every β-leg of each vertex will be attached to a β-end of a propagator. This allows some immediate conclusions. First, one notices that the arrowhead directions of propagator-ends and vertex-legs being connected in this way will always agree. There will thus be a consistent and continuous time-flow throughout the complete diagram. Second, the fact that there are no propagators or vertices with only incoming arrows implies that the time-flow has no sinks. In appendix C.1 we demonstrate how these two properties can be used to derive the following Causality rule A diagram is causally forbidden and vanishes identically if it has only incoming arrows on its outer legs or contains a subdiagram which does so.
If the system of interest is statistically homogeneous, the macroscopic perturbation theory has some additional very convenient properties. First of all, we recall from (2.37) that the free cumulants conserve spatial Fourier modes in this case. This property is translated to the propagator and vertices,

33)
V β···βf ···f (1, . . . , n β , 1 , . . . , n f ) ∝ (2π) 3 δ d k 1 + · · · + k n β + k 1 + · · · + k n f , The physical interpretation of this rule is straightforward. The only information carried by a diagram with only one external leg is the value of the mean phase-space density. But a homogeneous background density can not affect the dynamics of the particles, as these depend on potential gradients only.
The two Feynman rules will enable us to drastically reduce the number of contributions to the macroscopic perturbative expansion that have to be taken into account. The Homogeneity rule further implies the applicability of the conventional loop expansion scheme when considering statistically homogeneous systems, as this rule ensures that there will always only be a finite number of non-vanishing diagrams with a given number of loops. The complete procedure for calculating the L-loop expression for an interacting n f -point phase-space density cumulant of a statistically homogeneous system within the macroscopic perturbation theory is then as follows: 1. Draw n f points and label these with the Fourier space arguments 1 to n f . .

(3.36)
Having proven that for a statistically homogeneous system it is possible to systematically order the macroscopic perturbative contributions by their number of loops, it remains to be shown that this is also a physically motivated expansion scheme. Unfortunately, there is no obvious expansion parameter which would allow us to estimate the magnitude of the L-loop order contribution simply by the L-th power of said parameter. However, we can nevertheless gain some intuition for the relevance of different perturbative terms by taking a closer look at the structure of the macroscopic vertices.
From (3.19) and item 1 of section 2.3 we know that a vertex V β···βf ···f with n β β-legs and n f f -legs describes the free n β -point phase-space density cumulant and its n f -th order response to interactions between the particles. Such a vertex introduces two distinct types of nonlinearities into the perturbative calculation: On the one hand, it acts as an n β -point correlated source of the interaction potential, which is nonlinear in the initial phase-space density correlations if the initial particle momenta are correlated and n β ≥ 2 [9]. On the other hand, the vertex describes a nonlinear response of the phase-space density to the particle interactions if n f ≥ 2. The Homogeneity rule further implies that only vertices with at least three legs contribute, n β + n f ≥ 3, meaning that every single vertex will introduce at least one of those two types of nonlinearities.
Crucially, fixing the number of loops limits the maximal number of vertices appearing in a perturbative contribution as well as the maximal number of β-and f -legs that any of those vertices can have. The loop order thus characterises the overall degree of nonlinearity caused by interactions that are taken into account. Therefore, an expansion scheme in loop orders indeed seems to be a plausible perturbative approach for probing the nonlinear dynamics of an interacting system of classical particles. Further investigations of the validity and convergence of this scheme are still required, though, and will be subject of future work.

Tree-level calculation of cosmic structure formation
Having discussed the general formalism of RKFT, we will now use it specifically to describe the growth of cosmic structures in a standard ΛCDM cosmology. In this work, we consider the whole cosmic matter content to be made up of a single species of identical dark matter particles, and compute the macroscopic propagator to obtain the tree-level result of G f f . Loop corrections as well as the evolution of a joint system of dark and baryonic matter will be investigated in separate papers.

Newtonian dynamics of dark matter particles
First, we need to specify our choice of coordinates. Instead of the cosmic time t we use as our time coordinate, following [16], with D + being the usual linear growth factor. Note that this definition implies η i := η(t i ) = 0. The spatial coordinate is chosen to be comoving with the homogenous expansion of the background space-time, q := r/a, where r denotes the physical coordinate and a the cosmological scale factor. As momentum variable we use p := d q/dη, i. e. the comoving velocity with respect to the coordinate time η.
In appendix D we show that the resulting Newtonian equations of motion for the phase-space trajectory of a dark matter particle read with the dimensionless matter density parameter Ω m defined in (D.8), the growth function f + := d ln D + /d ln a and the rescaled Newtonian gravitational potentialṼ satisfying the Poisson equation Due to our choice of coordinates,ρ = Φ ρ is the mean comoving number density of dark matter particles. We additionally used Ω m /f 2 + ≈ 1, which is a very good approximation during the matter-and Λ-dominated cosmological epochs that we are interested in [1,11].
The free equations of motion, obtained by settingṼ = 0, can be written as (4.5) Following [14], the corresponding single-particle retarded Green's function G defined in (2.7) is then given by and thus has the components

Furthermore, the Poisson equation (4.4) is solved bỹ
with the single-particle gravitational potential reading v(k) = − 3 2 1 ρ k 2 (4.12) in Fourier space. We will use this potential in the definition (2.17) of the interaction matrix element σ f B . Now all that is left to do is to specify the initial conditions. We choose to fix the initial time to some instant early in the matter-dominated epoch when the cosmic density and momentum fields were still well-described by Gaussian random fields. For this case, an exact expression for the initial phase-space probability distribution P i was derived in [7], depending only on the initial density contrast power spectrum P (i) General expressions for the free collective-field cumulants resulting from this choice of initial conditions have been derived in [9], and in appendix A we list these for the 1-and 2-point cumulants used in this work. Inserting (4.7)-(4.10) and (4.12) into these expressions yields G (0) where we neglected contributions due to shot noise by taking the thermodynamic limit, as discussed in item 3 of section 2.3. This is an excellent approximation on all scales relevant for cosmic structure formation, as the number of dark matter particles contained in any volume of interest is huge. We also already imposed the constraints on k 2 and l 2 set by the delta functions, and defined the short-hand notations  for the time-dependencies introduced by the Green's function. Furthermore, σ 2 p denotes the variance of the initial momentum field, 18) and C 2 describes the 2-point correlations of the phase-space density emerging from the free-streaming of particles. Expanding the latter to first order in P As pointed out in [9], the full nonlinear expression for C 2 given in (A.8) has to be evaluated numerically.

Analytic large-scale limit
To determine the macroscopic propagator for cosmic structure formation, we first insert the expression (4.15) for G This immediately fixes the off-diagonal elements of the macroscopic propagator (3.22). In this limit, the damping factor in the expression (4.16) for G (0) f f simplifies as well, and restricting C 2 to its linear part becomes a very good approximation, as the nonlinear contributions are expected to be sub-dominant on large scales. Hence, the large-scale limit of G (ls) as shown in detail in appendix B. From (4.23) we can infer the large-scale solution of the tree-level density contrast power spectrum, , (4.25) using (2.24), (4.13) and (4.1). We see that we are precisely recovering the well-known linear growth of structures on the largest scales. It can not be stressed enough that at no point in the derivation did we have to employ the Zel'dovich approximation or assume Eulerian fluid dynamics. Instead, this result was derived completely analytically from Newtonian particle dynamics alone. To our knowledge, it is the first time this has been achieved, making it one of the main results of our paper. If we compare (4.25) with the large-scale limit of the freely evolved power spectrum, we can see that the growth of the latter is drastically slower and even bounded from above for late times. The reason for this is that, relative to the expanding space-time, freely evolving particles slow down in comoving coordinates since their initial momentum falls behind the cosmic expansion [14]. We conclude that the gravitational particle interactions resummed in the macroscopic propagator, as discussed in section 3.2, precisely compensate the lack of large-scale structure growth caused by this deceleration.

Numerical solution on all scales
Going beyond the large-scale limit requires a numerical evaluation of the macroscopic propagator. We have shown in appendix B how this can be reduced to solving a simple matrix equation and computing one-and two-dimensional integrals over time arguments. Performing these steps is in itself computationally cheap and numerically stable.
Although the numerical evaluation of the full nonlinear expression for the free 2-point cumulant G (0) f f (1, 2) entering this computation is generally more challenging, we are able to compute the full tree-level power spectrum of the density contrast, in a numerically fast and stable manner, exploiting that G (0) f f (1, 2) is evaluated at l 1 = l 2 = 0 in this case. Only when computing the full tree-level power spectrum of the momentum density, which requires to evaluate G (0) f f (1, 2) also for non-vanishing l 1 and l 2 , our numerical implementation is currently not sufficiently stable. Computing P (∆) π will thus be the subject of future work.
For our analysis, we set the initial time to a redshift of z = 1100, corresponding approximately to the time of CMB decoupling, and assume the initial power spectrum P (i) δ to be given by a BBKS spectrum [18] with spectral index n s = 1 normalized such that σ 8 = 0.8 today. The resulting full tree-level spectrum P (∆) δ is shown in figure 2 for the four exemplary redshifts 500, 200, 20 and 0. We compare it to the freely evolved spectrum, (4.34) as well as the power spectrum obtained from linear Eulerian perturbation theory (EPT), where the latter corresponds to using our large-scale limit (4.25) for all wavenumbers. To To reduce the dynamic range, all spectra are divided by the initial spectrum. The RKFT tree-level result follows the linear growth on large scales, drops below it on scales smaller than the particles' mean free-streaming length, and eventually approaches the freely evolved spectrum. The wavenumber k fs associated with the mean free-streaming length scale is defined in (4.36) and marked with a vertical dotted line for each redshift.
focus on the deviations of the full tree-level result from its large-scale limit, all spectra are divided by the initial spectrum P (i) δ and only wavenumbers k ≥ 10 h Mpc −1 are shown. While the amplitudes of the spectra of course decrease with redshift, we observe the same qualitative behaviour for all redshifts: The full tree-level result follows the linear EPT spectrum on small wavenumbers, but it drops and eventually approaches the freely evolved spectrum when going to higher wavenumbers. The overall shape of the tree-level spectrum thus interpolates between linear growth on large scales and free evolution of structures on small scales. We further see that the wavenumber where the tree-level spectrum starts to drop below the linear EPT spectrum increases with redshift, while the wavenumber where it starts to follow the freely evolved spectrum decreases.
To understand the origin of this behaviour, we have to take a closer look at the influence of the particles' initial momentum variance σ 2 p . Since the initial 1-particle momentum distribution function is given by a Maxwell-Boltzmann distribution [9], the initial variance σ 2 p is related to the mean initial particle momentum viap (i) = 8/π σ p . Using (2.6), the mean distance travelled by a free-streaming particle since the initial time is then obtained by multiplyingp (i) with the position-momentum component of the single-particle Green's function, Taking the inverse yields an estimate for the wavenumber associated with this mean freestreaming length, where we inserted (4.8).
For our choice of initial conditions we find σ p ≈ 7.16 · 10 −3 h −1 Mpc, yielding values of k fs ≈ 139, 79, 51 and 45 h Mpc −1 at redshifts z = 500, 200, 20 and 0, respectively. We marked the corresponding values by vertical dotted lines in figure 2, finding that they match the wavenumbers above which the tree-level power spectrum P (∆) δ drops below the linear EPT spectrum very well. This shows that the gravitational interactions resummed by the macroscopic propagator can not fully counteract the dissolution of structures caused by the particles' free-streaming.
The reason for this is that the propagator captures the full nonlinear free-streaming evolution of the phase-space density, but only its linear response to the gravitational interactions. On large scales this reproduces the linear EPT result, as these scales are dominated by contributions linear in the initial power spectrum. On scales smaller than the particles' mean free-streaming length, though, contributions of higher order in the initial spectrum become relevant and suppress the linear growth of structures.
This behaviour is reminiscent of the small-scale suppression of the power spectrum found in Zel'dovich dynamics [1,8] or Renormalized Perturbation Theory (RPT) [13,19]. Unlike those, however, we have seen in figure 2 that the RKFT tree-level result never drops below the structure growth found for freely evolving particles. This can be explained by the fact that the Neumann series expression (3.26) of the statistical propagator ∆ f f (1, 2) explicitly contains the free cumulant G Altogether, the full tree-level result for the density contrast power spectrum is found to capture the linear effects introduced by gravitational interactions in a way that is consistent with the underlying free Newtonian particle dynamics. Of course, nonlinear effects of the gravitational interactions will modify the small-scale behaviour found here, as these scales lie far within the nonlinear regime of cosmic structure formation. In upcoming papers we will investigate these effects by including loop corrections.

Summary and conclusions
Building on our previous work in [7][8][9], we developed an exact reformulation of the KFT partition function as a path integral over the macroscopic phase-space density field and a macroscopic auxiliary field only, while preserving all information on the microscopic particle dynamics. This reformulation gave rise to a new macroscopic perturbative expansion that resums an infinite subset of contributions appearing in the original microscopic perturbative series introduced in [7]. We further introduced a diagrammatic language that allows to systematically compute perturbative contributions to cumulants of the phase-space density within this new resummed KFT (RKFT) framework, following a simple set of Feynman rules.
Using this framework to describe the growth of cosmic structures in a standard ΛCDM cosmology, we calculated the tree-level results for the power spectra of the density contrast and the momentum density of dark matter particles following Newtonian dynamics. In the limit of large scales, these results precisely recovered the well-known linear growth of structures. To our knowledge, this is the first analytic derivation of the emergence of linear large-scale growth from Newtonian particle dynamics alone. At no point in the derivation did we employ the Zel'dovich approximation or assume Eulerian fluid dynamics.
On scales smaller than the particles' mean free-streaming length, the tree-level result for the density contrast power spectrum drops below linear growth, which is reminiscent of the small-scale suppression of the power spectrum found in Zel'dovich dynamics and RPT. But unlike those, our result never drops below the power spectrum resulting from freely streaming particles. This shows that the lowest-order perturbative result within RKFT is able to take into account the linear effects introduced by the gravitational interactions between particles in a way that is consistent with the underlying free Newtonian dynamics.
The small-scale behaviour of the tree-level power spectrum of the momentum density will be investigated in future work once we have implemented the challenging numerical evaluation of the full nonlinear expression for the free 2-point phase-space density cumulant in a sufficiently stable manner. In upcoming papers we will also discuss the computation of loop corrections within the macroscopic perturbation theory and extend the RKFT framework to the treatment of a joint system of dark and baryonic matter.
form a Gaussian random field with zero mean and covariance matrix .
Homogeneity and isotropy imply that the covariance matrix can only depend on the spatial distance | q 12 | := | q 1 − q 2 | between two points. For the non-vanishing 1-and 2-point cumulants used in this work these general expressions read Here, the vectors L q,r and L p,r encode the phase shift of the Fourier transformed phase-space density caused by the free particle motion from time t to t r , L q,r (t) := k r g qq (t r , t) + l r g pq (t r , t) , L q,r := L q,r (t i ) , If t = t i , we omit writing the time-dependence. Furthermore, σ 2 p is the initial 1-point momentum variance, and the function C 2 describes the contribution to the 2-point phase-space density cumulant emerging from the correlations between 2 freely evolving particles, In the case of cosmic structure formation, the initial momentum field is irrotational and related to the density contrast field via the continuity equation, which allows us to express all components of the initial covariance matrix in terms of the initial density contrast power spectrum P (i) δ defined in (4.13), as was derived in [7]. With these, (A.7) and (A.8) become and where we expanded C 2 only to first order in P (i) δ . Evaluating the full nonlinear expression for C 2 requires to perform the Fourier transform in (A.8) numerically.

B Computation of the macroscopic propagator
The causal macroscopic propagators ∆ R and ∆ A are given by the functional inverse (3.23) which is defined as the solution of the integral equation with the identity 2-point function I introduced in (3.20). The physical systems we are most interested in are statistically homogeneous and have a free Hamiltonian that only depends on the particle momenta but not their positions. The latter implies g qq (t, t ) = θ(t − t ) and g pq (t, t ) = 0, as shown in [7]. In this case we can conclude from (A.3) and (A.5) that G (0) f F can be written as where we introduced the reduced cumulantG . This reduces (B.1) to the following integral equation for which can be solved in two steps. First, we solve it in the case l 1 = 0, which can be done independently for all possible values of k 1 . Afterwards we insert the result ∆ R k 1 , 0; t 1 , t 2 of this into the right-hand side of (B.5) and perform the time integral to obtain the full l 1 -dependent solution∆ R k 1 , l 1 ; t 1 , t 2 .

IfG
(0) f F k 1 , 0; t 1 , t 2 depends on t 1 and t 2 only in terms of their difference t 12 := t 1 − t 2 , then according to (3.25) the same should hold for∆ R k 1 , 0; t 1 , t 2 , and the time integral in (B.6) becomes a convolution, with t1 2 := t1 − t 2 . In this case we can turn the integral equation into an algebraic equation by means of a Laplace transform with respect to t 12 , where z denotes the complex frequency conjugate to t 12 . Bringing∆ R to one side of the equation and performing an inverse Laplace transform then yields the solution of (B.6), Generally, however, the time-dependence ofG (0) f F will not be that simple and we have to determine∆ R numerically. This can be done by discretizing the time arguments in (B.6) and solving the resulting matrix equation for each k 1 -value of interest. Note that, due to their retarded causal structure,G (0) f F and∆ R become lower triangular matrices after discretization. Hence, this matrix equation can be solved with minimal computational effort via forward substitution.
Once the causal propagators are known we can insert them into the relation (3.22) for the complete macroscopic propagator, which immediately fixes the off-diagonal components. The computation of the remaining statistical propagator then reduces to performing the following time integrals, ∆ f f (1, 2) = d1 d2 ∆ R (1,1) G (0) f f (−1,2) ∆ A (−2, 2) (B.10) causal structures of the propagators (3.27) and vertices (3.28) imply that the diagram can only be non-vanishing if the time arguments along the loop satisfy t in 1 ≤ t out 1 ≤ t in 2 ≤ t out 2 ≤ · · · ≤ t in 1 . A closer inspection allows to tighten this restriction even further, though. For this, we first insert the proportionality relation (2.33) of the free collective-field cumulants into the expressions (3.19) and (3.24) for the vertices and the causal propagators, to infer Using this result in figure 3 tells us that the l-argument at every outgoing β-leg of a vertex through which the loop passes is set to zero, which means that we actually have to apply the stricter case of the vertices' causal restriction (3.28). As a consequence, the considered type of diagram could in fact only be non-vanishing if the time arguments along the loop satisfied t in 1 < t out 1 ≤ t in 2 < t out 2 ≤ · · · ≤ t in 1 . This, however, is always a contradiction and thus the diagram vanishes identically, concluding the proof of the Causality rule.

C.2 Homogeneity rule
Let the system under consideration be statistically homogeneous, and consider a diagram that contains a so-called tadpole subdiagram, i. e. a subdiagram which is connected to the rest of the diagram solely via a single propagator. We first use that due to the conservation of spatial Fourier modes, (3.33) and (3.34), such a tadpole subdiagram satisfies where the hatched circle represents a subdiagram with no further outer legs than the one attached to the propagator on the left. Note that we do not have to consider tadpole diagrams using the other two possible propagators since they would vanish identically according to the Causality rule. Then, we insert the proportionality relation (2.33) of the free collective-field cumulants into the definition (3.19) of the vertices, to obtain V β···βf ···f (1, . . . , n β , 1 , . . . , n f ) ∝ k r v(k r , t r ) ∀ r ∈ {1 , . . . , n f } . (C.5) If we now attach the tadpole diagram to the f -leg of some vertex, (C.4) and (C.5) together imply that the resulting diagram vanishes identically if k v(k, t) = 0 at k = 0. Fourier transforming this condition to real space, 0 ! = k v(k, t) k=0 = −i d 3 q ∇ q v(q, t) , (C.6) reveals that it is always satisfied, as the integrand on the right hand side switches sign under a sign flip of q. Physically, this corresponds to the fact that the volume integral of a particle's central force field over the whole space vanishes. This concludes the proof of the Homogeneity rule.
Note that for the gravitational potential (4.12) the condition (C.6) apparently does not hold in Fourier space, even though it is always satisfied in real space. The reason for this is that the Fourier transform of an 1/q potential is actually only properly defined as the limit of a large-scale regularisation, v(k) ∝ lim (C.7) Using this proper definition, the condition (C.6) is indeed satisfied for the Fourier transformed gravitational potential. This is related to the so-called "Jeans swindle" [21]. In practice, however, one rarely encounters situations where using (C.7) makes a difference. In particular, this is true for all calculations in section 4, which is why we directly set k c to zero in (4.12). Nevertheless, one should always keep (C.7) in mind.

D Equations of motion for point particles in an expanding space-time
The Lagrangian of a single classical point particle of mass m in an expanding space-time, expressed in comoving coordinates q = r/a and cosmic time t, is given by with the Newtonian gravitational potential V satisfying the Poisson equation as derived in [14,22]. Here, ρ m denotes the comoving mass density of the total cosmic matter content, andρ m is its mean value, which is constant in time.
Following the steps detailed in [14], we transform the Lagrangian to the new time coordinate η(t) := log D + (t)/D + (t i ) and deduce the resulting Hamiltonian equations of motion, and used Hf + =Ḋ + /D + as well as the fact that D + solves the linearised density perturbation evolution equationD see [1]. Here, is the dimensionless matter density parameter. The rescaled potentialṼ satisfies the modified Poisson equation (4.4), where we additionally assumed the whole matter content to be made up of point particles of mass m, such that ρ m = mΦ ρ .