\section*{Abstract}
{\bf 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.}

\vspace{10pt}
\noindent\rule{\textwidth}{1pt}
\tableofcontents\thispagestyle{fancy}
\noindent\rule{\textwidth}{1pt}
\vspace{10pt}

\section{Introduction}
\label{sec:1}

Kinetic Field Theory (KFT) describes ensembles of classical particles in and out of equilibrium \cite{2010PhRvE..81f1102M, 2012JSP...149..643D, 2013JSP...152..159D}. It is based upon the Martin-Siggia-Rose approach to classical statistical systems \cite{1973PhRvA...8..423M} and has been adapted to cosmological initial conditions and to the expanding cosmological background in previous papers \cite{2016NJPh...18d3020B, 2019AdP...531..446B}. Its central mathematical object is a generating functional encapsulating the statistical properties of the initial state, the Green's function or propagator of the equations of motion, and the particle-particle interactions. These interactions are described by an exponential operator acting on the free generating functional. In the conventional approach to statistical field theories, this operator is expanded into a Taylor series, leading to a systematic approach to perturbation theory in terms of Feynman diagrams.

In this paper, we show that the interaction operator can instead be approximated as an averaged interaction term using a mean-field approach. This can be done in such a way that its action on the generating functional can be separated from the integration over the initial phase-space distribution. This results in a numerical, time and scale-dependent factor multiplying the free generating functional. Averaging over a pair of density factors then leads to an approximate, but closed and analytic expression for the non-linear power spectrum of cosmic density fluctuations.

The mean-field approximation introduces two parameters, the non-linear scale and an effective viscosity reducing the velocity variance after shell-crossing.
Both of them can be calibrated from within KFT itself. With these parameters self-calibrated in this way, our mean-field approximation to the non-linear power spectrum agrees with results from numerical simulations with a relative deviation of typically $\lesssim 10\%$ up to $k\approx 10\,h\,\mathrm{Mpc}^{-1}$ at redshift $z = 0$. Alternatively, these parameters can be optimised to further improve the agreement between the non-linear power spectra from our mean-field approximation and from numerical simulations. Doing so, the relative deviation to numerical results for $\Lambda$CDM can be lowered to $\lesssim 5\,\%$ in the same range of wave numbers.

In Sect.~\ref{sec:2}, we discuss the trajectories of Hamiltonian particles in the expanding cosmic space-time. Using results derived in detail in Appendix~\ref{app:A}, we show that the conventional Zel'dovich approximation 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.~\ref{sec:3}, we briefly review the calculation of power spectra from kinetic field theory. In Sect.~\ref{sec:4}, we develop our mean-field approach to the particle-particle interaction term, and we summarise and discuss our results in Sect.~\ref{sec:5}.

We use the convention
\begin{equation}
  \mathcal{F}\left[f\right] =: \tilde f\left(\vec k\,\right) =
  \int_qf\left(\vec q\,\right)\E^{-\I\vec k\cdot\vec q}\;,\quad
  \mathcal{F}^{-1}\left[\tilde f\,\right] = f\left(\vec q\,\right) =
  \int_k\tilde f\left(\vec k\,\right)\E^{\I\vec k\cdot\vec q}
\label{eq:1}
\end{equation}
for the Fourier transform $\mathcal{F}$ and its inverse $\mathcal{F}^{-1}$, with the short-hand notations
\begin{equation}
  \int_q := \int\D^3q\;,\quad\int_k := \int\frac{\D^3k}{(2\pi)^3}\;.
\label{eq:2}
\end{equation}
Where needed, we adopt a (spatially flat) $\Lambda$CDM cosmological model with matter-density parameter $\Omega_\mathrm{m0} = 0.3$.

\section{Particle dynamics}
\label{sec:2}

\subsection{Particle trajectories}
\label{sec:2.1}

