Kinetic field theory: Non-linear cosmic power spectra in the mean-field approximation

We use the recently developed Kinetic Field Theory (KFT) for cosmic structure formation to show how non-linear power spectra for cosmic density fluctuations can be calculated in a mean-field approximation to the particle interactions. Our main result is a simple, closed and analytic, approximate expression for this power spectrum. This expression has two parameters characterising non-linear structure growth which can be calibrated within KFT itself. Using this self-calibration, the non-linear power spectrum agrees with results obtained from numerical simulations to within typically $\lesssim10\,\%$ up to wave numbers $k\lesssim10\,h\,\mathrm{Mpc}^{-1}$ at redshift $z = 0$. Adjusting the two parameters to optimise agreement with numerical simulations, the relative difference to numerical results shrinks to typically $\lesssim 5\,\%$. As part of the derivation of our mean-field approximation, we show that the effective interaction potential between dark-matter particles relative to Zel'dovich trajectories is sourced by non-linear cosmic density fluctuations only, and is approximately of Yukawa rather than Newtonian shape.

for the inertial motion of particles implies that the remaining particle-interaction potential is sourced only by the non-linearly evolved density contrast, which leads to an approximately Yukawa-shaped cut-off of the Newtonian gravitational potential. In Sect. 3, we briefly review the calculation of power spectra from kinetic field theory. In Sect. 4, we develop our mean-field approach to the particle-particle interaction term, and we summarise and discuss our results in Sect. 5.
We use the convention for the Fourier transform F and its inverse F −1 , with the short-hand notations Where needed, we adopt a (spatially flat) ΛCDM cosmological model with matter-density parameter Ω m0 = 0.3.

Particle dynamics 2.1 Particle trajectories
On the expanding spatial background of a Friedmann-Lemaître model universe with scale factor a, we introduce comoving coordinates q and use the linear growth factor D + as the time coordinate t. We set initial conditions at some time t i in the distant past when matter just began dominating the dynamics of the cosmic expansion. For later convenience, we set both the scale factor a and the growth factor D + to unity at the initial time such that t = D + − 1 and t i = 0 initially. In these coordinates, the Hamiltonian of point particles is with an effective, dimension-less, time-dependent particle mass m given by m = a 3 D + (a)E(a) , E(a) := H(a) H i (4) in terms of the derivative D + (a) = dD + (a)/da of the linear growth factor and the expansion function E(a), defined as the Hubble function H(a) divided by the Hubble constant H i at the initial time. The particle mass m used to be called g in our earlier papers on KFT, but we change the notation here to acquire a more intuitive meaning. Like a and D + , the expansion function E is supposed to be normalised to unity at the initial time such that m(t i ) = 1. In an Einstein-de Sitter universe, D + = a and E = a −3/2 , thus m = a 3/2 = (1 + t) 3/2 . The potential ϕ in (3) satisfies the Poisson equation sourced by the density contrast δ relative to the background density [7]. The Laplacian in (5) acts with respect to the comoving coordinates q.
Combining the phase-space coordinates of an individual particle j into a vector x j := ( q j ,˙ q j ) , the solution of the Hamiltonian equation of motion beginning at x (i) j = ( q (i) j , p (i) j ) at t i can be written in the form with the matrix-valued propagator G(t, t ) and the effective force f j (t ) on particle j, introduced and derived in Appendix A. Instead of individual particles, we consider canonical ensembles of N 1 classical particles j at positions q j with velocities˙ q j and introduce the tensor x = x j ⊗ e j to bundle the phase-space coordinates of the entire ensemble, where e j is the j-th Cartesian unit vector in N dimensions. This allows us to write the entire bundle of all particle trajectories in the compact form with G(t, t ) := G(t, t ) ⊗ 1 N and f (t ) := f j (t ) ⊗ e j . The inertial trajectories are x 0 = G(t, 0)x (i) , and y are the deviations therefrom caused by the effective force f .

