Cosmic Structure Formation with Kinetic Field Theory

TL;DR

Kinetic Field Theory (KFT) offers a comprehensive statistical framework for modeling cosmic structure formation, incorporating non-linear dynamics, initial conditions, and modified gravity, with applications to both cosmology and atomic physics.

Contribution

This paper extends KFT to cosmological contexts, develops new analytical methods for non-linear structures, and applies the theory to fluids, dark matter, modified gravity, and atomic systems.

Findings

Derived an analytic, parameter-free equation for the non-linear cosmic power spectrum.

Developed a resummation scheme for fluid and dark matter mixtures.

Applied KFT to cold Rydberg atoms, analyzing their spatial correlation functions.

Abstract

Kinetic Field Theory (KFT) is a statistical field theory for an ensemble of point-like classical particles in or out of equilibrium. We review its application to cosmological structure formation. Beginning with the construction of the generating functional of the theory, we describe in detail how the theory needs to be adapted to reflect the expanding spatial background and the homogeneous and isotropic, correlated initial conditions for cosmic structures. Based on the generating functional, we develop three main approaches to non-linear, late-time cosmic structures, which rest either on the Taylor expansion of an interaction operator, suitable averaging procedures for the interaction term, or a resummation of perturbation terms. We show how an analytic, parameter-free equation for the non-linear cosmic power spectrum can be derived. We explain how the theory can be used to derive the density profile of gravitationally bound structures and use it to derive power spectra of cosmic velocity densities. We further clarify how KFT relates to the BBGKY hierarchy. We then proceed to apply kinetic field theory to fluids, introduce a reformulation of KFT in terms of macroscopic quantities which leads to a resummation scheme, and use this to describe mixtures of gas and dark matter. We discuss how KFT can be applied to study cosmic structure formation with modified theories of gravity. As an example for an application to a non-cosmological particle ensemble, we show results on the spatial correlation function of cold Rydberg atoms derived from KFT.

Figures

Equations

$$\delta (a)={\delta }_0{D}_{+}(a),$$

$$Z=\int \mathcal{D}\phi P(\phi ).$$

$$\mathbf{x}=x_j\otimes e_j,$$

$$⟨\mathbf{x},\mathbf{y}⟩=\left(x_i\cdot y_j\right)\left(e_i\cdot e_j\right)=x_j\cdot y_j.$$

$$Z=\int \mathrm{d}\mathbf{x}P(\mathbf{x}).$$

$$Z\to Z[\mathbf{J}]=\int \mathrm{d}\mathbf{x}P(\mathbf{x})\exp \left\{\mathrm{i}{\int }_0^{\infty }\mathrm{d}{t}^{\prime }⟨\mathbf{J}(t^{\prime }),\mathbf{x}(t^{\prime })⟩\right\},$$

$$⟨x_j(t)⟩=-\mathrm{i}{\frac{\delta }{\delta J_j(t)}Z[\mathbf{J}]|}_{\mathbf{J}=0}.$$

$$\mathbf{J}=J_j\otimes e_j=\left(\begin{array}{c}{J}_{q_j}\\ J_{p_j}\end{array}\right)\otimes e_j.$$

$$P(\mathbf{x})=\int \mathrm{d}{\mathbf{x}}^{(\mathrm{i})}P\left(\mathbf{x}|{\mathbf{x}}^{(\mathrm{i})}\right)P\left({\mathbf{x}}^{(\mathrm{i})}\right).$$

$$P\left(\mathbf{x}|{\mathbf{x}}^{(\mathrm{i})}\right)={\delta }_{\mathrm{D}}\left[\mathbf{x}-{\Phi }_{\mathrm{cl}}\left({\mathbf{x}}^{(\mathrm{i})}\right)\right],$$

$$E\left(\mathbf{x},{\mathbf{x}}^{(\mathrm{i})}\right)=0,$$

$$P\left(\mathbf{x}|{\mathbf{x}}^{(\mathrm{i})}\right)={\delta }_{\mathrm{D}}\left[E\left(\mathbf{x},{\mathbf{x}}^{(\mathrm{i})}\right)\right]$$

$$\dot{x}-\mathcal{J}{\nabla }_{x}H=0,$$

$$\mathcal{J}:=\left(\begin{array}{cc}0& {\mathcal{I}}_3\\ -{\mathcal{I}}_3& 0\end{array}\right)$$

$$\dot{x}-\mathcal{J}{\nabla }_{x}H=K,$$

$$\bar{x}(t)=G(t,0)x^{(\mathrm{i})}+{\int }_0^t\mathrm{d}{t}^{\prime }G(t,t^{\prime })K(t^{\prime }).$$

$$K(t^{\prime })=\left(\begin{array}{c}0\\ -\nabla V(t^{\prime })\end{array}\right).$$

$$G(t,t^{\prime }):=\left(\begin{array}{cc}{g}_{qq}(t,t^{\prime }){\mathcal{I}}_3& g_{qp}(t,t^{\prime }){\mathcal{I}}_3\\ 0& g_{pp}(t,t^{\prime }){\mathcal{I}}_3\end{array}\right).$$

$$E\left(\mathbf{x},{\mathbf{x}}^{(\mathrm{i})}\right)=\mathbf{x}(t)-\bar{\mathbf{x}}(t)=0.$$

$$Z[\mathbf{J},\mathbf{K}]=\int \mathrm{d}\Gamma \exp \left\{\mathrm{i}{\int }_0^{\infty }\mathrm{d}{t}^{\prime }⟨\mathbf{J}(t^{\prime }),\bar{\mathbf{x}}(t^{\prime })⟩\right\},$$

$$\mathrm{d}\Gamma :=\mathrm{d}{\mathbf{x}}^{(\mathrm{i})}P\left({\mathbf{x}}^{(\mathrm{i})}\right).$$

$$\rho (q,t_1)=\sum _{j=1}^N{\delta }_{\mathrm{D}}\left(q-q_j(t_1)\right).$$

$$\tilde{f}(k)=\int {\mathrm{d}}^{3}qf(q){\mathrm{e}}^{-\mathrm{i}kq},f(q)=\int \frac{{\mathrm{d}}^{3}k}{(2\pi {)}^3}\tilde{f}(k){\mathrm{e}}^{\mathrm{i}kq},$$

$$\tilde{\rho }(k_1,t_1)=:\tilde{\rho }(1)=\sum _{j=1}^N{\mathrm{e}}^{-\mathrm{i}{k}_1{q}_j(t_1)},$$

$${\int }_q:=\int {\mathrm{d}}^{3}q,{\int }_k:=\int \frac{{\mathrm{d}}^{3}k}{(2\pi {)}^3}.$$

$$\hat{\rho }(1)=\sum _{j=1}^N{\hat{\rho }}_j(1),$$

$${\hat{\rho }}_j(1):=\exp \left(-k_1\frac{\delta }{\delta J_{q_j}(t_1)}\right).$$

$$\mathbf{J}\to \mathbf{J}-\sum _{j=1}^r{\delta }_{\mathrm{D}}\left(t^{\prime }-t_j\right)k_j\cdot \left(\begin{array}{c}1\\ 0\end{array}\right)\otimes e_j.$$

$$Z[\mathbf{L}]=\int \mathrm{d}\Gamma {\mathrm{e}}^{\mathrm{i}⟨{\mathbf{L}}_q,\mathbf{q}⟩+\mathrm{i}⟨{\mathbf{L}}_p,\mathbf{p}⟩+\mathrm{i}{S}_{\mathrm{I}}}$$

$${\mathbf{L}}_q$$

Full text

We explain how the theory can be used to derive the density profile of gravitationally bound structures and use it to derive power spectra of cosmic velocity densities. We further clarify how KFT relates to the BBGKY hierarchy. We then proceed to apply kinetic field theory to fluids, introduce a reformulation of KFT in terms of macroscopic quantities which leads to a resummation scheme, and use this to describe mixtures of gas and dark matter. We discuss how KFT can be applied to study cosmic structure formation with modified theories of gravity. As an example for an application to a non-cosmological particle ensemble, we show results on the spatial correlation function of cold Rydberg atoms derived from KFT.

\shortabstract

1 Introduction

In our cosmic neighbourhood, we see ourselves surrounded by rich and pronounced, large-scale structures. Considering individual galaxies as the smallest constituents of these structures, we see large, almost empty regions, so-called voids, surrounded by filaments which intersect in knots. The typical size of the voids is of order $10\,h^{-1}\,\mathrm{Mpc}$222Astronomical units of length: $1\,\mathrm{Mpc}=3.1\cdot 10^{24}\,\mathrm{cm}$; $h$ is the Hubble constant in units of $100\,\mathrm{km\,s^{-1}\,Mpc^{-1}}$. The knots form the so-called galaxy clusters. In addition to this final, present-day state of cosmic structures, we can also observe seed structures in the very early universe. According to the well-established standard model of cosmology (see [1] for a review), the universe originated in a hot, dense state called the Big Bang. The electromagnetic heat radiation left over from this event was predicted in [2], observed in [3] and confirmed to have a near-perfect black-body spectrum with a temperature of $2.7\,\mathrm{K}$ in [4]. It forms the so-called cosmic microwave background (CMB), which was released when the cosmic plasma recombined, approximately 400,000 years after the Big Bang.

The CMB is almost perfectly isotropic. At closer inspection, it reveals temperature fluctuations with a relative amplitude of $\approx 10^{-5}$. We have good reasons to believe that these temperature fluctuations trace the cosmic structures present very shortly after the Big Bang. Notwithstanding the important question as to how these structures originated, we consider them as reflecting the initial state of the evolved cosmic structures we see today. Figure 1 shows both, the initial state as reflected by the CMB temperature fluctuations as observed by the Planck satellite[5], and the final state as shown by a galaxy survey conducted at $2\,\mu\mathrm{m}$ wavelength (the 2MASS survey, [6]).

It is for several reasons most important for cosmology to understand the evolution, the amplitude and the morphology of these structures. In the cosmological standard model, reconciling the amplitudes of the initial and the final state is possible only if we assume that by far most of the matter in the universe does not partake in the electromagnetic interaction, and is therefore called dark matter [7]. If this was not the case, the amplitude of the CMB temperature fluctuations would have to be two orders of magnitude larger. Within the limits of even the most precise measurements, the initial state can be considered as a Gaussian random field [8]. The formation of filamentary structures from a Gaussian random field has been shown to be a natural consequence of gravitational collapse of overdense regions against the cosmic expansion [9, 10].

It is quite straightforward to analyze the early stages of cosmic structure formation analytically as long as the relative density fluctuations are small, more precisely as long as the so-called density contrast $\delta$ (cf. Eq. 45) is less than unity. The standard procedure begins with the equations of ideal hydrodynamics and evaluates them on an expanding background to show that the density contrast grows in place proportional to the so-called linear growth factor $D_{+}$,

[TABLE]

with $a$ being the cosmic scale factor describing the expansion of the universe (see Sect. 3 below), and $\delta_{0}$ the density contrast at an arbitrary reference time when $D_{+}$ is set to unity. As long as the growth remains linear, structures do not change their size against the expanding background. A Fourier analysis in wave numbers co-moving with the mean cosmic expansion shows that the Fourier modes of the density-contrast field neither couple to each other nor change their wave number with time. Large-scale density-fluctuation modes are still linear today.

Despite formidable efforts and successes, it has turned out to be quite hard to follow cosmic structure formation into the non-linear regime with analytic means. See [11] for an excellent review on standard perturbation theory, [12, 13, 14, 15, 16, 17, 18, 19, 20, 21] as examples for several innovative approaches, [22, 23, 24, 25, 26, 27, 28, 29, 30] as examples for Lagrangian perturbation theory and [31, 32, 33] for recent developments towards an effective field theory of cosmic structure formation. One of the main reasons for this is that the fluid description of the predominantly dark matter must ultimately fail where and when the matter flow converges. When streams meet in a fluid, a discontinuity or a shock forms, keeping the velocity field of the flow unique. This is a consequence of the central assumption of (ideal) hydrodynamics that the mean-free path of the fluid particles is negligibly small. When dark-matter streams meet, however, multiple streams will form that simply cross each other, leading to a multi-valued velocity field. This notorious shell-crossing problem of cosmic structure formation is most likely the most important reason hampering progress in analytic theories of cosmic structure formation.

Of course, one can resort to fully numerical simulations. Corresponding algorithms have reached an impressive level of sophistication and have delivered overwhelmingly detailed results that, by and large, agree very well with observations. Yet, good reasons remain for trying to understand the formation of late-time, non-linear cosmic structures on an analytical basis. Conceptually the most important of these reasons is that numerical simulations strictly speaking do not explain the properties of cosmic structures, even though they succeed in reproducing them in astounding detail. Of particular importance are universal features of cosmic structures, such as the density profiles of gravitationally-bound objects. Analytic theories, on the other hand, could be expected to trace physical properties of non-linear structures back to their foundations in fundamental physics.

A second convincing reason in favour of an analytic theory is its flexibility against changing assumptions. Running sufficiently detailed and well-resolved numerical simulations is computationally expensive, and scanning wider parameter ranges is thus often forbiddingly time consuming. A third good reason is that shot noise is inevitable in numerical simulations, caused by the finite number of particles in any simulation volume. An analytic theory could be expected to allow taking the thermodynamic limit. Higher-order measures for the statistics of the cosmic structures, such as three- or four-point correlation functions which are indispensable to quantify non-linear and in particular non-Gaussian structures, could then be calculated without uncertainties due to shot noise and the finite number of samples that can be drawn from any numerical simulation.

