Kinetic Field Theory: Exact free evolution of Gaussian phase-space correlations

In recent work we developed a description of cosmic large scale structure formation in terms of non-equilibrium ensembles of classical particles, with time evolution obtained in the framework of a statistical field theory. In these works, the initial Gaussian correlations between particles have so far been treated perturbatively or restricted to pure momentum correlations. Here we treat the correlations between all phase-space coordinates exactly by adopting a diagrammatic language for the different forms of correlations, directly inspired by the Mayer cluster expansion. We will demonstrate that explicit expressions for phase-space density cumulants of arbitrary $n$-point order, which fully capture the non-linear coupling of free streaming kinematics due to initial correlations, can be obtained from a simple set of Feynman rules. These cumulants will be the foundation for further investigations of interacting perturbation theory.


Introduction
Based on the works of Das & Mazenko [1,2,3,4] we have developed a variant of kinetic theory [5,6,7] we have dubbed 'Kinetic Field Theory' (KFT). Our primary goal is to find a new approach to cosmic structure formation which is free from the limitations of established methods like Eulerian standard perturbation theory (SPT). For a discussion of the advantages we refer to our previous work and announce an upcoming paper comparing KFT and SPT directly in more detail.
The central idea of KFT is to use the path integral approach to classical mechanics (cf. [8,9,10]) in order to directly encode the microscopic phase-space dynamics of individual particles in a generating functional. Collective macroscopic fields are then constructed from the microscopic information at the time of interest. We adopted this idea in order to describe the time evolution of the non-equilibrium statistics of a system of gravitationally interacting particles whose initial positions and momenta are correlated in such a way that they sample Gaussian density and velocity fields.
In [5] we employed perturbation theory in the interaction potential relative to a modified version of Zel'dovich trajectories [11]. This reproduced the non-linear growth of the CDM density fluctuation power spectrum known from N-body simulations over a wide range of scales remarkably well considering we included perturbative corrections only to first order. In this first work the initial correlations were also perturbatively expanded into orders of the initial density fluctuation power spectrum. This allowed straightforward numerical evaluation but restricted the validity of even the free theory to certain time and length scales and made the formalism appear less elegant than it actually is.
Consequently, we layed out a more general approach in [6] focusing on an ensemble of particles correlated in their momenta only. By taking these correlations into account exactly and artificially reducing damping effects (see [12] for a motivation of this workaround) we could show that a substantial part of the late-time non-linear growth behaviour of the density fluctuation power spectrum is actually already encoded in the initial correlations. This further corroborates the hope that low-order perturbation theory relative to Zel'dovich-type trajectories will yield accurate results for the evolution of density statistics.
Restricting initial correlations to the momenta only is an accurate description of cosmic structure formation at late times when using a particle propagator which is not bound from above like in the case of Zel'dovich-type trajectories, because momentum correlations scale with higher polynomial orders of the propagator compared to density correlations. However, when using a strict separation between non-interacting kinematic motion and interactions the particle propagator will be bound from above [11] due to cosmic expansion. Considering this case will be necessary for a direct comparison of the non-interacting evolution of density statistics of both KFT and SPT.
More importantly, we will need the cumulants including the exact initial correlations between all phase-space coordinates as the building blocks for an alternative formulation of perturbation theory built around macroscopic fields while still retaining the underlying particle dynamics. This alternative approach will also be developed in a separate paper. It will naturally lead to a partial resummation of the canonical perturbation theory in [5]. In order for this resummed theory to reproduce the familiar linear growth behaviour of the matter power spectrum known from SPT, one needs to include density and momentum correlations on equal footing.
In the main body of this work we concentrate on the conceptual steps leading to our central result, i. e. the general form of non-interacting collective field cumulants in terms of diagrams evaluated according to Feynman rules. Detailed calculations and derivations can be found in the Appendices. Since our results only depend on very general properties of the force acting between particles, our results in the present work should in principle have a much broader range of applicability than just cosmic structure formation. We will thus endeavour to keep the mathematical formulation of our results as general as possible, even though we will most often try to make connections to cosmology when commenting on them.
In Section 2 we give a short recapitulation of the KFT approach by setting up the generating functional for a particle ensemble described initially by a Gaussian random field of density and momentum. The focus of Section 3 lies on re-ordering the free generating functional in a way that reflects the particle-based nature of the theory more clearly. The main result will be a grand canonical version of the generating functional in the form of a sum over connected diagrams representing initial particle correlations ordered by the number of correlated particles. In Sections 4 and 5 we explain the calculation of exact and explicit expressions for non-interacting n-point cumulants of the collective fields. We extend the results of [5] from the spatial density ρ to the full phase-space density f . Their calculation essentially reduces to simple Feynman rules for the initial correlation diagrams. Finally, Section 6 will contain a summary of the paper and present our conclusions.

Notation
We will follow the notation of [5] for the microscopic phase-space coordinates of a large set of N particles confined to a volume V. Individual particles are enumerated with particle labels j = 1, . . . , N and have positions and momenta denoted as d-dimensional vectors q j and p j . We combine those into a 2d-dimensional phase-space coordinate vector as For all N particles we bundle these vectors with the help of the tensor product q = q j ⊗ e j , p = p j ⊗ e j , x = x j ⊗ e j , (1.2) where e j is the Cartesian base vector in N dimensions with entries ( e j ) i = δ i j . The Einstein summation convention is implied unless explicitly stated otherwise or obvious from context. The bold vectors from (1.2) follow the rules of matrix multiplication inducing a scalar product a · b = a b = a j ⊗ e j b k ⊗ e k = a j b k δ jk = a j · b j . (1.3) For functions of time we extend the dot product to include an implied time integration as (1.4) We will also encounter tuples of fields defined on the domains of time and space. For two tuples A, B, the dot product additionally includes an implied integration over Fourier space as which we extend in an analogous fashion to matrix products. We have abbreviated field arguments as field labels (±r) (± s r , t r ) = (± k r , ± l r , t r ), where s is Fourier-conjugate to x, k to q and l to p. We also denote integration over phase-space and its conjugate Fourier space by (1.6)