Effective gravitational potential
The potential of the effective force acting relative to Zel'dovich trajectories is given by Eq.
where ψ is the potential of the curl-free initial velocity field. With the Poisson equations (5) for ϕ and ∇ 2 ψ = −δ (i) for the initial velocity potential ψ, the Poisson equation for φ is where δ (lin) = D + δ (i) is the linearly growing density contrast. The potential φ describing the gravitational interaction relative to the inertial particle trajectories x 0 (t) is thus sourced exclusively by the non-linearly evolved contribution to the density fluctuations. The effective gravitational force mediated by φ is therefore confined to small scales. This has important consequences for kinetic theory. If we describe inertial particle orbits with the Zel'dovich propagator, their mutual interaction must be modified in such a way that only the non-linear density contrast contributes to the force between them. This reflects the fact that the Zel'dovich propagator already takes the large-scale part of the gravitational interaction between the particles into account. Since the Zel'dovich trajectories reflect effective inertial motion with respect to the time coordinate t = D + − 1, and since the gravitational force caused by the linear density contrast relative to these trajectories needs to vanish, only the deviation of the density contrast from its linear value can be the source of the effective gravitational interaction. On large scales, where the density contrast keeps growing linearly for all cosmologically relevant times, and where the Zel'dovich approximation describes the particle motion accurately, the effective force must vanish. On small scales, where non-linear structures build up, the effective gravitational interaction must set in as non-linear density contrasts develop.

Shape of the effective potential
The effective gravitational potential between particles following Zel'dovich trajectories must thus deviate from the Newtonian form in such a scale-dependent way that the force tends to zero for k k 0 and approaches the Newtonian form for k k 0 , with the wave number k 0 set by the time-dependent boundary between linear and non-linear scales.
For quantifying how the potential needs to be modified, we search for a particle-particle interaction potential v such that the collective potential φ = nδ * v of the particle ensemble with the mean number density n, i.e. the convolution of the particle number-density fluctuation nδ with the potential v, satisfies the Poisson equation (9). Using the Fourier convolution theorem, the Fourier transform of this Poisson equation readsδṽ Multiplying this equation once withδ, once withδ (lin) , and taking the ensemble average gives introducing the power spectrum P δ of the density contrast. Going from (10) to (11), we have implicitly assumed that the form ofṽ is independent of the density contrast, which should be a good approximation, but does not generally need to be the case. Eliminating δδ (lin) between these equations leads to a quadratic equation forṽ whose only meaningful solution is The function f v (k) turns to zero for wave numbers k small enough to fall into the linear regime, and to unity for k large enough to be deeply in the non-linear regime. Power spectra obtained from numerical simulations [8][9][10] suggest that the function represents the transition from large to small scales reasonably well, with k 0 quantifying the wave number above which non-linear evolution begins to dominate (see Fig. 1). With (12), this results in the effective potentialṽ which is of Yukawa rather than Newtonian form. Note that the expression (12) for the particle-particle interaction potential v is statistical, i.e. it depends on the spatial correlations within the particle ensemble.
The essential result of this discussion is thus that the effective interaction between particles following Zel'dovich trajectories is mediated by an approximately Yukawa-like rather than a Newtonian potential. It is important for our purposes to note that the scale k 0 can be determined from KFT itself in a way to be described in Sect. 4.3. The right panel in Fig. 1 shows the time-dependent Yukawa scale k 0 determined in this way.  δ /P δ ) 1/2 , scaled to unity at k 1, compared to the fitting function f v (k) from Eq. (13). The linear power spectrum was taken from [11], the non-linear from [9]. Right: The Yukawa scale k 0 as a function of the redshift z.

Power spectra from KFT
We briefly review in this section the KFT approach to cosmic structure formation. For further detail, we refer the reader to [5], [12], and the review [6].

Generating functional
The central mathematical object of KFT is its generating functional Z. Like a partition sum in thermodynamics, it is a phase-space integral over the probability P(x) for the phase-space positions x to be occupied. Splitting P(x) into a probability P(x (i) ) for an initial state times a transition probability P(x|x (i) ) from the initial to the final state, further introducing a generator field to allow extracting moments of particle positions and momenta via functional derivatives with respect to J later, and introducing the particle trajectories x(t) from (7), leads to the generating functional (16) derived in [12], with the dot denoting the time-integrated scalar product According to (7), the phase in (16) can be split into a free and an interacting part, Since the contribution S I [J] := iJ · y to the phase depends on all relative positions of the correlated particle ensemble, it seems impossible to evaluate the generating functional Z[J] analytically. A systematic approach to perturbation theory may begin with converting the interaction term into an operator acting on the free generating functional Z 0 defined in (20) below, followed by Taylor-expanding this operator. We have previously shown that even the first order of this perturbative approach leads to non-linear density-fluctuation power spectra close to results from numerical simulations [5], and we will further analyse KFT perturbation theory in future papers. Here, we follow a different path.
The essential purpose of this paper is to find a suitable average for the interacting part, which would allow us to write the generating functional (18) as Before we proceed to construct and analyse such an average, we briefly review how power spectra are derived in KFT from the generating functional Z[J].