Kinetic field theory is based on field-theoretical descriptions of classical particles [34, 35] and was developed mainly for studying glasses, fluctuation-dissipation theorems and the ergodic-non-ergodic transition [36, 37, 38, 39]. We have adapted it to cosmology and applied it to various aspects of cosmological structure formation [40, 41, 42, 43, 44], but also to completely different kinds of physical systems whose elementary constituents can be described as classical degrees of freedom. For cosmology, the KFT approach has several decisive advantages, the most important of which is that the motion of the classical particles is described in phase space by Hamilton’s equations. Since trajectories in phase space do not cross, this approach avoids the shell-crossing problem by construction. KFT incorporates the dynamics of the particles into a generating functional. Structurally, this generating functional resembles the generating functionals of statistical quantum field theories. One important simplification of the application to classical particles is due to the symplectic nature of Hamilton’s equations, which gives rise to Liouville’s theorem. Thus, by construction, KFT is built on a diffeomorphic map with unit functional determinant of an initial phase-space configuration to any later time.

We begin this review in Sect. 2 with the foundations of the theory and the construction of its central object, i.e. its generating functional. We specialize it to cosmology in Sect. 3 and summarize several further developments in Sect. 4, mainly offering different ways of taking particle interactions into account in approximate or averaged ways. Here, we derive a closed, analytic, non-perturbative, parameter-free equation for the non-linear density-fluctuation power spectrum which agrees very well with the results from numerical simulations up to wave numbers of $k\approx 10\,h^{-1}\,\mathrm{Mpc}$. In Sect. 5, we review how KFT can be used to explain the density profiles of dark-matter haloes. Section 6 presents the derivation of velocity power spectra. In Sect. 7, we establish the connection between KFT and the BBGKY hierarchy of kinetic theory, offering a closure condition for the hierarchy. We discuss in Sect. 8 how fluids can be incorporated into KFT. In Sect. 9, we give an exact re-formulation of KFT purely in terms of macroscopic fields and show how this leads to an efficient resummation of perturbation terms. The fluid description and the macroscopic re-formulation are then combined in Sect. 10 to describe mixtures of dark matter and gas with KFT. We outline in Sect. 11 how KFT can be extended to be applied within modified theories of gravity, and illustrate this outline with results for one specific example. Finally, in Sect. 12, we give an outlook on how to apply KFT to ensembles of cold Rydberg atoms as one example for a completely different physical system. We summarize this review in Sect. 13.

2 Kinetic Field Theory

2.1 Generating Functional

Kinetic field theory is a statistical field theory for classical particle ensembles in or out of equilibrium. Its central mathematical object is a generating functional $Z$. In close analogy to the partition sum of equilibrium thermodynamics, this generating functional integrates the probability distribution $P(\varphi)$ for system states $\varphi$ over the state space,

[TABLE]

It is generally a path integral if the states space is a function space of field configurations.

For classical (canonical) ensembles consisting of $N$ point particles, we can immediately specify the generating functional further. The state space is then the phase space $\Gamma$ of the $N$ particles whose trajectories $(q_{j},p_{j})=:x_{j}$ with $1\leq j\leq N$ are tuples of position $q_{j}$ and momentum $p_{j}$. To allow a compact notation, we bundle the phase-space trajectories $x_{j}$ into a tensorial object

[TABLE]

where summation over $j$ is implied and the vector $e_{j}$ has components $(e_{j})_{i}=\delta_{ij}$, $1\leq i\leq N$. For such tensors, we define the scalar product

[TABLE]

Point particles are described by Dirac delta distributions instead of smooth fields. The path integral in (2) then turns into an ordinary integral over the $N$-particle phase space,

[TABLE]

We introduce a generator field $\boldsymbol{J}$ conjugate to the phase-space trajectories $\boldsymbol{x}(t)$ into the generating functional $Z$,

[TABLE]

such that the functional derivative of $Z[\boldsymbol{J}]$ with respect to the component $J_{j}(t)$ of the generator field returns the average position $\langle x_{j}(t)\rangle$ of particle $j$ at time $t$,

[TABLE]

Like $x$, the generator field has components conjugate to position $q$ and momentum $p$,

[TABLE]

We further split the probability $P(\boldsymbol{x})$ for the state $\boldsymbol{x}$ to be occupied into a probability $P(\boldsymbol{x}^{\mathrm{(i)}})$ for the particle ensemble to occupy an initial state $\boldsymbol{x}^{\mathrm{(i)}}$ at time $t=0$, and the conditional probability $P(\boldsymbol{x}|\boldsymbol{x}^{\mathrm{(i)}})$ for the ensemble to move from there to the time-evolved state $\boldsymbol{x}$,

[TABLE]

For particles on classical trajectories, the transition probability $P(\boldsymbol{x}|\boldsymbol{x}^{\mathrm{(i)}})$ must be a functional delta distribution of the classical equation of motion,

[TABLE]

where $\Phi_{\mathrm{cl}}(\boldsymbol{x}^{\mathrm{(i)}})$ denotes the classical (Hamiltonian) flow on the phase space, beginning at the initial phase-space points $\boldsymbol{x}^{\mathrm{(i)}}$ of the particle ensemble. Writing the equation of motion in the form $E(\boldsymbol{x})=0$, the Hamiltonian flow consists of all solutions of this equation for initial points within a certain domain of phase space, as illustrated in Figure 2.

Denoting the solution of the classical equations of motion beginning at $\boldsymbol{x}^{\mathrm{(i)}}$ as

[TABLE]

we can write the transition probability as

[TABLE]

with a functional delta distribution singling out the deterministic trajectories solving the equations of motion.

For classical point particles, the equation of motion for a single particle is the Hamiltonian equation,

[TABLE]

where $H$ is the Hamiltonian function on the phase space and

[TABLE]

is the usual symplectic matrix, with $\mathcal{I}_{n}$ denoting the unit matrix in $n$ dimensions. Let us first assume the particles to be non-interacting, in which case the Hamiltonian equations are linear and admit the construction of a Green’s function $G(t,t^{\prime})$. To prepare for the inclusion of interactions, we augment this free equation of motion by an inhomogeneity or source term $K$,

[TABLE]

and write its solution $\bar{x}(t)$ beginning at $x^{\mathrm{(i)}}$ in terms of the Green’s function as

[TABLE]

We shall henceforth assume that the inhomogeneity is caused by an interaction with the potential $V$. Then, we can write

[TABLE]

The Green’s function $G$ itself is a $6\times 6$ matrix which we write as

[TABLE]

Its component functions $g_{qq}(t,t^{\prime})$, $g_{qp}(t,t^{\prime})$, and $g_{pp}(t,t^{\prime})$ will be specified later in a cosmological context. For simplicity, we abbreviate $g_{qp}(t,0)=:g_{qp}(t)$.

For the entire particle ensemble, we introduce the $N$-particle Green’s function $\boldsymbol{G}(t,t^{\prime}):=G(t,t^{\prime})\otimes\mathcal{I}_{N}$ and the $N$-component source field $\boldsymbol{K}:=K_{j}\otimes e_{j}$. Then, the solutions of the equations of motion for the particle ensemble are

[TABLE]

Inserting this expression into (12), the result into (9) and the probability distribution $P(\boldsymbol{x}^{\mathrm{(i)}})$ into the generating functional $Z[\boldsymbol{J}]$ from (6) gives

[TABLE]

where we have introduced the initial phase-space measure

[TABLE]

2.2 Density Operator

The standard application of KFT in this review will concern the calculation of density power spectra. We shall thus proceed to define a density operator and study its action on the generating functional (20). The number density of the particles at time $t_{1}$ is

[TABLE]

In a Fourier representation, defined by the convention

[TABLE]

the density mode with the wave number $k_{1}$ is

[TABLE]

where we have introduced the conventional short-hand notation $(k_{j},t_{j})=:(j)$. We shall abbreviate the integral expressions in (23) as

[TABLE]

Replacing the particle position $q_{j}(t_{1})$ by a functional derivative with respect to the generator-field component $J_{q_{j}}(t_{1})$, we find the density operator

[TABLE]

composed of the one-particle density operators

[TABLE]

For indistinguishable particles, which we shall henceforth assume, we can set the particle index $j$ to an arbitrary value without loss of generality. Since the operator (27) is an exponential of a derivative with respect to the generator field, it corresponds to a finite translation of the generator field. After applying $r\geq 1$ of these operators, the generator field is translated by

[TABLE]

Setting the generator field $\boldsymbol{J}$ to zero after application of the density operator turns the generating functional into

[TABLE]

with the components

[TABLE]

of the translations and the interaction term

[TABLE]

In (29), $\boldsymbol{q}$ and $\boldsymbol{p}$ in the exponent are the initial particle positions and momenta,

[TABLE]

The generating functional (20) and the result (29) after applying $r$ one-particle density operators and setting the source field $\boldsymbol{J}$ to zero are fully general in the sense that neither the initial particle distribution in phase space, nor the Green’s function, nor the interaction potential have been specified yet. So far, we have just assumed that a Green’s function and an interaction potential exist.

2.3 Interaction Operator

We suppose that the potential $V$ is a linear superposition of contributions by individual particles,

[TABLE]

where we have identified the density $\rho$ from (22). The potential gradient at position $q_{i}$ is

[TABLE]

where we have introduced the response field

[TABLE]

after partial integration. With (33), we can write this result as

[TABLE]

Applying the convolution theorem, this expression turns into

[TABLE]

where $(-1^{\prime})=:(-k_{1}^{\prime},t_{1}^{\prime})$. The Fourier transform of the density $\rho$ is expressed by the operator $\hat{\rho}$ from (26), while we can assign the operator

[TABLE]

to the response field.

With the definition $\mathrm{d}1^{\prime}:=\mathrm{d}t_{1}^{\prime}\mathrm{d}^{3}k_{1}^{\prime}/(2\pi)^{3}$, this allows us to write the interaction term $S_{\mathrm{I}}$ as the interaction operator

[TABLE]

acting on the free generating functional,

[TABLE]

with

[TABLE]

One approach to perturbation theory can now begin with a Taylor expansion of exponential factor in terms of the interaction operator. Anticipating the specialisation of KFT to cosmology detailed in Sect. 3, Fig. 3 shows cosmological results obtained in first-order perturbation theory [42]. A summary of this calculation is presented in Appendix A.

3 Specialisation to Cosmology

We shall now specify the theory to its cosmological application. Where not stated otherwise, we use the cosmological parameters $\Omega_{\mathrm{m}}=0.3$, $\Omega_{\Lambda}=0.7$, $\Omega_{\mathrm{b}}=0.04$ for the density parameters of matter, the cosmological constant, and baryonic matter, respectively; further $h=0.7$ for the dimension-less Hubble constant and $\sigma_{8}=0.8$ for the normalisation of the power spectrum. Moreover, we choose the form of the cold-dark matter transfer function provided by [10] to specify the initial power spectrum.

3.1 Equation of Motion for Point Particles

The cosmological standard model is built upon the theory of general relativity and two symmetry assumptions. They assert that there exists a mean flow in the Universe with respect to which all sufficiently averaged observable properties appear spatially isotropic and homogeneous. This mean flow can then only either expand or shrink isotropically and must thus be characterised by a scale factor $a(t)$, depending only on cosmic time $t$. The dynamics of $a(t)$ is determined by Einstein’s field equations which, under the symmetry assumptions made, simplify to Friedmann’s equations. The relative cosmic expansion rate is the Hubble function $H(a)=\dot{a}/a$. Usually, $a$ is normalised to unity today, but it is more convenient for our purposes to normalise $a$ to unity at the initial time $t=0$ such that $a(0)=:a_{\mathrm{i}}=1$.

The existence of the mean flow suggests to introduce coordinates $q$ comoving with the flow, defined in terms of the physical coordinates $r$ via $r=aq$. The Lagrange function of point particles of mass $m$, expressed in terms of the comoving coordinates, is

[TABLE]

(see, e.g. [50]), where $\Phi$ is the gravitational potential sourced by the fluctuations of the matter density $\rho$ around its time-dependent mean $\bar{\rho}$ via the Poisson equation

[TABLE]

The mean cosmic density after the radiation-dominated epoch is

[TABLE]

where $H_{\mathrm{i}}$ and $\Omega_{\mathrm{mi}}$ are the Hubble function and the cosmic matter-density parameter at the initial time. If we set this time early in the matter-dominated epoch, as we intend to do throughout, we can safely set $\Omega_{\mathrm{mi}}=1$ and thus bring the Poisson equation (43) into the form

[TABLE]

with the density fluctuations now being expressed by the dimension-less density contrast $\delta$.

Linear perturbation theory shows that the density fluctuations $\delta$ grow proportional to the so-called linear growth factor $D_{+}$. It is convenient to replace the cosmic time $t$ by the growth factor,

[TABLE]

where the time $t=0$ is now set to some initial time which is late enough for matter to dominate over radiation, but early enough for density fluctuations to be very small. A suitable choice is the time when the cosmic microwave background was released. In this new time coordinate, the Lagrange function reads

[TABLE]

(see e.g. [40]), with the potential $\bar{\Phi}$ satisfying the Poisson equation

[TABLE]

Here, $g(t)$ is defined to be the function

[TABLE]