On the expanding spatial background of a Friedmann-Lemaître model universe with scale factor $a$, we introduce comoving coordinates $\vec q$ and use the linear growth factor $D_+$ as the time coordinate $t$. We set initial conditions at some time $t_\mathrm{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_\mathrm{i} = 0$ initially. In these coordinates, the Hamiltonian of point particles is
\begin{equation}
  \mathcal{H} = \frac{\vec p^{\,2}}{2m}+m\varphi
\label{eq:3}
\end{equation}
with an effective, dimension-less, time-dependent particle mass $m$ given by
\begin{equation}
  m = a^3D_+'(a)E(a)\;,\quad
  E(a) := \frac{H(a)}{H_\mathrm{i}}
\label{eq:4}
\end{equation}
in terms of the derivative $D_+'(a) = \D D_+(a)/\D a$ of the linear growth factor and the expansion function $E(a)$, defined as the Hubble function $H(a)$ divided by the Hubble constant $H_\mathrm{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_\mathrm{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 $\varphi$ in (\ref{eq:3}) satisfies the Poisson equation
\begin{equation}
  \vec\nabla^2\varphi = A_\varphi\delta\quad\mbox{with}\quad
  A_\varphi := \frac{3}{2}\,\Omega_\mathrm{m}^\mathrm{(i)}\,\frac{a}{m^2}
\label{eq:5}
\end{equation}
sourced by the density contrast $\delta$ relative to the background density \cite{2015PhRvD..91h3524B}. The Laplacian in (\ref{eq:5}) acts with respect to the comoving coordinates $\vec q$.

Combining the phase-space coordinates of an individual particle $j$ into a vector $\vec x_j := (\vec q_j, \dot{\vec q}_j\,)^\top$, the solution of the Hamiltonian equation of motion beginning at $\vec x_j^{\,\mathrm{(i)}} = (\vec q_j^{\,\mathrm{(i)}}, \vec p_j^{\,\mathrm{(i)}})$ at $t_\mathrm{i}$ can be written in the form
\begin{equation}
  \vec x_j(t) = G(t,0)\vec x_j^{\,\mathrm{(i)}}+
  \int_0^t\D t'\,G(t,t')\cvector{0\\ \vec f_j(t')}
\label{eq:6}
\end{equation}
with the matrix-valued propagator $G(t,t')$ and the effective force $\vec f_j(t')$ on particle $j$, introduced and derived in Appendix~\ref{app:A}.

Instead of individual particles, we consider canonical ensembles of $N\gg1$ classical particles $j$ at positions $\vec q_j$ with velocities $\dot{\vec q}_j$ and introduce the tensor $\tens x = \vec x_j\otimes\vec e_j$ to bundle the phase-space coordinates of the entire ensemble, where $\vec 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
\begin{equation}
  \tens x(t) = \tens G(t,0)\tens x^\mathrm{(i)}+
    \int_0^t\D t'\,\tens G(t,t')\cvector{0\\ \tens f(t')} =:
    \tens x_0(t)+\tens y(t)
\label{eq:7}
\end{equation}
with $\tens G(t,t') := G(t,t')\otimes\id{N}$ and $\tens f(t') := \vec f_j(t')\otimes\vec e_j$. The inertial trajectories are $\tens x_0 = \tens G(t,0)\tens x^\mathrm{(i)}$, and $\tens y$ are the deviations therefrom caused by the effective force $\tens f$.

\subsection{Effective gravitational potential}
\label{sec:2.2}

The potential of the effective force acting relative to Zel'dovich trajectories is given by Eq.~(\ref{eq:73}) of Appendix~\ref{app:A},
\begin{equation}
  \phi = \varphi+A_\varphi D_+\psi\;,
\label{eq:8}
\end{equation}
where $\psi$ is the potential of the curl-free initial velocity field. With the Poisson equations (\ref{eq:5}) for $\varphi$ and $\vec\nabla^2\psi = -\delta^\mathrm{(i)}$ for the initial velocity potential $\psi$, the Poisson equation for $\phi$ is
\begin{equation}
  \vec\nabla^2\phi = A_\varphi\left(\delta-D_+\delta^\mathrm{(i)}\right) =
  A_\varphi\left(\delta-\delta^{\mathrm{(lin)}}\right)\;,
\label{eq:9}
\end{equation}
where $\delta^\mathrm{(lin)} = D_+\delta^\mathrm{(i)}$ is the linearly growing density contrast. The potential $\phi$ describing the gravitational interaction relative to the inertial particle trajectories $\tens x_0(t)$ is thus sourced exclusively by the non-linearly evolved contribution to the density fluctuations. The effective gravitational force mediated by $\phi$ 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 \emph{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.

\subsection{Shape of the effective potential}
\label{sec:2.3}

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\ll k_0$ and approaches the Newtonian form for $k\gg 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 $\phi = n\delta\ast v$ of the particle ensemble with the mean number density $n$, i.e.\ the convolution of the particle number-density fluctuation $n\delta$ with the potential $v$, satisfies the Poisson equation (\ref{eq:9}). Using the Fourier convolution theorem, the Fourier transform of this Poisson equation reads
\begin{equation}
  \tilde\delta\,\tilde v =
  -\frac{A_\varphi}{nk^2}\left(\tilde\delta-\tilde\delta^\mathrm{(lin)}\right)\;.
\label{eq:10}
\end{equation}
Multiplying this equation once with $\tilde\delta$, once with $\tilde\delta^\mathrm{(lin)}$, and taking the ensemble average gives
\begin{align}
  P_\delta\tilde v &\approx -\frac{A_\varphi}{nk^2}\left(
    P_\delta-\left<\tilde\delta\tilde\delta^\mathrm{(lin)}\right>
  \right)\quad\mbox{and}\nonumber\\
  \left<\tilde\delta\tilde\delta^\mathrm{(lin)}\right>\tilde v &\approx
  -\frac{A_\varphi}{nk^2}\left(
    \left<\tilde\delta\tilde\delta^\mathrm{(lin)}\right>-
    P_\delta^\mathrm{(lin)}
  \right)\;,
\label{eq:11}
\end{align}
introducing the power spectrum $P_\delta$ of the density contrast. Going from (\ref{eq:10}) to (\ref{eq:11}), we have implicitly assumed that the form of $\tilde v$ is independent of the density contrast, which should be a good approximation, but does not generally need to be the case. Eliminating $\langle\tilde\delta\tilde\delta^\mathrm{(lin)}\rangle$ between these equations leads to a quadratic equation for $\tilde v$ whose only meaningful solution is
\begin{equation}
  \tilde v = -\frac{A_\varphi}{nk^2}\,f_v(k)\quad\mbox{with}\quad
  f_v(k) = 1-\left(\frac{P_\delta^\mathrm{(lin)}}{P_\delta}\right)^{1/2}\;.
\label{eq:12}
\end{equation}
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 \cite{2003MNRAS.341.1311S, 2016MNRAS.459.1468M, 2015MNRAS.454.1958M} suggest that the function
\begin{equation}
  f_v(k) = \frac{k^2}{k_0^2+k^2}
\label{eq:13}
\end{equation}
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.\ \ref{fig:1}). With (\ref{eq:12}), this results in the effective potential
\begin{equation}
  \tilde v = -\frac{A_\varphi}{n\left(k_0^2+k^2\right)}\;,
\label{eq:14}
\end{equation}
which is of Yukawa rather than Newtonian form.

\begin{figure}[ht]
  \includegraphics[width=0.49\hsize]{fig_1a}\hfill
  \includegraphics[width=0.49\hsize]{fig_1b}
\caption{\emph{Left}: The ratio of power spectra $1-(P_\delta^\mathrm{(lin)}/P_\delta)^{1/2}$, scaled to unity at $k\gg1$, compared to the fitting function $f_v(k)$ from Eq.\ (\ref{eq:13}). The linear power spectrum was taken from \cite{1986ApJ...304...15B}, the non-linear from \cite{2016MNRAS.459.1468M}. \emph{Right}: The Yukawa scale $k_0$ as a function of the redshift $z$.}
\label{fig:1}
\end{figure}

Note that the expression (\ref{eq: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.~\ref{sec:4.3}. The right panel in Fig.~\ref{fig:1} shows the time-dependent Yukawa scale $k_0$ determined in this way.

\section{Power spectra from KFT}
\label{sec:3}

We briefly review in this section the KFT approach to cosmic structure formation. For further detail, we refer the reader to \cite{2016NJPh...18d3020B}, \cite{2017NJPh...19h3001B}, and the review \cite{2019AdP...531..446B}.

\subsection{Generating functional}
\label{sec:3.1}

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(\tens x)$ for the phase-space positions $\tens x$ to be occupied. Splitting $P(\tens x)$ into a probability $P(\tens x^\mathrm{(i)})$ for an initial state times a transition probability $P(\tens x|\tens x^\mathrm{(i)})$ from the initial to the final state, further introducing a generator field
\begin{equation}
  \tens J := \cvector{\vec J_{q_i}(t)\\ \vec J_{p_i}(t)}\otimes\vec e_i
\label{eq:15}
\end{equation}
to allow extracting moments of particle positions and momenta via functional derivatives with respect to $\tens J$ later, and introducing the particle trajectories $\tens x(t)$ from (\ref{eq:7}), leads to the generating functional
\begin{equation}
  Z[\tens J] = \int\D\Gamma\,\E^{\I\tens J\cdot\tens x}
\label{eq:16}
\end{equation}
derived in \cite{2017NJPh...19h3001B}, with the dot denoting the time-integrated scalar product
\begin{equation}
  \tens A\cdot\tens B = \int\D t\,\left<\vec A_i,\vec B_i\right>\;.
\label{eq:17}
\end{equation}
According to (\ref{eq:7}), the phase in (\ref{eq:16}) can be split into a free and an interacting part,
\begin{equation}
  Z[\tens J] = \int\D\Gamma\,
  \E^{\I\tens J\cdot(\tens x_0+\tens y)}\;.
\label{eq:18}
\end{equation}

Since the contribution $S_\mathrm{I}[\tens J] := \I\tens J\cdot\tens y$ to the phase depends on all relative positions of the correlated particle ensemble, it seems impossible to evaluate the generating functional $Z[\tens 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 (\ref{eq: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 \cite{2016NJPh...18d3020B}, 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,
\begin{equation}
  S_\mathrm{I}[\tens J] \to \left<S_\mathrm{I}[\tens J]\right>\;,
\label{eq:19}
\end{equation}
which would allow us to write the generating functional (\ref{eq:18}) as
\begin{equation}
  Z[\tens J] \approx \E^{\left<S_\mathrm{I}[\tens J]\right>}
  \int\D\Gamma\,\E^{\I\tens J\cdot\tens x_0} =:
  \E^{\left<S_\mathrm{I}[\tens J]\right>}Z_0[\tens J]\;.
\label{eq:20}
\end{equation}
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[\tens J]$.

\subsection{Density cumulants}
\label{sec:3.2}

The particle number density is a sum of delta distributions centered on the particle positions at time $t$ or, in a Fourier representation,
\begin{equation}
  \tilde\rho\left(\vec k, t\right) =
  \sum_{i=1}^N\E^{-\I\vec k\cdot\vec q_i(t)} =:
  \sum_{i=1}^N\tilde\rho_i\left(\vec k, t\right)\;.
\label{eq:21}
\end{equation}
Replacing the particle position $\vec q_i$ in each one-particle density contribution $\tilde\rho_i(\vec k, t)$ by a functional derivative with respect to $\vec J_{q_i}(t)$, we obtain the one- and $N$-particle density operators
\begin{equation}
  \hat\rho_i\left(\vec k, t\right) =
  \exp\left(-\I\vec k\cdot\frac{\delta}{\I\delta\vec J_{q_i}(t)}\right)\;,\quad
  \hat\rho\left(\vec k, t\right) =
  \sum_{i=1}^N\hat\rho_i\left(\vec k, t\right)\;.
\label{eq:22}
\end{equation}
Since the density operators are exponentials of derivatives with respect to components of $\tens J$, they generate translations of the generating functional, $Z[\tens J]\to Z[\tens J+\tens L]$.

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
\begin{equation}
  G_{\rho\rho}(1,2) =
  \sum_{i\ne j=1}^N\hat\rho_i(1)\hat\rho_j(2)\,Z[\tens J] =
  \sum_{i\ne j=1}^NZ[\tens J+\tens L]\;,
\label{eq:23}
\end{equation}
with the corresponding shift
\begin{equation}
  \tens L = -\delta_\mathrm{D}(t'-t)\cvector{1\\0}
  \left(\vec k_1\otimes\vec e_i+\vec k_2\otimes\vec e_j\right)\;.
\label{eq:24}
\end{equation}
The short-hand notation $(n)$ for the arguments in (\ref{eq:23}) indicates the Fourier-space position $(\vec k_n, t_n)$ at time $t_n$. The generator field $\tens 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 (\ref{eq:23}) gives the same result as for any particle pair arbitrarily labelled as $i, j = 1, 2$, and the cumulant becomes
\begin{equation}
  G_{\rho\rho}(1,2) = N(N-1)\,Z[\tens L] \approx N(N-1)\,
  \E^{\left<S_\mathrm{I}[\tens L]\right>}Z_0[\tens L]\;.
\label{eq:25}
\end{equation}

With initial conditions appropriate for the early universe, the free generating functional $Z_0[\tens J]$ after applying two density operators $\hat\rho_1(1)$ and $\hat\rho_2(2)$ can be written as
\begin{equation}
  Z_0[\tens L] =
  (2\pi)^3\delta_\mathrm{D}\left(\vec k_1+\vec k_2\right)
  V^{-2}\E^{-Q_\mathrm{D}}\,\mathcal{P}(k_1) =:
  (2\pi)^3\delta_\mathrm{D}\left(\vec k_1+\vec k_2\right)
  V^{-2}\bar{\mathcal{P}}(k_1)
\label{eq:26}
\end{equation}
with the damping term
\begin{equation}
  Q_\mathrm{D} := \frac{k^2}{3}\left(\sigma_1t\right)^2
\label{eq:27}
\end{equation}
and the freely evolved, non-linear power spectrum
\begin{equation}
  \mathcal{P}(k_1) := \int_q\left[
    \E^{-t^2k_1^2a_\parallel(q,\mu)}-1
  \right]\E^{\I\vec k_1\cdot\vec q}\;;
\label{eq:28}
\end{equation}
see \cite{2017NJPh...19h3001B} for the derivation. The function $a_\parallel$ appearing here is the auto-correlation function of momentum components parallel to the wave vector $\vec k_1$,
\begin{equation}
  a_\parallel =
  \mu^2\xi_\psi''(q)+\left(1-\mu^2\right)\frac{\xi_\psi'(q)}{q}\;,
\label{eq:29}
\end{equation}
where $\mu$ is the direction cosine of $\vec q$ relative to $\vec k_1$. The function
\begin{equation}
  \xi_\psi(q) = \frac{1}{2\pi^2}
  \int_0^\infty\frac{\D k}{k^2}P_\delta^\mathrm{(i)}(k)j_0(kq)
\label{eq:30}
\end{equation}
containing the spherical Bessel function $j_0$ is the auto-correlation function of the initial velocity potential $\psi$, determined by the initial density-fluctuation power spectrum $P_\delta^\mathrm{(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
\begin{equation}
  \sigma_n^2 := \int_kk^{2(n-2)}P_\delta^\mathrm{(i)}(k) =
  \frac{1}{2\pi^2}\int_0^\infty\D k\,k^{2n-2}P_\delta^\mathrm{(i)}(k)
\label{eq:31}
\end{equation}
of the density-fluctuation power spectrum and note that
\begin{equation}
  \lim_{q\to0}a_\parallel = -\frac{\sigma_1^2}{3}\;.
\label{eq:32}
\end{equation}
For small arguments of the first exponential in (\ref{eq:28}), $\mathcal{P}$ turns into the linearly evolved power spectrum,
\begin{equation}
  \mathcal{P}(k) \to (1+t)^2P_\delta^\mathrm{(i)}(k) = P_\delta^\mathrm{(lin)}(k)\;,
\label{eq:33}
\end{equation}
as shown in \cite{2017NJPh...19h3001B}.

\subsection{Damping and interaction}
\label{sec:3.3}

The damping term $Q_\mathrm{D}$ in (\ref{eq:26}) requires a careful discussion. Its definition (\ref{eq:27}) in conjunction with the velocity dispersion $\sigma_1^2$ from (\ref{eq:31}) shows that it arises because particles stream freely with an average velocity $\sigma_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 (\ref{eq: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 spectrum $\bar{\mathcal{P}}$ \emph{combined} with the damping term follows the linearly evolved power spectrum $P_\delta^\mathrm{(lin)}$ on large scales, drops below it on non-linear scales, but turns towards an asymptotic behaviour $\propto k^{-3}$ as $k$ increases further \cite{2017NJPh...19h3001B}.

Combining (\ref{eq:25}) and (\ref{eq:26}) to
\begin{equation}
  G_{\rho\rho}(1,2) \approx
  n^2\E^{\left<S_\mathrm{I}[\tens L]\right>-Q_\mathrm{D}}\,
  \mathcal{P}\left(k_1\right)
\label{eq:34}
\end{equation}
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 $\langle S_\mathrm{I}\rangle$ and $Q_\mathrm{D}$, leads to non-linear structure growth.

\section{Mean-field approach to non-linear power spectra}
\label{sec:4}

Based on these arguments, we now pursue the following approach. We wish to represent the particle interactions by a suitably averaged interaction term $\langle S_\mathrm{I}[\tens L]\rangle$. For simplicity, we further wish to approximate the damped, freely evolving power spectrum $\bar{\mathcal{P}}$ by the linearly evolving power spectrum $P_\delta^\mathrm{(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
\begin{equation}
  P^\mathrm{(nl)}_\delta(k) \approx \E^{\langle S_\mathrm{I}\rangle(k)}
  P^\mathrm{(lin)}_\delta(k)
\label{eq:35}
\end{equation}
for the non-linear power spectrum. Our main result will be a specific and simple equation for $\langle S_\mathrm{I}\rangle(k)$ reproducing numerically derived, non-linear power spectra remarkably well.

\subsection{Averaged particle-particle force}
\label{sec:4.1}

With $\tens y$ from (\ref{eq:7}) and $\tens L$ from (\ref{eq:24}), we have
\begin{equation}
  S_\mathrm{I}[\tens L] = \I\tens L\cdot\tens y =
  -\I\int_0^t\D t'\,(t-t')\left[
    \vec k_1\cdot\vec f_1(t')+\vec k_2\cdot\vec f_2(t')
  \right]
\label{eq:36}
\end{equation}
where $\vec 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 $\vec f_\mathrm{p}$ and the particle number density $\rho$, we can write the force $\vec f_i$ on particle $i$ as
\begin{equation}
  \vec f_i = \int_{q_1}\int_{q_2}
  \rho_i\left(\vec q_1\right)\vec f_\mathrm{p}\left(\vec q_1-\vec q_2\right)
  \rho\left(\vec q_2\right)\;.
\label{eq:37}
\end{equation}

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 $\langle\vec f_i\,\rangle$ on the wave vector $\vec k$, we project out the contribution by the mode $\vec k$ of the density field, writing
\begin{equation}
  \left<\vec f_i\,\right>\left(\vec k\,\right) = \int_{q_1}\int_{q_2}
  \vec f_\mathrm{p}\left(\vec q_1-\vec q_2\right)\left<
    \rho_i\left(\vec q_1\right)\rho\left(\vec q_2\right)
  \right>\E^{-\I\vec k\cdot(\vec q_1-\vec q_2)}\;.
\label{eq:38}
\end{equation}
The average over the product of densities introduces the correlation function $\xi(|\vec q_1-\vec q_2|)$ of the density field,
\begin{equation}
  \left<\rho_i\left(\vec q_1\right)\rho\left(\vec q_2\right)\right> =
  \frac{1}{N}\left<\rho\left(\vec q_1\right)\rho\left(\vec q_2\right)\right> =
  \frac{n^2}{N}\left[1+\xi\left(\left|\vec q_1-\vec q_2\right|\right)\right]\;.
\label{eq:39}
\end{equation}
We keep only the connected part of the correlation expressed by $\xi$ alone in (\ref{eq:39}) because the disconnected part cannot contribute in a homogeneous random field. Owing to homogeneity, the integrand in (\ref{eq:38}) depends only on the difference $\vec q_1-\vec q_2 =: \vec q$ of position vectors. We can thus integrate over $\vec q_1$, resulting in a factor $V$, and obtain
\begin{equation}
  \left<\vec f_i\,\right>\left(\vec k\,\right) = n\int_q
  \vec f_\mathrm{p}\left(\vec q\,\right)\xi(q)
  \E^{-\I\vec k\cdot\vec q}\;.
\label{eq:40}
\end{equation}
Since the remaining expression is the Fourier transform of a product, the scale-dependent, mean force term is the convolution of the particle-particle force $\tilde f_\mathrm{p}$ in Fourier space and the power spectrum $P_\delta$ of the particle distribution as the Fourier transform of the correlation function,
\begin{equation}
  \left<\vec f_i\,\right>\left(\vec k\,\right) =
  n\left(\tilde f_\mathrm{p}*P_\delta\right)\left(\vec k\,\right)\;.
\label{eq:41}
\end{equation}
Combining (\ref{eq:41}) with (\ref{eq:36}) and using Newton's third axiom in the form $\langle\vec f_i\,\rangle(-\vec k\,) = -\langle\vec f_i\,\rangle(\vec k\,)$, we thus find the expression
\begin{equation}
  \left<S_I\right>\left(\vec k\,\right) =
  -2\I n\vec k\cdot\int_0^t\D t'\left(t-t'\right)
  \left(\tilde f_\mathrm{p}*P_\delta\right)\left(\vec k\,\right)
\label{eq:42}
\end{equation}
for the averaged, scale-dependent, interaction term.

\subsection{Damping within Burgers' approximation}
\label{sec:4.2}

We now need to specify the power spectrum $P_\delta$ to be inserted into (\ref{eq:41}) for evaluating the average force term $\langle\vec f_i\,\rangle$. Following (\ref{eq:38}), evaluating the mean interaction term with the density correlation function suggests replacing $P_\delta = \bar{\mathcal{P}}$. Applying the same linear approximation leading from (\ref{eq:34}) to (\ref{eq:35}), we would then arrive at $P_\delta \approx P_\delta^\mathrm{(lin)}$. However, since the linear power spectrum in (\ref{eq: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 \cite{1940PKNAW..43....2B, 1972QAMat..30..195B, 1974nlde.book.....B}.

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,
\begin{equation}
  \frac{\D\dot{\vec q}}{\D t} = 0 \quad\to\quad
  \frac{\D\dot{\vec q}}{\D t} = \nu\vec\nabla^2\dot{\vec q}\;,
\label{eq:43}
\end{equation}
where $\nu$ is a viscosity parameter with the dimension of a squared length. A natural choice for the length scale $\nu^{1/2}$ is the non-linear radius $r_\mathrm{nl}$ defined by
\begin{equation}
  \sigma_{r_\mathrm{nl}}^2 = \left.
    \int_kP_\delta(k)W_R^2(k)
  \right\vert_{R = r_\mathrm{nl}} = 1
\label{eq:44}
\end{equation}
evaluated with the linear power spectrum from \cite{1986ApJ...304...15B} and a top-hat window function $W_R(k)$ at the present cosmic time. We set $\nu = r_\mathrm{nl}^2 = \mathrm{const}$ here for simplicity, but note that $\nu$ could also be generalised to become time-dependent.

Burgers' equation can be solved by a Hopf-Cole transformation \cite{1950CPApM...3..201H, 1951QAMat...9..225C}, which results in
\begin{equation}
  \dot{\vec q} = -2\nu\vec\nabla\ln U\;,
\label{eq:45}
\end{equation}
where $U$ is an exponential velocity potential given by the convolution
\begin{equation}
  U = \mathcal{N}_{\sqrt{2\nu t}}*\exp\left(-\frac{\psi}{2\nu}\right)
\label{eq:46}
\end{equation}
of the scaled, exponentiated initial velocity potential $\psi$ with a normal distribution $\mathcal{N}$ of width $(2\nu t)^{1/2}$. To linear order in $\psi$, (\ref{eq:45}) implies the velocity dispersion
\begin{equation}
  \sigma_v^2(t) = \int_k\frac{P_\delta(k)}{k^2}\exp\left(-2k^2\nu t\right)\;.
\label{eq:47}
\end{equation}
This expression clearly shows the effect of Burgers' approximation: small-scale modes with $k\gtrsim(2\nu 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 \cite{1989MNRAS.236..385G, 1990MNRAS.247..260W, 1992ApJ...393..437K, 2015JCAP...01..014R}. We take this reduced amount of particle motion into account by replacing the damping term $Q_D$ from (\ref{eq:27}) by
\begin{equation}
  \bar Q_D = k^2\lambda^2 \quad\mbox{with the damping scale}\quad
  \lambda(t) = \int_0^t\D t'\sigma_v(t')\;.
\label{eq:48}
\end{equation}
An excellent fit to $\lambda(t)$ is
\begin{equation}
  \lambda(t) \approx \frac{t}{\sqrt{1+t/\tau}}\quad\mbox{with}\quad
  \tau \approx 24.17\;.
\label{eq:49}
\end{equation}
The form of this fit expresses the transition from ballistic to diffusive particle motion.

\subsection{Averaged interaction term}
\label{sec:4.3}

Accordingly, we evaluate the averaged force term (\ref{eq:41}) as
\begin{equation}
  \left<\vec f_i\,\right>\left(\vec k\,\right) =
  n\left(\tilde f_\mathrm{p}*\bar P_\delta\right)\left(k\right)
  \quad\mbox{with}\quad
  \bar P_\delta(k) = \left(1+\bar Q_\mathrm{D}\right)^{-1}
  P_\delta^\mathrm{(lin)}(k)\;.
\label{eq:50}
\end{equation}
Inserting the averaged force term (\ref{eq:50}) into (\ref{eq:42}) results in the averaged interaction term
\begin{equation}
  \left<S_\mathrm{I}\right>(k) = -2\I n\vec k\cdot
  \int_0^t\D t'\left(t-t'\right)\left(
    \tilde f_\mathrm{p}*\bar P_\delta
  \right)(k)\;.
\label{eq:51}
\end{equation}
We finally need to evaluate the convolution of the particle-particle force term $\tilde f_\mathrm{p}$ with the damped power spectrum $\bar P$. This is done in Appendix~\ref{app:B} and results in the average interaction term
\begin{equation}
  \left<S_\mathrm{I}\right>(k) = 2\int_0^t\D t'\,(t-t')\frac{\dot m}{m}
  \left[
    D_+\sigma_J^2-
    \frac{1}{m}\int_0^{t'}\D\bar t\,
    \dot mD_+\sigma_J^2
  \right]\;,
\label{eq:52}
\end{equation}
where $\sigma_J^2$ is the moment (\ref{eq:83}) of the damped power spectrum. This interaction term $\langle S_\mathrm{I}\rangle$ in the mean-field approximation is shown in the right panel of Fig.\,\ref{fig:B} 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 (\ref{eq:14}) still needs to be set. KFT itself now suggests the following procedure. According to (\ref{eq:35}), the ratio between the linearly and non-linearly evolved power spectra is estimated by the exponential of the averaged interaction term. Combining (\ref{eq:35}) with (\ref{eq:12}), the factor $f_v$ from (\ref{eq:12}) is approximated by
\begin{equation}
  f_v(k) \approx 1-\E^{-\langle S_\mathrm{I}\rangle/2}\;.
\label{eq:53}
\end{equation}
We can thus determine a first estimate for $k_0$ by calculating $\langle S_\mathrm{I}\rangle$ with $k_0 = 0$ and fitting $f_v(k)$ from (\ref{eq:53}) with the functional form (\ref{eq:13}). With the resulting value of $k_0$, an updated estimate for $\langle S_\mathrm{I}\rangle$ 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.~\ref{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.

\subsection{Non-linear density-fluctuation power spectrum}
\label{sec:4.4}

The mean-field averaged interaction term (\ref{eq:52}), inserted into (\ref{eq:35}), is our approximate expression for the non-linear density-fluctuation power spectrum. Evaluating the averaged interaction term $\langle S_\mathrm{I}\rangle(k)$ from (\ref{eq:52}) numerically, determining the viscosity parameter $\nu$ and the Yukawa scale $k_0$ from KFT itself as described, and further assuming a spatially-flat $\Lambda$CDM model universe with matter-density parameter $\Omega_\mathrm{m0} = 0.3$, leads to the result shown in Fig.\,\ref{fig:2}.

\begin{figure}[ht]
  \includegraphics[width=\hsize]{fig_3}
\caption{Analytic and numerical power spectra at redshift $z = 0$. The linear power spectrum is taken from \cite{1986ApJ...304...15B}, the numerical spectrum from \cite{2016MNRAS.459.1468M}. The flat lower panel shows the relative deviation between the analytic and the numerical power spectra.}
\label{fig:2}
\end{figure}

In this Figure, our analytic approximation (\ref{eq:35}) is compared at redshift $z = 0$ to the description by \cite{2016MNRAS.459.1468M} of the power spectrum obtained from numerical simulations. As can be seen there, the analytic power spectrum agrees within typically $\lesssim10\%$ with the numerical expectation up to wave numbers of $k\lesssim 10\,h\,\mathrm{Mpc}^{-1}$. With the parameters $\nu$ and $k_0$ self-calibrated from within KFT, the analytic expression (\ref{eq: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 \cite{1986ApJ...304...15B}. Expression (\ref{eq: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 $\nu$ 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 (\ref{eq:35}) by minimising the squared difference between them. Doing so, using the power spectrum from \cite{2016MNRAS.459.1468M} as a reference, results in the power spectrum shown in Fig.~\ref{fig:3}.

\begin{figure}[ht]
  \includegraphics[width=\hsize]{fig_4}
\caption{Like Fig.~\ref{fig:2}, but with the parameters $\nu$ and $k_0$ modified to optimise the agreement with the numerical results by \cite{2016MNRAS.459.1468M}.}
\label{fig:3}
\end{figure}

The power spectrum shown there is obtained by reducing the Yukawa scale by $13\%$ and increasing the displacement (\ref{eq:48}) by $6\%$. The relative deviation of the mean-field approximated, analytic power spectrum from its numerical counterpart is now lowered to typically $\lesssim 5\,\%$ up to $k\lesssim10\,h\,\mathrm{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.

\section{Summary and conclusion}
\label{sec:5}

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 (\ref{eq:35}) for the non-linear power spectrum $P_\delta^\mathrm{(nl)}$ and the expression (\ref{eq: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 (\ref{eq:38}) and (\ref{eq: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\lesssim10\,h\,\mathrm{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 (\ref{eq: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 $\nu$, 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 (\ref{eq:44}).

We thus have two parameters, $k_0$ and $\nu$, 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, non-linear power spectrum shown in Fig.~\ref{fig:2}, which already agrees well with numerical results. This agreement can further be improved by slightly adjusting both parameters, as shown in Fig.~\ref{fig:3}.

The initial state of the microscopic degrees of freedom is fully determined by the linear density-fluctuation 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 dark-matter models, a different cosmological background, or alternative gravity theories. Non-linear cosmic power spectra such as these shown in Figs.~\ref{fig:2} and \ref{fig: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.

\section*{Acknowledgements}

We gratefully acknowledge fruitful discussions with many colleagues, most notably the always very helpful discussions with Manfred Salmhofer.

\paragraph{Funding information}

This work was supported in part by Deutsche Forschungsgemeinschaft (DFG) under Germany's Excellence Strategy EXC-2181/1 - 390900948 (the Heidelberg STRUCTURES Excellence Cluster), by the Heidelberg Graduate School of Physics (HGSFP), by the Centre for Quantum Dynamics at Heidelberg University, and by a Technion Fellowship.

\begin{appendix}

\section{Particle Trajectories}
\label{app:A}

\subsection{Particle mass}

Beginning with the effective particle mass $m$ defined in (\ref{eq:4}), the time derivative of $m$ is
\begin{equation}
  \dot m = \frac{\D m}{\D D_+} = \frac{m'}{D_+'}
\label{eq:54}
\end{equation}
by definition of the time $t = D_+-1$. The prime denotes differentiation with respect to the scale factor $a$. We can insert (\ref{eq:4}) once more into (\ref{eq:54}) to arrive at
\begin{equation}
  \dot m = \frac{a^3E}{D_+'}\left[
    D_+''+\left(\frac{3}{a}+\frac{E'}{E}\right)D_+'
  \right]\;.
\label{eq:55}
\end{equation}
The growth factor $D_+$ solves the linear growth equation. When transformed to the scale factor $a$ as an independent variable, this reads
\begin{equation}
  D_+''+\left(\frac{3}{a}+\frac{E'}{E}\right)D_+' =
  \frac{3}{2}\frac{\Omega_\mathrm{m}}{a^2}D_+\;.
\label{eq:56}
\end{equation}
Expressing the matter-density parameter $\Omega_\mathrm{m}$ by its value $\Omega_\mathrm{m}^\mathrm{(i)}$ at the initial time,
\begin{equation}
  \Omega_\mathrm{m} = \frac{\Omega_\mathrm{m}^\mathrm{(i)}}{a^3E^2}\;.
\label{eq:57}
\end{equation}
If we specify the initial conditions early in the matter-dominated epoch, we may further approximate $\Omega_\mathrm{m}^\mathrm{(i)} \approx 1$. We thus find the expression
\begin{equation}
  \dot m = \frac{3}{2}\,\Omega_\mathrm{m}^\mathrm{(i)}\,\frac{aD_+}{m}
\label{eq:58}
\end{equation}
for the time derivative of the effective particle mass $m$. Inserting finally the potential amplitude $A_\varphi$ from the Poisson equation (\ref{eq:5}), we arrive at the time derivative
\begin{equation}
  \dot m = mA_\varphi D_+ = mA_\varphi(t+1)
\label{eq:59}
\end{equation}
of the effective particle mass $m$.

\subsection{Solution of the Hamiltonian equations of motion}

The Hamiltonian equations of motion for the phase-space point $\vec x = (\vec q, \vec p\,)^\top$ can be written in the form
\begin{equation}
  \dot{\vec x} = A(t)\vec x-\cvector{0\\ m\vec\nabla\varphi}
  \quad\mbox{with}\quad
  A(t) := \matrix{cc}{0_3 & m^{-1}\id{3} \\ 0_3 & 0_3}\;.
\label{eq:60}
\end{equation}
Notice that $A(t)$ is a $6\times6$ matrix, with $0_3$ and $\id{3}$ representing the zero and unit matrices in three dimensions, respectively. The homogeneous equation $\dot{\vec x} = A(t)\vec x$ is solved by $\vec x_\mathrm{h} = \exp[\bar A(t,0)]\,\vec x_0$, with
\begin{equation}
  \bar A(t,t') := \int_{t'}^t\D\bar t\,A(\bar t) =
  \matrix{cc}{0_3 & g_\mathrm{H}(t,t')\,\id{3} \\ 0_3 & 0_3}\;,\quad
  g_\mathrm{H}(t,t') := \int_{t'}^t\frac{\D\bar t}{m(\bar t)}\;.
\label{eq:61}
\end{equation}
Since $\bar A$ is nilpotent, $\bar A^2(t,t') = 0_6$, the homogeneous solution shrinks to
\begin{equation}
  \vec x_\mathrm{h}(t) = \left(1+\bar A(t,0)\right)\vec x_0\;.
\label{eq:62}
\end{equation}
By variation of the constant vector $\vec x_0$, the inhomogeneous equation of motion (\ref{eq:60}) leads to
\begin{equation}
  \vec x_0(t) = \vec x^{\,\mathrm{(i)}}-
  \int_0^t\D t'\left(1+\bar A(t',0)\right)\cvector{0\\ m\vec\nabla\varphi}
\label{eq:63}
\end{equation}
and thus to the solution
\begin{equation}
  \vec x(t) = \left(1+\bar A(t,0)\right)\vec x^{\,\mathrm{(i)}}-
  \int_0^t\D t'\left(1+\bar A(t,t')\right)\cvector{0\\ m\vec\nabla\varphi}
\label{eq:64}
\end{equation}
for phase-space trajectories beginning at $\vec x^{\,\mathrm{(i)}} = (\vec q^{\,\mathrm{(i)}},\vec p^{\,\mathrm{(i)}})^\top$ at $t = 0$. The spatial trajectories are accordingly
\begin{equation}
  \vec q(t) =
  \underbrace
    {\vec q^{\,\mathrm{(i)}}+g_\mathrm{H}(t,0)\vec p^{\,\mathrm{(i)}}}_
    {=:\,\vec q_0(t)}-
  \int_0^t\D t'\,g_\mathrm{H}(t,t')m\vec\nabla\varphi =:
  \vec q_0(t)+\vec y_q(t)\;.
\label{eq:65}
\end{equation}

\subsection{Reference Trajectories}

It is often convenient in cosmology to replace the force-free trajectories $\vec q_0(t) = \vec q^{\,\mathrm{(i)}}+g_\mathrm{H}(t,0)\vec p^{\,\mathrm{(i)}}$ by the trajectories postulated in the Zel'dovich approximation,
\begin{equation}
  \vec q_0(t) \to
  \vec q^{\,\mathrm{(i)}}+t\vec p^{\,\mathrm{(i)}}\;,
\label{eq:66}
\end{equation}
exchanging the Hamiltonian propagator $g_\mathrm{H}(t,t')$ for $(t-t')$. The deviation $\vec y_q(t)$ defined in (\ref{eq:65}) as the difference between the actual and these reference trajectories is then determined by
\begin{equation}
  \vec y_q(t) = -\int_0^t\D t'\,\left[
    g_\mathrm{H}(t,t')m\vec\nabla\varphi+\left(
      1-\dot g_\mathrm{H}(t',0)
    \right)\vec p^{\,\mathrm{(i)}}
  \right]\;.
\label{eq:67}
\end{equation}
Implicitly defining an amplitude $A_\mathrm{p}(t')$ by
\begin{equation}
  \int_0^t\D t'\,\left(
      1-\dot g_\mathrm{H}(t',0)
  \right) \stackrel{!}{=} \int_0^t\D t'\,g_\mathrm{H}(t,t')A_\mathrm{p}(t')\;,
\label{eq:68}
\end{equation}
we can write (\ref{eq:67}) in the form
\begin{equation}
  \vec y_q(t) = -\int_0^t\D t'\,g_\mathrm{H}(t,t')
  \left(m\vec\nabla\varphi+A_\mathrm{p}(t')\vec p^\mathrm{\,(i)}\right)\;.
\label{eq:69}
\end{equation}
Since the initial momentum is the gradient of a potential, $\vec p^{\,\mathrm{(i)}} = \vec\nabla\psi$, (\ref{eq:69}) suggests defining an effective potential
\begin{equation}
  \phi := \varphi+\frac{A_\mathrm{p}}{m}\psi
\label{eq:70}
\end{equation}
such that
\begin{equation}
  \vec y_q(t) = -\int_0^t\D t'\,g_\mathrm{H}(t,t')m\vec\nabla\phi\;.
\label{eq:71}
\end{equation}

Differentiating (\ref{eq:68}) twice with respect to the time $t$ and using $\dot g_\mathrm{H}(t,t') = m^{-1}(t)$ gives
\begin{equation}
  A_\mathrm{p}(t) = \dot m = mA_\varphi D_+
\label{eq:72}
\end{equation}
with (\ref{eq:59}), and thus the trajectories
\begin{equation}
  \vec q(t) = \vec q^{\,\mathrm{(i)}}+t\vec p^{\,\mathrm{(i)}}-
  \int_0^t\D t'\,g_\mathrm{H}(t,t')\,m\vec\nabla\phi\;,\quad
  \phi = \varphi+A_\varphi D_+\psi\;.
\label{eq:73}
\end{equation}
Continutity demands that the initial velocity potential satisfies the Poisson equation $\vec\nabla^2\psi = -\delta^\mathrm{(i)}$.

\subsection{Unifying Propagators}

It is often convenient to replace the propagator $g_\mathrm{H}(t,t')$ in (\ref{eq:73}) also by the time difference $t-t'$. For doing so, we implicitly introduce an effective force $\vec f$ demanding
\begin{equation}
  \int_0^t\D t'\,g_\mathrm{H}(t,t')m\vec\nabla\phi \stackrel{!}{=}
  -\int_0^t\D t'\,\left(t-t'\right)\,\vec f(t')\;,
\label{eq:74}
\end{equation}
and solve for $\vec f(t)$. Differentiating (\ref{eq:74}) twice with respect to $t$, we find
\begin{equation}
  \vec f(t) = -\vec\nabla\phi+
  \frac{\dot m}{m^2}\int_0^t\D t'\,m\vec\nabla\phi \;.
\label{eq:75}
\end{equation}
With this expression for $\vec f(t)$, we can write the spatial trajectories as
\begin{equation}
  \vec q(t) = \vec q^{\,\mathrm{(i)}}+t\vec p^{\,\mathrm{(i)}}+
  \int_0^t\D t'\,\left(t-t'\right)\vec f(t')\;.
\label{eq:76}
\end{equation}
Finally replacing $\vec x = (\vec q,\vec p\,)^\top$ by $(\vec q,\dot{\vec q}\,)^\top$, we can bring the solution $\vec x(t)$ of the equations of motion into the form
\begin{equation}
  \vec x(t) = G(t,0)\vec x^{\,\mathrm{(i)}}+
  \int_0^t\D t'\,G(t,t')\cvector{0\\ \vec f(t')}
\label{eq:77}
\end{equation}
with the $6\times6$ matrix
\begin{equation}
  G(t,t') = \matrix{cc}{\id{3} & (t-t')\,\id{3} \\ 0_3 & \id{3}}\;.
\label{eq:78}
\end{equation}

\section{Convolving the particle-particle force with the power spectrum}
\label{app:B}

Based on the result (\ref{eq:75}) for the effective force, we begin with the expression
\begin{equation}
  \vec f_\mathrm{p} = -\vec\nabla v+
  \frac{\dot m}{m^2}\int_0^t\D t'\,m\vec\nabla v
\label{eq:79}
\end{equation}
for the particle-particle force represented by the potential $v$. The averaged interaction term (\ref{eq:51}) requires projecting the Fourier-transformed particle-particle force convolved with the power spectrum onto the wave vector $\vec k$. With this in mind, we first convolve the Fourier-transformed potential gradient $\widetilde{\nabla v} = \I\vec k\tilde v$ with the damped power spectrum. Inserting $\tilde v$ from (\ref{eq:14}) then leads to the intermediate equation
\begin{equation}
  \vec k\cdot\left(\widetilde{\nabla v}\ast\bar P_\delta\right)(k) =
  -\frac{\I A_\varphi}{n}\int_{k'}
    \frac
      {\vec k\cdot\left(\vec k-\vec k'\right)}
      {k_0^2+\left(\vec k-\vec k'\right)^2}
    \bar P_\delta\left(k'\right)\;.
\label{eq:80}
\end{equation}
We introduce spherical polar coordinates to evaluate the integral in (\ref{eq:80}) and turn the coordinate system such that $\vec k$ points into the direction of the polar axis. Defining the cosine $\mu$ of the polar angle between $\vec k$ and $\vec k'$, further substituting $\kappa := k'/k$ and $\kappa_0 := k_0/k$ then leads to
\begin{equation}
  \vec k\cdot\left(\widetilde{\nabla v}\ast\bar P_\delta\right)(k) =
  -\frac{\I A_\varphi}{n}\frac{k^3}{(2\pi)^2}
  \int_0^\infty\D\kappa\,\kappa^2\,\bar P_\delta\left(k\kappa, t'\right)\,
  J(\kappa, \kappa_0)
\label{eq:81}
\end{equation}
with
\begin{equation}
  J(\kappa, \kappa_0) := \int_{-1}^1\D\mu\,
  \frac{1-\kappa\mu}{1+\kappa_0^2+\kappa^2-2k\kappa\mu} =
  1+\frac{1-\kappa^2-\kappa_0^2}{4\kappa}
  \ln\frac{\kappa_0^2+(1+\kappa)^2}{\kappa_0^2+(1-\kappa)^2}\;.
\label{eq:82}
\end{equation}

\begin{figure}[ht]
  \includegraphics[width=0.49\hsize]{fig_Ba}
  \includegraphics[width=0.49\hsize]{fig_Bb}
\caption{\emph{Left panel}: The function $J(\kappa, \kappa_0)$ appearing in the averaged force term $\langle\tilde f_{12}\rangle$ is shown here for five values of $\kappa_0$. While the asymptotic behaviour of $J(\kappa,\kappa_0)\propto\kappa^{-2}$ for $\kappa\gg1$ is unaffected by $\kappa_0$, $J(\kappa, \kappa_0)$ approaches the value $2/(1+\kappa_0^2)$ for $\kappa\ll1$. \emph{Right panel}: Average interaction term $\langle S_\mathrm{I}\rangle(k)$ as a function of the wave number $k$.}
\label{fig:B}
\end{figure}

The left panel of Fig.~\ref{fig:B} shows the function $J(\kappa, \kappa_0)$ for different values of $\kappa_0$. It falls off $\propto\kappa^{-2}$ for $\kappa\gg1$ and tends towards the constant $2/(1+\kappa_0^2)$ for $\kappa\ll1$. The function $J(\kappa, \kappa_0)$ can thus be seen as a filter function for the power spectrum. We introduce the moment
\begin{equation}
  \sigma^2_J := \frac{k^3}{(2\pi)^2}
  \int_0^\infty\D\kappa\,\kappa^2\,\bar P_\delta^\mathrm{(i)}(k\kappa)\,
  J(\kappa, \kappa_0)
\label{eq:83}
\end{equation}
of the damped \emph{initial} power spectrum $\bar P_\delta^\mathrm{(i)}$, filtered with the function $J$, and bring (\ref{eq:81}) into the form
\begin{equation}
  \vec k\cdot\left(\widetilde{\nabla v}\ast\bar P\right)\left(\vec k\,\right) =
  -\frac{\I A_\varphi}{n}\,D_+^2\sigma_J^2 =
  -\frac{\I}{n}\frac{\dot m}{m}D_+\sigma_J^2\;.
\label{eq:84}
\end{equation}
In the last step, we have used the time derivative $\dot m$ from (\ref{eq:59}) in Appendix~\ref{app:A}. According to (\ref{eq:79}), and using (\ref{eq:84}), the convolved force $\tilde f_\mathrm{p}\ast\bar P$ projected on the wave vector $\vec k$ is
\begin{equation}
  \vec k\cdot\left(\tilde f_\mathrm{p}\ast\bar P\right) =
  \frac{\I}{n}\frac{\dot m}{m}\left[
    D_+\sigma_J^2-\frac{1}{m}\int_0^t\D t'\,\dot mD_+\sigma_J^2
  \right]\;.
\label{eq:85}
\end{equation}
With this result, we return to the averaged interaction term (\ref{eq:51}), finding the averaged interaction term (\ref{eq:52}).

\end{appendix}