Kinetic field theory
We will only briefly go through the main aspects of KFT, for a more detailed introduction see [5] or alternatively the original works by Das & Mazenko in [1,2,3,4]. As already mentioned in the introduction, the central idea of the KFT approach is to encapsulate both the dynamics and the initial statistics of a non-equilibrium N-particle system in a generating functional. Time evolution is accounted for by functionally integrating over all possible trajectories but forcing them to the classical solution of the equations of motion by means of a functional Dirac delta distribution. Stochasticity is introduced into the systems by averaging over the initial conditions with some initial phase-space probability distribution P(x (i) ). The canonical generating functional thus reads where E[x(t)] = 0 is the equation of motion for the N-particle system. We have rewritten the functional Dirac delta distribution as a Fourier transform by introducing new auxiliary fields χ(t). They play the role of the functional Fourier conjugates to the phase-space coordinates x(t). We will deal with systems with Hamiltonian equations of motion for which we assume the form i. e. the matrices F j encapsulate the free part of the equation of motion for individual particles and V( q, t) is the interaction potential. In order to separate free motion from interactions we first introduce generator fields J j , K j paired with x j , χ j respectively. Under the path integral, we can then replace .
Any function of these operators may now be pulled in front of the integrals. We restrict ourselves to systems of N identical particles in the absence of external forces and assume that the interaction potential between particles will only depend on the distance between their respective positions at the same point in time, implying that V may be written as a superposition of N single-particle potentials v. This naturally leads to the introduction [5,7] of two fundamental collective fields Φ(r) = (Φ f (r), Φ B (r)) . The Klimontovich phase-space density is then defined as It encodes the instantaneous occupation of a phase-space state x by the particles of the ensemble, which in the conjugate Fourier space turns into a superposition of phase-factors. According to (2.3), its operator version iŝ (1) The deviation of all individual particles from their inertial trajectories due to the interaction with all other particles is encoded in the response field Again using (2.3) we obtain its operator version The interaction potential is encapsulated in the matrix It will be advantageous to introduce source fields H f , H B paired with the collective fields Φ f , Φ B . The appropriate source term in Fourier space is given by where α = f, B and the hattedĤ operators allow to generate arbitrary factors of Φ by repeated application. Following [5,7] the generating functional may then be written as where we abbreviated dΓ i = dx (i) P(x (i) ). With the generating functional in this form, cumulants of the collective fields are generated by taking appropriate functional derivatives and then setting all sources to zero, (2.11) While the particle interactions will not be important for the remainder of this work, we include them in the appropriate places in order to show that they do not interfere with any of our calculations. In the case of (2.10) it is important to note that we were able to pull the interaction operator in front of the integration over the initial conditions which means the latter can be performed before any interactions need to be considered.
We will mostly concern ourselves with the free generating functional Z C,0 [J, K]. As shown in [5], the path integrals defining it can be executed once the Green's function G of the free equation of motion of a single particle is known. One finds where the solutionx(t) to the free and linear equations of motion is given bȳ 13) and the non-interacting particle propagator has the general form 14) The appearance of the additional source term K in the free solution allows for the generation of deviations from the inertial motion of particles. Since the auxiliary fields χ j appear in the response field (2.6) and are replaced according to (2.7), we now explicitly see how deviations from trajectories couple to interactions with the potential and thus form the actual perturbed physical quantity in KFT perturbation theory.

Initial phase-space probability distribution
Both the theory of inflation and observations of the cosmic microwave background imply that the initial cosmic density ρ (i) ( q) and momentum P (i) ( q) fields at the beginning of the matter dominated epoch should together constitute a Gaussian random field. Apart from its importance in cosmology, this initial condition is also an instructive example for the simplest possibility to include correlations, since any Gaussian field is completely defined by its mean and its covariance matrix. Due to the cosmological principle it must be statistically homogeneous and isotropic, implying where the average is taken over the Gaussian distribution of the field values. The homogeneous mean particle density is given byρ (i) = N/V and local fluctuations around it are described by the density contrast δ (i) . After centering the field appropriately, it is completely fixed by the N-point covariance matrix C, which we decompose into two-particle submatrices according to Homogeneity and isotropy require that any C jk may only depend on the distance | q (i) Isotropy furthermore implies that a rest frame of the matter distribution must exist in which C δ j p k = 0 for all j = k and C p j p k is a diagonal matrix with all equal entries. If this was not the case, there would be a preferred direction in the initial momentum field violating isotropy. Homogeneity implies that the self-correlation matrix of individual particles must be spatially constant and have the diagonal form As was shown in [5], we can apply Poisson sampling to individual field configurations and then average over their Gaussian distribution to obtain the initial phase-space distribution of our particle ensemble as where we introduced C pp = C p j p k ⊗ E jk . While the particle momentum correlations inherit the Gaussian nature of the momentum density field, the density and density-momentum correlations are contained in the polynomial expression The label set {n} contains all particle labels 1, . . . , N. The primed set {n} excludes the labels i, j chosen in the preceding sum over ordered pairs. This also applies to all higher order terms. The prime on the sum over ordered tupels of ordered pairs indicates that no label may be shared between the pairs in any individual term of the sum. This scheme is repeated up to the terms which have N/2 factors of C δ i δ j .

The initial correlation operator
The next step consists of inserting the explicit form of the initial phase-space distribution (2.18) into the free generating functional defined in (2.10). It will be advantageous to rewrite the result such that the effect of initial correlations becomes easier to understand. For this we need to split the covariance matrix of momentum auto-correlations as We define the uncorrelated Gaussian distribution specified by σ 2 p as and assemble all other correlations into the operator After some simple calculations detailed in Appendix A, we can write the free generating functional using the above definitions as Note that if we turn off all initial correlations thenĈ tot → 1 and we recover the generating functional of a non-interacting gas with momentum dispersion σ 2 p and uniform distribution in space. The total correlation operatorĈ tot induces shifts of the initial phase-space positions of the particles relative to this reference system. ‡ This is very similar to the standard scheme for the initial setup of correlations in cosmological N-body simulations [13], where the Zeldovich approximation is applied to a regular grid or uncorrelated glass state.
This shows that while the free dynamics of the particles can be obtained exactly, initial Gaussian correlations introduce another operatorĈ tot which couples the degrees of freedom of different particles, in addition to the interaction operator. Thus, even in the free regime we will have to perform an expansion in the number of particles being correlated. However, we will later find that statistical homogeneity truncates this expansion at the number n f of Φ f -fields in the cumulant.