normalised to unity at $t=0$ or $a=a_{\mathrm{i}}$. Notice that the potential $\bar{\Phi}$ defined in (48) now has the dimension of a squared length. The equation of motion reads

[TABLE]

with the potential $\varphi$ now obeying the Poisson equation

[TABLE]

The usual Legendre transform of the Lagrange function (47) leads to the Hamiltonian

[TABLE]

and thus to the Newtonian equation of motion

[TABLE]

which is solved by the Green’s function

[TABLE]

It is important to note that $\bar{g}_{qp}(t,t^{\prime})$ is limited from above because of the cosmic expansion. Consequently, the inertial trajectories described by this Green’s function deviate strongly from the true, fully interacting trajectories. It will thus be advantageous to find a replacement for $\bar{g}_{qp}$ that already captures part of the gravitational interaction. An example for such a replacement is given by the Zel’dovich approximation [51], but we prefer a slightly more general approach here.

3.2 Particle Trajectories and Effective Force

We wish to solve the equation of motion (50) with an expression of the form

[TABLE]

such that the propagator $g_{qp}(t,t^{\prime})$ can play the role of the $q$-$p$ component of the Green’s function $G(t,t^{\prime})$. Taking the first two time derivatives of (55) gives

[TABLE]

For (56) to agree with the equation of motion (50), the effective force $f(t)$ has to be appropriately adapted once $g_{qp}$ has been chosen. A frequent and convenient choice in cosmology is the Zel’dovich approximation [51], which describes particle trajectories as inertial in the time coordinate $t=D_{+}$,

[TABLE]

Inserting this particular choice into (56) and comparing to (50) immediately results in

[TABLE]

This is thus the effective force term in the Zel’dovich approximation. It is chosen such that the total initial force, i.e. the sum of the potential gradient and the velocity-dependent contribution, vanishes initially, corresponding to the initial inertial motion of the particles.

Another choice for the propagator $g_{qp}$ begins with defining the effective force term such that only its velocity-dependent contribution is initially absent,

[TABLE]

The remaining initial force then causes the particles to lag behind the inertial motion. This choice has been suggested to avoid part of the re-expansion in the Zel’dovich approximation of cosmic structures after their formation [40]. Relative to this effective force, the equation of motion (50) reads

[TABLE]

which implies the homogeneous solution $\dot{q}(t)=\dot{q}(0)\mathrm{e}^{h}$ and thus the propagator

[TABLE]

At late times, $h\to-1$, and the propagator (61) approaches $\mathrm{e}^{-1}$ times the Zel’dovich propagator.

Since the peculiar velocity $\dot{q}$ itself depends on the time-integrated force $f(t)$, (58) and (59) are integral equations for the force. After transforming them to differential equations, they can be solved by variation of constants. The solution of (59), which we will use here together with the propagator (61), turns out to be

[TABLE]

We shall use this expression later as a starting point for an averaging and an approximation scheme.

3.3 Initial Conditions

Supported by observations of the temperature fluctuations in the cosmic microwave background, we can assume that the density fluctuations in the early universe can be characterised as a Gaussian random field (for recent empirical support of this assumption, see [8]). By the Helmholtz theorem, the peculiar velocity field can be decomposed into a curl and a gradient. By the cosmic expansion and angular-momentum conservation, the curl component will quickly decay, so we can model the initial peculiar velocity as the gradient of a velocity potential $\Psi$. Then, by the continuity equation, the density contrast $\delta$ has to satisfy the Poisson equation $\delta=-\nabla^{2}\Psi$. Thus, given the velocity potential, both density contrast and peculiar velocity will be determined. The velocity potential also has to be a Gaussian random field. If the density contrast is statistically characterised by its power spectrum $P_{\delta}(k)$, the velocity potential has the power spectrum

[TABLE]

according to the Poisson equation. With the peculiar velocity and the density contrast both being Gaussian random fields derived from the velocity potential, either of the power spectra $P_{\Psi}(k)$ or $P_{\delta}(k)$ completely determines their statistical properties. Note that we assume irrotationality of the velocity field and thus absence of vorticity at the initial time only. During the subsequent nonlinear evolution of the particles vorticity can and will nevertheless be created.

To specify the generating functional of KFT, we need the probability distribution $P(\boldsymbol{x}^{\mathrm{(i)}})$ for the initial phase-space coordinates of the point particles of our ensemble.

Drawing particle positions randomly by a Poisson process with a probability proportional to the density contrast $\delta=-\nabla^{2}\Psi$, and assigning momenta to these particles proportional to $\nabla\Psi$, the probability distribution $P(\boldsymbol{x}^{\mathrm{(i)}})=P(\boldsymbol{q},\boldsymbol{p})$ turns out to be

[TABLE]

where $C_{pp}$ is the momentum-correlation matrix depending on the particle positions [42]. For late-time cosmological applications using the unbound Zel’dovich propagator (57) or its improvement (61), the factor $\mathcal{C}(\boldsymbol{p})$, which is a polynomial in the momenta, can safely be set to unity. We shall see in Sect. 9 how $\mathcal{C}(\boldsymbol{p})$ can be fully taken into account.

The momentum-correlation matrix is given by

[TABLE]

where $j\neq k$ and $E_{jk}=e_{j}\otimes e_{k}$. The variance $\sigma_{1}^{2}$ is defined by

[TABLE]

The first term on the right-hand side of (65) is the momentum dispersion caused by the momenta being drawn from a Gaussian random field. The second term on the right-hand side is the correlation between the momenta of different particles, separated by $q_{jk}:=|q_{j}-q_{k}|$. Since we assume that the initial density fluctuations are a statistically isotropic and homogeneous random field, the momentum correlations can only depend on the absolute value of the relative particle separation.

4 Further Developments

After these preparations, we return to some further developments simplifying the evaluation of the free generating functional (41) and the interaction term (31).

4.1 Factorisation of the Free Generating Functional

With the Gaussian initial distribution (64) of particle positions and momenta, the momenta can immediately be integrated out in the free generating functional (41), leading to

[TABLE]

With the momentum-correlation matrix depending only on the relative particle separations, this remaining expression can be fully factorised. In the simplest case of a two-point function, leaving out an irrelevant delta distribution, the result is

[TABLE]

where the non-linearly evolved power spectrum $\mathcal{P}(k_{1},t)$ is given by

[TABLE]

with the expressions

[TABLE]

appearing in the exponentials. The quadratic form $Q_{0}$ independent of $q$ appearing in the exponential damping term derives from the momentum dispersion.

The function $a_{\parallel}(q)$ appearing in this expression is the correlation function of those momentum components parallel to the line connecting the positions of the two momenta being correlated. In terms of the power spectrum $P_{\delta}(k)$ of the density fluctuations, it is given by

[TABLE]

with

[TABLE]

Here, the spherical Bessel functions $j_{1,2}(kq)$ appear, and $\mu$ is the cosine of the angle enclosed by $k_{1}$ and $q$. Notice that, since $j_{1}(x)\approx x/3$ and $j_{2}(x)\approx x^{2}/5$ for $x\ll 1$, the functions $a_{1,2}(q)$ and $a_{\parallel}(q)$ have the limits

[TABLE]

for $q\to 0$, and thus $Q\to Q_{0}$ in the same limit.

If the exponential in (69) can be approximated by its Taylor expansion to first order, the power spectrum $\mathcal{P}$ evolves linearly,

[TABLE]

4.2 Averaged Interaction Term

We now return to the interaction term (31) and evaluate it for a two-point function at equal times, $t_{1}=t_{2}=t$, taking the constraint $k_{1}+k_{2}=0$ into account which is enforced by the delta distribution in (68). We then have

[TABLE]

where $f_{i}$ is the force (62) acting on particle $i$. Dropping the time argument for brevity, we bring the force terms into the form

[TABLE]

where $f_{ij}$ is the force on particle $i$ due to particle $j$, and we have used $f_{12}-f_{21}=2f_{12}$ by Newton’s third law. If we can neglect three-point correlations for now, the second term on the right-hand side of (76) can be neglected in an isotropic random field since the forces exerted by particles $j$ with $j\geq 3$ on particles $1$ and $2$ will vanish on average because there is no preferred direction they could point to. We can then simplify the interaction term to

[TABLE]

containing the projection of the force $f_{12}$ between particles $1$ and $2$ on the wave vector $k_{1}$.

The Fourier transform of the potential $v$ of a unit point mass satisfying the Poisson equation (51) is

[TABLE]

thus the gradient of the potential of particle $2$ at the position of particle $1$ is

[TABLE]

A seemingly radical approximation of the interaction term (77) consists in replacing the potential gradient $\nabla_{1}\varphi_{2}$ between the particles $1$ and $2$ by a suitable average. Simply averaging $\nabla_{1}\varphi_{2}$ over all particle pairs would return zero because of the statistical isotropy of the particle distribution. It is thus important to realise that the interaction term (77) contains the projection of the potential gradient on the wave vector $k_{1}$ of the density mode considered. We thus need to calculate the projection $k_{1}\cdot\nabla_{1}\varphi_{2}$ in a suitable average which does not vanish for particles correlated with the density mode $k_{1}$.

The probability $P_{12}$ for finding particles $1$ and $2$ at a separation $q_{12}:=|q_{1}-q_{2}|$ has a spatially independent mean contribution $\bar{\rho}$ and a conditional contribution due to the particle correlations, expressed by the spatial correlation function $\xi(q_{12})$ of the particles,

[TABLE]

For uncorrelated particles in a homogeneous random field, the direction of $\nabla_{1}\varphi_{2}$ is random with respect to $k_{1}$, thus its contribution to the average must vanish, while the correlated part remains. We thus weigh the potential gradient $\nabla_{1}\varphi_{2}$ with the conditional probability $\xi(q_{12})$ and take the Fourier transform of the result to obtain the average Fourier component of the potential gradient at wave number $k_{1}$,

[TABLE]

By means of the Fourier convolution theorem, we can then write the Fourier transform of the averaged potential gradient as a convolution of the Fourier transforms of the potential gradient itself with that of the correlation function, i.e. with the power spectrum. For the latter, we take the linearly time-evolved initial density-fluctuation power spectrum $P_{\delta}(k)$, damped on the free-streaming scale. Thus, we insert

[TABLE]

for the Fourier transform of the correlation function and obtain

[TABLE]

The angular integral remaining in (83) can be carried out and results in

[TABLE]

where $\hat{k}_{1}$ is the unit vector in direction of $k_{1}$, and $y:=\kappa/k_{1}$. We truncate the integration at $y=1$ ($\kappa=k_{1}$) to suppress modes on scales smaller than the density fluctuation considered. The function $J(y)$ is

[TABLE]

Inserting the averaged potential gradient (83) into the force term (62), we obtain a scale-dependent, mean force term $\langle f_{12}\rangle(k_{1},t)$ weighed by the correlation function between correlated particle pairs, which defines the mean interaction term via (77),

[TABLE]

which by construction does not depend on the particle positions or momenta any more.

Replacing $S_{\mathrm{I}}$ by this average $\langle S_{\mathrm{I}}\rangle$, we can thus pull the interaction term in front of the integral (29) and use the form (68) of the free generating functional to obtain the closed expression

[TABLE]

for the non-linearly evolved, density-fluctuation power spectrum. Note that this expression is non-perturbative and parameter-free. Comparing the perturbative first-order result shown in Figure 3 with the power spectrum obtained from (87) shown in Figure 4 clearly demonstrates the significant improvement.

The non-linear power spectrum in the mean-field approximation shown in this Figure falls below the numerically simulated spectrum for $k\gtrsim 10\,h\,\mathrm{Mpc}^{-1}$. This does not indicate a conceptual breakdown of KFT, but rather a limitation of the mean-field approximation of the force term. The scale of $\approx 10\,h\,\mathrm{Mpc}^{-1}$ at $z=0$ is thus not of fundamental importance, but depends on the normalisation of the power spectrum.

4.3 Interaction in the Born Approximation

Let us now return to the force term (58) with the potential gradient expressed by its Fourier transform (79). If we wish to avoid averaging over particle positions, as we did before, the essential difficulty is that the actual particle positions $q_{1,2}(t)$ appear in the Fourier phase. This difficulty can be by-passed by approximating the particle trajectories by their unperturbed or inertial trajectories,

[TABLE]

where $q_{j}$ and $p_{j}$ without a time argument are meant to indicate initial particle positions and momenta. We abbreviate $q_{ij}=q_{i}-q_{j}$, $p_{ij}=p_{i}-p_{j}$ and introduce $K=g_{qp}(t)k$. Exchanging the order of integration over $k$ and $t$, we find

[TABLE]

The second term in brackets in (89),

[TABLE]

with $\kappa:=k\cdot p_{ij}$, is easily evaluated by numerically integrating its real and imaginary parts. The potential of the force term $f_{ij}$ from (89) thus has the effective, $k$-dependent amplitude

[TABLE]

whose real and imaginary parts we show in Figs. 5 and 6. Figure 5 shows the dependence on the scale factor for three different choices of $\kappa$, while Fig. 6 displays the dependence on $\kappa$ for $a=1$.

These Figures show several interesting properties of the effective potential amplitude. First, for large values of $\kappa$, i.e. for small scales or large initial momenta, the real part of $V_{\mathrm{eff}}$ begins oscillating at late times. This indicates that the displacement of the particles becomes comparable or larger than the structures they belong to, thus reducing the effective potential amplitude at late times: structures on such scales are then not built up any more, but particles oscillate within them. Alternatively, these oscillations can be interpreted as the effect of an oscillating effective propagator $g_{qp}$, which mimicks the motion of particles in gravitationally-bound structures.

