Kinetic field theory: effects of momentum correlations on the cosmic density-fluctuation power spectrum

The central object of the theory is the free generating functional ${Z}_{0}[{\boldsymbol{J}},{\boldsymbol{K}}]$ with generator fields ${\boldsymbol{J}}$ and ${\boldsymbol{K}}$ coupled to the phase-space coordinates ${\boldsymbol{x}}=({\boldsymbol{q}},{\boldsymbol{p}})$ and the one-particle response fields, respectively. The bold-faced symbols denote tensors bundling contributions from all N particles in the ensemble. Let ${\vec{x}}_{j}=({\vec{q}}_{j},{\vec{p}}_{j})$ be the phase-space coordinates of particle j, and ${\vec{e}}_{j}$ an N-dimensional unit vector with components ${({\vec{e}}_{j})}_{k}={\delta }_{{jk}}$, then

$$\begin{eqnarray}&&{\boldsymbol{x}}:= {\vec{x}}_{j}\otimes {\vec{e}}_{j}.\end{eqnarray} \tag{ 1 }$$

We define a scalar product between two such tensors by

$$\begin{eqnarray}&&\langle {\boldsymbol{a}},{\boldsymbol{b}}\rangle =\langle {\vec{a}}_{j}\otimes {\vec{e}}_{j},{\vec{b}}_{k}\otimes {\vec{e}}_{k}\rangle ={\vec{a}}_{j}\cdot {\vec{b}}_{j},\end{eqnarray} \tag{ 2 }$$

with Einstein's summation convention implied.

This generating functional is the integral

$$\begin{eqnarray}&&{Z}_{0}[{\boldsymbol{J}},{\boldsymbol{K}}]=\int {\rm{d}}{\rm{\Gamma }}\,\exp \left({\rm{i}}{\int }_{0}^{\infty }{\rm{d}}t\langle {\boldsymbol{J}},\bar{{\boldsymbol{x}}}\rangle \right)\end{eqnarray} \tag{ 3 }$$

over the N-particle phase space at the initial time ${t}_{0}=0$. In terms of the Green's function ${ \mathcal G }$ of the free equations of motion, the particle trajectories $\bar{{\boldsymbol{x}}}$ in phase space are

$$\begin{eqnarray}&&\bar{{\boldsymbol{x}}}(t)={ \mathcal G }(t,0){{\boldsymbol{x}}}^{({\rm{i}})}-{\int }_{0}^{t}{\rm{d}}t^{\prime} \,{ \mathcal G }(t,t^{\prime} ){\boldsymbol{K}}(t^{\prime} ).\end{eqnarray} \tag{ 4 }$$

The phase-space measure in (3) is

$$\begin{eqnarray}&&{\rm{d}}{\rm{\Gamma }}=P({{\boldsymbol{q}}}^{({\rm{i}})},{{\boldsymbol{p}}}^{({\rm{i}})}){\rm{d}}{{\boldsymbol{q}}}^{({\rm{i}})}{\rm{d}}{{\boldsymbol{p}}}^{({\rm{i}})}\end{eqnarray} \tag{ 5 }$$

with a suitable probability distribution P to be adapted to the initial conditions at hand. For N particles initially correlated in phase space, the probability distribution $P({{\boldsymbol{q}}}^{({\rm{i}})},{{\boldsymbol{p}}}^{({\rm{i}})})$ will be given in (21) below.

More explicitly, the N-particle Green's function ${ \mathcal G }$ is the tensor product

$$\begin{eqnarray}&&{ \mathcal G }(t,t^{\prime} )=G(t,t^{\prime} )\otimes {{ \mathcal I }}_{N}\end{eqnarray} \tag{ 6 }$$

of the matrix-valued, one-particle Green's function

$$\begin{eqnarray}G(t,t^{\prime} )=\left(\begin{array}{cc}{g}_{{qq}}(t,t^{\prime} ) & {g}_{{qp}}(t,t^{\prime} )\\ 0 & {g}_{{pp}}(t,t^{\prime} )\end{array}\right)\end{eqnarray} \tag{ 7 }$$

and the unit matrix ${{ \mathcal I }}_{N}$ in N dimensions.

Particle interactions are included by applying an interaction operator to the free generating functional, producing the full generating functional

$$\begin{eqnarray}&&Z[{\boldsymbol{J}},{\boldsymbol{K}}]=\exp ({\rm{i}}{\hat{S}}_{{\rm{I}}}){Z}_{0}[{\boldsymbol{J}},{\boldsymbol{K}}],\end{eqnarray} \tag{ 8 }$$

with

$$\begin{eqnarray}&&{\hat{S}}_{{\rm{I}}}=-\int {\rm{d}}1\,\hat{B}(-1)v(1)\hat{\rho }(1).\end{eqnarray} \tag{ 9 }$$

Here, v is the interaction potential, and $\hat{\rho }$ and $\hat{B}$ are the density- and response-field operators, respectively. The many-particle response field B describes how a particle ensemble responds to a change in the phase-space coordinates of one of its particles. The arguments abbreviate $1:= ({t}_{1},{\vec{k}}_{1})$ and $-1:= ({t}_{1},-{\vec{k}}_{1})$. In a Fourier-space representation, these operators are sums over one-particle operators

$$\begin{eqnarray}{\hat{\rho }}_{j}(1) & = & \exp \left(-{\rm{i}}{\vec{k}}_{1}\cdot \displaystyle \frac{\delta }{{\rm{i}}\delta {\vec{J}}_{{q}_{j}}(1)}\right),\\ {\hat{B}}_{j}(1) & = & \left({\rm{i}}{\vec{k}}_{1}\cdot \displaystyle \frac{\delta }{{\rm{i}}\delta {\vec{K}}_{{p}_{j}}(1)}\right){\hat{\rho }}_{j}(1)=: {\hat{b}}_{j}(1){\hat{\rho }}_{j}(1)\end{eqnarray} \tag{ 10 }$$

and the integral in (9) is taken over $1:= ({t}_{1},{\vec{k}}_{1})$. The response-field operator $\hat{B}$ thus contains a density operator $\hat{\rho }$.

Correlators of order n in the density, say, are obtained by applying n density operators to $Z[{\boldsymbol{J}},{\boldsymbol{K}}]$,

$$\begin{eqnarray}&&{G}_{\rho \ldots \rho }(1\ldots n)=\hat{\rho }(1)\cdots \hat{\rho }(n)Z[{\boldsymbol{J}},{\boldsymbol{K}}].\end{eqnarray} \tag{ 11 }$$

As usual in statistical field theory, each of the generator fields ${\boldsymbol{J}}$ and ${\boldsymbol{K}}$ is set to zero once all functional derivatives with respect to ${\boldsymbol{J}}$ or ${\boldsymbol{K}}$ have been applied.

An approach to a perturbative evaluation of (11) begins with expanding the exponential interaction operator $\exp ({\rm{i}}{\hat{S}}_{{\rm{I}}})$ into a power series, introducing two density operators $\hat{\rho }$ and one response-field operator $\hat{b}$ per power of ${\hat{S}}_{{\rm{I}}}$. Thus, for an nth order correlator with mth order particle interaction, we need to evaluate one-particle expressions of the form

$$\begin{eqnarray}&&{\hat{\rho }}_{{j}_{1}}(1)\cdots {\hat{\rho }}_{{j}_{r}}(r){{Z}_{0}[{\boldsymbol{J}},{\boldsymbol{K}}]| }_{{\boldsymbol{J}}=0}={Z}_{0}[{\boldsymbol{L}},{\boldsymbol{K}}]\end{eqnarray} \tag{ 12 }$$

with $r=n+2m$, and with the particle indices ${j}_{s}=1\ldots N$. The density operators thus replace the generator field ${\boldsymbol{J}}$ by the shift tensor

$$\begin{eqnarray}&&{\boldsymbol{L}}:= -\displaystyle \sum _{s=1}^{r}\left(\begin{array}{c}{\vec{k}}_{s}\\ 0\end{array}\right){\delta }_{{\rm{D}}}(t-{t}_{s})\otimes {\vec{e}}_{{j}_{s}}.\end{eqnarray} \tag{ 13 }$$

We define the position and momentum components of the shift tensor ${\boldsymbol{L}}$ out by the projections

$$\begin{eqnarray}{\vec{L}}_{{q}_{j}}(c) & := & {\displaystyle \int }_{0}^{\infty }{\rm{d}}t\left\langle {\boldsymbol{L}}(t),G(t,{t}_{c})\left(\begin{array}{c}{{ \mathcal I }}_{3}\\ 0\end{array}\right)\,\otimes \,{\vec{e}}_{j}\right\rangle ,\\ {\vec{L}}_{{p}_{j}}(c) & := & {\displaystyle \int }_{0}^{\infty }{\rm{d}}t\left\langle {\boldsymbol{L}}(t),G(t,{t}_{c})\left(\begin{array}{c}0\\ {{ \mathcal I }}_{3}\end{array}\right)\,\otimes \,{\vec{e}}_{j}\right\rangle \end{eqnarray} \tag{ 14 }$$

and abbreviate

$$\begin{eqnarray}&&{\vec{L}}_{{q}_{j}}:= {\vec{L}}_{{q}_{j}}(0),\quad {\vec{L}}_{{p}_{j}}:= {\vec{L}}_{{p}_{j}}(0)\end{eqnarray} \tag{ 15 }$$

with ${t}_{0}=0$. The symbol ${{ \mathcal I }}_{3}$ denotes the unit matrix in three-dimensions.