Cluster expansion
In the caseĈ tot → 1, the free generating functional (3.4) can be factorised into oneparticle contributions, i. e. all integrations over phase-space coordinates and auxiliary fields of different particles can be separated from one another. Since the collective fields Φ f (2.4) ‡ Note that we could also describe initial correlations as a shift relative to a uniform system at rest if we had not split off the one-particle momentum variance in (3.1). and Φ B (2.6) can also be separated into one-particle contributions, the free evolution of its cumulants can then be obtained by calculating the generating functional of only a single particle. Starting from (2.10) we apply e iH·Φ to (3.4) and then use the appropriate one-particle version of (3.2) to define this for any particle j as the trace operation Tr j The free generating functional would just be the product over N of these identical traces. In the caseĈ tot 1, the correlation operator couples particles so that integrals over their respective coordinates can no longer be separated, where the correlation operator still acts to the left on the exponential contained in the trace. Any factorisation must now be performed in terms of 'clusters' which are formed from a subset of all particles due to their pairwise correlations in some term ofĈ tot . This concept is essentially the same as the Mayer cluster expansion [14], where the configurational part of the partition sum of an interacting gas is expanded in the number of particles taking part in interactions.
There exists a very illustrative visualisation of this scheme in terms of Feynman-like diagrams which we adapt from [15]. A detailed derivation of the form of diagram lines and the rules constraining their combination into diagrams representing terms ofĈ tot can be found in Appendix B. The four fundamental line types representing correlation operators between particle dots i and j are All these lines depend on the relative coordinate distance q (i) i j through the correlation functions. The diagrams which can be constructed from these four line types are constrained by three rules, which are derived from the form of the initial phase-space distribution (2.18).
No self-correlation rule: Subdiagrams of the form for any kind of correlation line are forbidden.
No 2-particle loop rule: Any pair of particles can be connected directly by at most one line. This means that the subdiagram for any combination of correlation lines is forbidden.
Note that this does not exclude loops involving three particles or more. The third rule is the δ-Line Rule: No particle may have more than one solid δ-line attached to it, thus all subdiagrams which include are forbidden.
The correlation operator (3.3) can then be obtained by summing over all ordered n-tuples of particle labels, where 2 ≤ n ≤ N, and for each of them summing all diagrams one can draw under the above rules into a fixed graphical arrangement of particle labels representing their ordering.We write this in a schematic way aŝ To illustrate this, consider the contribution from three particle diagrams representing pure momentum correlations, which according to the above rules are We call any connected subdiagram of particles an -cluster, i. e. those subdiagrams through which one can trace a continuous path. Starting at the four-particle level, disconnected diagrams will appear inĈ tot . Factorising their contribution to the generating functional into connected diagrams, i. e. clusters, is the central idea of the cluster expansion. Consider for example 1 3 For any diagram one can count the number m of clusters of size , where in this case we find m 1 = N − 11, m 2 = 2, m 3 = 1, m 4 = 1 and zero for all larger clusters. One collects all diagrams with the same set {m } of numbers into a cluster configuration. Since there are only N particles any {m } must obey the constraint As we show in Appendix B, the sum over all diagrams in a cluster configuration factorises into products of the free -particle cluster generating functional of connected correlations Here Σ ( ) CDiag (1, . . . , ) is the sum of all connected -particle diagrams obtained with the same logic as Σ (n) Diag . Due to the integrations contained in the trace operator this quantity has the same value for any set of particles chosen from the overall N particles, thus we define it in terms of a representative set of particles labeled 1, . . . , . We then obtain the free generating functional by summing over all possible cluster configurations where the asterisk in the sum represents the constraint (3.11).

Grand canonical generating functional
In order to write down a grand canonical generating functional for our system we need a probability distribution P(N) for the number of particles, which replaces the constraint on the particle number N. Finding P(N) for a general non-equilibrium system is a formidable task because it may in principle depend on the complete instantaneous microscopic phase-space configuration x(t). However, we are considering a statistically homogeneous and isotropic system. In the absence of external forces and if the interaction potential also respects these symmetries, they are conserved in time. If we thus can derive P(N) at the initial time t i by employing only these two symmetries, it must also be valid at any later time.
Consider the usual setup where we embed our grand canonical system with volume V GC into a much larger canonical system with particle number N C and volume V C . Due to statistical homogeneity all choices for the location of the embedded subvolume V GC are equivalent. Homogeneity also implies that the probability for one of the N C particles to be in the subvolume V GC must be given by p = V GC /V C . The number of particles inside the subvolume is then binomially distributed, i. e. N GC ∼ B(N C , p). Taking the limit N C , V C → ∞ with the mean particle densityρ N C /V C =ρ (i) = const. and p → 0, the binomial distribution of N GC turns into a Poisson distribution. Concentrating on the grand canonical subsystem and dropping all suffixes, i. e. V GC = V and N GC = N, we thus find (3.14) The grand canonical generating functional is now obtained by averaging the canonical one with (3.14). Using (3.6) and that the interaction operator in (2.10) is independent of particle number, we find The modified trace operatorTr amounts to replacing 1/V →ρ in (3.5). From here on we will always understand the modified trace to be used and drop the tilde. The prefactor e − N can either be absorbed into the normalisation of the path integrals or dropped altogether since we only consider functional derivatives of ln Z GC later on. Using (3.13) and according to [15] we can execute the sum over particle numbers to obtain the free grand canonical generating functional as In the first step we use that first summing over all cluster configurations with the N-particle constraint {m } * and then summing over all particle numbers is identical to summing over all cluster configurations without this constraint. In the second line we reorder the sum in terms of cluster size rather than cluster configurations and then use the series definition of the exponential function.
Readers familiar with quantum field theory will recognize (3.16) as the usual relation between disconnected and connected Feynman diagrams, i. e. the sum of all disconnected diagrams is given by the exponential of the sum of all connected diagrams. The only difference is that diagrams are now physically interpreted as representing initial phase-space correlations. We will in fact see shortly that they can be evaluated by applying simple rules very similar to those of traditional Feynman diagrams.