Second, the imaginary part of the potential amplitude is typically small at early times and for small values of $\kappa$, but becomes large at late times for large $\kappa$. This indicates damping by the transport of particles and the associated mixing of phases of particle trajectories. This damping justifies the exponential cut-off in the power spectrum (82) which we have used for calculating the mean force term before.

It remains to be seen in comparison with numerical simulations how useful this variant of the Born approximation will be for cosmological structure formation. The agreement between numerical simulations and the non-linear power spectrum (87) including the averaged interaction, however, suggests that the Born-approximated interaction term may return similarly or even more accurate results.

5 Halo profiles in KFT

The formation of large-scale, dark-matter structures in our KFT formulation as well as in numerical simulations is governed by a few physical properties only: the interaction potential between dark matter particles, the equations of motion governing particle trajectories, the underlying cosmological background model, i.e. the expansion history of the model universe, and the initial density and momentum correlations set, as is most commonly believed, by quantum fluctuations during inflation.

It would appear that in numerically simulated as well as in observed, gravitationally-bound dark-matter structures, these physical properties or their interplay lead to a particular shape of the radial density profile [53, 54, 55] that is valid for halos on scales ranging from dwarf galaxies to massive galaxy clusters. Since the KFT formalism allows to easily change all the ingredients for large-scale structure formation, it is a natural playground for investigating the origin of the dark-matter, halo-density profile.

We found in [43] that taking into account initial momentum correlations leads to a characteristic deformation of the non-linear power spectrum on scales of the order of $k\approx 0.3\,h\,\mathrm{Mpc}^{-1}$ but not at wave numbers large enough for dark-matter halos to appear. Furthermore, we have shown that the full initial density correlations play a role for the shape of the density-fluctuation power spectrum at early times, but are washed out at late times when $a\to 1$, and therefore are unlikely responsible for the shape of density profiles of highly non-linear objects [44].

In fact, it is the Newtonian gravitational potential that is often charged with being responsible for the shape and universality property of dark matter halos. We have therefore used the KFT formalism to investigate whether the shape of the interaction potential is responsible for the shape of the non-linear power spectrum on small scales, $k\geq 1\,h\,\mathrm{Mpc}^{-1}$, where contributions to the non-linear power spectrum from inner structures of dark-matter halos begin to dominate.

According to the widely used halo model, the non-linear power spectrum is a convolution of Fourier-transformed Navarro-Frenk-White (NFW) profiles [53], weighed with the mass function (for an extensive review see [56]). Since the non-linear power spectrum is reproduced with KFT at least up to $k\approx 10\,h\,\mathrm{Mpc}^{-1}$, un-weighing with the mass function leads to the density profile of dark-matter halos.

The non-linear contributions to the power spectrum obtained from KFT approximately correspond to the one-halo term of the halo model. Thus, leaving the relative abundances of haloes with different mass unchanged, we can study with KFT how different choices for the gravitational potential affect the density-profile shape.

We point out that, while the halo model must be provided with the explicit form of the halo density profile, the KFT approach does not need any input of this kind. The non-linear density power spectrum and the halo density profile that is produced in the KFT approach are thus the result of particle dynamics including particle interactions and correlations.

To account for the fact that a general power-law potential leads to a linear growth factor that depends on the wave number $k$, the Poisson equation has to be modified accordingly. This is achieved through an algebraic, ad-hoc modification of Poisson’s equation in Fourier space,

[TABLE]

with $\tilde{\Phi}$ the gravitational potential, $G$ the gravitational constant and $\bar{\rho}$ the mean matter density. The modifying function $f(k)$ is the inverse of the Fourier-transformed particle interaction potential and is, in that sense, fixed. For the Newtonian gravitational interaction, it is simply $f(k)=k^{2}$.

A first analysis of the small-scale behaviour of the non-linear contributions to the power spectrum under different gravitational laws using the perturbative approach described in Sect. 2.3 yields the results shown in Fig. 7 for interaction potentials of the form

[TABLE]

where the parameter $\varepsilon$ is introduced for technical reasons in order to perform the Fourier transforms of the interaction potential for all $n$. It can be interpreted as a smoothing scale. We set $\varepsilon=10^{-4}$, which corresponds to scales much smaller than the scales we are currently interested in. We have verified that a small enough $\varepsilon$ does not influence the power spectrum on the much larger scales that we are considering.

Particle interactions are included up to first order in the interaction operator (39) in the same way as for the results shown in Fig. 3 since the agreement with non-linear power spectra from Cosmic Emulator [46, 47, 48, 49] that uses state-of-the-art $N$-body simulations is already very good.

The initial correlations appear in a quadratic form in an exponential in the KFT formalism (cf. 64) and must therefore be Taylor-expanded before the integration over initial particle positions and momenta can be performed in (40). Initial correlations are included up to second order here. Including the full hierarchy of initial momentum correlations, as shown in Sect. 4.1, proves to be much more cumbersome in the perturbative approach compared to the non-perturbative schemes of Sects. 3.2 and 3.3. Including the full hierarchy of initial momentum correlations in this context would not alter the results shown here, as they only affect the shape of the power spectrum on intermediate scales ($k\approx 0.3\,h\,\mathrm{Mpc^{-1}}$) as shown in [43].

The amplitudes of the interaction potential (93) and of the initial density-fluctuation power spectrum are chosen such that the amplitudes of the non-linear contributions to the power spectrum for any potential match that of the Newtonian case, and that the amplitude of linear growth on large scales, i.e. $k\ll 1\,h\,\mathrm{Mpc}^{-1}$, also matches the one observed today.

The linear power spectrum can be obtained, as usual, by multiplying the initial power spectrum $P^{\mathrm{(i)}}_{\delta}$ with the appropriate linear growth factor corresponding to the interaction potential being used,

[TABLE]

It is clear from this expression that the linear growth for non-Newtonian potentials will be different from the Newtonian case due to the wave-number dependence of the growth factor. Therefore, the linear growth together with the non-linear growth for non-Newtonian gravity will produce power spectra that look quite different from the Newtonian case. However, this is irrelevant for this analysis, because only those contributions to the density-fluctuation power spectrum which are due to the inner structure, i.e. density profile, of dark matter halos are of interest here.

Switching to the picture of the halo model, the argument becomes clearer: While the linear power spectrum, corresponding to the two-halo term, describes the correlations between two different dark matter halos, the one-halo term describes the contributions from the inner structures of halos, which are highly non-linear. Therefore, also in KFT the information on the density profile of dark matter halos will be contained exclusively in the non-linear contributions.

Figure 7 clearly shows that, for smaller scales, i.e. $k>1\,h\,\mathrm{Mpc}^{-1}$, the non-linear contributions are almost insensitive to the slope of the interaction potential. In the gray-shaded area in Fig. 7, there is some sensitivity to the potential slope, however these scales are too large to be relevant for the density profiles of even the largest dark matter halos.

Since the slopes of the non-linear contributions remain the same even for strongly varying interactions laws, it appears that different choices for the gravitational potential would not affect the density-profile shape of dark matter halos.

This result is surprising and should still be perceived with some caution. Although the non-linear power spectrum for Newtonian gravity obtained with the first-order perturbation term in KFT agrees well with predictions from numerical simulations as shown in Fig. 3, it is still uncertain whether the perturbation series converges regardless of the interaction potential. If the perturbation series does converge, for which there are good reasons, then higher-order terms should yield ever smaller corrections to the non-linear power spectrum.

However, for short-ranged interaction potentials, most of the structure is expected to be accumulated on very small scales, and higher-order interaction terms should become more dominant. Hence, the perturbative ansatz is expected to perform worse and may eventually break down. To test whether this will indeed be the case, higher-order terms and a check of the convergence of the perturbation series are needed.

We emphasise that the above results were interpreted assuming that the relative abundances of haloes with different mass remains unchanged. This assumption was made because we also assumed that only non-linear contributions corresponding to the one-halo term were considered. In actuality, this does not need to be the case. We could and should as well have contributions from the two-halo term in the non-linear part of the power spectrum. This poses a considerable problem for the analysis using the perturbative approach as well as for the same analysis employing the non-perturbative approaches presented in Sects. 3.2 and 3.3. To consider the fully non-linear density-fluctuation power spectrum which corresponds to the one-halo and two-halo terms combined, the mass function in the halo model would have had to be adjusted accordingly, i.e. more low-mass halos and less high-mass halos for short ranged interaction potentials. In short, the halo mass function has to be known for each particle interaction potential that is considered. It is only then possible to make conclusive statements about the halo density profiles.

6 Momentum-Density Correlations

6.1 Operator Expression

Another important and conceptually straightforward application of KFT is the calculation of momentum power spectra [57]. These are hardly accessible with the standard analytic approaches to cosmic structure formation. The momentum field could naively be constructed as

[TABLE]

with the velocities $u_{j}(t_{1})$ of the particles $j$, assuming all particle masses to be equal. This momentum field lacks any spatial dependence, since the microscopic degrees of freedom are the phase-space coordinates of all particles. In order to retain any spatial information, we have to impose that each particle can contribute to the momentum at a position $q$ if and only if it is at this position,

[TABLE]

In the last step, we have identified the Dirac delta distribution with the density of the $j$-th particle at position $q$ and time $t_{1}$. Thus, the new field $\Pi$ is a momentum density.

Fourier-transforming $\Pi(q,t_{1})$ and replacing the phase-space coordinates of particle $j$ by functional derivatives with respect to the corresponding source fields, we find the momentum-density operator

[TABLE]

composed of the one-particle operators

[TABLE]

with $\hat{\rho}(1)$ from (27). At this point, we note that components of the momentum-density field are available by specifying the operator to

[TABLE]

where $\alpha=(1,2,3)$ enumerates the Cartesian vector components.

The application of $r\geq 1$ of the operators (98) to the generating functional translates the generator field, as discussed in Sect. 2.2, and pulls down the momentum trajectories from the phase factor in (6). Having applied these operators and setting the generator field $\mathbf{J}$ to zero afterwards, we arrive at the expression

[TABLE]

where we have used the definitions (30) and (31).

In order to facilitate results obtained for density power spectra from Sect. 4.1, we employ the Fourier transform of the potential gradient (37) and replace the initial momentum of particle $j$ by a partial derivative with respect to the shift vector,

[TABLE]

In this way, we can write the momentum operator as

[TABLE]

such that the previous expression may be written as

[TABLE]

with $Z[\boldsymbol{L}]$ given by (40). Thus, the calculation of $r$-point correlations of the momentum-density differs from the calculation of density correlations only by the application of the product of $\hat{p}_{j}$-operators. In this way, the factorisation of the free generating functional can be applied in this context as well.

6.2 Perturbation Theory and Diagrams

In order to include particle interactions, a perturbative approach is necessary. One such approach begins with the Taylor expansion of the exponential factor in $Z[\boldsymbol{L}]$ in terms of the interaction operator $\hat{S}_{\mathrm{I}}$. In this Section, we formulate rules allowing the systematic calculation of all terms in the perturbation series by employing diagrams representing these terms. These rules have been stated for pure density correlations in [43]. We extend them here to also include momentum-density correlations. We note that this is effectively a perturbative expansion in the interaction potential, and thus corrections to the momentum trajectories need to be accounted for at the appropriate order. Thus, the linear correction to a momentum-density correlation is given by the acceleration of one particle in (102), which is of zero-th order in the Taylor expansion of the exponential, plus the linear term of the perturbative expansion. At the end of this Section, we will give a simple example of the linear correction to a $\left\langle\rho\Pi\right\rangle$ power spectrum.

Formulating the rules below, we assume that an $n$-point correlation function is to be calculated at $m$-th perturbation order. The diagrammatic representation of the terms correcting for forces between the particles can then be constructed following these rules:

• (i)

• a.

Attach $s=n+2m$ wave vectors pointing outward (represented by arrows) to a circle marking the free generating functional $Z_{0}$, where $s$ is the total number of density, momentum-density, and response field operators.

• b.

Of these, distinguish $r$ by dashed arrows, representing the momentum-density operators $\hat{\Pi}_{j_{l}}$, from the $(s-r)$ solid arrows representing density fields associated with either $\rho_{j_{l}}$ or $B_{j_{l}}$.

• (ii)

The time ordering is counter-clockwise along the circle, with the latest time being at the top. Interactions are represented by two lines attached to $Z_{0}$ at the same point in time. If the correlator is simultaneous, the external wave vectors are also attached to the same point.

• (iii)

The interaction potential is represented by a circled $v$ connected to a pair of internal wave vectors, which are marked with a prime. If the potential is translation-invariant, the connected internal wave vectors point in opposite directions and have the same magnitude.

• (iv)

• a.

Each response field is marked by a circle segment starting at a negative internal wave vector and connecting two different wave vectors. The circle segments always end at a later time.

• b.

Distinguish dashed circle segments, representing a deviation of the actual from the freely evolved particle position and given by the application of $\hat{S}_{\mathrm{I}}$, from solid circle segments corresponding to a change of the particle momentum by the second term in (102). These lines can only connect a response field (solid arrow) with a momentum-density (dashed arrow).

• c.

At a dashed wave vector, either no or one solid circle segment must end, but arbitrarily many dashed circle segments may end at the same wave vector, internal or external.

• (v)

Each diagram is assigned a multiplicity counting the equivalent diagrams.