Density cumulants
The particle number density is a sum of delta distributions centered on the particle positions at time t or, in a Fourier representation,ρ Replacing the particle position q i in each one-particle density contributionρ i ( k, t) by a functional derivative with respect to J q i (t), we obtain the one-and N-particle density operatorŝ Since the density operators are exponentials of derivatives with respect to components of J, they generate translations of the generating functional, Density cumulants of order r are obtained by applying r density operators to the generating functional. For synchronous power spectra, i.e. cumulants of order r = 2 with t 1 = t 2 =: t, we have with the corresponding shift The short-hand notation (n) for the arguments in (23) indicates the Fourier-space position ( k n , t n ) at time t n . The generator field J can be set to zero once all density operators have acted on the generating functional. Since the particles of the ensemble are indistinguishable, each term under the sum in (23) gives the same result as for any particle pair arbitrarily labelled as i, j = 1, 2, and the cumulant becomes With initial conditions appropriate for the early universe, the free generating functional Z 0 [J] after applying two density operatorsρ 1 (1) andρ 2 (2) can be written as with the damping term and the freely evolved, non-linear power spectrum see [12] for the derivation. The function a appearing here is the auto-correlation function of momentum components parallel to the wave vector k 1 , where µ is the direction cosine of q relative to k 1 . The function containing the spherical Bessel function j 0 is the auto-correlation function of the initial velocity potential ψ, determined by the initial density-fluctuation power spectrum P (i) δ (k). The initial velocity field is supposed to be the gradient of a velocity potential because any initial curl would decay quickly due to cosmic expansion and angular-momentum conservation. We further define the moments of the density-fluctuation power spectrum and note that For small arguments of the first exponential in (28), P turns into the linearly evolved power spectrum, as shown in [12].

Damping and interaction
The damping term Q D in (26) requires a careful discussion. Its definition (27) in conjunction with the velocity dispersion σ 2 1 from (31) shows that it arises because particles stream freely with an average velocity σ 1 in the Zel'dovich approximation. We should emphasise that this velocity dispersion is not of thermal origin, but arises from drawing initial particle velocities from a velocity potential which is a homogeneous and isotropic Gaussian random field. Where this velocity field converges, structures form, but these structures are smoothed in a system of free particles once caustics have been formed and converging particle streams have crossed.
In (26), it appears as if the damping term would exponentially reduce the power. However, this is not the case. Rather, numerical integration and asymptotic analysis alike show that the freely-evolved power spectrumP combined with the damping term follows the linearly evolved power spectrum P (lin) δ on large scales, drops below it on non-linear scales, but turns towards an asymptotic behaviour ∝ k −3 as k increases further [12].
Combining (25) and (26) to shows that the gravitational interaction between the particles counteracts this characteristic reduction of the power. In fact, a major part of the interaction term is required for keeping structures in place once they have formed. Any surplus, i.e. any positive difference between S I and Q D , leads to non-linear structure growth.
4 Mean-field approach to non-linear power spectra Based on these arguments, we now pursue the following approach. We wish to represent the particle interactions by a suitably averaged interaction term S I [L] . For simplicity, we further wish to approximate the damped, freely evolving power spectrumP by the linearly evolving power spectrum P (lin) δ . As just discussed, this implies that we are ignoring the reduction of power by the velocity variance after shell crossing. Since the interaction term counter-acts the damping, we then also need to reduce the interaction term on non-linear scales. We will do so by using Burgers' approximation. The result will be the simple expression for the non-linear power spectrum. Our main result will be a specific and simple equation for S I (k) reproducing numerically derived, non-linear power spectra remarkably well.

