Kinetic Field Theory: Cosmic Structure Formation and Fluctuation-Dissipation Relations

Building upon the recently developed formalism of Kinetic Field Theory (KFT) for cosmic structure formation by Bartelmann et al., we investigate a kinematic relationship between diffusion and gravitational interactions in cosmic structure formation. In the first part of this work we explain how the process of structure formation in KFT can be separated into three processes, particle diffusion, the accumulation of structure due to initial momentum correlations and interactions relative to the inertial motion of particles. We study these processes by examining the time derivative of the non-linear density power spectrum in the Born approximation. We observe that diffusion and accumulation are delicately balanced because of the Gaussian form of the initial conditions, and that the net diffusion, resulting from adding these two counteracting contributions, approaches the contributions from the interactions in amplitude over time. This hints at a kinematic relation between diffusion and interactions in KFT. Indeed, in the second part, we show that the response of the system to arbitrary gradient forces is directly related to the evolution of particle diffusion in the form of kinematic fluctuation-dissipation relations (FDRs). This result is independent of the interaction potential. We show that this relationship roots in a time-reversal symmetry of the underlying generating functional. Furthermore, our studies demonstrate how FDRs originating from purely kinematic arguments can be used in theories far from equilibrium.


Introduction
Recently, Bartelmann et al. [1,2,3] developed a Kinetic Field Theory (KFT) to treat cosmic structure formation based on methods introduced first by Das and Mazenko in [4,5,6,7] and structurally similar to non-equilibrium quantum field theory. This formalism mirrors the approach of N-body simulations following particles in phase space and, thus, avoiding difficulties with shell crossing ubiquitous in conventional approaches to cosmic structure formation.
The canonical, N-particle ensemble considered is initially correlated in phase space and subject to the Hamiltonian equations of motion. The central object of the formalism is a generating functional containing the complete statistical information about the initial conditions and the propagation of the particles. Correlators, e.g. the density power spectrum, can be extracted from the generating functional using functional derivatives. Gravitational interactions between particles are treated perturbatively using a response function in the spirit of Martin-Siggia-Rose theory [8,9] or can be approximated in the Born approximation [3].
In [1] it was demonstrated that already at first order in the interactions the non-linear power spectrum is in good agreement with N-body simulations down to remarkably small scales. In [2] we have shown that our formalism allows to take the full non-linear coupling of free-streaming trajectories due to initial momentum correlations into account and that the free generating functional factorizes into a single numerically tractable integral of standard form. In a separate analysis [3] we show that averaging the interactions in the Born approximation allows for a computation of the non-linear power spectrum which is in remarkable agreement with N-body simulations with relative differences being of order ≈ 15 % up to a wave number of k ≤ 10 h Mpc −1 for a scale factor of a = 1. With the present analysis, we wish to prove that diffusion and gravitational interactions are kinematically related in KFT.
As a first step we consider the time evolution of the density power spectrum approximating gravitational interactions in the Born approximation as in [3]. One major advantage of KFT over N-body simulations is that we can compute analytic expressions for density correlation functions and, in this way, the study of time derivatives enables us to separate three fundamental processes in structure formation, particle diffusion, the accumulation of structure due to the initial momentum correlations and gravitational interactions. Our analysis shows that diffusion and accumulation are delicately balanced, demonstrating the eminent role of Gaussian initial conditions. We observe that the resulting net diffusion seems to be closely related to the contributions from interactions. This suggests that the processes of diffusion and interactions are kinematically related to each other. We discuss the reliability of the Born approximation for our purposes by comparing the time derivative of the non-linear power spectrum in the Born approximations with N-body simulations.
Motivated by this result, we show in the second part that the time evolution of particle diffusion is related to the response of the system to an arbitrary gradient force by kinematic fluctuation-dissipation relations (FDRs). This relationship is a consequence of a timereversal symmetry of the generating functional which respects the Gaussian form of the initial conditions. Although these FDRs are purely kinematic statements, our analysis shows that they can give insight into processes far from equilibrium. To our knowledge, this is a novel application of FDRs.
This article begins with an introduction into KFT in Section 2, which summarizes the main results from the previous works [1,2,3] and introduces our notation and methods. We then study in Section 3 the time derivative of the non-linear density power spectrum and discuss FDRs in Section 4. Finally we conclude in Section 5 and give an outlook into future applications and relevance of FDRs for cosmic structure formation.

Kinetic Field Theory for Cosmic Structure Formation
In this section we review the formalism of the kinetic field theory recently developed in [1] and continued in [2] and [3]. This serves as an introduction to our notation and methods.