The subsequent application of a single, one-particle response-field operator ${\hat{b}}_{{j}_{c}}(c)$ to ${Z}_{0}[{\boldsymbol{L}},{\boldsymbol{K}}]$, taken at ${\boldsymbol{K}}=0$, simply returns a response-field pre-factor ${b}_{{j}_{c}}(c)$,

$$\begin{eqnarray}&&{\hat{b}}_{{j}_{c}}(c){{Z}_{0}[{\boldsymbol{L}},{\boldsymbol{K}}]| }_{{\boldsymbol{K}}=0}={b}_{{j}_{c}}(c){Z}_{0}[{\boldsymbol{L}},0],\end{eqnarray} \tag{ 16 }$$

given by

$$\begin{eqnarray}&&{b}_{{j}_{c}}(c)=-{\rm{i}}{\vec{k}}_{c}\cdot {\vec{L}}_{{p}_{{j}_{c}}}(c)\end{eqnarray} \tag{ 17 }$$

in terms of ${\vec{L}}_{{p}_{{j}_{c}}}$ as defined in (14). Inserting (14) results in

$$\begin{eqnarray}&&{b}_{{j}_{c}}(c)={\rm{i}}\displaystyle \sum _{s=1}^{r}({\vec{k}}_{c}\cdot {\vec{k}}_{s}){g}_{{qp}}({t}_{s},{t}_{c}){\delta }_{{j}_{c}{j}_{s}},\end{eqnarray} \tag{ 18 }$$

which has two important consequences for our later considerations. First, the position-momentum component ${g}_{{qp}}({t}_{s},{t}_{c})$ appears here, evaluated at the two times t_s and t_c. Due to causality, the particle ensemble can only respond to causes prior to the response, expressed by ${g}_{{qp}}({t}_{s},{t}_{c})=0$ for ${t}_{c}\geqslant {t}_{s}$. Thus, only response-field factors with ${t}_{c}\lt {t}_{s}$ do not vanish. Second, the Kronecker symbol ${\delta }_{{j}_{c}{j}_{s}}$ identifies two particle indices. We shall return to these two properties of the response-field factors later.

Our next goal will now be to evaluate the free generating functional ${Z}_{0}[{\boldsymbol{L}},0]$ after application of an arbitrary number r of density operators. Constructing the phase-space probability distribution $P({\boldsymbol{q}},{\boldsymbol{p}})$ in [6], we have assumed that a statistically homogeneous and isotropic, Gaussian random velocity potential ψ exists such that the momentum $\vec{p}$ at an arbitrary position is its gradient,

$$\begin{eqnarray}&&\vec{p}=\vec{{\rm{\nabla }}}\psi .\end{eqnarray} \tag{ 19 }$$

Continuity then demands that the density contrast δ is its negative Laplacian,

$$\begin{eqnarray}&&\delta =-{\vec{{\rm{\nabla }}}}^{2}\psi .\end{eqnarray} \tag{ 20 }$$

Then, the density-fluctuation power spectrum ${P}_{\delta }(k)$ specifies the initial phase-space probability distribution $P({\boldsymbol{q}},{\boldsymbol{p}})$ completely. As we have shown in [6], it is given by

$$\begin{eqnarray}&&P({\boldsymbol{q}},{\boldsymbol{p}})=\displaystyle \frac{{V}^{-N}}{\sqrt{{(2\pi )}^{3N}\det {C}_{{pp}}}}\,{ \mathcal C }({\boldsymbol{p}})\exp \left(-\displaystyle \frac{1}{2}{{\boldsymbol{p}}}^{\top }{C}_{{pp}}^{-1}{\boldsymbol{p}}\right),\end{eqnarray} \tag{ 21 }$$

where C_pp is the momentum-correlation matrix to be defined in (28) and discussed in the following section. The correlation operator ${ \mathcal C }({\boldsymbol{p}})$ appearing here and defined in [6] can safely be approximated by unity,

$$\begin{eqnarray}&&{ \mathcal C }({\boldsymbol{p}})\approx 1,\end{eqnarray} \tag{ 22 }$$

for correlators evaluated at sufficiently late times if the q–p-component of the Green's function is unbounded. In the cosmological application we are aiming at here, sufficiently late means that the cosmological scale factor a needs to be much larger than the scale factor a_i corresponding to the time when the phase-space distribution of the particles is initially set, $a\gg {a}_{{\rm{i}}}$. Adopting ${a}_{{\rm{i}}}\approx {10}^{-3}$ according to the release of the cosmic microwave background, $a\gt 0.01$ seems safe for approximation (22) to hold.

With the probability distribution (21), the integrations over the momenta in (3) can be carried out straightforwardly. Then, after applying an arbitrary number r of density operators and setting the generator fields to zero, the free generating functional of our microscopic, non-equilibrium, statistical field theory for canonical ensembles of N classical particles enclosed by the volume V can be written in the form

$$\begin{eqnarray}&&{Z}_{0}[{\boldsymbol{L}},0]={V}^{-N}\int {\rm{d}}{\boldsymbol{q}}\exp \left(-\displaystyle \frac{1}{2}{{\boldsymbol{L}}}_{p}^{\top }{C}_{{pp}}{{\boldsymbol{L}}}_{p}+{\rm{i}}\langle {{\boldsymbol{L}}}_{q},{\boldsymbol{q}}\rangle \right),\end{eqnarray} \tag{ 23 }$$

valid at all sufficiently late times. We note that this expression for ${Z}_{0}[{\boldsymbol{L}},0]$ needs to be summed over all different particle configurations, expressed by the indices ${j}_{1},\,\ldots ,\,{j}_{r}$ appearing in (12). The integral in (23) is carried out over all particle positions ${\boldsymbol{q}}$, and ${{\boldsymbol{L}}}_{q}$ and ${{\boldsymbol{L}}}_{p}$ are the position and momentum shift tensors resulting from applying the density operators, with components defined in (14).

If only density operators are applied, the shift tensors ${{\boldsymbol{L}}}_{q}$ and ${{\boldsymbol{L}}}_{p}$ will be sums over as many terms as density operators have been applied, with each term representing the contribution of a particle to the density. Since the possible later application of response-field operators will identify two particles per operation, the number of particles involved will be lowered by one for each response-field operator applied. For each particle, one pair of shift vectors $({\vec{L}}_{{q}_{j}},{\vec{L}}_{{p}_{j}})$ will remain. Thus, if r density operators and m response-field operators have been applied, the number of different particles involved will be $l=r-m=n+m$.

The momentum-correlation matrix is

$$\begin{eqnarray}&&{C}_{{pp}}=\displaystyle \frac{{\sigma }_{1}^{2}}{3}{{ \mathcal I }}_{3}\otimes {{ \mathcal I }}_{N}+\displaystyle \sum _{j\ne k}{C}_{{p}_{j}{p}_{k}}\otimes {E}_{{jk}},\end{eqnarray} \tag{ 24 }$$

where the matrix E_jk singles out the particles j and k from the ensemble of N particles,

$$\begin{eqnarray}&&{E}_{{jk}}:= {\vec{e}}_{j}\otimes {\vec{e}}_{k}.\end{eqnarray} \tag{ 25 }$$

The amplitude ${\sigma }_{1}^{2}$ is defined as the first moment of the power spectrum P_ψ of the velocity potential ψ. Generally, the moments ${\sigma }_{n}^{2}$ are defined by

$$\begin{eqnarray}&&{\sigma }_{n}^{2}:= {\int }_{k}{k}^{2n}\,{P}_{\psi }(k),\end{eqnarray} \tag{ 26 }$$

and the velocity-potential power spectrum is related to the density-fluctuation power spectrum P_δ by

$$\begin{eqnarray}&&{k}^{4}{P}_{\psi }(k)={P}_{\delta }(k)\end{eqnarray} \tag{ 27 }$$

due to (20).

By definition, the momentum-correlation matrix is given by

$$\begin{eqnarray}&&{C}_{{p}_{j}{p}_{k}}={\int }_{k}(\vec{k}\otimes \vec{k}\,){P}_{\psi }(k){{\rm{e}}}^{{\rm{i}}\vec{k}\cdot {\vec{q}}_{{jk}}},\end{eqnarray} \tag{ 28 }$$

where ${\vec{q}}_{{jk}}$ is the separation vector between particles j and k. Expression (28) is equivalent to

$$\begin{eqnarray}&&{C}_{{p}_{j}{p}_{k}}=-(\vec{{\rm{\nabla }}}\otimes \vec{{\rm{\nabla }}}){\int }_{k}{P}_{\psi }(k){{\rm{e}}}^{{\rm{i}}\vec{k}\cdot {\vec{q}}_{{jk}}}=-(\vec{{\rm{\nabla }}}\otimes \vec{{\rm{\nabla }}}){\xi }_{\psi }({q}_{{jk}}),\end{eqnarray} \tag{ 29 }$$

where ${\xi }_{\psi }(q)$ is the correlation function of the velocity potential, taken at distance q and normalised to ${\sigma }_{1}^{2}$. Due to isotropy, ${\xi }_{\psi }$ can only depend on q, but not on the direction of $\vec{q}$.

Accordingly, the tensor of second derivatives $\vec{{\rm{\nabla }}}\otimes \vec{{\rm{\nabla }}}$ needs to be expressed in terms of derivatives with respect to the particle separation q. Let $\hat{q}$ be the unit vector in the direction of $\vec{q}$, and by its means define the projectors

$$\begin{eqnarray}&&{\tilde{\pi }}_{\parallel }:= \hat{q}\otimes \hat{q},\quad {\tilde{\pi }}_{\perp }:= {{ \mathcal I }}_{3}-{\tilde{\pi }}_{\parallel }\end{eqnarray} \tag{ 30 }$$

parallel and perpendicular to $\hat{q}$. Then,