Averaged particle-particle force
With y from (7) and L from (24), we have where f i (t) is the effective force on an arbitrary particle i. Note again that the time t here is not the cosmological time, but defined by the linear growth factor D + . In terms of an effective particle-particle interaction force f p and the particle number density ρ, we can write the force f i on particle i as We now average this force term over particle ensembles drawn from a statistically homogeneous, correlated random density field. Since we wish to retain the dependence of the resulting, averaged, effective force term f i on the wave vector k, we project out the contribution by the mode k of the density field, writing The average over the product of densities introduces the correlation function ξ(| q 1 − q 2 |) of the density field, We keep only the connected part of the correlation expressed by ξ alone in (39) because the disconnected part cannot contribute in a homogeneous random field. Owing to homogeneity, the integrand in (38) depends only on the difference q 1 − q 2 =: q of position vectors. We can thus integrate over q 1 , resulting in a factor V, and obtain Since the remaining expression is the Fourier transform of a product, the scale-dependent, mean force term is the convolution of the particle-particle forcef p in Fourier space and the power spectrum P δ of the particle distribution as the Fourier transform of the correlation function, Combining (41) with (36) and using Newton's third axiom in the form f i (− k ) = − f i ( k ), we thus find the expression for the averaged, scale-dependent, interaction term.

Damping within Burgers' approximation
We now need to specify the power spectrum P δ to be inserted into (41) for evaluating the average force term f i . Following (38), evaluating the mean interaction term with the density correlation function suggests replacing P δ =P. Applying the same linear approximation leading from (34) to (35), we would then arrive at P δ ≈ P (lin) δ . However, since the linear power spectrum in (35) ignores damping and thus overestimates the power on non-linear scales, we need to reduce the average force term on these scales appropriately. The amount of this reduction can be effectively estimated by Burgers' approximation [13][14][15].
Burgers' approximation changes the inertial motion of particles in the Zel'dovich time coordinate in a way derived from the Navier-Stokes equation of hydrodynamics, where ν is a viscosity parameter with the dimension of a squared length. A natural choice for the length scale ν 1/2 is the non-linear radius r nl defined by evaluated with the linear power spectrum from [11] and a top-hat window function W R (k) at the present cosmic time. We set ν = r 2 nl = const here for simplicity, but note that ν could also be generalised to become time-dependent.
Burgers' equation can be solved by a Hopf-Cole transformation [16,17], which results iṅ where U is an exponential velocity potential given by the convolution of the scaled, exponentiated initial velocity potential ψ with a normal distribution N of width (2νt) 1/2 . To linear order in ψ, (45) implies the velocity dispersion This expression clearly shows the effect of Burgers' approximation: small-scale modes with k (2νt) −1/2 are removed from the velocity field, slowing down the particles on such scales and thus reducing the re-expansion of structures after stream crossing [18][19][20][21]. We take this reduced amount of particle motion into account by replacing the damping term Q D from (27) bȳ An excellent fit to λ(t) is The form of this fit expresses the transition from ballistic to diffusive particle motion.

Averaged interaction term
Accordingly, we evaluate the averaged force term (41) as Inserting the averaged force term (50) into (42) results in the averaged interaction term We finally need to evaluate the convolution of the particle-particle force termf p with the damped power spectrumP. This is done in Appendix B and results in the average interaction term where σ 2 J is the moment (83) of the damped power spectrum. This interaction term S I in the meanfield approximation is shown in the right panel of Fig. 4 as a function of wave number k for redshift z = 0. As it has to be, the averaged interaction term is dimension-less.
The scale k 0 of the Fourier-transformed Yukawa-like potential (14) still needs to be set. KFT itself now suggests the following procedure. According to (35), the ratio between the linearly and non-linearly evolved power spectra is estimated by the exponential of the averaged interaction term. Combining (35) with (12), the factor f v from (12) is approximated by We can thus determine a first estimate for k 0 by calculating S I with k 0 = 0 and fitting f v (k) from (53) with the functional form (13). With the resulting value of k 0 , an updated estimate for S I can be calculated, and so forth. This iteration quickly converges; it turns out that even one step suffices. The scale k 0 determined from KFT in this way is shown in the right panel of Fig. 1 as a function of redshift z. For simplicity, we adopt the constant value of k 0 at z = 0 and ignore its time dependence.