General form of cumulants
Grand canonical cumulants are defined in the same manner as in (2.11). In the non-interacting regime the cumulant generating functional reads This simple sum makes calculating free cumulants a lot more comfortable when compared to the canonical ensemble. It also shows that the notion of connected particle correlation diagrams directly translates into connected correlations of collective fields and that this notion is conserved under free motion. From the definition of the one-particle trace (3.5) it is trivial to see that where bold vectors are now defined for the set of representative particles,Φ ( ) is the corresponding collective field operator andx is the free solution (2.13) for the phase-space trajectories of all particles. A general free cumulant can be written as where general free -particle cumulants were defined in the last line. The product of collective field operators in (4.3) can be written out into a sum of n terms with the help of (2.5) and (2.7), each of them being a product of n single particle operators. In every such term each field label 1, . . . , n appears exactly once. We will use the notion that some particle j carries such a label if the corresponding operatorΦ (1) f j (1) is present. Statistical homogeneity enforces certain constraints on the two-particle correlations represented by diagram lines. In Appendix C we show how these lead to the Homogeneity Restriction: Such terms in the operator productΦ ( ) α 1 (1) . . .Φ ( ) α n (n) vanish in which at least one of the particles j ∈ {1, . . . , } carries no field label.
The physical picture behind this is best understood by considering the simple example of the two-particle contribution to the two-point Φ f cumulant. If we ask for the Fourier mode contribution of the first particle a to both Φ f fields, i. e. the particle carries both field labels, the second particle b can only influence the result due to its initial correlation with a. But since we also need to integrate the initial position of b over all of space, these correlations must cancel out in a statistically homogeneous situation, leading to a factor of zero.
If we only consider pure phase-space density cumulants G (0) f (1)··· f (n) this has the important consequence that because with each field label appearing exactly once in every term of the operator product, each term will have at least one particle not carrying a label if n < , so they all vanish. We will extend this statement to mixed cumulants including response fields in section 5. Statistical homogeneity thus naturally truncates the sum over particle numbers in (4.3) at the number of phase-density fields Φ f present in the cumulant. This has a very important implication. Since the diagrammatic rules of section 3.2 impose a finite limit on the possible number of diagrams for a fixed number of particles, we can write down exact and explicit expressions accessible to numerical evaluation for any given non-interacting n-point cumulant. In a separate paper, we will argue that one cannot do this in the Eulerian standard perturbation theory of cosmic structure formation.
We can now change our perspective and only specify how field labels are grouped into a collection of sets without specifying which particle carries which set. We call such a collection a label grouping. It turns out that summing over all possibilities of assigning sets to particles merely cancels the factor of 1/ ! in (4.3). For any grouping, we may thus freely choose one arbitrary assignment of label sets to particles by enumerating the collection of sets as I j , j = 1, . . . , . A general phase-space density cumulant is then simply a sum over all possible groupings The quantityΣ ( ) CDiag is defined as the Fourier transformed sum of connected -particle correlation diagrams The uncorrelated Gaussian momentum distribution P σ 2 p ( p (i) ) in (4.2) has turned into a Gaussian damping factor e −Q(I 1 ,...,I ) exp which has been discussed at length in [6,12]. The effects of applying collective field operators are parametrised in terms of phase shift vectors L q,I j (t) = r ∈ I j L q,r (t) r ∈ I j k r g qq (t r , t) + l r g pq (t r , t) , (4.8) r ∈ I j k r g qp (t r , t) + l r g pp (t r , t) . (4.9) These are the generalisations of the shift vectors of our previous works [5,6] to the complete phase-space. They encode how the phases of the Fourier factors contributing to the phasespace density field labels r = ( k r , l r , t r ) change due to the free motion from time t to time t r of those particles which carry these field labels. Causality, i. e. t r ≥ t, is encoded in the propagators.

Feynman rules
The final step is to evaluate the Fourier transformΣ ( ) CDiag (4.6) of the connected correlation diagrams. This is easily done be following a set of rules that are very similar to Feynman diagram rules in quantum field theory. These rules are derived in Appendix D. The basic idea is to extend the Fourier transforms in (4.6) from coordinates q (i) j centered on the -th coordinate to all relative coordinates q (i) jk . The Fourier transformed correlation lines in the diagrams are then defined as where j < k in all cases by definition. Along each line there is thus a flow of a general Fourier momentum K (i) jk from smaller to larger labels. The algorithm for evaluating (4.6) for a given grouping {I 1 , . . . , I } of field labels is then as follows.
(i) Choose a fixed graphical arrangement of dots labeled j = 1, . . . , representing the particles carrying the label sets I j . These are the vertices of the diagrams.
(ii) With this fixed set of vertices, draw all possible diagrams using the lines in (4.10) while obeying the rules detailed in 3.2. Repeat the next three steps for all diagrams.
(iii) For any closed loop in a diagram, assign one of the lines in the loop connecting vertices labeled j < k a loop momentum k (i) jk and introduce an integral k (i) jk .
(iv) Pick a vertex j which has only one line left with undetermined momentum and use the conservation law to fix the momentum of this line. Incoming momenta from vertices with smaller labels are counted positive while outgoing momenta to vertices with larger labels are counted negative. Momenta associated with vertices not connected to vertex j are zero.
(v) Consecutively go through all other vertices to fix the remaining undetermined momenta by repeatedly applying the previous step.
One may always use that the leading Dirac delta distribution in (4.5) ensures the global conservation of the position phase-shift vectors L q,I 1 + . . . + L q,I = 0. Loop momenta k (i) jk may be shifted freely in position Fourier space and general momenta K (i) jk can be reversed in sign because due to statistical isotropy the correlation functions only depend on the absolute value | q (i) jk | of the separation vectors. These rules are best understood by considering some easy examples. The most trivial are two-particle cumulants. A general grouping of the field labels 1, . . . , n is given by I 1 , I 2 . The sum of all possible two-particle diagrams is Since there are no loops we directly go to step (iv) of our algorithm and pick vertex 1.
According to our convention, we have to count the general momentum K (i) 12 as outgoing leading to the conservation law L q,I 1 − K (i) 12 = 0 ⇒ K (i) 12 = L q,I 1 . (4.13) We thus find that the sum of correlation diagrams evaluates tõ L q,I 1 (4.14) According to (4.5) the two-particle contribution to the n-point phase-space density cumulant is then given by In order to demonstrate the determination of Fourier momenta when loops are involved we consider a simple one-loop diagram of four particles.
In the first diagram we have chosen to assign the loop momentum to the line connecting vertices 1 and 3. We then go through the vertices 1, 2 and 3 consecutively and for each of them fix the undetermined K (i) i j by using the conservation law (4.11).