$$\begin{eqnarray}&&\vec{{\rm{\nabla }}}\otimes \vec{{\rm{\nabla }}}={\tilde{\pi }}_{\parallel }\,\displaystyle \frac{{{\rm{d}}}^{2}}{{\rm{d}}{q}^{2}}+{\tilde{\pi }}_{\perp }\,{q}^{-1}\displaystyle \frac{{\rm{d}}}{{\rm{d}}q}\end{eqnarray} \tag{ 31 }$$

and

$$\begin{eqnarray}&&{C}_{{p}_{j}{p}_{k}}=-{\tilde{\pi }}_{\parallel }\,{\xi }_{\psi }^{\prime\prime }({q}_{{jk}})-{\tilde{\pi }}_{\perp }\,\displaystyle \frac{{\xi }_{\psi }^{\prime }({q}_{{jk}})}{{q}_{{jk}}},\end{eqnarray} \tag{ 32 }$$

with the prime denoting the derivative with respect to the argument. The correlation function ${\xi }_{\psi }$ of the velocity potential and its derivatives ${\xi }_{\psi }^{\prime }$ and ${\xi }_{\psi }^{\prime\prime }$ are worked out in appendix A together with accurate fit formulae convenient for fast numerical evaluations.

The quadratic form

$$\begin{eqnarray}&&Q:= {{\boldsymbol{L}}}_{p}^{\top }{C}_{{pp}}{{\boldsymbol{L}}}_{p}\end{eqnarray} \tag{ 33 }$$

remaining in (23) splits into two terms by inserting (24),

$$\begin{eqnarray}&&Q=\displaystyle \frac{{\sigma }_{1}^{2}}{3}\displaystyle \sum _{j}{\vec{L}}_{{p}_{j}}^{2}+\displaystyle \sum _{j\ne i}{\vec{L}}_{{p}_{i}}^{\top }{C}_{{p}_{i}{p}_{j}}{\vec{L}}_{{p}_{j}}.\end{eqnarray} \tag{ 34 }$$

Replacing the sum of squares by a squared sum, we can write instead

$$\begin{eqnarray}&&Q={Q}_{0}-{Q}_{{\rm{D}}}+\displaystyle \sum _{j\ne i}{\vec{L}}_{{p}_{i}}^{\top }{C}_{{p}_{i}{p}_{j}}{\vec{L}}_{{p}_{j}}\end{eqnarray} \tag{ 35 }$$

with the damping terms

$$\begin{eqnarray}&&{Q}_{0}:= \displaystyle \frac{{\sigma }_{1}^{2}}{3}{\left(\displaystyle \sum _{j}{\vec{L}}_{{p}_{j}}\right)}^{2},\quad {Q}_{{\rm{D}}}:= \displaystyle \frac{{\sigma }_{1}^{2}}{3}\displaystyle \sum _{j\ne k}{\vec{L}}_{{p}_{j}}\cdot {\vec{L}}_{{p}_{k}}.\end{eqnarray} \tag{ 36 }$$

We shall see below that Q₀ will vanish identically in important cases and that Q_D has an intuitive and important effect on the time evolution of the density-fluctuation power spectrum. We shall refer to the Q₀ and Q_D as dispersion and diffusion terms, respectively.

We now turn to factorising the generating functional in the form (23), which is a lengthy procedure detailed in appendix B. The essential ideas are that, in a statistically homogeneous field, only relative particle coordinates ${\vec{q}}_{j}-{\vec{q}}_{i}$ must matter, and that all these coordinate differences must be statistically indistinguishable.

The final result of the calculations presented in appendix B is that the free generating functional (23) for a shift tensor ${\boldsymbol{L}}$ with contributions from l different particles can be completely factorised,

$$\begin{eqnarray}&&{Z}_{0}[{\boldsymbol{L}},0]={V}^{-l}{(2\pi )}^{3}{\delta }_{{\rm{D}}}\left(\displaystyle \sum _{j=1}^{l}{\vec{L}}_{{q}_{j}}\right){{\rm{e}}}^{-({Q}_{0}-{Q}_{{\rm{D}}})/2}\displaystyle \prod _{2\leqslant b\lt a}^{l}{\int }_{{k}_{{ab}}}\displaystyle \prod _{1\leqslant k\lt j}^{l}({{\rm{\Delta }}}_{{jk}}+{{ \mathcal P }}_{{jk}}).\end{eqnarray} \tag{ 37 }$$

The index pairs (a, b) and (j, k) are defined in (B.5) and (B.8), respectively. The function ${{ \mathcal P }}_{{jk}}$ appearing in each of the factors in (37),

$$\begin{eqnarray}&&{{ \mathcal P }}_{{jk}}({k}_{{jk}},\tau )={\int }_{q}\{{{\rm{e}}}^{{g}_{{qp}}^{2}(\tau ,0){k}_{{jk}}^{2}({a}_{\parallel }{\lambda }_{{jk}}^{\parallel }+{a}_{\perp }{\lambda }_{{jk}}^{\perp })}-1\}{{\rm{e}}}^{{\rm{i}}{\vec{k}}_{{jk}}\cdot \vec{q}},\end{eqnarray} \tag{ 38 }$$

is a nonlinearly (and non-trivially) time-evolved density-fluctuation power spectrum, as we shall discuss in detail in the next section. The expression ${{\rm{\Delta }}}_{{jk}}$ abbreviates

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{jk}}:= {(2\pi )}^{3}{\delta }_{{\rm{D}}}({\vec{k}}_{{jk}}).\end{eqnarray} \tag{ 39 }$$

The wave vectors ${\vec{k}}_{{jk}}$ are defined in (B.15), with the indices (j, k) given in (B.5) and the indices (a, b) defined in (B.8). We shall call the wave vectors ${\vec{k}}_{j1}$ external because they contain the shift vectors ${\vec{L}}_{{q}_{j}}$, and the remaining wave vectors ${\vec{k}}_{{jk}}$ with $k\geqslant 2$ internal because they can entirely be integrated out. The quantities ${\lambda }_{{jk}}^{\parallel ,\perp }$ are defined by

$$\begin{eqnarray}&&{\lambda }_{{jk}}^{\parallel }=\displaystyle \frac{{\vec{L}}_{{p}_{j}}^{\top }{\pi }_{{jk}}^{\parallel }{\vec{L}}_{{p}_{k}}}{{g}_{{qp}}^{2}(\tau ,0){k}_{{jk}}^{2}}\quad {\rm{and}}\quad {\lambda }_{{jk}}^{\perp }=\displaystyle \frac{{\vec{L}}_{{p}_{j}}^{\top }{\pi }_{{jk}}^{\perp }{\vec{L}}_{{p}_{k}}}{{g}_{{qp}}^{2}(\tau ,0){k}_{{jk}}^{2}}\end{eqnarray} \tag{ 40 }$$

according to (B.30), with the projectors ${\pi }_{{jk}}^{\parallel }$ and ${\pi }_{{jk}}^{\perp }$ being defined with respect to ${\vec{k}}_{{jk}}$. Let $\hat{k}$ be the unit vector in the direction of ${\vec{k}}_{{jk}}$, then

$$\begin{eqnarray}&&{\pi }_{{jk}}^{\parallel }=\hat{k}\otimes \hat{k},\quad {\pi }_{{jk}}^{\perp }={{ \mathcal I }}_{3}-{\pi }_{{jk}}^{\parallel }.\end{eqnarray} \tag{ 41 }$$

The factorisation (37) of the free generating functional and the expression (38) for the time-evolved density-fluctuation power spectrum are the first main results of our paper. The power spectrum ${{ \mathcal P }}_{{jk}}$ is particularly relevant for cosmological structure formation.

The main goal of this section is to develop a systematic approach to evaluating (37). It begins by realising that the product over the index pair (j, k) in (37) can be expanded into a sum ordered by an increasing power of power-spectrum factors. Since there are

$$\begin{eqnarray}&&F=\displaystyle \frac{l(l-1)}{2}\end{eqnarray} \tag{ 52 }$$

factors in this product, we can write formally

$$\begin{eqnarray}&&\displaystyle \prod _{j\gt k}({{\rm{\Delta }}}_{{jk}}+{{ \mathcal P }}_{{jk}})=\displaystyle \sum _{f=0}^{F}\displaystyle \left(\genfrac{}{}{0em}{}{F}{f}\right)\,{{\rm{\Delta }}}^{F-f}\,{{ \mathcal P }}^{f}.\end{eqnarray} \tag{ 53 }$$

If the power-spectrum factors ${ \mathcal P }$ are small enough, this sum can be truncated at low powers f.

According to (37), all internal wave vectors ${\vec{k}}_{{ab}}={\vec{k}}_{{ab}}^{\prime }$ are then to be integrated over. For evaluating the result of this sequence of integrals applied to the sum terms in (53), it is important to see how many of these integrals can trivially be carried out due to the Δ factors appearing in these sum terms. For f = 0, for example, the remaining integrals set all of the internal wave vectors to zero, leaving ${\vec{k}}_{j1}={\vec{L}}_{{q}_{j}}$ and

$$\begin{eqnarray}&&\displaystyle \prod _{a\gt b}{\int }_{{k}_{{ab}}}\displaystyle \prod _{j\gt k}{{\rm{\Delta }}}_{{jk}}=\displaystyle \prod _{j=2}^{l}{\rm{\Delta }}({\vec{L}}_{{q}_{j}}).\end{eqnarray} \tag{ 54 }$$

This is a pure shot-noise term. Identifying and counting the power-spectrum factors with external or internal wave numbers thus allows to simplify the remaining expressions substantially.