To illustrate the construction of the diagrams and the physical meaning of the terms, we consider the linear-order correction to an equal-time 2-point cross-correlation function of a density $\rho$ and a momentum-density field $\Pi$. The respective diagrams are shown in Fig. 8.

By rules (i) and (ii), we attach one solid density line and one dashed momentum-density line at the same point to the circle representing the free generating functional for an equal-time correlation. Rule (iii) connects the internal wave vectors by the interaction potential. Following rule (iv), we construct all possible particle identifications. The first two diagrams represent the linear term of the Taylor expansion

[TABLE]

Without loss of generality, we here enumerate the particles in clockwise order starting at the top left. This term accounts for the displacement of a particle from its position according to inertial motion at the time of evaluation. The third diagram represents the deviation of the particle’s momentum due to the interaction with other particles. The term is

[TABLE]

The translation tensor is given by (30), and the Kronecker symbols identify an internal with an external particle.

At quadratic order in the interaction potential, considerably more diagrams contribute because two interactions of particles occur in the evolution of the system. Three scenarios can be distinguished: (1) one external particle undergoes two interactions, (2) each external particle undergoes one interaction, (3) one external particle scatters from a previously scattered internal particle.

6.3 Momentum-Density Power Spectra

In a similar fashion as for the density power spectra (69) and (87), we can write down a closed expression for the momentum-density power spectrum. We focus here on one scalar function that can be constructed from the $\left\langle\Pi\times\Pi\right\rangle$ correlation tensor. More specifically, we consider the projection of $\Pi$ perpendicular to the mode $k$. This power spectrum is needed in order to calculate the amplitude of secondary temperature fluctuations in the cosmic microwave background caused by Thomson scattering off free electrons moving with the bulk of structures. This is called the kinetic Sunyaev-Zel’dovich effect [58].

The result is

[TABLE]

with the definitions for $Q_{0}$, $Q$ and $\left\langle S_{I}(t)\right\rangle$ of equations (70) and (86).

The function $a_{\bot}(q)$ is the correlation function of the initial momentum components perpendicular to the line connecting the positions of the two momenta being correlated. It can also be expressed by the initial power spectrum $P_{\delta}(k)$ of density fluctuations and is given by

[TABLE]

using the definitions for $a_{1}(q)$ and $a_{2}(q)$ in (72).

Taking the limit of expression (106) for large scales or small wave numbers, $k\rightarrow 0$, we recover the standard result of [59] for the slope of the power spectrum,

[TABLE]

In Fig. 9, we compare our results with those obtained by [60] who used Eulerian standard perturbation theory. The authors of [60] calculate the power spectrum of the momentum-density by separating it into unconnected and connected terms using Wick’s theorem and evaluating these at one-loop order. This captures the interactions only partially, while the KFT approach predicts a higher amplitude when using the averaged interaction term (86). In addition, our results suggest that the spectrum increases proportional to $k^{2}$ at wave numbers $\lesssim 10^{-2}\,h/\mathrm{Mpc}$, while the SPT results fall below this slope approximately at wavenumbers above $10^{-3}\,h/\mathrm{Mpc}$.

7 The BBGKY hierarchy in KFT

To relate KFT to conventional kinetic theory, we now want to extract evolution equations for the phase-space density and its higher-order correlators. We will first define a phase-space density operator and generalise the generating functional as well as the interaction term from Sect. 2 to cover all of phase space.

In analogy to the number density (22) and its operator (27), we introduce the phase-space density

[TABLE]

and its operator expression in a six-dimensional Fourier space

[TABLE]

where $x=(q,p)$ is again the six-dimensional phase-space position with its Fourier conjugate $s=(k,\ell)$.

In this Section, we will always consider equal-time, $r$-point phase-space density cumulants with a Taylor expansion of the interaction up to $N$-th order,

[TABLE]

where $\tilde{f}(1)=\tilde{f}(t,s_{1})$, $\tilde{f}(2)=\tilde{f}(t,s_{2})$ and so on, and $Z_{0}[\boldsymbol{L}]$ denotes the free generating functional with $r$ phase-space density operators having been applied and the source fields set to zero. Similar to (20), the free generating functional is then

[TABLE]

with the shift operators

[TABLE]

The interaction (39) can then be extended to cover all of phase space,

[TABLE]

Here, the potential $\tilde{v}(1^{\prime})$ from (39) is extended to read

[TABLE]

the integration measure is $\mathrm{d}1^{\prime}=\mathrm{d}t_{1}^{\prime}\mathrm{d}^{6}s^{\prime}_{1}/(2\pi)^{6}$, and the phase-space response-field operator is defined by

[TABLE]

Before computing the evolution equations of phase-space density cumulants, we need to introduce the phase-space density current

[TABLE]

with its Fourier-space operator

[TABLE]

Applying a phase-space density-current operator $\widehat{fu}$ together with $r-1$ phase-space density operators $f$ to the free generating functional and setting the source fields to zero results in

[TABLE]

The time evolution of the phase space correlator (111) can now be found by direct computation,

[TABLE]

In the first term, the time derivative acts on the generating functional

[TABLE]

while it acts on the time-integral boundary in the interaction in the second term,

[TABLE]

where we used $g_{qp}(t,t)=0$ and set $t_{1}^{\prime}=t$.

Inserting the derivatives back into (120), we can use (119) in the first term, and (116) in the second term. Identifying applied operators in the form of the free generating functional and collecting them into operators finally results in

[TABLE]

where $\tilde{f}(j-1^{\prime}):=\tilde{f}(t,s_{j}-s_{1}^{\prime})$.

This shows that the time evolution of the $r$-point phase-space correlator is divided into two effects: The first is the convection of the phase-space density described by the first term. The second term contains all possible interactions with an additional point in phase space. Here, $r+1$-point phase-space correlators appear. Consequently, evolution equations for $r+1$-point phase space density correlators are needed in order to solve the $r$-point evolution equation. These again contain convection terms as well as $r+2$-point phase space density correlators. For the latter, evolution equations are needed that contain $r+3$-point correlators, and so on. An (almost) infinite hierarchy of coupled differential equations unfolds, known as the BBGKY hierarchy. A Fourier transform back into configuration space shows that (123) are the exact terms of the conventional BBGKY hierarchy [41].

A truncation criterion of the BBGKY hierarchy emerges here from the order of the Taylor expansion of the interaction: In the evolution equation for the $r$-point correlator, the $r+1$-point correlators appear only in the $(N-1)$-st order in the interaction, reduced by one. This is repeated in the time evolution of the two $r+1$-point correlators: the $r+2$ point correlator is again reduced by one order in the interaction and so on for each further step up the hierarchy. Once the [math]-th order is reached, only convection terms and no higher-order correlators appear in the evolution equation, and the hierarchy ends.

8 Fluids in KFT

8.1 Introduction

Cosmic structure formation is governed by both dark and baryonic matter. Therefore, it is crucial to extend KFT to capture the physics of both particle types. Unlike dark matter, baryonic interactions are not limited to gravity only. This causes the cosmic matter power spectrum to attain additional features such as baryon acoustic oscillations (BAO) as well as disturbances in the spectrum at small scales due to e.g. pressure and baryonic cooling. We describe here one way towards including mixtures of dark and baryonic matter into KFT.

8.2 Extension of the generating functional

In principle, fluid dynamics should be captured by KFT if microscopic gas-particle interactions were taken into account. However, it is a regime hard to reach in any expansion of the interaction. To acquire access to fluid dynamics nonetheless, we implement a fluid model based on an idea close to the conventional approach to hydrodynamics: We introduce mesoscopic particles similar to conventional fluid elements and demand a hydrodynamical scale hierarchy. Then, on the scale of the mesoscopic particles, local thermal equilibrium can be assumed to have been established. The properties of the mesoscopic particles are then described by state variables such as pressure and energy density. At the same time, we demand that the mesoscopic particles are much smaller than any scale we are interested in and can be modelled as point-like.

For simplicity, we assume an isotropic fluid which is fully described by a spatially dependent velocity field $u(q)$, a pressure field $P(q)$ and an energy-density field $\epsilon(q)$. In addition, we assume that the pressure is only caused by the random velocities of the microscopic particles. Then, both the pressure and the internal energy-density are proportional to the enthalpy-density $h$,

[TABLE]

In order to sample the fluid, the mesoscopic particles need to capture these properties: they need to contain information about their position, the momentum, and the enthalpy at that position. Hence, we endow each particle with these three degrees of freedom,

[TABLE]

and describe their non-interacting dynamics by ballistic motion, conserving momentum and enthalpy. These dynamics are expressed by the Green’s function

[TABLE]

From here on, a free generating functional for KFT can be constructed in analogy to (20)

[TABLE]

where the tilde marks the free generating functional for the fluid model, $J_{j}=(J_{q_{j}},J_{p_{j}},J_{\mathcal{H}_{j}})$ and

[TABLE]

8.3 Interaction operator

In addition, acceleration due to pressure gradients and pressure-volume work must be taken into account. For viscous hydrodynamics, diffusion of energy and momentum must be added. Within KFT, these effects are included via interactions between the mesoscopic particles. The exact form of the interaction operator is found via an approach similar to smooth-particle hydrodynamics (SPH) [61]. For the purposes of this review, we only sketch the derivation of the interaction operator for pressure-volume work. A thorough derivation including all terms of ideal and viscous fluid dynamics can be found in [62].

The pressure field split into particle contributions reads

[TABLE]

where $\mathcal{H}_{i}$ is the enthalpy and $q_{i}$ the position of the $i$-th particle.

Following the Euler equation, the time evolution of the momentum field contains the term

[TABLE]

where we inserted the discretised pressure field from (129).

The momentum change at position $q_{1}$ is associated to particles according to their (spatial) contribution at that position. To this end, we weigh the momentum change at $q_{1}$ with a Dirac delta distribution around the $j$-th particle’s position and integrate over the entire space,

[TABLE]

For the interaction operator, this object must be expressed by functional derivatives acting on the free generating functional. However, this is not possible for the inverse density. We handle this complication by approximating the inverse densities in a Taylor series around the mean density $\bar{\rho}$ of the ensemble. As a first step, we truncate the approximation already at [math]-th order,

[TABLE]

For the pressure-volume work, a similar term for the enthalpy change of the $j$-th particle, $\dot{\mathcal{H}}_{j}$, can be derived in an analogous calculation. From the momentum and enthalpy changes, the interaction operator for an ideal, isotropic fluid can then be constructed,

[TABLE]

with $\mathrm{d}1=\mathrm{d}t\mathrm{d}^{3}k_{1}$ and the one-particle operators

[TABLE]

For an ensemble characterized by the free generating functional (127) and the interaction operator (133), macroscopic evolution equations for density, momentum density, and energy density can be extracted. As shown in [62], these are indeed the continuity, Euler, and energy-conservation equations of ideal hydrodynamics. If an interaction operator for diffusive effects is added, the Navier-Stokes and energy conservation of viscous hydrodynamics are recovered.

9 Macroscopic Reformulation of KFT

From the perspective of most quantum and statistical field theories, the generating functional (5) of KFT is rather unusual in the sense that the path integral is expressed in terms of microscopic degrees of freedom even though we are actually interested in macroscopic fields like the density. As described before, the decisive advantage of this approach is the simplicity of the equations of motion for the microscopic degrees of freedom.

However, to facilitate the application of established perturbative as well as non-perturbative field-theoretical techniques, we would like to reformulate the KFT partition function as a path integral over macroscopic fields. It turns out that this will lead to a resummation of an infinite subset of terms appearing in the microscopic perturbative expansion in the interaction operator. This will allow us to treat dark matter particles in terms of their fundamental Newtonian dynamics rather than the improved Zel’dovich dynamics. The presentation here is based on the more detailed and general derivation in [63] which uses the full phase-space density rather than the spatial density to preserve momentum information.

9.1 Macroscopic action

We begin by using (9) and (12) to write the partition function as

[TABLE]

where we represented the delta distribution in terms of a functional Fourier integral with respect to an auxiliary field $\boldsymbol{\chi}$ with components \chi_{j}=\bigl{(}\chi_{q_{j}},\chi_{p_{j}}\bigr{)}. We further introduce the combined microscopic field $\boldsymbol{\psi}:=(\boldsymbol{x},\boldsymbol{\chi})$ and its action

[TABLE]

which splits into a part $S_{\psi,0}$ describing the free motion of the particles and a part $S_{\psi,\mathrm{I}}$ describing their interactions. Using (15) and (17), the latter can be expressed as

[TABLE]

Inserting (37) and defining the dressed response field $\mathcal{F}$ via its Fourier transform,

[TABLE]

allows us to write the interaction term as

[TABLE]

where we also introduced the dot product as a short-hand notation for integrating over field arguments.

We now replace the explicitly $\psi$-dependent field $\rho$ with a new formally $\psi$-independent field $\phi_{\rho}$, using a functional delta distribution,

[TABLE]

This way, the new field $\phi_{\rho}$ effectively still carries all the information contained in $\rho$. Most importantly, $\phi_{\rho}$- and $\rho$-correlation functions are identical. However, to emphasise their different origin we will deliberately call $\phi_{\rho}$ the macroscopic and $\rho$ the collective density field in the following.

Similar to (92), we now express the delta distribution as a functional Fourier transform with respect to a new macroscopic auxiliary field $\phi_{\beta}$,

[TABLE]