Initially correlated Hamiltonian system
We study a Hamiltonian system of N classical particles with identical mass, which we set to unity for simplicity. The individual particles are described by their phase-space coordinates x j = ( q j , p j ) ⊤ , where the index j = 1, . . . , N denotes the particle number. Introducing the N-dimensional unit vector e j in j-direction, we collect the N phase-space coordinates x j into a phase-space tensor: where summation over j is implied. In the following, bold-faced symbols always denote tensors combining contributions from all N particles. We define a scalar product between two such tenors by: The unit-mass particles are subject to the Hamiltonian equations of motion, which we sometimes write schematically as E(x) = 0. Using the linear growth factor D + − D (i) + as time coordinate, the Hamiltonian of our system in expanding space is given by (see [10]): where g is normalised to 1 at the initial time, a is the cosmological scale factor, H is the Hubble function and v is the Newtonian potential satisfying the Poisson equation with density contrast δ. We assume the particles to be initially correlated in phase space. Every realization x (i) of the initial conditions has a probability described by a phase-space distribution P(x (i) ). Under the standard cosmological assumptions that the initial velocity is the gradient of a Gaussian random velocity field p = ∇ψ and that the initial particle distribution obeys the continuity equation δ = − ∇ 2 ψ, the initial phase-space distribution is completely determined by the initial power spectrum and has the Gaussian form (see [1] for a careful derivation): The factor C(p (i) ), which additionally appears in [1] and describes initial density correlations, is here assumed to be unity. This is a reasonable assumption if we are interested in the late evolution of cosmic structures where the initial density correlations become subdominant compared to initial momentum correlations. The initial momentum correlations are given by the initial density power-spectrum P δ (k): which defines the ( j, k)-component of the matrix C pp appearing in (5). The C −1 p j p k in (5) are then defined as the ( j, k)-component of the inverse matrix C −1 pp . Furthermore, we define σ 2 1 3 ½ 3 ≔ C p j p j as the initial momentum variance.
We can already see from equations (5) and (6) that our theory will contain two competing processes. On the one side there is particle diffusion due to the initial momentum variance σ 2 1 /3. This diffusion should not be confused with thermal diffusion, but should rather be seen as an ensemble effect: The momentum of each particle seen on its own has the variance σ 2 1 when averaging over all realisations of the initial conditions. Every single realisation of the initial conditions has a completely deterministic velocity field without any local variance, i.e. there is no 'thermal' component. Furthermore, the conditional probability C p j p k takes into account that the momenta of any two particles are not independent of each other, but depend on their initial distance | q (i) j − q (i) k | leading to the accumulation of structure already in the free theory.

Generating functional
It was shown in [1] that the entire statistical information on the system is encoded in a generating functional in the spirit of Martin-Siggia-Rose (MSR) theory [8,9]: where we defined the MSR-action This generating functional averages over all possible initial configurations according to the phase-space measure dΓ = dx (i) P(x (i) ), which represents an ensemble average. We introduced an MSR-field χ as the Fourier conjugate to the equations of motion (integration over χ gives a functional delta distribution ensuring that the equations of motion with an auxiliary source field K are satisfied). The auxiliary generator fields K and J are introduced to allow computing correlation functions through functional derivatives of the generating functional: This gives a physical meaning to the field χ as a measure for the response of the system to an external force K. We define the free generating functional Z 0 [J, K] by replacing the full equations of motion by the free equations of motion in the MSR-action (8).
If we denote the solution to the equations of motion with auxiliary source (E + K = 0) byx, or equivalentlyx 0 in the free case, the generating functionals take the form: Usually one is interested in collective observables like the particle density rather than in all the microscopic degrees of freedom described by x. The statistical information about the density in Fourier space can be extracted from the generating functional with the density operator:Φ whereΦ ρ j is a one-particle operator. In the following we abbreviate the arguments of Fourierspace operators by 1 ≔ t 1 , k 1 and −1 ≔ t 1 , − k 1 . We also introduce a collective response field B combining the microscopic response field χ with the gradient of a density field in Fourier space (i k 1 ρ(1)), thus describing the reaction of the system to a gradient force: The one-particle operator for this field is then given by: These collective-field operators enable us to write the full generating functional in terms of the free generating functional: where v is the Fourier-transform of the interaction potential.