Let us illustrate this simplification with one more example. For f = 1, the single power-spectrum factor can depend either on an external or an internal wave number. If it is external, the Δ factors set all internal wave numbers to zero as for f = 0 before, leaving again the external wave numbers ${\vec{k}}_{j1}={\vec{L}}_{{q}_{j}}$. All of these wave vectors ${\vec{k}}_{j1}$ except one are also set to zero by the remaining Δ factors. Without loss of generality, we can label the one non-zero wave vector by the index j = 2 because the particles are indistinguishable. Then, the power-spectrum factor ${ \mathcal P }({\vec{L}}_{{q}_{2}})$ appears besides the remaining Δ factors. If, however, the single power-spectrum factor depends on an internal wave number, all except one internal wave numbers are integrated out. Again without loss of generality, we can label the only remaining internal wave number ${\vec{k}}_{32}={\vec{k}}_{32}^{\prime }$. The external wave vectors are then

$$\begin{eqnarray}&&{\vec{k}}_{21}={\vec{L}}_{{q}_{2}}+{\vec{k}}_{32}^{\prime },\quad {\vec{k}}_{31}={\vec{L}}_{{q}_{3}}-{\vec{k}}_{32}^{\prime }\end{eqnarray} \tag{ 55 }$$

and ${\vec{k}}_{j1}={\vec{L}}_{{q}_{j}}$ for $j\gt 3$. All external wave vectors appear as arguments of Δ factors in this case. The remaining integration over ${\vec{k}}_{32}^{\prime }$ sets ${\vec{k}}_{32}^{\prime }={\vec{L}}_{{q}_{3}}$ by means of the factor ${{\rm{\Delta }}}_{31}$ and leaves ${\vec{k}}_{21}={\vec{L}}_{{q}_{2}}+{\vec{L}}_{{q}_{3}}$.

There are $l-1$ external and $F-l+1$ internal wave vectors. Thus, there are $l-1$ possibilities for choosing an external and $F-l+1$ possibilities for choosing an internal wave number, allowing us to write

$$\begin{eqnarray}&&\displaystyle \prod _{a\gt b}{\int }_{{k}_{{ab}}}\displaystyle \left(\genfrac{}{}{0em}{}{F}{1}\right)\,{ \mathcal P }{{\rm{\Delta }}}^{F-1}=\displaystyle \left(\genfrac{}{}{0em}{}{l-1}{1}\right)\,{{ \mathcal P }}_{21}({\vec{L}}_{{q}_{2}})\displaystyle \prod _{j=3}^{l}{{\rm{\Delta }}}_{j1}+\displaystyle \left(\genfrac{}{}{0em}{}{F-l+1}{1}\right)\,{{ \mathcal P }}_{32}({\vec{L}}_{{q}_{3}}){{\rm{\Delta }}}_{21}({\vec{L}}_{{q}_{1}})\displaystyle \prod _{j=4}^{l}{{\rm{\Delta }}}_{j1},\end{eqnarray} \tag{ 56 }$$

where the overall Dirac delta distribution

$$\begin{eqnarray}&&{\delta }_{{\rm{D}}}\left(\displaystyle \sum _{j=1}^{l}{\vec{L}}_{{q}_{j}}\right)\end{eqnarray} \tag{ 57 }$$

in (37) allows us to replace ${{\rm{\Delta }}}_{21}({\vec{L}}_{{q}_{2}}+{\vec{L}}_{{q}_{3}})$ by ${{\rm{\Delta }}}_{21}({\vec{L}}_{{q}_{1}})$.

From these examples for f = 0 and f = 1, a straightforward scheme emerges for evaluating the sum terms in (53) and the subsequent integrations over all internal wave vectors:

• (i)  
  Expand each term into a sum ordered by the number of power-spectrum factors depending on external wave numbers,

  $$\begin{eqnarray}&&\displaystyle \left(\genfrac{}{}{0em}{}{F}{f}\right){{\rm{\Delta }}}^{F-f}{{ \mathcal P }}^{f}=\displaystyle \sum _{e=0}^{f}\displaystyle \left(\genfrac{}{}{0em}{}{l-1}{e}\right)\displaystyle \left(\genfrac{}{}{0em}{}{F-l+1}{f-e}\right)\,{{ \mathcal P }}_{\mathrm{ext}}^{e}{{ \mathcal P }}_{\mathrm{int}}^{f-e}{{\rm{\Delta }}}_{\mathrm{ext}}^{l-e-1}{{\rm{\Delta }}}_{\mathrm{int}}^{F-l+1-f+e},\end{eqnarray} \tag{ 58 }$$

  where the subscripts 'ext' and 'int' indicate that respective factors depend on external or internal wave vectors.
• (ii)  
  Use the Δ factors depending on internal wave vectors to set $(F-l+1-f+e)$ of the internal wave vectors to zero. Label the remaining $(f-e)$ internal wave vectors beginning with ${\vec{k}}_{32}$ and determine the external wave vectors ${\vec{k}}_{j1}$ according to (B.15).
• (iii)  
  Use the Δ factors depending on external wave vectors and the integrals over the remaining internal wave vectors to eliminate as many of the internal wave numbers as possible from the external wave vectors ${\vec{k}}_{j1}$ and from the arguments of the power-spectrum factors depending on internal wave numbers.

This procedure lends itself to be evaluated by a symbolic computer code returning the terms appearing in the sum (58), given l and f. We are now in the process of developing such a code, aiming at possibly completely automatising the evaluation of perturbation terms.

As often in perturbation theories, diagrams are useful to keep an overview of the terms involved. We shall now develop suitable diagrams for the perturbative approach to our non-equilibrium field theory for correlated classical particle ensembles.

We return to the free generating functional ${Z}_{0}[{\boldsymbol{L}},0]$, evaluated at a specific shift tensor ${\boldsymbol{L}}$, with the generator field ${\boldsymbol{K}}$ set to zero. The position and momentum components of ${\boldsymbol{L}}$ are given in (14) or, after inserting (13) there, by

$$\begin{eqnarray}{\vec{L}}_{{q}_{j}} & = & -\displaystyle \sum _{s=1}^{r}{\vec{k}}_{s}\,{\delta }_{{{jj}}_{s}},\\ {\vec{L}}_{{p}_{j}} & = & -\displaystyle \sum _{s=1}^{r}{\vec{k}}_{s}\,{g}_{{qp}}({t}_{s},0){\delta }_{{{jj}}_{s}}.\end{eqnarray} \tag{ 59 }$$

In these expressions, the index $1\leqslant j\leqslant l$ identifies the l particles whose positions are being correlated, while the indices j_s assign particles to wave vectors with field labels s. The Kronecker symbols appear in (59) because only such phase-space positions contribute to the density which are occupied by particles. Given a specific set of particle indices $\{{j}_{1},\,\ldots ,\,{j}_{r}\}$, we write

$$\begin{eqnarray}&&{Z}_{0}[\{{j}_{1},\,\ldots ,\,{j}_{r}\}]:= {Z}_{0}[{\boldsymbol{L}},0]\end{eqnarray} \tag{ 60 }$$

to denote the particle indices explicitly. As indicated above, the response-field factors given by (18) have two crucial properties affecting the selection of terms that can or cannot contribute to the perturbation series: the Kronecker symbol appearing there identifies the two particles with indices j_s and j_c and thus assigns the same particle to the wave vectors labelled by c and s. Furthermore, since the propagator ${g}_{{qp}}({t}_{s},{t}_{c})$ vanishes if ${t}_{s}\leqslant {t}_{c}$, only such terms can contribute for which ${t}_{s}\gt {t}_{c}$.

Applying r one-particle density and m one-particle response-field operators to the free generating functional thus leaves us with the expression

$$\begin{eqnarray}&&\displaystyle \prod _{c=r-m+1}^{r}{b}_{{j}_{c}}(c){Z}_{0}[\{{j}_{1},\,\ldots ,\,{j}_{r}\}]=\displaystyle \prod _{c=r-m+1}^{r}\left({\rm{i}}\displaystyle \sum _{s=1}^{r}({\vec{k}}_{c}\cdot {\vec{k}}_{s}){g}_{{qp}}({t}_{s},{t}_{c}){\delta }_{{j}_{c}{j}_{s}}\right){Z}_{0}[\{{j}_{1},\,\ldots ,\,{j}_{r}\}].\end{eqnarray} \tag{ 61 }$$

Several aspects are important to note at this point. First, each response-field factor with index c identifies two particles with the indices j_c and ${j}_{{s}_{c}}$. Second, no terms can contribute to any response-field factor which belong to the same times because ${g}_{{qp}}(t,t)=0$. Therefore, only such particles may be identified which are assigned to wave vectors at different times. Fourth, since ${g}_{{qp}}({t}_{1},{t}_{2})=0$ for ${t}_{2}\geqslant {t}_{1}$, any response-field factor implies a time ordering in the sense ${t}_{1}\gt {t}_{2}$.

Beginning with a set of particle indices $\{{j}_{1},\,\ldots ,\,{j}_{r}\}$, we thus proceed as follows to evaluate the terms in (61): we first identify as many particle pairs as response fields occur, i.e. we identify m pairs of particle indices, taking care that particle pairs must not identify positions with equal times and that a time-ordering is involved in all response-field factors. The remaining r − m particles form a reduced set $\{j\}^{\prime}$ of particle indices. All possible reduced sets $\{j\}^{\prime}$ must finally be summed over.

Let us illustrate this procedure with the simple example r = 4 and m = 1, corresponding to the terms contributing to a two-point density correlator calculated at first-order in the particle interaction. Since the interaction is instantaneous, ${t}_{3}={t}_{4}$, and if the correlator is chosen to be simultaneous, ${t}_{1}={t}_{2}$.

In this case, (61) simplifies to the two possible terms