Pulling all $\psi$-independent parts to the front, we find that the remaining microscopic part of the path integral assumes the form of the free generating functional of $\rho$- and $\mathcal{F}$-correlators,

[TABLE]

with $\phi_{\beta}$ and $\phi_{\rho}$ playing the roles of the source fields for the collective fields $\rho$ and $\mathcal{F}$, respectively.

Finally, defining the combined macroscopic field $\phi:=(\phi_{\rho},\phi_{\beta})$ and the free generating functional of collective-field cumulants $W^{\rho,\mathcal{F}}_{0}:=\ln Z^{\rho,\mathcal{F}}_{0}$, we arrive at the result

[TABLE]

with the macroscopic action

[TABLE]

We emphasise that this reformulation is exact and hence (143) still contains the complete information on the microscopic dynamics, even though $S_{\phi}$ does not depend on $\psi$ any more. The microscopic information is now encoded in the free generating functional W^{\rho,\mathcal{F}}_{0}\bigl{[}\phi_{\beta},\phi_{\rho}\bigr{]} and thus, by means of a functional Taylor expansion, in the free collective-field cumulants

[TABLE]

9.2 Including density-density and density-momentum correlations

Since we will use this macroscopic reformulation of KFT in Sect. 9.4 to treat Newtonian dynamics, whose Green’s function (54) is limited from above, we have to take the complete expression for the initial probability distribution (64) into account, which includes all correlations between the initial phase-space coordinates of particles. This means that the lengthy polynomial $\mathcal{C}(\boldsymbol{p})$ introduced in (64) must be included. This makes the explicit calculation of the above cumulants anything but straightforward. By adopting the essential ideas of the so-called Mayer cluster expansion [64], we have however been able to condense this process into a small set of rules for Feynman-like diagrams. Here, we will restrict ourselves to reviewing the recipe for cumulants involving the particle number density $\rho$, while the complete technical derivation for cumulants involving the phase-space density $f$ can be found in [44].

The general form of the diagrams in this approach is very simple: particles are represented by dots, and correlations between them by different types of connecting lines. Since we are only interested in connected contributions to correlations, we need to consider connected diagrams only, i.e. such through which we can trace a continuous path. A general free cumulant as in (145) is then ordered in terms of the number of particles $\ell$ being correlated in such a connected way,

[TABLE]

Note that, as a consequence of statistical homogeneity, the sum over particle numbers truncates at the number of density fields in the cumulant. Combined with the diagram rules, this ensures that all cumulants have a finite number of terms, i.e. KFT can exactly describe the highly non-linear effects of initial correlations on the free-streaming evolution with a finite number of explicit expressions.

We begin with the calculation of pure density cumulants since mixed cumulants involving the dressed response field $\mathcal{F}$ are derived from the former by the insertion of simple response factors. We first need the crucial concept of a field-label grouping. For any fixed number of particles $\ell$, we group the field labels $(1,\ldots,n_{\rho},1^{\prime},\ldots,n_{\mathcal{F}}^{\prime})$ of the cumulant into $j=1,\ldots,\ell$ non-empty sets $\mathrm{I}_{j}$. Any such collection of sets is called a field-label grouping. We define a variant of the phase-translation vectors (30) in terms of these groupings as

[TABLE]

A general $\ell$-particle density cumulant is then given as a sum over all possible field-label groupings,

[TABLE]

where we introduced the generalisation of the damping amplitude $Q_{0}$ from (70) as

[TABLE]

The term $\tilde{\Sigma}_{\mathrm{CDiag}}^{(\ell)}(\mathrm{I}_{1},\ldots,\mathrm{I}_{\ell})$ represents the sum over all possible connected $\ell$-particle correlation diagrams. The possible diagram line types are

[TABLE]

with the phase $\theta:=\mathrm{i}\theta_{q}+\theta_{p}$ and

[TABLE]

In addition to the momentum correlations introduced in Sect. (3.3), $C_{\delta_{j}\delta_{k}}$ represents density autocorrelations, and $C_{\delta_{j}p_{k}}$ density-momentum correlations between the initial positions of particles $j$ and $k$. Due to statistical homogeneity, the above Fourier integrals only need to be performed over the relative initial coordinates $q^{\,(\mathrm{i})}_{jk}=q^{\,(\mathrm{i})}_{j}-q^{\,(\mathrm{i})}_{k}$ of particles, where $j<k$ in all cases by convention. All translation vectors are evaluated at the initial time. One may draw these lines between particle dots according to the following rules:

• •

No self-correlations in the form of subdiagrams may occur for any kind of correlation line.

• •

Any pair of particles can be connected directly by at most one line, i.e. no subdiagrams may occur for any combination of correlation lines.

• •

No particle may have more than one solid $\delta$-line attached to it, i.e. no subdiagrams my occur which include .

By convention, there is a flow of a general Fourier momentum $\mathcal{K}^{\,(\mathrm{i})}_{jk}$ along each line from smaller to larger labels. The algorithm of Feynman rules for evaluating $\tilde{\Sigma}_{\mathrm{CDiag}}^{(\ell)}(\mathrm{I}_{1},\ldots,\mathrm{I}_{\ell})$ for a given grouping of field labels is then as follows.

1. (1)

Choose a fixed graphical arrangement of $\ell$ dots labeled $j=1,\ldots,\ell$, representing the particles carrying the label sets $\mathrm{I}_{j}$. These are the vertices of the diagrams. 2. (2)

Given this fixed set of vertices, draw all possible diagrams using the lines in (151), subject to the above rules. Repeat the next three steps for all these diagrams. 3. (3)

For any closed loop in a diagram, assign one of the lines in the loop connecting vertices labeled $j<k$ a loop momentum ${k}^{\,(\mathrm{i})}_{jk}$ and introduce an integral $\int_{k^{(\mathrm{i})}_{jk}}$. 4. (4)

Pick a vertex $j$ which has only one line left with undetermined momentum and use the conservation law

[TABLE]

to fix the momentum of this line. Incoming momenta from vertices with smaller labels are counted positive while momenta outgoing to vertices with larger labels are counted negative. Momenta associated with vertices not connected to vertex $j$ are zero. 5. (5)

Consecutively go through all other vertices to fix the remaining undetermined momenta by repeatedly applying the previous step.

Once a pure density $n$-point cumulant is known any mixed $(n=n_{\rho}+n^{\prime}_{\mathcal{F}})$-point cumulant is simply obtained by first inserting $n^{\prime}_{\mathcal{F}}$ response factors and interaction potentials as

[TABLE]

The response factors are defined simply as

[TABLE]

where $\mathrm{I}(r)$ is the set of field labels containing $r$. The retarded nature of the particle propagators contained in these response factors also leads to the general property that only such field-label groupings $\{\mathrm{I}_{1},\ldots,\mathrm{I}_{\ell}\}$ give non-vanishing contributions to mixed cumulants (154) which have the following property: For all $\mathrm{I}_{j},j=1,\ldots,\ell$, there must be at least one $\rho$-label $r\in\mathrm{I}_{j}$ such that $t_{r}\geq t_{u^{\prime}}$ for all $\mathcal{F}$-labels $u^{\prime}\in\mathrm{I}_{j}$. This can be understood as the causal flow of responses having to terminate at some density field label $r$ at later time since there is no instant response of the spatial density to forces acting on particles.

9.3 Macroscopic perturbation theory

The reformulated path integral (143) allows us to set up a new perturbative approach to KFT following the standard procedure familiar from quantum and statistical field theory, i.e. in terms of propagators and vertices. For this purpose, we first split the action (144) into a propagator part $S_{\Delta}$, collecting all terms quadratic in $\phi$, and a vertex part $S_{\mathcal{V}}$, containing the remaining terms,

[TABLE]

with the definitions

[TABLE]

introducing the inverse macroscopic propagator $\Delta^{-1}$ and the macroscopic $(n_{\beta}+n_{\rho})$-point vertices $\mathcal{V}_{\beta\dotsm\beta\rho\dotsm\rho}$.

We further define the macroscopic generating functional $Z^{\phi}$ by introducing a source field {M=\bigl{(}M_{\rho},M_{\beta}\bigr{)}} for the combined macroscopic field $\phi$ into the partition function,

[TABLE]

Then, the vertex part of the action can be pulled in front of the path integral by replacing its $\phi$-dependence with functional derivatives with respect to $M$, \hat{S}_{\mathcal{V}}:=S_{\mathcal{V}}\bigr{|}_{\phi\rightarrow\frac{\delta}{\mathrm{i}\delta M}}, acting on the remaining path integral,

[TABLE]

Expanding the first exponential in (160) in powers of the vertices now gives rise to a new perturbative approach that we will refer to as the macroscopic perturbation theory.

Within this approach, the interacting macroscopic-field cumulants are obtained by taking appropriate functional derivatives of the macroscopic cumulant-generating functional $W^{\phi}[M]:=\ln{Z^{\phi}[M]}$. In particular, the lowest perturbative order of the 2-point density cumulant $G_{\rho\rho}$ is given by the $\rho\rho$-component of the macroscopic propagator,

[TABLE]

Adopting the usual field-theoretical language, we will refer to this as the tree-level expression for $G_{\rho\rho}$. For further detail on the general properties of the macroscopic perturbation theory, we refer the reader to [63], where we introduce a Feynman-diagram language to compute the higher-order perturbative contributions – i.e. the so-called loop-contributions – systematically.

To find expressions for the inverse propagator $\Delta^{-1}$ and the vertex term $\mathcal{V}_{\beta\dotsm\beta\rho\dotsm\rho}$, we insert the functional Taylor expansion of $W^{\rho,\mathcal{F}}_{0}$ into the macroscopic action (144) and identify terms with (157) and (158), respectively,

[TABLE]

Here, $\mathcal{I}$ denotes the identity 2-point function,

[TABLE]

The propagator $\Delta$ is then obtained by a combined matrix and functional inversion of (162), defined via the following matrix integral equation,

[TABLE]

The matrix part of this inversion can be performed immediately and yields

[TABLE]

where we defined the retarded and advanced (macroscopic) propagators

[TABLE]

containing the remaining functional inverses. They describe the linear response of the density at time $t_{1}$ to a perturbation of the interacting system at time $t_{2}$.

Consequently, $\Delta_{\rho\rho}=\Delta_{\mathrm{R}}\cdot G^{(0)}_{\rho\rho}\cdot\Delta_{\mathrm{A}}$ is the 2-point density cumulant resulting from the linear response of the interacting system to the correlated motion of the free-streaming particles. Generally, $\Delta_{\rho\rho}$ will nevertheless be nonlinear in the initial density correlation because the free-streaming of particles in $G^{(0)}_{\rho\rho}$ itself builds up those nonlinearities as long as the initial particle momenta are correlated[44].

In [63], we describe how $\Delta_{\mathrm{R}}$, $\Delta_{\mathrm{A}}$ and $\Delta_{\rho\rho}$ can be computed explicitly for a given physical system. But even without specifying the system, it is always possible to express (167) formally in terms of a Neumann series by expanding the functional inverse in orders of $\mathrm{i}G^{(0)}_{\rho\mathcal{F}}$,

[TABLE]

From the definition of $\mathcal{F}$ in (138), it follows that $G^{(0)}_{\rho\mathcal{F}}\propto v$. Therefore, $\Delta_{\mathrm{R}}$, $\Delta_{\mathrm{A}}$ and $\Delta_{\rho\rho}$ contain terms of arbitrarily high order in the interaction potential $v$. This demonstrates that the lowest perturbative order within the macroscopic approach already captures effects which could only be described at infinitely high order within the microscopic perturbation theory.

A closer analysis reveals that, in a statistically homogeneous system, the transition from the micro- to the macroscopic formulation entails a resummation of all contributions which do not lead to any mode-coupling beyond the one already introduced by the free evolution. Due to this property, we will refer to the macroscopic reformulation of KFT as resummed KFT (RKFT).

9.4 Resumming Newtonian dynamics

Using RKFT enables us to treat dark-matter particles in terms of their fundamental Newtonian dynamics rather than the improved Zel’dovich dynamics. Since the inertial motion of particles on Zel’dovich-type trajectories already contains part of the gravitational interaction, even a first-order calculation in the microscopic perturbation theory captures the nonlinear evolution of the power spectrum over a wide range of scales remarkably well, as shown in Fig. 3 and discussed in detail in [42].

This is not the case for Newtonian dynamics, and any expansion to finite order in the Newtonian gravitational potential does not even reproduce the linear growth of structures correctly. Working with an averaged interaction term or the Born approximation, as described in Sects. 4.2 and 4.3, does not improve the situation either. This is because both approaches require deviations from the inertial trajectories caused by the (effective) forces acting between them to be small, which is not the case for the truly non-interacting inertial trajectories in Newtonian dynamics.

Using RKFT, however, allows to overcome this limitation. In Fig. 10, we show the tree-level RKFT result for the density contrast power spectrum,

[TABLE]

evolved from the time of CMB decoupling until today assuming Newtonian particle dynamics, and compare it to the power spectrum resulting from Newtonian free-streaming as well as the power spectrum obtained from linear Eulerian perturbation theory (EPT).

At small wavenumbers, $P_{\delta}^{(\Delta)}$ follows the linear EPT spectrum, before it starts to fall below the linear prediction at wavenumbers $k\gtrsim 20\,h\,\mathrm{Mpc}^{-1}$. When going to even higher wavenumbers it eventually approaches and follows the free-streaming spectrum for $k\gtrsim 10^{4}\,h\,\mathrm{Mpc}^{-1}$.