Derived cumulants
Preserving the full microscopic phase-space statistics by working in terms of the phase-space density Φ f has the benefit that at any n-point level, once the corresponding Φ f -cumulant (4.5) has been calculated, one can directly derive cumulants between any observables that can be defined in terms of Φ f . Momentum moments are of special interest since they lead to particle density and velocity density statistics. These are defined as integrals over real momentum space whereas the cumulants (4.5) are calculated in Fourier momentum space. We thus need to transform the Fourier momentum of the external label under consideration back into real momentum, multiply by the appropriate factors of momentum and integrate over momentum space. For some field label r this leads to the operator for the m-th order moment, where i a ∈ {1, . . . , d} ∀a ∈ {1, . . . , m}. The n-point particle density cumulant is particularly straightforward since it is the zeroth order moment in all Φ f entries, (4.20) As an explicit example, let us consider the two-point particle density cumulant G (0) ρ(1)ρ (2) for a 3D-system of particles on straight trajectories with g qq = 1 and g pq = 0. We set l 1 = l 2 = 0 according to (4.20). This reduces the phase-shift vectors (4.8) and (4.9) to L q,r (t i ) = k r and L p,r (t i ) = g qp (t r , t i ) k r . According to the homogeneity restriction the only non shot-noise contribution is the two-particle term. There is only one possible grouping of the field labels 1, 2 into two non-empty sets. We choose to enumerate them as I 1 = {1}, I 2 = {2}. Setting t i = 0 and t 1 = t 2 = t, we use the general result (4.15) to find We want to emphasize at this point that besides the assumption of statistical homogeneity no approximation or perturbative expansion has gone into deriving this explicit expression. The underlying microscopic physics driving the time evolution of the initially Gaussian density field encompasses the entirety of its free streaming kinematics including all effects that must be considered non-linear at the level of macroscopic fields. Since velocity density statistics will be derived from the same fundamental Φ f -cumulants this feature will apply to them as well.

Mixed cumulants and causality
So far we only considered the explicit calculation of Φ f -cumulants. However, when calculating perturbative corrections in the particle interaction we also need to use mixed cumulants of both Φ f and the response field Φ B . According to (2.7) the single particleΦ (1) B joperator is the product of the response factor operatorb j and the phase-space density operator Φ (1) f j . In any term of the operator productΦ ( ) α 1 (1) . . .Φ ( ) α n (n) of (4.3), we may thus push all density operators to the right and execute them first.
Theb j contain only derivatives w. r. t. to the K p j sources which do not couple to the initial phase-space coordinates and may be pulled in front of the integral over them. Mixed cumulants can thus always be derived from pure Φ f cumulants by multiplying evaluated terms in the operator product with the appropriate response factors b j resulting from application of operatorsb j . By combining (C.6) with (2.7) we find these response factors to be given by The physical interpretation is straightforward. The spatial gradient in (2.7) has been Fourier transformed into k u . Multiplication by a potential v(k u ) and integration over k u will generate a force onto particle j leading to an instantaneous momentum change δ p at time t u . The corresponding response of the Fourier phase contribution of particle j to the phase-space density at the field labels in I j is then evolved under free motion forward in time with the particle propagator G contained in the momentum phase shift vector. In order to write down a general mixed cumulant we denote the response factors as b I(r) b j (r), with I j = I(r) the set containing the field label r. A mixed cumulant is then given by We now take a closer look at the causal structure of these response factors. We first recognise that in any perturbative calculation aiming at the statistics of physical observables, i. e. cumulants of Φ f and integrals thereof, every Φ B -field will always appear in combination with the σ-matrix (2.8) of particle interactions. For such calculations it is thus advantageous to introduce a new collective field directly in Fourier space as Its physical content is just what we described following (5.1), i. e. it gives the linear response of the system per phase-space density Φ f at the Fourier phase-space point r and as such must always be contracted with another Φ f field appearing in the perturbation series, e. g.
. Since the σ-matrix only depends on field labels, it does not interfere with any of our previous calculations and we can obtain cumulants involving Φ F as compound cumulants by contracting mixed cumulants with the σ-matrix, Since we assume particle interactions to act instantaneously we have σ f B (r , −r) ∝ δ d (t r − tr), so we can always set t r = tr from here on and consider Φ B -labels and Φ Flabels equivalent for causal considerations. Since we also describe the interaction by a purely space-dependent potential we obtain the factor δ d ( l r ) δ d ( l r ) in the σ(r, r )-matrix (2.8). If we combine this with g qp (t, t) = 0, we see that for compound cumulants (5.4) we only ever need to consider effective b j factors without instantaneous response, Let us generalise our arguments following (5.1). We can imagine the individual phaseshift factors as they appear in (5.5) as arrows connecting dots representing points in time, where the direction of the arrow is given by the time ordering of the propagator. For a compound cumulant (5.4), let the subset {ū 1 , . . . ,ū m } ⊂ {1, . . . ,n B } be carried by particle j, i. e. {ū 1 , . . . ,ū m } ⊂ I j . The product of response factors b j is then If we write this out into a sum, each individual term is a 'time flow' diagram, where arrows represent how the effects of interactions on particle j at some instance of time 'flow' with time to some other field label. This other field label may in turn also be an instance where the particle interacts and again the effect needs to flow to a field with another label. The sum of all these diagrams thus represents all possible time orderings of interaction events for particle j given the set I j of field labels it carries.
As an example consider a particle j carrying the set I j = {ū 1 , . . . ,ū 5 , r}, i. e. carrying five Φ B -labels and the Φ f -label r without loss of generality. One possible flow diagram would bē The only ab initio topological restriction on these time flow diagrams resulting from their definition (5.6) is that all Φ B -labelsū 1 , . . . ,ū m ∈ I j appear exactly once as a dot in every diagram, have exactly one line line pointing away from them and that any Φ f -label r ∈ I j may only appear as an end-point of time flow. In Appendix E we show the additional restriction that due to causality there may be no closed loops in these diagrams. If we combine these restrictions we see that all possible diagrams must be composed of tree-like structures with an unambigious direction of time flow terminating at some Φ f -label r (compare also with the diagrams in [6]). Consider again the above example and assume e. g. t r < tū 5 , then the diagram immediately vanishes sinceū 5 r ∝ Θ(t r − tū 5 ). Also, if there is no Φ f -label, time must flow from theū 5 -dot to some otherū i -dot which creates a closed loop.
For the compound cumulants these topological restrictions translate into the Causality Restriction: Only such field label groupings {I 1 , . . . , I } give non-vanishing contributions to compound cumulants (5.4), which have the following property: For all I j , j = 1, . . . , there must be at least one Φ f -label r ∈ I j such that t r ≥ tū = t u for all Φ B -labelsū ∈ I j .
This has the direct consequence that . . , n B with t u > t r ∀r ∈ 1, . . . , n f . (5.7) Even more importantly, we can combine both the causality and the homogeneity restriction to obtain a truncation criterion for compound cumulants Adding any number of Φ B -labels to the case of (4.4) will still leave particles carrying no labels or only Φ B -labels, so again all terms vanish. This shows that not only can we write down exact and explicit cumulants in the non-interacting regime, but also that at any fixed order of perturbation theory in the interaction we can account for all effects of initial correlations coupling freely streaming particles with a finite number of perturbative terms.