Non-linear density-fluctuation power spectrum
The mean-field averaged interaction term (52), inserted into (35), is our approximate expression for the non-linear density-fluctuation power spectrum. Evaluating the averaged interaction term S I (k) from (52) numerically, determining the viscosity parameter ν and the Yukawa scale k 0 from KFT itself as described, and further assuming a spatially-flat ΛCDM model universe with matter-density parameter Ω m0 = 0.3, leads to the result shown in Fig. 2.   [11], the numerical spectrum from [9]. The flat lower panel shows the relative deviation between the analytic and the numerical power spectra.
In this Figure, our analytic approximation (35) is compared at redshift z = 0 to the description by [9] of the power spectrum obtained from numerical simulations. As can be seen there, the analytic power spectrum agrees within typically 10% with the numerical expectation up to wave numbers of k 10 h Mpc −1 . With the parameters ν and k 0 self-calibrated from within KFT, the analytic expression (35) has no adjustable parameters since the Yukawa scale k 0 is set by KFT itself, and the viscosity is set to the square of the non-linear scale determined by the linear power spectrum from [11]. Expression (35) is also non-perturbative in the sense that the original exponential interaction operator is not expanded into a power series, but averaged in a mean-field approach.
We should emphasise that the derivation of the mean-field approximation is mathematically not fully rigorous, but in several steps guided by some intuition and the principles of statistical field theory. These are that we evaluate the mean interaction term with the linearly evolved power spectrum, reduce its damping by means of Burgers equation, and approximate the gravitational potential of the particles by a Yukawa form. Nonetheless, the agreement between our mean-field approximated analytic and the numerical results well into the non-linear regime of cosmic structure formation suggests that the microscopic approach of kinetic field theory, combined with a suitable choice for the inertial reference motion and adapting the effective force between particles to this reference motion, captures essential aspects of the physics of large-scale cosmic structure formation. The notorious shell-crossing problem does not occur in this approach, which is the main reason for the possibility to extend it far into the non-linear regime.
Instead of self-calibrating the two parameters ν and k 0 from within KFT, they can be considered as free parameters of the theory and chosen to optimise the agreement between numerical power spectra and the mean-field expression (35) by minimising the squared difference between them. Doing so, using the power spectrum from [9] as a reference, results in the power spectrum shown in Fig. 3.   Fig. 2, but with the parameters ν and k 0 modified to optimise the agreement with the numerical results by [9].
The power spectrum shown there is obtained by reducing the Yukawa scale by 13% and increasing the displacement (48) by 6%. The relative deviation of the mean-field approximated, analytic power spectrum from its numerical counterpart is now lowered to typically 5 % up to k 10 h Mpc −1 . For even smaller scales, the analytic power spectrum falls below the numerical expectation because then the mean-field approximation of the interaction term is no longer strong enough.

Summary and conclusion
The kinetic field theory for classical particle ensembles encapsulates the statistical information on the initial state and the propagator for the equations of motion in a generating functional which is closely analogous to the canonical or grand-canonical partition sum in thermodynamics. This generating functional evolves in time. Statistical macroscopic information is obtained from it by applying suitable operators. In this paper, we have used this approach to derive an analytic expression for the non-linear power spectrum of cosmic density fluctuations in a mean-field approximation of the particle-particle interaction term. Our main results are the closed, analytic approximation (35) for the non-linear power spectrum P (nl) δ and the expression (52) for the mean-field averaged interaction term. We have derived this form of the interaction term from KFT, averaging it over particle ensembles drawn from a statistically homogeneous, Gaussian random density field as shown in (38) and (42). This derivation is not mathematically rigorous because we have bypassed several complications and subtleties for the sake of simplicity. A rigorous assessment of the mean-field approximation needs to be based on systematic perturbation theory and will be worked out in a forthcoming paper. Nonetheless, the agreement with numerical results up to wave numbers k 10 h Mpc −1 is very good and encouraging.
It is important for our result that the microscopic, Hamiltonian equations of motion allow the introduction of an inertial motion with respect to the time t = D + − 1, corresponding to the celebrated Zel'dovich approximation, which captures the linear evolution of cosmic structures on large scales. Linearly growing, large-scale density fluctuations must then not exert any gravitational force on the inertial particle trajectories with respect to this time coordinate. This requires us to replace the Newtonian gravitational potential by an approximately Yukawa-shaped gravitational potential which ensures that only small-scale, non-linearly growing modes contribute to the particle-particle interaction. The Yukawa scale k 0 can be determined from kinetic field theory itself in a quickly converging iteration. We emphasise that the Yukawa shape is suggestive, but approximate and has no fundamental justification yet.
Our aim expressed in (35) to capture the non-linear evolution of the density-fluctuation power spectrum simply by a multiplicative, exponential interaction term applied to the linearly evolved power spectrum requires us to damp part of the interaction term on small scales. We do so by means of the damping term naturally appearing in the mean-field expression for the interaction term via KFT, but lowering the damping scale in a way derived from Burgers' equation. This introduces a viscosity parameter ν, which is the square of a length scale characterising non-linear structures. A natural choice for this length scale is the non-linear radius defined in (44).
We thus have two parameters, k 0 and ν, which can either be set by KFT itself or seen as free parameters. Self-calibrating both parameters with KFT leads to the mean-field approximated, nonlinear power spectrum shown in Fig. 2, which already agrees well with numerical results. This agreement can further be improved by slightly adjusting both parameters, as shown in Fig. 3.
The initial state of the microscopic degrees of freedom is fully determined by the linear densityfluctuation power spectrum at the initial time, which can (and should) be set as early as the onset of the matter-dominated epoch. We have used the cold-dark matter power spectrum here. Since we use the growth factor D + of linear density fluctuations as a time coordinate, the cosmological framework model enters only through the relation between redshift or scale factor and time, and through the time dependence of the effective particle mass. It can thus easily be generalised towards alternative darkmatter models, a different cosmological background, or alternative gravity theories. Non-linear cosmic power spectra such as these shown in Figs. 2 and 3 can be calculated within seconds on conventional laptops.
The approach followed in this paper can be improved in several ways. So far, we substantially simplified our mean-field approximation scheme, and we have modelled the particle-particle interaction potential by a Yukawa form for intuitive simplicity. The detailed form of the interaction potential could, however, also be derived from kinetic field theory itself; this would just cause the calculation of the mean interaction term to become more involved. The results shown should be seen as a further step towards a systematic interpretation of non-linear cosmic structures in terms of fundamental physics.