First of all, this shows that the partial resummation of Newtonian gravitational interactions captured by the macroscopic propagator allows to precisely recover the linear growth of structures on the largest scales even though the structure growth resulting from non-interacting Newtonian dynamics is orders of magnitude weaker. At the same time, the resummation does not compensate the small-scale damping effect caused by the particles’ initial momentum variance $\sigma_{1}^{2}$. This behaviour is reminiscent of the small-scale decay of the power spectrum found in Renormalized Perturbation Theory (RPT) [16, 17].

However, we can also infer that RKFT does not enhance this damping effect either, as $P_{\delta}^{(\Delta)}$ never drops below the free-streaming spectrum. This differs from the power spectra found in RPT or Zel’dovich dynamics, which overpredict the dissolution of small-scale structures[16, 17, 11, 43].

Altogether, the tree-level RKFT result for the density contrast power spectrum is thus found to capture the linear effects introduced by gravitational interactions in a way that is fully consistent with the underlying non-interacting Newtonian particle dynamics. Higher-order perturbative corrections within RKFT are currently investigated.

10 Mixtures of gas and dark matter in KFT

Baryonic matter with its additional non-gravitational interactions has a profoundly different effect on large-scale cosmic structures than dark matter. To gain a complete understanding of structure growth, it is therefore necessary to investigate how a mixture of dark and baryonic matter co-evolves through the cosmic history.

To describe such a mixture, we extend the macroscopic reformulation of KFT to include two different particle species. The long term goal hereby is to use this approach to study the non-linear evolution of structures which contain both baryonic and dark matter. In this formalism, the collisionless character of dark matter is kept while at the same time the fluid properties of baryonic matter needs to be taken into account.

We demonstrate the functionality of this formalism on structure growth in the linear regime. Our aim hereby is a proof of concept showing that KFT is able to correctly describe the co-evolution of mixtures. The foundations presented here can be used in the future to investigate nonlinear structure growth.

Furthermore, we demonstrate that our formalism can qualitatively reproduce some of the key effects baryons have on linear structure formation. Firstly, the suppression of structures on small scales due to the repulsive effects of baryonic pressure is demonstrated. Secondly, by also taking interactions between photons and baryonic matter into account, essential features of baryon-acoustic oscillations are reproduced.

10.1 KFT with two particle types

To describe the interactions of two different particle types in the framework of KFT we expand the macroscopically reformulated generating functional (160) to a system of two particle species in a straightforward fashion by replacing the macroscopic fields $\phi_{\rho}$ and $\phi_{\beta}$ as well as the collective dressed response field $\mathcal{F}$ by two-component vectors containing the respective field contributions for baryons and dark matter,

[TABLE]

Proceeding in a completely analogous fashion to Sect. 9, we obtain a macroscopic propagator that takes exactly the same form as in the single-particle species case (166). The only difference being that all 2-point functions appearing in the propagator are now replaced by $2\times 2$ matrices,

[TABLE]

The individual components of $G^{(0)}_{\rho\rho}$ describe the auto- and cross-correlations between the density fields of baryonic and dark matter, while the components of $G^{(0)}_{\rho\mathcal{F}}$ describe the response of these density fields to the interactions between particles of either the same or different species. These interactions are accordingly collected in a symmetric $2\times 2$ potential matrix,

[TABLE]

To be able to calculate the power spectrum for matter mixtures, we still need to find the appropriate interaction potential $v^{\mathrm{b}\mathrm{b}}$, which is the only entry containing not only gravitational interactions.

10.2 Baryon interactions

Unlike collisionless dark matter, baryonic matter is well-described by an ideal gas, and hence its dynamics are governed by the ideal fluid dynamics discussed in Sect. 8. Accordingly, we treat baryons in KFT as mesoscopic particles characterised by a position, a momentum and an enthalpy.

However, to facilitate the incorporation of mesoscopic particles into the RKFT framework, we approximate the gas as being isothermal. That is, we approximate the baryon temperature as being spatially homogeneous throughout the Universe, following the evolution of the mean gas temperature $T(t)$. In this case, every mesoscopic particle has the same enthalpy

[TABLE]

where $m^{\mathrm{b}}$ is the mass of a mesoscopic particle, $\gamma$ is the adiabatic index of the gas and $c_{\mathrm{s}}(T)$ is its sound velocity. This effectively reduces the dynamical degrees of freedom to just a position and a momentum, allowing us to treat baryons like microscopic particles experiencing a gravitational potential

[TABLE]

as well as an additional effective repulsive interaction potential modelling the effects of pressure

[TABLE]

where $\bar{\rho}^{\mathrm{b}}$ is the mean number density of the mesoscopic particles [65]. The different components of the potential matrix (172) are then given by $v^{\mathrm{dd}}=v^{\mathrm{bd}}=v^{\mathrm{db}}=v_{\mathrm{g}}$ and $v^{\mathrm{bb}}=v_{\mathrm{g}}+v_{\mathrm{p}}$. To arrive at the expressions (175) and (176), we have set the masses of both dark matter and mesoscopic particles equal. We can do this without loss of generality since we are dealing with effective particles. We are now able to obtain quantitative results for the evolution of mixed dark and isothermal baryonic matter.

10.3 Results at late times

Their density fluctuation power spectra are evolved from the epoch when matter and radiation decoupled at a redshift of $z=1100$ to a redshift of $z=100$. We choose a relatively high redshift for the final time since we currently perform the calculations only in the linear order of the resummed theory. Furthermore, the beginning of the reionization era makes it difficult to rigorously track the mean temperature evolution of the gas.

We use the same initial BBKS power spectra [10] for both dark and baryonic matter. However, in Fig. 12, we set the initial baryonic power spectrum to zero on wavenumbers $k\geq 1\,h\,\mathrm{Mpc}^{-1}$ as a crude approximation for the Silk damping, which suppressed baryonic structures at small scales during radiation decoupling [66].

We plot the results in Figs. 11 and 12, where we have divided the obtained power spectra by the power spectrum for pure dark matter to show the relative deviations caused by the incorporation of baryonic matter.

The results show that the repulsive pressure interactions suppress structure growth at small scales. The fact that the dark-matter power spectrum is also suppressed shows that RKFT is able to capture the interplay between the baryonic pressure interactions and the gravitational interactions between the two particle species. This also allows to describe the build-up of baryonic structures on scales where the initial baryonic power spectrum was set to zero, as can be seen in Fig. 12. Even the subsequent suppression of the dark-matter spectrum caused by the pressure of these new baryonic structures is captured.

The next steps towards gaining further insight into baryonic effects on structure formation after radiation decoupling are to drop the isothermal approximation and to develop effective potentials capturing additional physical properties of baryonic matter, e.g. radiative cooling effects.

10.4 Mesoscopic particles with photons

To investigate structure formation before radiation decoupling at $z\approx 1100$, it is necessary to also take into account that the pressure forces acting on the baryons are mainly due to photons. Therefore, we modify the set-up of the mesoscopic particle by constructing it in such a way that it contains both baryons and photons (cf. Fig. 13).

When investigating structure growth between matter-radiation equality and photon decoupling, we assume that the mass of the mesoscopic particle and hence its gravitational potential are dominated by its baryonic content. Hence, the gravitational potential of the particle is the same as in (175). However, the enthalpy of the particle is dominated by the photons since they are much more abundant – the ratio in the number densities being given by $\eta\approx 6\cdot 10^{-10}$ – as well as due to their ultrarelativistic nature.

The resulting pressure potential in a baryon-photon fluid is given by

[TABLE]

with $k_{\mathrm{B}}$ being the Boltzmann constant, $m_{\mathrm{P}}$ the proton mass, $T_{\mathrm{CMB}}$ the current temperature of the CMB, $H_{0}$ the current value of the Hubble function and $\Omega_{\mathrm{m},0}$ the current dimensionless matter density parameter.

Assuming an initial BBKS power spectrum, the resulting matter power spectrum can be seen in Fig. 14. It is compared to the power spectrum obtained using the Eisenstein-Hu transfer function [67].

We clearly see that our model qualitatively reproduces features of BAOs close to their positions in the semi-analytic Eisenstein-Hu model. The main difference between the two models is a stronger suppression of structures above a wavenumber of $k\approx 10^{-2}\,h\,\mathrm{Mpc}^{-1}$.

This suppression is a consequence of assuming a simplified initial BBKS spectrum, which does not take the effect of photons on the transfer function into account. In the future, we will correct for this to improve the quantitative agreement between our results and Eisenstein-Hu. Beyond this, the mesoscopic particles need to be modified to also take into account the dissolution of structures at small scales due to Silk damping. This can be achieved by treating the enthalpy as an additional dynamical property of the particles.

11 Structure formation with modified theories of gravity

General relativity successfully describes the physics of gravitational interactions on a vast range of scales under the assumption that the total matter and energy budget of the Universe is to be augmented by the presence of three unknown ingredients, namely dark matter, dark energy in form of a cosmological constant, and the inflaton field. It is understood that the theory can only be used as an effective field theory valid up to the Planck scale, which is not renormalizable. Even though it is observationally in good agreement with cosmological data, general relativity suffers from severe and persistent theoretical obstacles, one of them being the cosmological constant itself. These have motivated the construction of alternative theories of gravity, both in the large- and in the small-scale regimes.

In the same spirit as general relativity, one can construct other effective field theories for the underlying gravitational interactions. Keeping the fundamental properties of general relativity intact (for instance unitarity, locality and invariance under Lorentz symmetry), any modification of gravity will inevitably introduce new propagating degrees of freedom. They usually come in form of additional scalar, vector, or tensor fields. Therefore, the majority of alternative theories of gravity can be classified by means of three important classes (see [68] for a recent review on alternative theories of gravity)

• •

scalar-tensor theories: quintessence, Brans-Dicke, $k$-essence, Galileon, Horndeski, beyond-Horndeski, and DHOST-type theories;

• •

vector-tensor theories: vector-Galileon, generalized Proca theories, beyond-generalized Proca, or multi-Proca theories;

• •

tensor-tensor theories: massive gravity, bigravity, or multi-gravity theories;

If present, these additional degrees of freedom weaken gravity on cosmological scales. They either form a condensate whose energy density sources self-acceleration of the cosmic expansion, or they contribute as a condensate whose energy density compensates the cosmological constant. The latter is used to solve the cosmological-constant problem, while the former is useful for applications to dark energy and inflation. The presence of some of these degrees of freedom, for instance cosmic vector fields, could even give rise to a violation of the cosmological principle. Quite generically, these alternative theories of gravity will modify

• •

the background evolution by altering the Hubble function $H$;

• •

perturbations by enforcing altered structure formation due to: an altered time sequence of gravitational clustering; the evolution of peculiar velocities and the number density of collapsed objects; modifications in the gauge-invariant matter density contrast and its relation to the gravitational potential; changes in the gravitational slip parameter, the effective gravitational potential, and the growth rate;

These alternative theories come with free parameters or even with free functions. Restricting their allowed parameter space is an indispensable task in order to test them against general relativity. At the background level, geometrical probes measuring the angular-diameter distance as a function of redshift (CMB and BAO) and the distance-redshift relation of supernovae can be used to constrain the allowed field space. Concerning linear perturbations, the CMB temperature anisotropies provide a supplementary tool to probe the evolution of the inhomogeneous perturbations. However, these probes are not sufficiently powerful to disentangle certain degeneracies between different alternative gravity theories. Hence, testing alternative theories of gravity on cosmological scales will require a thorough understanding of non-linear cosmic structure formation. Nevertheless, high-resolution simulations of cosmic structures in the non-linear regimes for the many existing alternative theories of gravity would not be affordable or almost impossible. Therefore, an essential task for future inference from cosmological data will be to develop analytical approaches to non-linear cosmic structure formation.

KFT offers a promising framework for this purpose. Modifications of gravity theories do not alter the KFT formalism itself. Many alternative theories of gravity can thus easily be incorporated into the fundamental KFT equations, and e.g. non-linear power spectra obtained from KFT can be directly compared to observations. The few necessary adaptations concern

• •

the background evolution, encoded in the Hubble function $H$;

• •

the Hamiltonian of the system, including additional degrees of freedom;

• •

a possible evolution of the effective gravitational coupling constant $G_{\mathrm{eff}}$;

• •

the linear growth factor $D_{+}$; and

• •

changes to the gravitational potential.

As a proof of concept, we show here results obtained within a specific dark-energy model taken from the class of generalized Proca theories [69, 70, 71, 72]. Concrete self-accelerating solutions have successfully been constructed [73, 74, 75] and compared with background probes [76]. The non-linear power spectrum obtained for this specific modification of gravity is shown in Fig. 15. One immediate observation is the increase of power on scales $k\sim 2\,h\,\mathrm{Mpc}^{-1}$, which would give rise to enhanced structure formation on these scales. Direct comparisons of such results with $N$-body simulations would be highly interesting.

12 Beyond cosmological applications

As stated at the beginning of this article, Kinetic Field Theory can be applied to any classical ensemble of particles in or out of equilibrium. To illustrate this, an interesting exercise is to consider an ensemble of Rydberg atoms out of equilibrium.