Conclusions
Following up on previous work in [5,6] we refined the treatment of non-interacting cumulants of collective fields in the framework of KFT. We extended the position-space density to the complete phase-space density and included all initial correlations of a Gaussian density and velocity field in an exact fashion. The latter was achieved by condensing the effects of the complicated initial phase-space distribution (2.18) into a coupling operator (3.3) between particles acting at the initial time. We saw that the idea of how KFT incorporates initial correlations is actually quite close to the initial setup scheme of cosmological N-body simulations in that it applies phase-space shifts to an uncorrelated system. We expressed this correlation operator in a diagrammatic language with rules derived from the form of the initial phase-space distribution (2.18). Using the principles of the Mayer cluster expansion we factorised the free generating functional into connected diagrams. Exploiting the statistical homogeneity and isotropy of our system we changed our description to a grand canonical ensemble which made the calculation of cumulants of the collective fields straightforward. We saw that the notion of connected correlations of collective fields is equivalent to connected particle diagrams in the initial correlations and that this notion is conserved under free evolution.
We then re-derived the general scheme of how to calculate such cumulants for the full phase-space density Φ f . As shown in Appendix B and Appendix D, the use of a diagrammatic approach allows us to work out all tedious combinatorics in a general way such that the calculation of any cumulant can be reduced to evaluating all possible correlation diagrams according to familiar Feynman-type rules without the need for any kind of symmetry factors. We saw that the combination of statistical homogeneity and causal consistency of interaction events impose restrictions on which kind of terms may contribute to given free and compound cumulants. These restrictions naturally truncate the number of particles one has to consider at the n f -point order of Φ f -fields in the cumulant. At any given order of perturbation theory in the interaction we can thus in principle obtain results for the time evolution of phase-space statistics which take the full non-linear coupling of free streaming kinematics due to initial correlations into account with only a finite number of explicit terms.
This result also has the profound consequence that in the non-interacting regime KFT generates exact time-evolved cumulants of a Gaussian density field. We gave the explicit result for the two-point density cumulant in (4.21), i. e. the density power spectrum which is of prime interest in cosmology. We will show in a separate paper that this is next to impossible to achieve in the Eulerian approach of SPT, which includes gravity at the linear level at the price of giving up knowledge of the exact free dynamics. The linear growth behaviour of cosmic density fluctuations derived in linear SPT is however easily recovered in a resummed form of KFT perturbation theory without losing any information about the free dynamics. This result will be developed in a separate paper, where our insights into the causal structure of compound cumulants obtained in the present work will prove beneficial for understanding the causal flow through entire terms in perturbative expansions.
We have not presented any numerical evaluation of quantities like (4.21) because implementing the necessary Fourier transforms (4.10) of the correlation lines is a quite formidable task due to very large dynamic range of scales that must be considered and is thus beyond the scope of this work. For numerical evaluations of the density fluctuation power spectrum restricted to initial momentum correlations we refer the reader to [6,12]. In future work, we will endeavour to extend these calculations to include all correlation types in order to test the free theory against non-interacting N-body simulations.
Nonetheless it is obvious that once (4.10) is precomputed for general arguments, the presented approach should be well suited for automated calculation of arbitrary free cumulants by a computer. The idea is the same as in [6] in that one specifies the desired cumulant, and a symbolic code generates the appropriate correlation diagrams which can then be evaluated according to the Feynman rules presented in section 4.2. This would then serve as a basis for perturbative and resummed treatments of the interacting system.
Beyond this application to cosmic structure formation our results should be applicable with little modification to any system which is statistically homogeneous and whose dynamics can be described in terms of an Hamiltonian with an interaction respecting homogeneity. Two examples would be the statistical evolution of both cold Rydberg gases and classical spin lattices.
Using (A.2) and (3.2) we may rewrite the exponential in (A.4) as We can rewrite factors ofJ p in terms of functional derivatives by using With the definition of the correlation operator (3.3) we finally rewrite the free generating functional as which is equivalent to (3.4) after rewriting the free solution in terms of path integrals.