$$\begin{eqnarray} & & {\rm{i}}({\vec{k}}_{4}\cdot {\vec{k}}_{1}){g}_{{qp}}({t}_{1},{t}_{4}){Z}_{0}[\{{j}_{1}={j}_{4},{j}_{2},{j}_{3}\}],\\ & & {\rm{i}}({\vec{k}}_{4}\cdot {\vec{k}}_{2}){g}_{{qp}}({t}_{2},{t}_{4}){Z}_{0}[\{{j}_{1},{j}_{2}={j}_{4},{j}_{3}\}].\end{eqnarray} \tag{ 62 }$$

Other terms would identify particles at the same time and thus return zero.

Combining (37) and (53) can profitably be represented by diagrams which greatly help constructing and ordering the terms appearing in the generating functional. The essential point of the diagrammatic representation which we are going to construct now is to systematically construct all wave vectors ${\vec{k}}_{{jk}}$ according to (B.15) which enter into the factors $({{ \mathcal P }}_{{jk}}+{{\rm{\Delta }}}_{{jk}})$ appearing in the generating functional.

The diagrams are constructed according to the following rules:

• (i)  
  Mark the free generating functional ${Z}_{0}[{\boldsymbol{J}},{\boldsymbol{K}}]$ by a circle. According to (10) and (12), each one-particle density operator ${\hat{\rho }}_{{j}_{s}}$ applied to Z₀ corresponds to a functional derivative with respect to a component ${\vec{J}}_{{q}_{{j}_{s}}}$ of the generator field ${\boldsymbol{J}}$. According to (12) and (13), each of these operations adds a wave vector ${\vec{k}}_{s}$ to the shift tensor ${\boldsymbol{L}}$ that we need to find for each term in a perturbation series.Thus, for each of the $r=n+2m$ density operators applied for a perturbation term appearing in an nth order correlator at mth order in the particle interaction, we attach a wave vector to the circle symbolising Z₀, pointing outward.
• (ii)  
  The times when these density operators act are represented by filled dots on the circumference of the Z₀ symbol. Each $\vec{k}$ vector thus begins at a filled dot. Since each interaction is instantaneous, all internal $\vec{k}$ vectors need to be pairwise attached to the same time. If the correlator to be calculated is simultaneous, the external $\vec{k}$ vectors are also attached to the same time.Thus, the internal wave vectors representing interactions are pairwise attached to the same points in time. Times are ordered counter-clockwise, with the latest times appearing on top. The external wave vectors appearing in the final correlator are attached to the same time if the correlator is simultaneous.
• (iii)  
  Between the internal wave vectors corresponding to the density and the response-field operator of an interaction term (9), a further factor appears representing the interaction potential v. If this potential depends only on the separation between the two particles interacting, it requires the two $\vec{k}$ vectors of the density and the response-field operator it connects to be equal in magnitude and opposite in sign.Internal wave vectors are marked with primes. To include the interaction potential into the diagram, we connect any two internal wave vectors with a circled v which enter into a single interaction term. For translation-invariant potentials v, the two internal wave vectors connected by a potential are equal and opposite.
• (iv)  
  According to (18), each response field identifies two particles at different times, i.e. it assigns the same particle to two different wave vectors at two different times. Adopting the convention in (9), we assign response fields to the negative internal wave vectors appearing in the interaction terms.Thus, response fields are represented by dashed circle segments between two different wave vectors attached to different times. Each response field must begin at a negative internal wave vector and must be connected to exactly one other wave vector, internal or external.
• (v)  
  Equivalent diagrams can appear multiple times. For example, in the diagram shown in figure 5 representing a contribution to the density power spectrum at first order in the particle interactions, the response-field arrow can end on ${\vec{k}}_{1}$ or ${\vec{k}}_{2}$. In a homogeneous random field, these two wave vectors are equivalent, and the diagram with the response-field arrow attached to the vector ${\vec{k}}_{2}$ corresponds to an identical perturbation term.Thus, diagrams are assigned a multiplicity corresponding to the number of equivalent configurations they express.

Zoom In Zoom Out Reset image size

Figure 5 shows the single diagram according to these rules representing the first-order interaction contribution to the second-order density power spectrum.

Diagrams for higher-order interactions or correlators of higher order are now easily constructed. To give an example, we show in figure 6 the four non-equivalent diagrams contributing to a second-order density correlator at second order in the particle interactions.

Zoom In Zoom Out Reset image size

Position and momentum shift vectors according to (59) can now be read off these diagrams as follows: randomly assign particle indices to the $\vec{k}$ vectors extending from Z₀, thereby assigning the same particle index to such $\vec{k}$ vectors connected by a dashed line representing a response field. These will be r − m indices in total, for which the numbers $1\ldots r-m$ can be chosen without loss of generality. Any permutation of these indices will result in an equivalent set of shift vectors.

For example, we can use the diagram in figure 5 to assign the particle indices $(1,2,3)$ clock-wise to the $\vec{k}$ vectors extending from Z₀, beginning with ${\vec{k}}_{1}$. This results in the assignment of particles $(1,2,3,1)$ to the wave vectors $({\vec{k}}_{1},{\vec{k}}_{2},{\vec{k}}_{1}^{\prime },-{\vec{k}}_{1}^{\prime })$ because ${\vec{k}}_{1}$ and $-{\vec{k}}_{1}^{\prime }$ are connected by a response field. With (59), this implies the position shift vectors

$$\begin{eqnarray}&&{\vec{L}}_{{q}_{1}}=-({\vec{k}}_{1}-{\vec{k}}_{1}^{\prime }),\quad {\vec{L}}_{{q}_{2}}=-{\vec{k}}_{2},\quad {\vec{L}}_{{q}_{3}}=-{\vec{k}}_{1}^{\prime }\end{eqnarray} \tag{ 63 }$$

and the momentum shift vectors

$$\begin{eqnarray}{\vec{L}}_{{p}_{1}} & = & -({g}_{{qp}}(t,0){\vec{k}}_{1}-{g}_{{qp}}(t^{\prime} ,0){\vec{k}}_{1}^{\prime }),\quad {\vec{L}}_{{p}_{2}}=-{g}_{{qp}}(t,0){\vec{k}}_{2},\\ {\vec{L}}_{{p}_{3}} & = & -{g}_{{qp}}(t^{\prime} ,0){\vec{k}}_{1}^{\prime },\end{eqnarray} \tag{ 64 }$$

with $t^{\prime}$ denoting the time of the interaction and $t\gt t^{\prime}$ the time where the correlator is to be evaluated. Any permutation of the particle indices would merely permute the labels on these shift vectors.

These diagrams allow a quick construction of all terms contributing to the generating functional Z for given orders of correlators and of particle interactions. From these diagrams, the shift vectors ${\vec{L}}_{q}$ and ${\vec{L}}_{p}$ as well as the response-field factors can be read off. They can then be inserted into the complete factorisation of the generating functional to evaluate the perturbation terms. The procedures involved can now be implemented in a symbolic computer code.

The potential correlation function ${\xi }_{\psi }(x)$ introduced in (29) is

$$\begin{eqnarray}&&{\xi }_{\psi }(q)=\displaystyle \frac{1}{2{\pi }^{2}}{\int }_{0}^{\infty }{k}^{2}{\rm{d}}k\,{P}_{\psi }(k){{j}}_{0}({kq}),\end{eqnarray} \tag{ A.1 }$$

where P_ψ the power spectrum of the velocity potential. With

$$\begin{eqnarray}&&{\partial }_{q}^{\alpha }{\xi }_{\psi }(q)=\displaystyle \frac{1}{2{\pi }^{2}}{\int }_{0}^{\infty }{k}^{2+\alpha }{\rm{d}}k\,{P}_{\psi }(k){{j}}_{0}^{(\alpha )}({kq})\end{eqnarray} \tag{ A.2 }$$

and using the recursion relations for the spherical Bessel functions and their derivatives, we find

$$\begin{eqnarray}&&{\xi }_{\psi }^{\prime }(q)=-\displaystyle \frac{1}{2{\pi }^{2}}{\int }_{0}^{\infty }\displaystyle \frac{{\rm{d}}k}{k}{P}_{\delta }(k){{j}}_{1}({kq}),\end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}&&{\xi }_{\psi }^{\prime\prime }(q)-\displaystyle \frac{{\xi }_{\psi }^{\prime }(q)}{q}=\displaystyle \frac{1}{2{\pi }^{2}}{\int }_{0}^{\infty }{\rm{d}}k\,{P}_{\delta }(k){{j}}_{2}({kq}),\end{eqnarray} \tag{ A.4 }$$

where the density-fluctuation power spectrum ${P}_{\delta }(k)={k}^{4}{P}_{\psi }(k)$ was introduced. The asymptotic behaviour for small arguments of the spherical Bessel functions of the first kind ensure that

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{q\to 0}\displaystyle \frac{{\xi }_{\psi }^{\prime }(q)}{q}=-\displaystyle \frac{{\sigma }_{1}^{2}}{3}=\mathop{\mathrm{lim}}\limits_{q\to 0}{\xi }_{\psi }^{\prime\prime }(q).\end{eqnarray} \tag{ A.5 }$$

The functions ${\xi }_{\psi }^{\prime }(q)/q$ and ${\xi }_{\psi }^{\prime\prime }(q)$ are shown in figure A1.

Zoom In Zoom Out Reset image size

In view of later numerical evaluations, it will be advantageous to introduce

$$\begin{eqnarray}&&{a}_{1}(q):= \displaystyle \frac{{\xi }_{\psi }^{\prime }(q)}{q},\quad {a}_{2}(q):= {\xi }_{\psi }^{\prime\prime }(q)-\displaystyle \frac{{\xi }_{\psi }^{\prime }(q)}{q}.\end{eqnarray} \tag{ A.6 }$$