Density correlators and interactions
The aim of KFT is to compute cosmological density correlators by applying the density operators (11) to the generating functional: The application of n single-particle density operators to Z[J, K] leads to the translation J → J + L, and thus: where the shift tensor L is given by: Density correlators are thus given by Z[L, 0]. The full generating functional is, however, not tractable in an analytic fashion. The result (14) allows for two different ways of computing density correlators in KFT. We can either apply the density operators to the free generating functional and include the interactions by expanding the exponential in (14), or we approximate the full solutionsx appropriately and work with the full generating functional. We briefly discuss both methods in the following.
If we want to work with the free generating functional we have to define at first what we mean by 'free'. If we decompose the Hamiltonian (3) into a free and an interacting part, H = H 0 + H I : then the solution to the free equation of motion with auxiliary source K is given by: where we neglected the K q j since they are not acted upon by the response operator (13). We also defined the propagator components ‡: Since t = D + − D (i) + this choice of 'free motion' is equivalent to the Zeldovich approximation and thus contains already part of the interactions. The 'remainder' of the interactions is introduced perturbatively via the interaction operator (14). An application of a one-particle response operatorb j r on the generating functional Z 0 [L, K] gives the response factor b j r : At order m in the interaction operator (14) we have to apply m response operators and m additional density operators. So the most general object we have to consider in this scope is the correlator of m response fields and l = m + n density fields computed from Z 0 [J, K]: It remains to compute the density correlator Z 0 [L, 0] from (10): where we introduced the spatial and momentum shift tensors: ‡ In this work, the letter G denotes correlation functions, see (15), as well as Green's functions, however their different index structure should prevent confusion. and we integrated over the initial momenta. The momentum correlations C p j p k depend only on the relative particle separations r jk = q (i) j − q (i) k . This allows us to write Z 0 [L, 0] in the form: where we introduced the one-and two-particle factors P 1 (L p j ) and P 2 (L) as As already discussed in Section 2.1 the free theory contains two competing processes, diffusion due to the initial momentum variance σ 2 1 /3 of every particle seen on its own and accumulation of structure due to the conditional probability C p j p k between the momenta of different particles. This becomes more explicit in (26), where the one-particle factors P 1 L p j describe the diffusion of particle j and the two-particle factor P 2 (L) describes the accumulation of structure due to the conditional probability C p j p k between two particles.
We now discuss an alternative approach presented in detail in [3] approximating the solutions to the full equations of motion in the Born approximation. The full equations of motion following from the Hamiltonian (3) are given by: We again want to write the solution of this equation in the following form: In [10] a propagator of the form: was proven to be particularly useful for the study of cosmic structure formation as it contains an even larger part of the interactions than the Zel'dovich propagator. Substituting the solution (29) with the improved Zel'dovich propagator (30) into the equation of motion (28) shows that the force kernel f j (t) in (29) has the form: Substituting the solution (29) into the full generating functional (10) gives the density correlations: It was shown in [3] that the function F(t, t ′ ) ≔ j F j (t, t ′ ) in the case of the power spectrum, i.e. n = 2, can be averaged and approximated in the spirit of the Born approximation: where with the initial density power spectrum P (i) δ . In [3] it is shown that an approximation of the full power spectrum in this approach is in remarkable agreement with the power spectrum from N-body simulations, with relative deviations being of order ≈ 15 % up to a wavenumber of k ≤ 10 h Mpc −1 .

Factorization
We have not explained yet how Z 0 [L, 0] can be treated. We have shown in [2] that, by introducing the internal wave vectors k ab for a > b and a = 3, . . . , l, the two-particle term can be factorized into a single, numerically tractable integral of standard form: where we defined the wave vectors k j1 for j = 2, . . . , l as: Furthermore, we introduced the Dirac-delta term: and the function The exponent inside the parentheses is a decomposed version of the quadratic form L ⊤ p j C p j p k L p k appearing in the two-particle term (27): where we defined λ jk and λ ⊥ jk implicitly and used the projectors parallel π jk ≔k jk ⊗k jk and perpendicular π ⊥ jk ≔ I 3 − π jk to the unit vectork jk in the direction of k jk .
The function P jk acquires an intuitive meaning when considering its limit for large scales or early times (g 2 qp (t, 0)k 2 jk ≪ 1). It was proven in [2] that, in this limit, P jk is linear in the initial power spectrum: For the free two-point function, i.e. the free power spectrum, λ 21 = −1 showing that, in this case, P 21 reduces to the linearly evolved power spectrum on large scales or early times. Thus, we can interpret P jk as a generalization of the linearly evolved power spectrum which takes the full non-linear coupling of free trajectories by initial momentum correlations into account. Crucially, the function P jk can be quickly evaluated numerically using e.g. a Levin collocation scheme [2,11,12,13].