Appendix B. Correlation diagrams and cluster expansion
We need to express the elementary building blocks of the correlation operatorĈ tot (3.3) in terms of line diagrams. It will be easier to define at first only three line types, which while different from those in (3.7) will have some notational overlap. OnceĈ tot is expressed in terms of these three elementary line types we redefine the meaning of diagram lines as a combination of the elementary lines types to arrive at (3.7). We begin by introducinĝ where Einstein summation is not implied. Note that the symmetry of the diagram reflects the symmetry of C p j p k . The above term corresponds to a Mayer function f i j = e −βv i j − 1 in the Mayer cluster expansion. We can now express the Gaussian inĈ tot as In the first line we used the symmetry of C p i p j to reduce the sum to ordered pairs of particles. When expanding out the product over them into sums we have to exclude identical pairs and prevent overcounting, hence the ordering i 1 < i 2 between pairs. We now want to rewrite this not as a sum over multiple ordered pairs but over ordered n-tupels. In the diagrammatic language this means we connect lines sharing a common particle label, i. e. those with j n = j m in the above case of (B.2). The fact that we only sum over ordered tupels of ordered pairs i n < i m , ∀n < m translates into the p-Line Rule: Any pair of particles may at most be connected by one line directly.
This means that the subdiagram is forbidden.
Writing down any n-tupel sum in (B.2) then amounts to drawing and summing up all diagrams possible under this rule for a fixed graphical arrangement of a general n-tupel of particle labels representing their ordering by ascending numerical value and then summing over all possible realisations of that tupel. This leads to where N is the number of particles in the system. Notice that since we removed the momentum self-correlation fromĈ tot in (3.1) no subdiagrams of the form appear in the above expression.
Turning to density auto-correlations and density-momentum cross-correlations we introduce diagrams We can write the operatorĈ inĈ tot following the same logic used in the step from (B.2) to (B.3), i. e. we multiply out all products over particle labels into sums of ordered tupels of labels by connecting lines at identical particle labels. However, we need to take into account that the form (2.19) ofĈ imposes certain restrictions on the possible diagrams.
• In any term of the ordered sums {{i 1 < j 1 },...,{i n < j n }} i 1 <...<in over the labels of the C δ i k δ j k any particle label may only appear once. We may thus not connect two -lines with one another.
• In the products {n} 1 + m nĈδn p m each label n pertaining to δ n only appears once.
We may thus not connect two -lines with their solid δ end.
• The above products do not share labels with the above sums, as indicated by the prime. This forbids connecting the solid δ end of a -line to a -line.
These three restrictions can be combined into the δ-line rule given in section 3.2, albeit that for now it holds for the three elementary line types. Together with the m n restriction in the second of the above points the δ-line rule also excludes any form of self-correlations, i. e. one-particle loops. TheĈ-operator may thus be written aŝ Diag | δδ,δp ( j 1 , . . . , j n ) The last step is to multiplyĈ (B.6) with the Gaussian exponential (B.3). Since there are no restrictions on connecting with either or , we simply obtainĈ tot as the sum over all diagrams that can be constructed from the three line types subject to both the δ-line and p-line rules. We illustrate this by writing it out for the two-particle contribution In this form of the diagrammatic language the complete three-particle contribution would already contain a plethora of diagrams with up to six lines. We can remedy this to some degree by recognizing the trivial fact that the diagrams in (B.7) are also the building blocks of all diagrams with three particles or more. By combining diagrams and thus reducing their number in (B.7) we effectively do so in all higher order terms as well. The form of (B.7) suggests the diagram redefinition The meaning of the line is unchanged. It is straightforward to see by checking all possible combinations that for these new diagrams the former δ-line rule and p-line rule combine into the same δ-line rule and a general rule forbidding two-particle loops. Together with the exclusion of one-particle loops these are the three diagram rules given in section 3.2. The operatorĈ tot is then expressed as in (3.8), where the diagram sum now includes all diagrams that can be drawn with the redefined line types under the three restriction rules.
In order to arrive at (3.13) we now need to reorderĈ tot and through it the free generating functional Z C,0 (3.4)   We now add all other possible connected diagrams for the three-particle cluster (1, 2, 3) to (B.10) while keeping the rest of the diagram fixed. Afterwards we do the same for the twoparticle cluster (4,5) and, using the definition of the connected -particle generating functional From this we directly see that the contribution of any particle assignment in a cluster configuration {m } to the free generating functional (3.6) may be written as The last step is to figure out how many particle assignments are actually present in (3.6) per given cluster configuration. By simply rearranging the particle order freely in e. g. [ (1,2,3), (4,5), (6), . . . , (N)], there are in principle N! possible assignments for any {m }. However, we defined the correlation operator (B.7) as a sum over ordered tupels of labels. Once we have fixed which particles are clustered together we must thus only include one assignment realising this clustering. We therefore divide out all N =1 ! m reorderings of particle labels inside the clusters and all N =1 m ! ways to exchange all particles between two clusters of equal size. This finally leads to the expression (3.13) for the free generating functional where the asterisk in the sum represents the constraint (3.11).

Appendix C. Operator evaluation and homogeneity restriction
In order to calculate general cumulants (4.3) we need to evaluate the collective field operatorŝ Φ ( ) α r (r) and the correlation operators (3.7) represented by the diagrams. We first introduce the unaveraged -particle free generating functional We write out the product of collective field operators in (4.3) into individual terms and in each of them pull all one-particle phase-space density operators to the right. As already shown in [5], the effect of applying such an operator iŝ This is easily extended to the case where the particle j carries some subset I j of the field labels 1, . . . , n by defining the following phase shift vector which allows us to write the effect of applying all single-particle phase-space density operators in some term of (4.3) as a shift of the J sourcẽ This shows that knowledge of L is equivalent to the explicit grouping of field labels into sets I j carried by particles j. Assume for now that all functional derivatives w. r. t. K p (t i ) have been executed, which leaves no operators remaining and we may therefore set all sources J and K to zero. For the unaveraged generating functional we then find where the one-particle phase shift vectors L q,I j , L p,I j have been defined in (4.8), (4.9). Both the single particle response operatorsb j and all correlation operators are defined purely in terms of operatorsˆ χ p j . For a given shift vector L they evaluate tô The physical interpretation ofˆ χ p j is thus that it generates the phase shift vector which when multiplied with some momentum perturbation ∆ p gives the linear response of the system in form of a Fourier phase shift. As argued in section 3.1, the initial correlations introduce such momentum perturbations relative to the uncorrelated system. The correlation operators represented by lines (3.7) in correlation diagrams thus evaluate to j kZ and analogously for the three other line types. We now understand Σ ( ) CDiag (L) as the sum over all -particle diagrams with lines evaluated in the sense of (C.7). A general -particle phasespace density cumulant is then written as Note that so far the sum over all possible L also includes all cases where label sets may be empty. We now show how the homogeneity restriction of section 4.1 comes about by picking any term, where without loss of generality particle j carries no field label, i. e. I j = ∅. Since particle j is involved in the cumulant it must appear in the diagrams of Σ ( ) CDiag (L) and be connected with one of the four line types shown in (3.7). If the particle is connected by more than one line, then the δ-line rule demands that at most one of the lines connect to j with a dashed p-side. But since I j = ∅ we have L p,I j = 0 and by combining (C.7) with (3.7) we see that in this case any corresponding correlation factor is zero. If the particle is connected to only one line, we need to distinguish three different cases to complete the proof: • The particle is connected to a dashed p-type line so the same argument as above applies.
• It is connected to a C δ j δ k line. Since I j = ∅ we also have L q,I j = 0 and according to (C.5) Thus the only remaining quantity in the term that depends on q (i) j is C δ j δ k . Using the definition of the density contrast gives • It is connected to the solid δ-side of aĈ δ j p k . The same arguments as in the previous case apply and we find The homogeneity restriction excludes all shift vectors L corresponding to terms with empty label sets from the sum in (C.8). We now make the change of perspective to label groupings where we only specify which labels are grouped into non-empty sets enumerated as I m without assigning them to particles j. We use the convention that writing the sets in the order (I 1 , . . . , I ) is equivalent to using the identity mapping between set labels m and particle labels j, i. e. particle j carries set I m with m = j. All other possible ! assignments of this specific grouping are then obtained by all permutations of the order of the sets. At this point the particle labels are essentially redundant information, so we may always understand labels j as pertaining to the sets of I j of the grouping.
The sum over L in (C.8) can then be written as summing over all possible label groupings and then summing each of them over all permutations of ordering the label sets. With the exception of Σ ( ) CDiag all other quantities in (C.8) are invariant under set order permutations. We now use the fact that Σ ( ) CDiag does not depend on the initial momenta of the particles to execute the integration over them. By combining the momentum part of (C.5) with (3.2), this leads to the Gaussian damping factor e −Q(I 1 ,...,I ) dp (i) P σ 2 where we factorised the uncorrelated Gaussian distribution into its individual particle contributions in the second step. At this point we may write the general phase-space density cumulant as (1) , . . . , I π( ) . (C.10)