For the ΛCDM density-fluctuation power spectrum, it turns out that these two functions a₁ and a₂ allow simple fits by rational functions,

$$\begin{eqnarray}{a}_{1} & = & {A}_{1}{\left(1+\displaystyle \frac{q}{{q}_{11}}\right)}^{-{\alpha }_{1}},\\ {a}_{2} & = & {A}_{2}{q}^{{\alpha }_{2}}{\left(1+{\left(\displaystyle \frac{q}{{q}_{21}}\right)}^{{\beta }_{2}}+{\left(\displaystyle \frac{q}{{q}_{22}}\right)}^{{\gamma }_{2}}\right)}^{-{\delta }_{2}},\end{eqnarray} \tag{ A.7 }$$

both of which fall like ${q}^{-2}$ asymptotically for q ≫ q₁ or q ≫ q₃. For a ΛCDM power spectrum according to [22] with cosmological parameters ${{\rm{\Omega }}}_{{\rm{m}}0}=0.3$ and ${{\rm{\Omega }}}_{{\rm{\Lambda }}0}=0.7$, we find the best-fitting parameters listed in table A1.

Table A1.  Parameters of the rational functions (A.7) fitting the functions a₁ and a₂ defined in (A.6) for a ΛCDM power spectrum according to [22] with cosmological parameters ${{\rm{\Omega }}}_{{\rm{m}}0}=0.3$ and ${{\rm{\Omega }}}_{{\rm{\Lambda }}0}=0.7$.

$\mathrm{ln}{A}_{1}$	−10.0039	$\mathrm{ln}{A}_{1}$	−11.2727
q11	61.6918	q21	0.3608
		q22	11.9575
${\alpha }_{1}$	2.0507	${\alpha }_{2}$	2.044 59
		${\beta }_{2}$	0.5645
		${\gamma }_{2}$	1.7061
		${\delta }_{2}$	2.3780

Figure A2 shows the functions ${a}_{\mathrm{1,2}}(q)$ together with the fits (A.7) specified by the parameters in table A1.

Zoom In Zoom Out Reset image size

Returning to (23), we introduce the coordinate differences with respect to the arbitrarily chosen particle 1,

$$\begin{eqnarray}&&{\vec{q}}_{j1}:= {\vec{q}}_{j}-{\vec{q}}_{1}\quad \forall \quad j=2\ldots l.\end{eqnarray} \tag{ B.1 }$$

Then, the scalar product of the positions ${\boldsymbol{q}}$ with the spatial shift tensor ${{\boldsymbol{L}}}_{q}$ can be written as

$$\begin{eqnarray}&&\langle {{\boldsymbol{L}}}_{q},{\boldsymbol{q}}\rangle ={\vec{L}}_{{q}_{1}}\cdot {\vec{q}}_{1}+{\vec{L}}_{{q}_{2}}\cdot ({\vec{q}}_{21}+{\vec{q}}_{1})\,+\,\ldots \,=\,\left(\displaystyle \sum _{j=1}^{l}{\vec{L}}_{{q}_{j}}\right){\vec{q}}_{1}+\displaystyle \sum _{j=2}^{l}{\vec{L}}_{{q}_{j}}\cdot {\vec{q}}_{j1}.\end{eqnarray} \tag{ B.2 }$$

Since the correlation matrix C_pp depends on coordinate differences only and not on absolute positions, the integration over ${\vec{q}}_{1}$ can be carried out. This results in

$$\begin{eqnarray}&&{Z}_{0}[{\boldsymbol{L}},0]=\displaystyle \frac{{(2\pi )}^{3}}{{V}^{l}}{\delta }_{{\rm{D}}}\left(\displaystyle \sum _{j=1}^{l}{\vec{L}}_{{q}_{j}}\right)\displaystyle \prod _{j=2}^{l}{\int }_{{q}_{j1}}\exp \left(-\displaystyle \frac{1}{2}{{\boldsymbol{L}}}_{p}^{\top }{C}_{{pp}}{{\boldsymbol{L}}}_{p}+{\rm{i}}\displaystyle \sum _{j=2}^{l}{\vec{L}}_{{q}_{j}}\cdot {\vec{q}}_{j1}\right).\end{eqnarray} \tag{ B.3 }$$

The momentum-correlation matrix depends on the absolute values of all pair-wise coordinate differences

$$\begin{eqnarray}&&{\vec{q}}_{{jk}}:= {\vec{q}}_{j}-{\vec{q}}_{k}\end{eqnarray} \tag{ B.4 }$$

with

$$\begin{eqnarray}&&k=1\ldots (l-1),\quad j=(k+1)\ldots l,\end{eqnarray} \tag{ B.5 }$$

not just on the coordinate differences ${\vec{q}}_{j1}$ with respect to particle number 1. Since there are l particles to consider in total, the number of coordinate differences that C_pp depends on is

$$\begin{eqnarray}&&{N}_{\mathrm{pairs}}=\displaystyle \frac{l(l-1)}{2}.\end{eqnarray} \tag{ B.6 }$$

Of these coordinate differences, $(l-1)$ are taken into account by the $(l-1)$ difference vectors ${\vec{q}}_{j1}$. We now extend the integration in (B.3) to all these coordinate differences by introducing the remaining difference vectors

$$\begin{eqnarray}&&{\vec{q}}_{{ab}}={\vec{q}}_{a}-{\vec{q}}_{b}={\vec{q}}_{a1}-{\vec{q}}_{b1}\end{eqnarray} \tag{ B.7 }$$

with

$$\begin{eqnarray}&&b=2\ldots (l-1),\quad a=(b+1)\ldots l.\end{eqnarray} \tag{ B.8 }$$

These contribute

$$\begin{eqnarray}&&\displaystyle \frac{l(l-1)}{2}-(l-1)=\displaystyle \frac{(l-1)(l-2)}{2}\end{eqnarray} \tag{ B.9 }$$

additional, dependent coordinates related by (B.7), which must be ensured by including appropriate delta distributions

$$\begin{eqnarray}&&{\delta }_{{\rm{D}}}({\vec{q}}_{{ab}}-{\vec{q}}_{a1}+{\vec{q}}_{b1})\end{eqnarray} \tag{ B.10 }$$

into the integrand, with a and b from (B.8). Thus, the free generating functional turns into

$$\begin{eqnarray}{Z}_{0}[{\boldsymbol{L}},0] & = & \displaystyle \frac{{(2\pi )}^{3}}{{V}^{l}}{\delta }_{{\rm{D}}}\left(\displaystyle \sum _{j=1}^{l}{\vec{L}}_{{q}_{j}}\right)\displaystyle \prod _{j\gt k}{\displaystyle \int }_{{q}_{{jk}}}\exp \left(-\displaystyle \frac{1}{2}{{\boldsymbol{L}}}_{p}^{\top }{C}_{{pp}}{{\boldsymbol{L}}}_{p}+{\rm{i}}\displaystyle \sum _{j=2}^{l}{\vec{L}}_{{q}_{j}}\cdot {\vec{q}}_{j1}\right)\\ & & \cdot \displaystyle \prod _{a\gt b}{\delta }_{{\rm{D}}}({\vec{q}}_{{ab}}-{\vec{q}}_{a1}+{\vec{q}}_{b1}).\end{eqnarray} \tag{ B.11 }$$

Replacing all of these delta distributions by their Fourier expansions,

$$\begin{eqnarray}&&{\delta }_{{\rm{D}}}({\vec{q}}_{{ab}}-{\vec{q}}_{a1}+{\vec{q}}_{b1})={\int }_{{k}_{{ab}}}{{\rm{e}}}^{{\rm{i}}{\vec{k}}_{{ab}}^{\prime }\cdot ({\vec{q}}_{{ab}}-{\vec{q}}_{a1}+{\vec{q}}_{b1})}\end{eqnarray} \tag{ B.12 }$$

with auxiliary wave vectors ${\vec{{k}_{ab}^{\prime }}}_{}^{}$ conjugate to the coordinate differences ${\vec{q}}_{{ab}}$, we arrive at the expression

$$\begin{eqnarray}{Z}_{0}[{\boldsymbol{L}},0] & = & \displaystyle \frac{{(2\pi )}^{3}}{{V}^{l}}{\delta }_{{\rm{D}}}\left(\displaystyle \sum _{j=1}^{l}{\vec{L}}_{{q}_{j}}\right)\displaystyle \prod _{a\gt b}{\displaystyle \int }_{{{k}_{ab}^{\prime }}_{}^{}}\displaystyle \prod _{j\gt k}{\displaystyle \int }_{{q}_{{jk}}}\\ & & \cdot \exp \left(-\displaystyle \frac{1}{2}{{\boldsymbol{L}}}_{p}^{\top }{C}_{{pp}}{{\boldsymbol{L}}}_{p}+{\rm{i}}\displaystyle \sum _{j=2}^{l}{\vec{L}}_{{q}_{j}}\cdot {\vec{q}}_{j1}+{\rm{i}}\displaystyle \sum _{a\gt b}{\vec{k}}_{{ab}}^{\prime }\cdot ({\vec{q}}_{{ab}}-{\vec{q}}_{a1}+{\vec{q}}_{b1})\right)\end{eqnarray} \tag{ B.13 }$$

for the free generating functional. Reordering terms, we can rewrite the phase factor as