Time derivatives of correlation functions
We have seen in the last Section that in the free theory we have two physical processes governing cosmic structure formation, i.e. diffusion due to the initial momentum variance σ 2 1 /3 and accumulation of structure due to the initial conditional probability C p j p k between the momenta of two different particles. Including gravitational interactions into the theory adds a third process to cosmic structure formation, viz. the mutual gravitational attraction between two particles. In this section we treat the interactions in the fashion of the Born approximation as studied in [3] and summarised in the last Section. So the non-linearly evolved density power spectrum is given by: In this Section, the propagator g qp (t, t ′ ) is always assumed to be the improved Zel'dovich propagator (30). We want to understand the effect that all three processes have on cosmic structure formation. Since they appear as factors in (42), this is best done by examining the time derivative of the non-linear power spectrum. The time derivative of the power spectrum in the Born approximation (42) is straightforward to compute: where the operators D (1) t , D (2) t and D I t are defined as time derivative operators acting only on the 1-particle factors, i.e. e Q D , the 2-particle factor, i.e. P 21 , and the interaction factor e − F (t) respectively.
The action of these time derivative operators on the non-linear power spectrum can be computed to be: where the derivative of the interaction factor is given by: Note that the derivative of the boundary of the integral vanishes as F(t, t) = 0 since g qp (t, t) = 0. The derivative of the integrand F(t, t ′ ) does not vanish and becomes: In most of the following discussion we neglect the k-independent part of this term as it simply describes an overall rescaling of the spectrum and we are mostly interested in k-dependent features. For this purpose, we introduce the notationsD (I) t andD (I) t as differential operators acting only on the k-dependent and k-independent part of F(t, t ′ ) respectively.
Similarly, the two-particle term (45) contains the linear evolution of the power spectrum as well as additional contributions due to the fact that we consider the full non-linear coupling of free trajectories by initial momentum correlations, see (41) and the discussion thereafter. In the following we are mostly interested in the latter part. Thus, we define: in this way we have subtracted the linear evolution times the Born factor e − F (t) .

Balance between diffusion and accumulation of structure
In a first step we want to examine the one-and two-particle contributions, (44) and (49), describing the diffusion and accumulation of structure due to the initial conditions. Both contributions are depicted in Fig. 1 as well as their sum and for comparison also the time derivative of the linearly evolved power spectrum g 2 qp (t, 0)P (i) δ . We divided the one-and twoparticle contributions by the Born factor e − F(t) in order to make the curves independent of any shortcomings of the Born approximation on small scales. In this way the results are exact at any scale and it makes sense to depict them up to a wave number of k = 100 h Mpc −1 .
We observe a remarkable balance between the one-and two-particle contributions on small scales as their sum is several orders of magnitude smaller than both terms individually. This balance originates from the Gaussian form of the cosmological initial conditions and is only possible because our formalism allows to take the full non-linear coupling of free trajectories due to the initial momentum correlations into account -accounting only for part of them would lead to a strong domination of diffusion. Any small violation of this balance would lead to either very strong diffusion, thus preventing any structure formation, or much stronger structure formation than currently observed (many orders of magnitude larger than the linear evolution). An interesting consecutive question for a future study is whether this observation could constrain initial non-Gaussianities.
The sum of the one-and two-particle contributions can be positive or negative, depending on the time we are looking at, see Fig. 2 for its behaviour at different times. While we can  see at early times that this sum leads to some structure formation on small scales, diffusion dominates at late times where it results in a net diffusion effect.