Appendix D. Derivation of Feynman rules
We first need to show that the sum over permutations of field-label sets in (C.10) cancels against the preceding factor of 1/ !. For this it is advantageous to go back to the definition (3.12) of W ( ) 0 and classify the connected diagrams in terms of an equivalence relation. For a given representative diagram of some fixed ordered n-tupel of particle labels { j 1 < . . . < j n } we define the sum over its permutation equivalence class as the sum of all diagrams that can be transformed into the representative under some permutation of these particle labels. For the first three diagrams in (3.9) this would be The second line shows that we could also define the sum over the equivalence class as the sum of the representative diagram over all possible permutations of particle labels modulo those permutations which leave the diagram invariant, i. e. those which define the symmetry of the diagram. This particular diagram has a symmetry factor of S = 2 which leads us to exclude the diagrams in order to avoid overcounting. In the original definition of the permutation equivalence class this overcounting was prevented by the fixed order of the particle labels. In general there are ! possible permutations of a given ordered -tuple of particle labels. Consider someparticle diagram representative R, then it must hold that ! = S(R) N [R] for the number N [R] of diagrams in its equivalence class.
Since the -particle trace is invariant under particle label permutations, all diagrams in a permutation equivalence class give the same contribution to W ( ) 0 . For each permutation equivalence class we may choose a representative, giving us a set R ( ) a , a = 1, . . . , M ( ) ∼ . In terms of these we can write the -cluster generating functional as . The second step uses the fact that summing the representative over all permutations results in the sum over the entire equivalence class, but without excluding symmetric terms like in (D.2). These cancel against the symmetry factors. The fact that due to statistical homogeneity initial correlations between particles only depend on their relative distance has been exploited in [6] to factorize correlator contributions into products of a general quantity P, which is essentially given by the Fourier transform of the -line. We repeat the same logic in the present case where we consider all types of correlations and not only momentum correlations. Without loss of generality we shift the origin of our coordinate system to the initial position q (i) leading to new coordinates Relative distance vectors are then given by where j, k < . (D.6) Since Σ ( ) CDiag only depends on the q (i) jk we can execute the integration over q (i) to find Inserting (D.4) and (D.7) into (C.10) then leads to the general form (4.5) for a phase-space density cumulant. By introducing appropriate Dirac delta distributions that enforce the constraints (D.6) and rewriting them in terms of their Fourier transforms we can extend the integration to all relative position vectors This introduces Fourier vectors K (i) jk conjugate to the q (i) jk with j < k which are defined as Applying these Fourier transforms to Σ ( ) CDiag then allows us to transform each correlation line independently which leads to the general definition of Fourier transformed correlation lines given in (4.10), where P p j p k corresponds to the P jk already discussed in [6].
One can easily check that the Fourier momenta (D.9) obey the conservation law (4.11). We first consider the case j .
The case j = is only slightly more complicated. Since K (i) i = 0 by definition we find The Feynman rules for assigning momenta follow from the fact that any open pair j k of particle vertices j < k in a diagram of Σ ( ) CDiag leads to a factor (D.12) We can first go through all open pairs where k and thus K (i) jk = k (i) jk . For these we can execute the appropriate k (i) jk from (D.8) to obtain a factor of unity and drop k (i) jk from all conservation laws. For any remaining open pairs j < k with k = we can choose and execute one of the remaining k (i) ab to replace k (i) ab in all conservation laws with the constraint obtained from solving the equation K (i) j = 0 for k (i) ab . Once all Dirac delta distributions of the form (D.12) have been exploited, some non-trivial k (i) ab -integrations will remain as the loop momenta of our Feynman-rules, i. e. their choice is equivalent to the choice of the k (i) ab -integrations we use to resolve the δ d ( K (i) j ). To see that this conforms to the familiar idea of loops in Feynman diagrams, consider the number of particle pairs N pairs , the number of lines P playing the role of Feynman-propagators, the number of particle vertices V = and the total number I of k (i) ab -integrals we had to introduce. If L is the number of non-trivial k (i) ab -integrals we find L = I − (N pairs − P) = (N pairs − ( − 1)) − (N pairs − P) = 1 + P − V , (D. 13) which is the familiar QFT relation between loops, propagators and vertices. We finally need to acknowledge that our choice of q (i) as the center of our new coordinate system was arbitrary. Any other choice would still lead to a set of generalised Fourier momenta like in (D.9) which obey the conservation law (4.11) and allow to integrate out momenta associated with open pairs. We therefore can always reconstruct the appropriate set of Fourier momenta directly from the conservation law by following the Feynman rules of section 4.2.

Appendix E. Time flow loops
We need to show that any loop-subdiagram in the time-flow diagram language of (5.6) leads to a factor of zero. In such a loop, each dot corresponding to a field labelū must be the origin of an arrow line, thus there must be a response factor b j (ū) andū is therefore a Φ B -label. In a compound cumulant (5.4) such labels are always paired with a σ-matrix element which forces lū = 0. The phase-shift vectors corresponding to the arrow lines in the loop thus reduce to tū tr = i L p,r (tū) · kū lr=0 = g qp (tr, tū) kr · kū ∝ Θ (tr − tū) . We also note that if we integrate out the momentum space information from the collective fields, i. e. we only work in terms of the spatial particle density Φ ρ as in [5], the shift vectors L p,r have the form (E.1) even before contracting mixed cumulants with the σ-matrix. In this case the causality restriction applies directly to the bare mixed cumulants.