$$\begin{eqnarray}&&\displaystyle \sum _{j=2}^{l}{\vec{L}}_{{q}_{j}}\cdot {\vec{q}}_{j1}+\displaystyle \sum _{a\gt b}{\vec{k}}_{{ab}}^{\prime }\cdot ({\vec{q}}_{{ab}}-{\vec{q}}_{a1}+{\vec{q}}_{b1})\\ &&\quad =\,\displaystyle \sum _{j=2}^{l}\left({\vec{L}}_{{q}_{j}}-\displaystyle \sum _{b=2}^{j-1}{\vec{k}}_{{jb}}^{\prime }+\displaystyle \sum _{a=j+1}^{l}{\vec{{k}_{aj}^{}}^{\prime} }_{}^{}\right)\cdot {\vec{q}}_{j1}+\displaystyle \sum _{a\gt b}{\vec{k}}_{{ab}}^{\prime }\cdot {\vec{q}}_{{ab}}=: \displaystyle \sum _{j\gt k}{\vec{k}}_{{jk}}\cdot {\vec{q}}_{{jk}}\end{eqnarray} \tag{ B.14 }$$

with

$$\begin{eqnarray}{\vec{k}}_{{jk}}:= \left\{\begin{array}{lll}{\vec{L}}_{{q}_{j}}-\displaystyle \sum _{b=2}^{j-1}{\vec{k}}_{{jb}}^{\prime }+\displaystyle \sum _{a=j+1}^{l}{\vec{k}}_{{aj}}^{\prime } & {\rm{for}} & \quad k=1,\quad j=2\ldots l\\ {\vec{k}}_{{jk}}^{\prime } & {\rm{for}} & \quad k=2\ldots (l-1),\quad j=(k+1)\,...\,l\end{array}\right.\end{eqnarray} \tag{ B.15 }$$

and $a,b$ from (B.8).

With these results and the definition of Q₀ and Q_D in (35), the generating functional (B.13) factorises completely in the integrations over the coordinate differences ${\vec{q}}_{{jk}}$,

$$\begin{eqnarray}&&{Z}_{0}[{\boldsymbol{L}},0]=\displaystyle \frac{{(2\pi )}^{3}}{{V}^{l}}{\delta }_{{\rm{D}}}\left(\displaystyle \sum _{j=1}^{l}{\vec{L}}_{{q}_{j}}\right){{\rm{e}}}^{-({Q}_{0}-{Q}_{{\rm{D}}})/2}\displaystyle \prod _{a\gt b}{\int }_{{k}_{{ab}}}\displaystyle \prod _{j\gt k}{I}_{{jk}}\end{eqnarray} \tag{ B.16 }$$

with

$$\begin{eqnarray}&&{I}_{{jk}}:= {\int }_{{q}_{{jk}}}{{\rm{e}}}^{-{\vec{L}}_{{p}_{j}}^{\top }{C}_{{p}_{j}{p}_{k}}{\vec{L}}_{{p}_{k}}+{\rm{i}}{\vec{k}}_{{jk}}\cdot {\vec{q}}_{{jk}}}\quad (j\ne k)\end{eqnarray} \tag{ B.17 }$$

since the correlation matrix ${C}_{{p}_{j}{p}_{k}}$ depends on the distances $| {\vec{q}}_{{jk}}|$ between the particles only. We can thus decompose the generating functional into independent factors for all particle pairs, which are then to be convolved in Fourier space by integrating over all auxiliary wave vectors ${\vec{k}}_{{ab}}^{\prime }$.

For two-point density correlations, $l=2+m$, hence the number of coordinate pairs is

$$\begin{eqnarray}&&{N}_{\mathrm{pair}}=\displaystyle \frac{(m+2)(m+1)}{2},\end{eqnarray} \tag{ B.18 }$$

of which

$$\begin{eqnarray}&&\displaystyle \frac{(l-1)(l-2)}{2}=\displaystyle \frac{m(m+1)}{2}\end{eqnarray} \tag{ B.19 }$$

are dependent.

For first-order perturbation theory of a two-point spectrum, we have l = 3, and we need to introduce one auxiliary wave vector ${\vec{k}}_{32}^{\prime }$. According to (B.15), we then have

$$\begin{eqnarray}&&{\vec{k}}_{21}={\vec{L}}_{{q}_{2}}+{\vec{k}}_{32}^{\prime },\quad {\vec{k}}_{31}={\vec{L}}_{{q}_{3}}-{\vec{k}}_{32}^{\prime },\quad {\vec{k}}_{32}={\vec{k}}_{32}^{\prime },\end{eqnarray} \tag{ B.20 }$$

and ${\vec{L}}_{{q}_{1}}=-({\vec{L}}_{{q}_{2}}+{\vec{L}}_{{q}_{3}})$ because of the δ distribution in (B.3).

For second-order perturbation theory of a two-point spectrum, l = 4, and three auxiliary wave vectors ${\vec{k}}_{32}^{\prime }$, ${\vec{k}}_{42}^{\prime }$ and ${\vec{k}}_{43}^{\prime }$ need to be introduced. Then,

$$\begin{eqnarray}{\vec{k}}_{21} & = & {\vec{L}}_{{q}_{2}}+{\vec{k}}_{32}^{\prime }+{\vec{k}}_{42}^{\prime },\\ {\vec{k}}_{31} & = & {\vec{L}}_{{q}_{3}}-{\vec{k}}_{32}^{\prime }+{\vec{k}}_{43}^{\prime },\\ {\vec{k}}_{41} & = & {\vec{L}}_{{q}_{4}}-{\vec{k}}_{42}^{\prime }-{\vec{k}}_{43}^{\prime },\\ {\vec{k}}_{32} & = & {\vec{k}}_{32}^{\prime },\quad {\vec{k}}_{42}={\vec{k}}_{42}^{\prime },\quad {\vec{k}}_{43}={\vec{k}}_{43}^{\prime },\end{eqnarray} \tag{ B.21 }$$

and ${\vec{L}}_{{q}_{1}}=-({\vec{L}}_{{q}_{2}}+{\vec{L}}_{{q}_{3}}+{\vec{L}}_{{q}_{4}})$.

We now need to evaluate the integrals I_jk defined in (B.17), which are all of the type

$$\begin{eqnarray}&&{I}_{21}:= {\int }_{q}\,{{\rm{e}}}^{-{\vec{L}}_{{p}_{2}}^{\top }{C}_{{p}_{2}{p}_{1}}{\vec{L}}_{{p}_{1}}+{\rm{i}}{\vec{k}}_{21}\cdot \vec{q}},\end{eqnarray} \tag{ B.22 }$$

with an independent wave vector ${\vec{k}}_{21}$ representing any of the vectors ${\vec{k}}_{{jk}}$ defined in (B.15).

The momentum-correlation matrix ${C}_{{p}_{2}{p}_{1}}$ defined in (32) depends on the particle separation q only. Repeating (32), we can write ${C}_{{p}_{2}{p}_{1}}$ as

$$\begin{eqnarray}&&{C}_{{p}_{2}{p}_{1}}=-{\tilde{\pi }}_{\parallel }\,{\xi }_{\psi }^{\prime\prime }(q)-{\tilde{\pi }}_{\perp }\displaystyle \frac{{\xi }_{\psi }^{\prime }(q)}{q}\end{eqnarray} \tag{ B.23 }$$

with the projectors ${\tilde{\pi }}_{\parallel }$ and ${\tilde{\pi }}_{\perp }$ defined in (30).

We now expand the projectors ${\tilde{\pi }}_{\parallel ,\perp }$ with respect to $\hat{q}$ into the respective projectors ${\pi }_{\parallel ,\perp }$ with respect to ${\vec{k}}_{21}$. Let $\hat{k}$ be the unit vector in the direction of ${\vec{k}}_{21}$, then

$$\begin{eqnarray}&&{\pi }_{21}^{\parallel }=\hat{k}\otimes \hat{k},\quad {\pi }_{21}^{\perp }={{ \mathcal I }}_{3}-{\pi }_{21}^{\parallel }.\end{eqnarray} \tag{ B.24 }$$

For doing so, we expand the Hessian of the potential-correlation function into the projectors ${\pi }_{21}^{\parallel }$ and ${\pi }_{21}^{\perp }$,

$$\begin{eqnarray}&&{D}^{2}{\xi }_{\psi }(q)={\tilde{\pi }}_{\parallel }\,{\xi }_{\psi }^{\prime\prime }(q)+{\tilde{\pi }}_{\perp }\displaystyle \frac{{\xi }_{\psi }^{\prime }(q)}{q}={a}_{\parallel }{\pi }_{21}^{\parallel }+{a}_{\perp }{\pi }_{21}^{\perp },\end{eqnarray} \tag{ B.25 }$$

multiply this equation by ${\pi }_{21}^{\parallel }$ and ${\pi }_{21}^{\perp }$ and take the trace of the resulting two equations to find

$$\begin{eqnarray}{a}_{\parallel } & = & {\xi }_{\psi }^{\prime\prime }(q)\mathrm{tr}{\tilde{\pi }}_{\parallel }{\pi }_{21}^{\parallel }+\displaystyle \frac{{\xi }_{\psi }^{\prime }(q)}{q}\mathrm{tr}{\tilde{\pi }}_{\perp }{\pi }_{21}^{\parallel },\\ 2{a}_{\perp } & = & {\xi }_{\psi }^{\prime\prime }(q)\mathrm{tr}{\tilde{\pi }}_{\parallel }{\pi }_{21}^{\perp }+\displaystyle \frac{{\xi }_{\psi }^{\prime }(q)}{q}\mathrm{tr}{\tilde{\pi }}_{\perp }{\pi }_{21}^{\perp }.\end{eqnarray} \tag{ B.26 }$$

Introducing the cosine $\mu := \hat{q}\cdot \hat{k}$ of the angle between $\vec{q}$ and ${\vec{k}}_{21}$, the traces on the right-hand sides are

$$\begin{eqnarray}&&\mathrm{tr}{\tilde{\pi }}_{\parallel }{\pi }_{21}^{\parallel }={\mu }^{2},\quad \mathrm{tr}{\tilde{\pi }}_{\perp }{\pi }_{21}^{\parallel }=1-{\mu }^{2}=\mathrm{tr}{\tilde{\pi }}_{\parallel }{\pi }_{21}^{\perp },\quad \mathrm{tr}{\tilde{\pi }}_{\perp }{\pi }_{21}^{\perp }=1+{\mu }^{2}.\end{eqnarray} \tag{ B.27 }$$

We thus have

$$\begin{eqnarray}{a}_{\parallel } & = & {\mu }^{2}{\xi }_{\psi }^{\prime\prime }(q)+(1-{\mu }^{2})\displaystyle \frac{{\xi }_{\psi }^{\prime }(q)}{q},\\ 2{a}_{\perp } & = & (1-{\mu }^{2}){\xi }_{\psi }^{\prime\prime }(q)+(1+{\mu }^{2})\displaystyle \frac{{\xi }_{\psi }^{\prime }(q)}{q},\end{eqnarray} \tag{ B.28 }$$

and the quadratic form in (B.22) turns into

$$\begin{eqnarray}&&{\vec{L}}_{{p}_{2}}^{\top }{C}_{{p}_{2}{p}_{1}}{\vec{L}}_{{p}_{1}}=-{\vec{L}}_{{p}_{2}}^{\top }{\pi }_{21}^{\parallel }{\vec{L}}_{{p}_{1}}\,{a}_{\parallel }-{\vec{L}}_{{p}_{2}}^{\top }{\pi }_{21}^{\perp }{\vec{L}}_{{p}_{1}}\,{a}_{\perp }.\end{eqnarray} \tag{ B.29 }$$

It is now convenient to quantify the remaining projections by ${\lambda }_{21}^{\parallel ,\perp }$ defined by

$$\begin{eqnarray}&&{\vec{L}}_{{p}_{2}}^{\top }{\pi }_{21}^{\parallel }{\vec{L}}_{{p}_{1}}={g}_{{qp}}^{2}(\tau ,0){k}_{21}^{2}{\lambda }_{21}^{\parallel },\quad {\vec{L}}_{{p}_{2}}^{\top }{\pi }_{21}^{\parallel }{\vec{L}}_{{p}_{1}}={g}_{{qp}}^{2}(\tau ,0){k}_{21}^{2}{\lambda }_{21}^{\perp }.\end{eqnarray} \tag{ B.30 }$$

With these definitions, we arrive at the form

$$\begin{eqnarray}&&{I}_{21}={\int }_{q}\,{{\rm{e}}}^{{g}_{{qp}}^{2}(\tau ,0){k}_{21}^{2}({a}_{\parallel }{\lambda }_{12}^{\parallel }+{a}_{\perp }{\lambda }_{21}^{\perp })+{\rm{i}}{\vec{k}}_{21}\cdot \vec{q}}\end{eqnarray} \tag{ B.31 }$$

for the integral to be solved.

This integral has an intuitive physical meaning. To clarify it, we first split off a delta distribution,

$$\begin{eqnarray}&&{I}_{21}={{\rm{\Delta }}}_{21}+{{ \mathcal P }}_{21}\end{eqnarray} \tag{ B.32 }$$

with

$$\begin{eqnarray}&&{{ \mathcal P }}_{21}:= {\int }_{q}\{{{\rm{e}}}^{{g}_{{qp}}^{2}(\tau ,0){k}_{21}^{2}({a}_{\parallel }{\lambda }_{12}^{\parallel }+{a}_{\perp }{\lambda }_{21}^{\perp })}-1\}{{\rm{e}}}^{{\rm{i}}{\vec{k}}_{21}\cdot \vec{q}}\end{eqnarray} \tag{ B.33 }$$

and

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{21}:= {(2\pi )}^{3}{\delta }_{{\rm{D}}}({\vec{k}}_{21}).\end{eqnarray} \tag{ B.34 }$$

Due to the large dynamic range of the argument of the exponential in (B.33) and the fast oscillations of the Fourier phase, the integration required to evaluate ${{ \mathcal P }}_{21}$ is numerically difficult in some regions of parameter space. We use Levin collocation [23–25] for a fast and reliable integration scheme.

The meaning of ${{ \mathcal P }}_{21}$ becomes best visible in the limit of early times. Then, the argument of the exponential in (B.33) is small, the exponential can be approximated by a first-order Taylor expansion, and we are left with

$$\begin{eqnarray}&&{{ \mathcal P }}_{21}\approx {g}_{{qp}}^{2}(\tau ,0){k}_{21}^{2}{\int }_{q}({a}_{\parallel }{\lambda }_{12}^{\parallel }+{a}_{\perp }{\lambda }_{21}^{\perp }){{\rm{e}}}^{{\rm{i}}{\vec{k}}_{21}\cdot \vec{q}}.\end{eqnarray} \tag{ B.35 }$$

To evaluate this expression, we define the integrals

$$\begin{eqnarray}&&{I}_{n}^{\alpha }(k):= 2\pi {\int }_{0}^{\infty }{\rm{d}}q\,{q}^{n}{\xi }_{\psi }^{(n)}(q){\int }_{-1}^{1}{\rm{d}}\mu \,{\mu }^{\alpha }{{\rm{e}}}^{{\rm{i}}{kq}\mu },\end{eqnarray} \tag{ B.36 }$$

where ${\xi }_{\psi }^{(n)}(q)$ is the nth derivative of the correlation function ${\xi }_{\psi }(q)$. Using

$$\begin{eqnarray}&&{\int }_{-1}^{1}{\rm{d}}\mu \,{\mu }^{\alpha }{{\rm{e}}}^{{\rm{i}}{kq}\mu }=\displaystyle \frac{1}{{({\rm{i}}k)}^{\alpha }}\,{\partial }_{q}^{\alpha }\,{\int }_{-1}^{1}{\rm{d}}\mu \,{{\rm{e}}}^{{\rm{i}}{kq}\mu }=\displaystyle \frac{2}{{({\rm{i}}k)}^{\alpha }}\,{\partial }_{q}^{\alpha }\,{{j}}_{0}({kq})\end{eqnarray} \tag{ B.37 }$$

and the recurrence relations for the derivatives of the spherical Bessel functions, we find

$$\begin{eqnarray}{I}_{1}^{0}(k) & = & 4\pi {\displaystyle \int }_{0}^{\infty }{\rm{d}}q\,{\xi }_{\psi }(q)({kq}{{j}}_{1}({kq})-{{j}}_{0}({kq})),\\ {I}_{2}^{0}(k) & = & 4\pi {\displaystyle \int }_{0}^{\infty }{\rm{d}}q\,{\xi }_{\psi }(q)[(2-{k}^{2}{q}^{2}){{j}}_{0}({kq})-2{kq}{{j}}_{1}({kq})],\\ {I}_{1}^{2}(k) & = & 4\pi {\displaystyle \int }_{0}^{\infty }{\rm{d}}q\,{\xi }_{\psi }(q)\left[{{j}}_{0}({kq})+{kq}\left(1-\displaystyle \frac{4}{{k}^{2}{q}^{2}}\right){{j}}_{1}({kq})\right],\\ {I}_{2}^{2}(k) & = & 4\pi {\displaystyle \int }_{0}^{\infty }{\rm{d}}q\,{\xi }_{\psi }(q)\left[(2-{k}^{2}{q}^{2}){{j}}_{0}({kq})-\displaystyle \frac{4}{{kq}}{{j}}_{1}({kq})\right].\end{eqnarray} \tag{ B.38 }$$

Using these results with ${a}_{\parallel }$ and ${a}_{\perp }$ from (B.28), we find immediately

$$\begin{eqnarray}&&{\int }_{q}{a}_{\parallel }{{\rm{e}}}^{{\rm{i}}\vec{k}\cdot \vec{q}}={I}_{2}^{2}(k)+{I}_{1}^{0}(k)-{I}_{1}^{2}(k)=-{k}^{2}{P}_{\psi }(k)\end{eqnarray} \tag{ B.39 }$$

and

$$\begin{eqnarray}&&2{\int }_{q}{a}_{\perp }{{\rm{e}}}^{{\rm{i}}\vec{k}\cdot \vec{q}}={I}_{2}^{0}(k)-{I}_{2}^{2}(k)+{I}_{1}^{0}(k)+{I}_{1}^{2}(k)=0.\end{eqnarray} \tag{ B.40 }$$

Therefore, at early times,

$$\begin{eqnarray}&&{{ \mathcal P }}_{21}\approx -{g}_{{qp}}^{2}(\tau ,0){\lambda }_{21}^{\parallel }\,{k}_{21}^{4}\,{P}_{\psi }({k}_{21})=-{g}_{{qp}}^{2}(\tau ,0){\lambda }_{21}^{\parallel }\,{P}_{\delta }({k}_{21}).\end{eqnarray} \tag{ B.41 }$$

This is the density-fluctuation power spectrum, linearly evolved with ${g}_{{qp}}^{2}(\tau ,0)$, times $-{\lambda }_{21}^{\parallel }$. As we shall see in the main text, ${\lambda }_{21}^{\parallel }=-1$ in the most straightforward applications. This shows that the function ${{ \mathcal P }}_{21}$ generalises the linearly evolved density-fluctuation power spectrum and thus that ${{ \mathcal P }}_{21}$ is a nonlinearly time-evolving density-fluctuation power spectrum. Its Fourier transform is the corresponding generalisation of the spatial density correlation function.