Interactions and net diffusion
We now want to include the contributions from interactions (46) into the picture. The time derivative of the power spectrum in the Born approximation, see (43), is given by the net diffusion, i.e. the sum of the one-and two-particle contributions, now including the Born factor, the linear evolution times the Born factor, i.e. the part of the two-particle contributions which we neglected so far, see (49), and the contributions from interactionsD (I) t P nl andD (I) t P nl . We depict these terms individually as well as their sum in Fig. 3. Also shown is the linear evolution for comparison. We observe that the time derivative of the k-independent part of the Born approximation is negligible on all scales. The time derivative of the linear power spectrum times the Born factor ensures that on large scales the linear evolution is reproduced and is also important for the non-linear evolution on small scales. The time derivative of the k-dependent part of the Born factor (D (I) t P nl ) follows the net diffusion closely on large scales, but has the opposite sign, on intermediate scalesD (I) t P nl is bigger than the net diffusion and on small scales these two contributions again come quite close to each other in amplitude.
Before we can draw any further conclusions, it is important to check the reliability of the Born approximation by comparing our results with the time derivative of the power spectrum   obtained from N-body simulations [14]. We depict the result from the Born approximation (43) in comparison to N-body simulations in Fig. 4 for three choices of the scale factor a = 1, 0.3, 0.1. In Fig. 5 we plot the relative difference between our results and N-body simulations for the same choices of the scale factor as well as the scales factors a = 0.9, 0.8 which we discuss in more detail later on. At first we discuss the results for the scale factor a = 1 in Figs. 4 and 5. When comparing the results in the Born approximation with simulations, we observe a qualitatively similar behaviour as in [3] where the non-linear power spectrum (42) was compared with simulations. The Born approximation is able to describe the non-linear structures up to a scale of k ∼ 5 h Mpc −1 to remarkable accuracy with the absolute value of the relative error being on average of the order of ∼ 15 % and never bigger than 25 %. On scales beyond k ∼ 5 h Mpc −1 the Born approximation falls strongly below the results from simulations.
Considering the evolution of the relative error in time we see that, as expected, the Born approximation predicts the early evolution (a = 0.1) reasonably well. However, the relative error is of order ∼ 30 % for scale factors between a = 0.3 and a = 0.9 and on large scales the error does not increase monotonically with redshift. The curves for a = 0.3 in Fig. 4 show that the Born approximation fails to predict the correct form for the onset of the non-linear structure at 1 h Mpc −1 k 5 h Mpc −1 . Notice that for k 5 h Mpc −1 and a = 0.3 the magnitude of the relative error decreases again. This contradicts the intuition that a perturbative theory of structure formation should be accurate on linear and mildly non- linear scales, but should fail on highly non-linear scales. However, the Born approximation is a non-perturbative, but approximate approach to cosmic structure formation and it is not clear whether this intuition is reasonable in this case. The Born approximation seems to be especially well suited for the study of the late time non-linear evolution. The time evolution of the relative error needs to be studied in detail together with the reliability and limitations of the Born approximation in a future analysis. For our purposes here it is merely important that the Born approximation at late times (a = 1, 0.9, 0.8) is in agreement with simulations at a level of 30 % on scales up to k ∼ 5 h Mpc −1 .
With these caveats in mind we want to analyse the processes of net diffusion D (1) t +D (2) t P nl and the k-dependent part of the interactionsD (I) t more closely. We plot the sum of both contributions relative to the net diffusion in Fig. 6. We observe that the two contributions become significantly closer in amplitude for late times. Notice that this tendency of the two terms to approach each other in amplitude is most notably visible on scales k < 5 h Mpc −1 where the Born approximation is in reasonably good agreement with simulations. On large scales the two contributions seem to follow a close relationship.
These observations would need to be considered unnatural unless there is some mechanism relating the amplitude of both processes to each other. This motivates our analysis in the next Section where we indeed find a relation between interactions and diffusion in terms of fluctuation-dissipation relations.

Fluctuation-Dissipation Relations
In this Section we prove a fundamental connection between diffusion and interactions in terms of fluctuation-dissipation relations. We show that FDRs in KFT relate the process of diffusion with the reaction of the free system to an arbitrary gradient force. In this Section we work in the 'free' theory, where the trajectories of free particles are given by (19) and (20) for K = 0, the propagator is given by the Zel'dovich propagator (21) and the reaction to the system to an arbitrary force K is computed via the response operator (13). This makes the relations independent of the interaction potential.

Density autocorrelation
In a first step we examine the single-particle density autocorrelation G (0) ρ j ρ j (12) in more detail. In this case the system is effectively a one-particle system and the initial probability distribution (5) reduces to a Maxwellian velocity distribution: Thus, the autocorrelation is purely described by the dissipative one-particle part P 1 (L p j ) and two-particle factors are absent: The Dirac δ-distribution reflects the homogeneity of the system and leads to the following form of the momentum shift vector L p 1 : which is invariant under time-translation in the case of Zel'dovich trajectories. Having determined L p j , the time derivative of the autocorrelation can be computed: where we used ∂ t 1 g qp (t 2 , t 1 ) = −1 for the Zel'dovich propagator. This expression for the time derivative of the autocorrelation becomes most interesting when comparing it with the response function: We can conclude: This relation between the response function and the time derivative of the autocorrelation has exactly the form of a fluctuation-dissipation relation.
We have multiplied the response field on the left-hand side of (55) with −i because this is the correctly normalized response field measuring the reaction of the system to two-particle interactions as can be seen in (14), where this factor also appears.
Fluctuation-dissipation relations are known to hold for many statistical systems [15]. For example, in classical linear response theory one considers a system in thermal equilibrium and examines the response of an observable O to a small perturbation away from equilibrium. This is measured by the response function χ(t, t ′ ). Using time-translation invariance in thermal equilibrium, the standard textbook result can be derived: This result is remarkable as it tells us that the response of the system to a small perturbation is the same as the evolution of equilibrium correlations. Examination of the equilibrium systems thus gives information on non-equilibrium physics. β is here the inverse temperature, so comparison of (55) and (56) shows that in KFT σ 2 1 /3 assumes the role of the temperature. The FDR in KFT (55) has the same mathematical form as the conventional FDR (56). However, there are some significant physical differences between the relations since KFT is a theory far from equilibrium. We discuss this in Section 4.5 in more detail, but at first we want to generalise the result (55).