In the past decade Rydberg atoms – atoms excited into a very high principal quantum number $n$ – have enjoyed increasing popularity in theoretical as well as experimental physics. Due to the large separation of the outer electron and the nucleus, Rydberg atoms have a strong dipole moment. Thus, the interaction of atoms in Rydberg states with each other even when they are separated by a microscopic distance is still strong, whereas the van der Waals forces between two ground-state atoms separated by a macroscopic distance would be negligible. In addition, Rydberg atoms have a long lifetime of $\sim 100\,\mu\mathrm{s}$. This makes them excellent candidates for the study of interactions in many-body systems [77, 78].

The two properties of Rydberg atoms that make them such interesting systems for KFT are their strong interactions and the Rydberg blockade [79, 80] illustrated in Fig. 16 and Fig. 17. The Rydberg blockade provides a natural (anti-)correlation function for these systems, as no two Rydberg atoms can be closer than two times the Rydberg blockade radius $R_{b}$ initially.

Strictly speaking, a system of Rydberg atoms has to be treated in the scope of quantum mechanics. However, if treated as a classical system with given anti-correlation function due to the Rydberg blockade radius, they can be useful to distinguish between quantum and classical correlations in strongly-correlated many-body systems. Finding the barrier at which quantum effects can no longer be ignored is a central task in the field of quantum mechanics and may help to understand the nature of correlations due to quantum-mechanical effects.

12.1 Rydberg systems and KFT

The initial correlation function is set up with the following picture in mind: $N$ ground-state atoms with a high packing fraction are simultaneously excited into the Rydberg state. Due to the Rydberg blockade, no two Rydberg atoms are closer than $2\,R_{b}$. Thus, the Rydberg atoms can initially well be approximated as hard spheres. The excitation is assumed to be a Gaussian random process, even though, in reality, Gaussianity is broken precisely because of the Rydberg blockade effect due to a strong interaction between the atoms. It stands to reason that, although Gaussianity is broken on very small scales around each Rydberg atom, it is again restored on scales outside the Rydberg blockade, where interactions between ground-state atoms are negligible and atoms are excited at random. Since the initial conditions are formulated at the scales of the Rydberg-blockade radius, and the volume of the Rydberg blockade is very small compared to the relevant scales, Gaussianity can be safely assumed. Thus, the initial distribution of Rydberg atoms will be described by a multivariate Gaussian with a correlation matrix determined by anti-correlations due to the Rydberg blockade. The initial correlations can thus be formulated just as in Sect. 3.3.

The system of $N$ Rydberg atoms is assumed to be homogeneous and isotropic, with the particle positions initially being correlated due to the Rydberg blockade, but momentum correlations being absent except for a momentum dispersion due to a physical temperature that will be set externally. In Fig. 18 the initial correlation function which is directly sampled from a molecular-dynamics (MD) simulation is shown. When treating Rydberg atoms as classical particles, their trajectories will be subject to Hamilton’s equations of motion.

The result for the free evolution, i.e. no interaction potential between particles, of the system is shown in Fig. 18. The predictions from KFT for the free evolution of particles agree with those from the MD simulations on all scales. The initial structures due to the Rydberg blockade are gradually washed out by thermal motion. The inclusion of an interaction potential between particles is now straightforward in the scope of the KFT formalism. For this system the resummation scheme introduced in Sect. 9 is the most efficient approach to compute the density-fluctuation power spectrum with linear effects. With the tree-level results from Sect. 9, it is already possible to capture all linear effects in the initial power spectrum to infinite order in the interaction potential.

13 Summary and conclusions

Kinetic field theory applies the principles and concepts of a statistical field theory to ensembles of classical particles in or out of equilibrium. The particle ensemble is characterized by the probability distribution of its phase-space coordinates and by the equations of motion. An exact and complete generating functional can be defined such that the statistical properties of the particle ensemble at any later time can then be derived by applying suitable functional derivatives to this generating functional. Up to this point, the theory is applicable to wide classes of classical particle ensembles, or, more generally, to ensembles of classical degrees of freedom.

The focus of this review has been the application of KFT to cosmology. In contrast to other analytic approaches to cosmic structure formation, KFT has the decisive advantage of avoiding the notorious shell-crossing problem by construction. Compared to the structurally similar statistical quantum field theories, the Hamiltonian equations of motion for classical particles are deterministic and, by their symplectic nature, have unit functional determinant. The initial state for the particles in cosmology is empirically well-defined as a statistically homogeneous and isotropic, Gaussian random field. Based on these conditions, applying KFT to cosmology is quite straightforward.

Particle interactions have to be approximated, though. We have described three ways for doing so; first, by expanding the interaction operator into a Taylor series, leading to the equivalent of Feynman rules; second, by applying Born’s approximation relative to suitably chosen inertial particle trajectories; and third, by a resummation scheme that emerges from a reformulation of KFT purely in terms of macroscopic fields. We have shown how non-linear power spectra for the cosmic density and velocity fields can be derived from KFT and how KFT can be used to explain the internal structure of gravitationally-bound, dark-matter structures.

We have further clarified the relation between KFT and the BBGKY hierarchy of more conventional kinetic theory. The application of KFT to fluids is possible by introducing mesoscopic pseudo-particles which, in addition to their microscopic interaction potential, acquire a repulsive potential for mimicking pressure forces. This approach then allows to describe mixtures between dark matter and gas, and we have shown how effects like the baryonic acoustic oscillations can be described in this way. As a further extension of KFT within cosmology which is quite straightforwardly possible, we have given a first example for its application to generalizations of general relativity [68]. Finally, we have shown an application of KFT to a completely different physical system, composed of cold Rydberg atoms, to illustrate the flexibility of the KFT approach.

KFT itself, but in particular its applications to cosmology, are in an early stage of development, and much work remains to be done. The main purpose of this review is three-fold: (1) It should demonstrate how natural the foundations of KFT are and how straightforward its generating functional can be derived. (2) It should sketch formal extensions of the theory, in particular in view of different ways to include particle-particle interactions and resummation schemes. (3) It should describe possibilities for manifold applications of KFT to cosmological problems, including the joint evolution of gas and dark matter, and give an outlook on applications of KFT to systems entirely unrelated to cosmology. Fundamental clarifications of statistical properties of self-gravitating systems are also in reach and have been initiated in terms of fluctuation-dissipation relations [81]. Apart from this, calculating higher-order spectra of cosmological quantities and quantifying the accuracy of KFT are among the next important steps to be undertaken.

Acknowledgements.

First and foremost, we wish to acknowledge the numerous and most valuable contributions of many students to the development of KFT and its applications to cosmology and Rydberg gases, in particular Arne Bahr, Daniel Berg, Moritz Beutel, Sebastian Bieringer, Marc Boucsein, Igor Bröckel, Robin Bühler, Frederick del Pozo, Johannes Dombrowski, Leander Fischer, Marlene Götz, Lukas Heizmann, Patrick Jentsch, Felix Kuhn, Laila Linke, Julius Mildenberger, Robert Ott, Martin Pauly, Christophe Pixius, Thimo Preis, Philipp Saake, André Salzinger, Jan Schneider, Alexander Schuckert, Johannes Schwinn, Christian Sorgenfrei, Marie Teich, and Stefan Zentarra. We are grateful to Manfred Salmhofer for his particularly precious help and advice, to Jürgen Berges, Jan Pawlowski, Matthias Weidemüller, Christof Wetterich, Stefan Hofmann, Adi Nusser, Andreas Wipf, Ruth Durrer, Norbert Straumann, Luca Amendola, and Björn M. Schäfer for clarifying discussions and support. Martin Feix and Sven Meyer contributed creative ideas. The Heidelberg Graduate School for Fundamental Physics, the Centre for Quantum Dynamics of Heidelberg University, and the Max Planck Institute for Astronomy provided substantial financial support, partly from funds of the Excellence Initiative of the German Federal and State Governments. This work was supported in part by the Transregional Collaborative Research Centre “The Dark Universe” of the German Science Foundation.

Appendix

Appendix A: Guide to calculating the power spectrum in the perturbative approach

This appendix is intended as a comprehensive guide to provide the reader with a concise picture of the steps needed to arrive at the power spectrum presented in Fig. 3 of this review. 333Note that all vectors in this appendix will be indicated by an arrow to help avoid possible confusion during implementation.

We begin with (40) and apply

1. (1)

a Taylor expansion of the exponential factor in terms of the interaction operator to first order; and 2. (2)

a Taylor expansion of the exponential factor in the initial conditions (64) in terms of the initial correlations up to second order.

Step 1 results in the following expression for the power spectrum

[TABLE]

To prepare for step 2 we split the momentum correlations (65) into a damping term

[TABLE]

and a term contaning momentum correlations between any two distinct particles,

[TABLE]

The Taylor expansion up to second order of step 2 then results in the following expressions,

[TABLE]

$Z_{0}^{(1)}[\boldsymbol{L},0]$ contains all terms up to first order in initial correlations and $Z_{0}^{(2)}[\boldsymbol{L},0]$ only the second order terms. We can then write the free generating functional as

[TABLE]

Note that we only consider higher-order momentum correlations here and ignore cross terms of the form $QC_{\delta_{j}\delta_{k}}$ and $QC_{\delta_{j}p_{k}}\vec{L}_{p_{k}}$ since terms containing momentum correlations dominate at late times due to the time dependence of the momentum propagator $g_{qp}(\tau):=g_{qp}(\tau,0)$. The initial momentum correlations and, of course, the interaction potential between different particles is responsible for the growth of structures while the damping term is responsible for dissolving those structures. If we take initial momentum correlations as well as interactions into account order by order, we must make sure that the damping term does not enter exponentially into the low-order contributions since this would lead to an overestimation of the damping effect. We shall therefore later approximate the damping term $\exp(-Q_{D}/2)$ consistently at one order less than the term $\exp(-Q/2)$. This implies that damping will only be included in terms of at least second-order in the initial momentum correlations.

To first and second order in the initial correlations, with the improved Zel’dovich propagator $g_{qp}(\tau)$ defined in (61) and suitably approximated damping terms, we thus find the following expressions for the power spectrum in our free theory

[TABLE]

where we have specified for clarity that the power spectra on the right-hand sides are the power spectra $P_{\delta}^{\mathrm{(i)}}$ characterising the initial particle distribution.

In our perturbative approach the interaction potential is given by

[TABLE]

where we have introduced an additional term as compared to (78) which corresponds to an adhesion term following the idea of the adhesion approximation [82, 83, 84]. This approximation was introduced in order to compensate the effect of free streaming in the Zel’dovich approximation once particle trajectories have crossed. The necessity for such a term in the perturbative scheme arises due to the fact, that to first order in the gravitational potential the overshooting of particle trajectories due to the improved Zel’dovich propagator is not adequately compensated, causing a loss of power on small scales, where structures are wiped out. This overshooting is compensated by going to second-order perturbation theory, but has not yet been implemented in our code. In our more advanced schemes to include particle interactions which we present in Sect. 4.2 and 4.3 the adhesion term becomes obsolete, since the choice for the free propagator describing inertial motion and the form of the force term are consistently linked.

In (184), $\bar{\nu}$ is an amplitude with the dimension of a length scale. As a suitable length scale, we choose the velocity dispersion $\sigma_{v}$, propagated to the time $\tau$ by the improved Zel’dovich propagator $g_{qp}(\tau)$,

[TABLE]

We approximate $\bar{\nu}$ by its late-time behaviour, setting $g_{qp}(\tau)\,\sigma_{v}\approx 2$.

With the interaction potential (184), the first-order perturbation contributions to the non-linear power spectrum are

[TABLE]

where we introduced the notations $T:=g_{qp}(\tau)\vec{k}$ and $T^{\prime}:=g_{qp}(\tau^{\prime})\vec{k}^{\prime}$ for brevity. The terms $Z_{\mathrm{A}}$, $Z_{\mathrm{B}}$, $Z_{\mathrm{C}}$, and $Z_{\mathrm{D}}$ have been derived in detail in [42] so that we merely state the results here,

[TABLE]

The terms $P_{\delta}^{(1)}$ and $\delta^{(1)}P_{\delta}^{(1)}$ that are proportional to the initial power spectrum must reproduce the linear growth of the initial power spectrum that can be given as $D^{2}_{+}P_{\delta}^{\mathrm{(i)}}(k)$. We therefore replace the terms linear in the initial power spectrum in (183) and (186) by $D^{2}_{+}P_{\delta}^{\mathrm{(i)}}(k)$. The final result for the non-linear density-fluctuation power spectrum to first order in the interaction potential and to second order in the initial correlations is then given by

[TABLE]

We still need to evaluate the integrals over $\tau^{\prime}$ and $k^{\prime}$ that are left in (191) which is done numerically.

The final ingredients still required are the cosmological model of our Universe and an initial density-fluctuation power spectrum. The cosmological parameters are needed for the computation of the expansion function of the Universe which, in turn, is needed for the growth factor $D_{+}(a)$ that serves as a time coordinate for us. We have, of course, now the complete freedom to take the cosmological model of our choosing; we can choose any initial power spectrum that we like, start our calculations at the appropriate scale factor $a$ at which our initial power spectrum is defined and can then compute the non-linear power spectrum (191) at a desired final scale factor. To arrive at Fig. 3 we made the following choices:

• •

$\Lambda$CDM Universe with $\Omega_{\mathrm{m}}=0.3$, $\Omega_{\Lambda}=0.7$, $\Omega_{\mathrm{b}}=0.04$, $h=0.7$, $\sigma_{8}=0.8$.

• •

The initial power spectrum was generated with CAMB [45] at redshift $z=1100$.

• •

Our final time is chosen to be today (corresponding to $z=0$).

There are no other parameters in KFT.