A Particle Trajectories
A.1 Particle mass Beginning with the effective particle mass m defined in (4), the time derivative of m iṡ by definition of the time t = D + − 1. The prime denotes differentiation with respect to the scale factor a. We can insert (4) once more into (54) to arrive aṫ The growth factor D + solves the linear growth equation. When transformed to the scale factor a as an independent variable, this reads Expressing the matter-density parameter Ω m by its value Ω (i) m at the initial time, If we specify the initial conditions early in the matter-dominated epoch, we may further approximate Ω (i) m ≈ 1. We thus find the expressionṁ = 3 2 for the time derivative of the effective particle mass m. Inserting finally the potential amplitude A ϕ from the Poisson equation (5), we arrive at the time derivativė of the effective particle mass m.

A.2 Solution of the Hamiltonian equations of motion
The Hamiltonian equations of motion for the phase-space point x = ( q, p ) can be written in the forṁ Notice that A(t) is a 6 × 6 matrix, with 0 3 and 1 3 representing the zero and unit matrices in three dimensions, respectively. The homogeneous equation˙ SinceĀ is nilpotent,Ā 2 (t, t ) = 0 6 , the homogeneous solution shrinks to x h (t) = 1 +Ā(t, 0) x 0 .
By variation of the constant vector x 0 , the inhomogeneous equation of motion (60) leads to and thus to the solution for phase-space trajectories beginning at x (i) = ( q (i) , p (i) ) at t = 0. The spatial trajectories are accordingly q(t) = q (i) + g H (t, 0) p (i) =: q 0 (t) − t 0 dt g H (t, t )m ∇ϕ =: q 0 (t) + y q (t) .

A.3 Reference Trajectories
It is often convenient in cosmology to replace the force-free trajectories q 0 (t) = q (i) + g H (t, 0) p (i) by the trajectories postulated in the Zel'dovich approximation, q 0 (t) → q (i) + t p (i) , exchanging the Hamiltonian propagator g H (t, t ) for (t − t ). The deviation y q (t) defined in (65) as the difference between the actual and these reference trajectories is then determined by Implicitly defining an amplitude A p (t ) by t 0 dt 1 −ġ H (t , 0) we can write (67) in the form Since the initial momentum is the gradient of a potential, p (i) = ∇ψ, (69) suggests defining an effective potential such that Differentiating (68) twice with respect to the time t and usingġ H (t, t ) = m −1 (t) gives A p (t) =ṁ = mA ϕ D + with (59), and thus the trajectories Continutity demands that the initial velocity potential satisfies the Poisson equation ∇ 2 ψ = −δ (i) .

A.4 Unifying Propagators
It is often convenient to replace the propagator g H (t, t ) in (73) also by the time difference t − t . For doing so, we implicitly introduce an effective force f demanding and solve for f (t). Differentiating (74) twice with respect to t, we find With this expression for f (t), we can write the spatial trajectories as q(t) = q (i) + t p (i) +