General density correlators
Time translation invariance (TTI) is an important requirement for the validity of FDRs. In the case of the autocorrelation, TTI follows from homogeneity allowing us to write L p j in the form (52). If the structure of the correlator is not that simple, time-translation invariance might not hold for all momentum shift vectors L p j . This leads to deviations from the FDR in its conventional form.
The most general density correlator Z 0 [L, 0] has an arbitrary number n j of density fields with particle index j and equivalently for any other particle index: We introduced the last term as an abbreviation which we use in the following. For this correlator, the momentum shift vector L p j has the form: where n j is the number of applied density operators with particle index j. We split L p j into a part respecting TTI and a part which breaks this invariance: where the time variable t 1 is arbitrary so far and the TTI breaking part of L p j is given by: We see that ∆ L p j vanishes for any autocorrelation, where n j is identical with the full number of applied density operators and statistical homogeneity thus ensures n j s=1 k s = 0. In the general case, however, ∆ L p j is non-zero and leads to a correction term in the fluctuation-dissipation relation. Analogous to the steps (53) to (55) we derive the relation: where we assumed t s > t 1 for all 1 < s ≤ n j and t 1 is now the time where we evaluate the response function. We see that breaking TTI leads to a violation of the FDR in its conventional form. The response factor on the left hand side measures the propagation of an inhomogeneity from t 1 to the times t s . However, P 1 L p j describes the diffusion from the initial time to the times t s . Only if the modes satisfy n j s=1 k s = 0, can we neglect the propagation of particle j from the initial time to t 1 . Thus, in general, we have to subtract the diffusion taking place between the initial time and t 1 , which is encoded in ∆ L p j (t 1 ).
In order to quantify the diffusion taking place between t 1 and the times t s , we introduce a covariant time derivative being invariant under time translation: where the operator D (1, j) t is defined as a time derivative acting only on the one-particle factor of particle j, i.e. P 1 L p j . In terms of this time derivative, the FDR for a general density correlator, assuming t s > t 1 for all 1 < s ≤ n j , has the form: We see that the diffusion of particle j and the reaction of particle j to an arbitrary external force are kinematically related and this result is not restricted to the autocorrelation (55), but is valid for an arbitrary correlator if we substitute the full time derivative by the covariant time derivative D (1, j) t describing the evolution of diffusion. In the following we show that this connection between the one-particle ensemble diffusion and the response function goes even deeper and is a consequence of the structure of the generating functional and the Gaussian form of the cosmological initial conditions.

Time-Reversal Symmetry
In statistical field theories, FDRs are typically connected to a time-reversal symmetry of the generating functional [16]. We will show that the same is the case for our kinetic field theory. This gives an alternative way of deriving the FDRs considered above and leads to an easy generalization to higher-order response functions, which will turn out to be related to higherorder time derivatives.
First of all we note that the formalism of KFT is originally designed for times later than the initial time t = 0. Hence it is not immediately clear what a time-reversal symmetry means within this formalism. However, in principle there is no need to restrict the formalism to positive times. The solutions to the equations of motion (19) and (20) are also valid for negative times. Thus, we can propagate the particles from their initial conditions backwards in time and calculate correlation functions at negative times by functional derivatives of the generating functional Z[J, K] with respect to the auxiliary fields at negative times: where t, t ′ > 0. Our aim is to find a transformation T of the microscopic fields x and χ which reverses the time coordinate and leaves the generating functional invariant. The time-reversed generating functional is defined as: and in this way contains the statistical information of the system with time-reversed dynamics.
The phase-factor with the auxiliary fields is not transformed. Thus, correlation functions are computed in the same way as before, but they are now evaluated with time-reversed dynamics. If the generating functional is invariant under the transformation T , any correlation function will be invariant as well. This can be shown by substituting x → T x and χ → T χ in the path integrals in (65) and using that a time-reversal symmetry has to be its own inverse: In this form, we see that functional derivatives of the time-reversed generating functional give time-reversed correlation functions, but with the original dynamics. Thus, invariance of the generating functional also shows invariance of correlation functions under time reversal.
We now aim at finding the explicit form of T . Since the density and response field operators, (11) and (13), act only on the spatial part of J and the momentum part of K, we can restrict ourselves here to the case where the momentum part of J and the spatial part of K vanish. Using time-reversal invariance of the Hamiltonian equations of motion, we prove in Appendix A that the following transformation is a time-reversal symmetry of the free generating functional: T : where the term c j [J, t] is a functional of J and is defined as: It describes the breaking of time-translation invariance similar to the term ∆ L p j in (60). Indeed, if we substitute J → L in (68), as is the case for density correlators, (60) and (68) become completely equivalent. The transformation laws for the microscopic degrees of freedom allow us to derive transformation laws for the macroscopic fields B and ρ. Defining1 ≔ −t 1 , k 1 , the timereversed density and response fields become: with The first term on the right-hand side has a simple interpretation if we consider this term to be embedded into a density correlation. We prove in Appendix B.2 that: where the dots denote an arbitrary number of further density fields. Thus, we can interpret the term ∆B j (1) as the evolution of time-reversed diffusion: In the following two paragraphs we want to build up some intuition for this term. A time-reversal symmetry relates the probability of a path and its time-reversed version. To illustrate this, we consider a single particle of the ensemble and choose one realization of its initial position q, but leave the initial momentum free. The trajectories of the different realizations of the momentum form a double cone as shown in Figure 7 if we propagate the particle also backward in time. At each time the probability to find the particle at some position q is given by a Gaussian with variance σ 2 1 /3g 2 qp (t, 0), which is symmetric in time t q t 1 q p(q, t 1 ) = p(q, −t 1 ) −t 1 Figure 7. Sketch of different realizations (dashed lines) of the trajectory of one particle with fixed initial position propagated also to negative times backwards from the initial conditions. The probability p(q, t 1 ) to find the particle at position q at time t 1 is Gaussian and symmetric under time reversal.
reversal. This symmetry is the reason why density correlations are invariant under time reversal, see (69). However, we see in Figure 7 that for positive times diffusion takes place, while for negative times the process of diffusion is reversed. The response field, being the functional Fourier conjugate of the equations of motion, see (8), is sensitive to the direction of time. Thus, the response field is not invariant under time reversal, but the time-reversed diffusion has to be taken into account in the transformation (70).
The form (70) of the symmetry together with (73) is reassuringly similar to the timereversal symmetries known from statistical field theory. For example Andreanov et al. [16] consider the Langevin dynamics: with potential V(X) and where the stochastic force η(t) has variance 2T . They find the timereversal symmetry: whereX takes on the role of the MSR-response function equivalent to χ. Comparing with (70) and (73), we see that the temperature plays the same role as the velocity variance σ 2 1 /3. The second term on the right-hand side in (76) describes diffusion analogous to (73).

Fluctuation-dissipation relations from the time-reversal symmetry
In statistical field theory, the time-reversal symmetry can be used to derive FDRs. For the Langevin dynamics (74) invariance of correlations under the symmetry (76) gives [16]: Due to causality the first term on the right-hand side has to vanish if t > t ′ and we arrive at the fluctuation-dissipation relation: In the following we aim at a similar derivation of the FDRs for KFT using the timereversal symmetry (70). Since the transformation T leaves correlation functions invariant, we can conclude for example: For t 2 > t 1 the left hand-side has to vanish due to causality. Using that the operator D (1, j) t 1 is equivalent to a regular time derivative for the auto-correlation, we arrive at the FDR (55): The time-reversal symmetry not only reproduces the FDR (55), but also gives rise to a whole hierarchy of relations between response functions of arbitrary order and time derivatives of correlation functions. Most generally, the time-reversal symmetry gives us: This proves that response functions of arbitrary order are related to diffusion.

Kinematic FDRs far from equilibrium
Since KFT is a theory far from equilibrium, there are some significant physical differences between FDRs in KFT and in thermal systems. The physical picture behind the FDR in a thermal system is the following. Particles in a thermal system necessarily interact with each other as this enables the system to distribute its total energy among the degrees of freedom according to equipartition. The interactions between particles lead to the diffusion of correlations, which is described by ∂ t O(t)O(t ′ ) eq in (56). If the system is slightly perturbed away from equilibrium, the energy of the system will be redistributed towards equipartition due to the interactions in the system. This will lead to a dissipation of the perturbation and a relaxation back to equilibrium, which is described by the left hand side of equation (56). Thus, the response of the system χ(t, t ′ ) and the diffusion ∂ t O(t)O(t ′ ) eq have the same physical origin, viz. the momentum transfer between the microscopic degrees of freedom.
The FDRs in KFT are derived within the free theory far away from equilibrium, thus the physical picture is different. Within a single realization of the initial conditions, the particle j has a constant momentum p (i) j . However, in the ensemble seen as a whole, the momentum of particle j is random with velocity variance σ 2 1 /3 due to the averaging over all realizations. The solutions to the free equations of motion with inhomogeneity K, see (19), describe the propagation of this random initial momentum in terms of the propagator g qp . However, the propagation of the inhomogeneity is described by g qp as well. The physical origin of the FDRs in KFT could thus be seen as a consequence of Newton's 2nd axiom: The inhomogeneity changes the momentum by an amount which then has to propagate like a momentum. The FDRs in KFT are thus purely kinematic arguments which makes them, however, not less valuable. In our derivation of the time-reversal symmetry we crucially used the Gaussian form of the initial conditions, so these kinematic FDRs seem to rely on the special form of the initial conditions.
In thermal systems, FDRs describe the linear response of the system to small departures from equilibrium. In KFT we are dealing with a system far from equilibrium and consider departures from the free theory. Keep in mind that the 'free' theory already contains some interactions since we use the Zel'dovich propagator, cf. (18). Through the time-reversal symmetry we were able to prove a whole hierarchy of higher-order FDRs which means that our interpretation of FDRs is not limited to linear departures from the free theory.
To our knowledge this type of kinematic FDRs far from equilibrium is a novel relation and is of significant interest because it can describe systems far away from equilibrium.

Conclusion and Outlook
In the context of the recently developed formalism of KFT [1,2] we have shown that the process of cosmic structure formation can be split into three processes: Particle diffusion due to the initial momentum variance σ 2 1 /3, accumulation of structure due to the initial conditional probability C p j p k between the momenta of two particles, and interactions relative to the inertial evolution.
We observed that the processes of diffusion and accumulation of structure are delicately balanced and for late times result in a net diffusion. The delicate balance is a consequence of the fact that our formalism allows to take the full non-linear coupling of free trajectories by initial momentum correlations into account and so far explicitly relies on the Gaussian form of the initial conditions.
Including the contributions from interactions in the Born approximation, we were able to compute the time derivative of the full non-linear density power spectrum and compare the Born approximation with simulations. We observed that the naive intuition, coming from perturbative approaches to cosmic structure formation, that the relative error of the approximation should be smaller at high redshifts and large scales, is not valid here as our approach is non-perturbative, but approximative. At low redshift, the relative error is of order 30 % even on highly non-linear scales up to k ∼ 5 h Mpc −1 . We postpone a thorough analysis of the time evolution of the relative error to a future study.
We observed a tendency of the net diffusion to approach the interaction contribution over time suggesting that diffusion and interaction are kinematically related to each other. The FDRs discussed in Section 4 show that the evolution of diffusion is indeed kinematically related to the reaction of the system to a gradient force. This result is independent of the form of the gravitational potential.
Although KFT describes a system far from equilibrium, we find that kinematic FDRs hold and in fact a whole hierarchy of higher order FDRs follows from a time-reversal symmetry of the generating functional. The Gaussian form of the initial conditions seems to be crucial for the validity of the FDRs as well as the form (19) of the trajectories showing that the initial momentum is propagated in the same way as an inhomogeneity.
Finally, we want to give a rather speculative outlook into future applications and relevance of the FDRs found in this work. It is yet unclear whether the FDRs are somehow related to the virial theorem, both being derived from kinematic arguments and relating the interactions and inertial movement in a system. The FDRs might, thus, be used to study virialization and the stable clustering regime. In fact our results in Fig. 6 indicate that the relative sum of net diffusion and interactions drops to zero on very small scales. If this remains true when taking post-Born approximations into account, which is necessary on these scales, this would imply that an equilibrium between clustering and diffusion is reached on very small scales. Furthermore, since the Gaussian form of the initial conditions plays an important role for the validity of the FDRs and the cancellation of the one-and two-particle contributions in Fig. 1, FDRs might be used in future studies to examine how non-Gaussianities effect the formation of a stable-clustering regime.
where c j (t) remains to be determined. In order to prove invariance of the generating functional we first have to check how the MSR-action, determining the dynamics of the system, changes under the transformation. After we have shown how the dynamics of the system transform, we use these dynamics to prove invariance of the generating functional and determine the appropriate form of the remaining free parameter c j (t) of the transformation.
× exp i dt J q j (t) q (i) j + g qp (t, 0) p (i) j − dt ′ G (adv) qp (t, t ′ ) K p j (t ′ ) = dΓ exp −i dt ′ K p j (t ′ ) · dt J q j (t)g qp (t, 0) exp i dt K p j (t) · c j (t) where we executed the derivative with respect to the initial momenta in the last step. Careful comparison with the non-transformed generating functional (A.15) shows that c j has to be a functional of the auxiliary field J and has the form: This finally proves the time-reversal symmetry (67) together with c j in the form (68).
where we used that p k (t 1 ) = p (i) k since no response fields are applied. In a final step we evaluate the integrals in the initial momenta: k 1 C −1 p j p k p k (t 1 )ρ j (1) . . . = V −N dq (i) k 1 C −1 p j p k ∂ i∂ L p k e − 1 2 L p l C p l pm L pm e i L qr q (i) r = iV −N dq (i) k 1 C −1 p j p k C p k p s L p s e − 1 2 L p l C p l pm L pm e i L qr q (i) where we used C −1 p j p k C p k p s = δ js ½ 3 .