Dark matter vorticity and velocity dispersion from truncated Dyson–Schwinger equations

. Large-scale structure formation is studied in a kinetic theory approach, extending the standard perfect pressureless fluid description for dark matter by including the velocity dispersion tensor as a dynamical degree of freedom. The evolution of power spectra for density, velocity and velocity dispersion degrees of freedom is investigated in a non-perturbative approximation scheme based on the Dyson–Schwinger equation. In particular, the generation of vorticity and velocity dispersion is studied and predictions for the corresponding power spectra are made, which qualitatively agree well with results obtained from 𝑁 -body simulations. It is found that velocity dispersion grows strongly due to non-linear effects and at late times its mean value seems to be largely independent of the initial conditions. By taking this into account, a rather realistic picture of non-linear large-scale structure formation can be obtained, albeit the numerical treatment remains challenging, especially for very cold dark matter models.


Introduction
One of the primary goals of contemporary cosmology is to describe the evolution of dark matter under the influence of gravity, in order to understand the observed large-scale structure of the Universe.As modern and future observations can probe increasingly smaller scales, there is a genuine need to understand the dynamics governing the relevant physical processes at these scales.
In kinetic theory, dark matter is described by a phase-space distribution function which evolves with Boltzmann's equation.At cosmological late times, microscopic interaction rates are negligible and the distribution function's evolution is sufficiently well described by the (collisionless) Vlasov equation [1].On scales much smaller than the Hubble horizon, gravity around an expanding, homogeneous and isotropic background cosmology is, in leading order approximation, described by Poisson's equation.The so-called Vlasov-Poisson system of equations is a closed system of non-linear integro-differential equations in seven variables and, as such, is rather difficult to solve in full generality.From an observational point of view, one is often interested in velocity moments or cumulants instead of the full phase-space distribution function.The Vlasov-Poisson system of equations can be cast into an infinite hierarchy of coupled evolution equations for the distribution function's moments or cumulants [2].
In order to obtain a closed and finite set of equations, a common approach is to truncate the distribution function's infinite cumulant expansion after the first order.This often used single-stream approximation models dark matter as a perfect pressureless fluid and is preserved by the Vlasov-Poisson system of equations, in the sense that no higher order cumulants are sourced if initially absent.However, it is an apparent self-consistency due to the phenomenon of shell-crossing.During gravitational collapse dark matter particles meet in position space and generate a non-trivial velocity dispersion tensor which indicates the break down of the single-stream approximation.While the perfect pressureless fluid description provides a surprisingly good account of early stage gravitational instabilities, one naturally longs for a description that is capable to capture the physics of a multi-stream flow.Ultimately, this relies on including some notion of velocity dispersion.To this end, various approaches have been developed, ranging from a description including the velocity dispersion tensor [3][4][5][6][7], and more recently even higher-order cumulants [8,9], through functional approaches, either over the phase-space distribution function [10] or the Lagrangian displacement field [11], to using the Schrödinger method to tackle the cumulant hierarchy beyond a truncation [12,13], to name only a few.
In this work a truncation of the cumulant expansion after the second order is pursued, which includes the velocity dispersion tensor as dynamical degrees of freedom.Within this truncation, the single-stream approximation's degenerate distribution function is smoothened out, which allows to capture the average motion of a multi-stream flow.This should be understood in terms of an effective macroscopic description which models dark matter in a fluid dynamical approach with some notion of 'warmness'.Moreover, the inclusion of velocity dispersion is naturally paired with non-trivial dynamics for the velocity's rotational mode.Vorticity is not sourced if primordially absent in the single-stream approximation, but is known to be generated at shell-crossing [14].
A common approach to study cosmic structure formation is standard perturbation theory.Formally assuming that the deviations from a homogeneous and isotropic background cosmology are small, the dynamics are solved around linear theory and non-linear corrections are organised in an expansion in the initial correlation of mass density fluctuations.For sufficiently early times and large scales, perturbation theory captures the relevant dynamics, whereas at later times and smaller scales the perturbative series badly fails to converge since the deviations from the background become large.In order to access the non-linear regime various approaches have been proposed, including partial resummations of the perturbation series [15][16][17][18][19][20][21][22][23][24][25], effective theories [26][27][28][29] or renormalisation group approaches [30][31][32][33][34][35].To this end, this work studies a non-perturbative approach based on the Dyson-Schwinger equations of the underlying statistical field theory.An approximation involving cosmic one-and twopoint correlation functions is studied for the cosmological fluid.However, this approach can easily be systematically extended as is done for example in quantum field theories.
The paper is structured as follows.In Section 2 the Vlasov-Poisson system of equations is presented and the distribution function's cumulant expansion reviewed.Section 3 revises the Martin-Siggia-Rose/Janssen-de Dominicis formalism for Langevin type dynamics with Gaussian random field initial conditions.The corresponding Dyson-Schwinger equations are derived and truncated in order to obtain a closed set of evolution equations for one-and twopoint correlation functions.Section 4 utilises the truncated Dyson-Schwinger equations in a variant of Hartree's approximation to evolve cosmic correlation functions.This is extended to a self-consistent Hartree-Fock approximation in Section 5, an approach which features genuinely non-perturbative physics.Some conclusions are drawn in Section 6. Appendix A provides details about the employed stochastic response field formalism and Appendix B collects explicit expressions for bare vertices.

Gravitational dynamics of dark matter and statistical description
The gravitational dynamics of an ensemble of non-relativistic and collisionless point particles of mass  is described by the Vlasov-Poisson system of equations [1], Here,  is conformal time,  are comoving coordinates and  =  d/d are the corresponding conjugate momenta. 1 The former two are related to cosmic time and proper coordinates by d = ( ) d and  = ( ) , respectively, where ( ) is the scale factor, parametrising the expansion of the spatial Universe.Further, H = d/d is the conformal Hubble rate and  m is the (time-dependent) dark matter density parameter.Finally, ( , , ) is the dark matter one-particle phase-space distribution function, here normalised to and ( , ) is the Newtonian gravitational potential. 2 When studying late-time cosmic structure formation, one is often less interested in the full phase-space distribution function but rather in velocity moments or cumulants thereof, the full set of which completely characterises ( , , ).The first few velocity moments are given by the (background-normalised) dark matter mass density, the momentum density, and the stress tensor, here expressed in terms of the density contrast field ( , ), quantifying the local deviation from the mean dark matter mass density, the velocity vector field   ( , ) and the velocity 1 Symbols with indices from the middle of the Latin alphabet refer to spatial vector components, while boldface symbols denote the corresponding vector.Indices appearing twice in a single term are summed over, where vectors and covectors are not distinguished and are both denoted with subscripted indices, as is common in flat space.Partial derivatives with respect to the th comoving coordinate component are abbreviated by   , while all other partial derivative are denoted by a suitably subscripted .Finally, function arguments are often suppressed for the sake of brevity. 2 The expectation value ⟨ ⋅ ⟩ is either taken as an ensemble average over cosmic histories with stochastic initial conditions or as a sample average over large spatial volumes in a single cosmic history.
dispersion tensor field   ( , ).The latter two are the first-and second-order velocity cumulants of the distribution function, respectively.
Instead of dealing with the momentum dependence of the distribution function, the Vlasov-Poisson system of equations (2.1) can be cast into an infinite tower of coupled evolution equations for the velocity moments or cumulants [2].For practical purposes, one often turns to approximations describing dark matter by a finite number of velocity cumulants to obtain a closed system of equations.The most prominent example is the so-called single-stream approximation, which models dark matter as a perfect pressureless fluid that is described by the density contrast field  and velocity field   only [36].
Although the single-stream approximation successfully describes the gravitational dynamics of dark matter at early times and large scales, it fails to do so at later times and smaller scales.Most notably, the single-stream approximation cannot account for the phenomenon of shell-crossing, when particle trajectories meet in position space, since higher-order velocity cumulants are naturally generated, which indicates the break-down of the single-stream approximation [14].
In the following, the description of dark matter is extended beyond the single-stream approximation by including the velocity dispersion tensor   as a dynamical degree of freedom to describe gravitational dynamics also after shell-crossing.On the level of the phase-space distribution function, this corresponds to smoothening out the degenerate single-stream approximation distribution function to a normal distribution, where the velocity dispersion tensor is the covariance matrix of the spatial momenta [7,35].Naturally, this approximation cannot describe shell-crossing microscopically, but it supports the average motion of a multistream fluid [3] and captures signatures associated to shell-crossing, as is shown a posteriori.While a truncation of the cumulant expansion after the second-order is not self-consistent, in the sense that higher-order velocity cumulants are dynamically generated even if initially absent [14], one can argue that this approximation could provide a viable description of virialised clumps of dark matter, which seem to be near to a non-relativistic thermal state, at least for simple halo models [37].
In the following, higher-order velocity cumulants are neglected cum grano salis and dark matter is described in terms of an effective theory with ten matter degrees of freedom.These are parametrised in terms of the density contrast field , the velocity vector field   and the (symmetric) velocity dispersion tensor field   .The relevant evolution equations are given by the continuity equation, the Cauchy momentum equation, and the velocity dispersion equation, which together with Poisson's equation, form a closed system of equations. 3n cosmology it is common to adopt a statistical description for dark matter, where the relevant degrees of freedom are described in terms of random fields.This is most often motivated by a very early inflationary epoch of the Universe, which predicts the initial state of the theory to be very well described by Gaussian random fields that are statistically homogeneous and isotropic in space [38].The fair sample hypothesis asserts that the employed statistical field theory can be understood as describing an ensemble of cosmic histories with stochastic initial conditions or a sample of large spatial volumes within a single cosmic history [39].
To investigate the statistical properties of cosmic large-scale structures, one is interested in correlation functions of the velocity moments and cumulants, or more generally the phase-space distribution function.Since the underlying random fields are homogeneous and isotropic in space, it is useful to decompose the matter degrees of freedom with respect to spatial translations and rotations, which is most conveniently done using the Fourier transform. 4 For the velocity vector field the decomposition can be parametrised as where i is the imaginary unit and   is the totally antisymmetric Levi-Civita symbol.The scalar velocity divergence mode () parametrises the conservative part of the velocity field, while the pseudovector vorticity mode   () quantifies the solenoidal part.Similarly, the velocity dispersion tensor field can be decomposed as where k  =   / and   is the Kronecker delta.The scalar trace mode () parametrises isotropic velocity dispersion, while the scalar mode (), the transverse vector mode   () and the transverse traceless tensor mode   () parametrise anisotropic velocity dispersion degrees of freedom.The decompositions (2.10) and (2.11) correspond to splitting the velocity and velocity dispersion degrees of freedom 3 → 2 + 1 and 6 → 2 + 2 + 1 + 1, respectively.Due to the statistical symmetries, the one-point correlation functions (mean fields) are given by ⟨( , )⟩ = 0 , ⟨  ( , )⟩ = 0 , ⟨  ( , )⟩ =   σ ( ) , (2.12) 4 The splitting into Fourier modes is a decomposition into irreducible unitary representations of the translation group and for any suitable function or distribution ℎ() the convention is employed, where integrals are abbreviated as The Fourier transform is denoted by the same symbol and distinguished by its argument if not clear from context.The three-dimensional Euclidean inner product is written as  ⋅  =     and the modulus of threevectors is denoted by the corresponding lightface symbol, such as the wave number  = ||.Finally, wave vector Dirac delta functions are often denoted by δ() = (2) where σ is a (time-dependent) isotropic velocity dispersion background. 5For thermally produced dark matter σ is related to the particle mass and decoupling temperature, but generally depends on the type of candidate considered [40].Two-point correlation functions (covariance functions) are of natural importance in cosmology due to the Gaussian nature of the initial state.In terms of the Fourier transform these can be quantified by their power spectral density, e.g. the matter power spectrum, ⟨( , ) ( ′ ,  ′ )⟩ = δ( +  ′ )   ( ,  ′ , ) . (2.13) Due to the statistical symmetries, only covariance functions of the same type of modes (scalar, transverse vector or transverse traceless tensor) are non-vanishing.As such, the transverse vector mode power spectra carry an additional projector P  () =   − k  k  , e.g. the vorticity power spectrum, and similar for the transverse traceless tensor mode power spectra.
Although late-time cosmic structure formation is naturally non-Gaussian since the evolution equations are non-linear, the power spectra are usually the objects of main interest in cosmology.However, higher-order correlation functions, such as the bi-and trispectrum, are also frequently studied and provide insights into deviations from Gaussian statistics of dark matter.

Functional approach
To study correlation functions, such as mean fields and covariance functions, it is convenient to work in a field-theoretic functional formulation [41][42][43], some aspects of which are recapitulated in Appendix A.
The focus of this work is to investigate the evolution of the four scalar and the vorticity vector mode, while the velocity dispersion vector and tensor modes are neglected here. 6The decomposed degrees of freedom are assembled in the field multiplet where the index  runs through the field content, is summed over if appearing twice in a single term and carries any additional tensorial substructure of the fields.In the following, it is useful to work in terms of the time evolution parameter  = ln( + ), where  + is the 5 The velocity dispersion mean field naturally feeds into Friedmann's equations that govern the evolution of the scale factor (), but the corresponding contribution is of the order σ / 2 , with  the speed of light, and is thus negligible for cold enough dark matter. 6The inclusion of velocity dispersion degrees of freedom allows to source the vorticity field non-linearly, whereas the vorticity field decays linearly in the single-stream approximation due to the Hubble expansion and cannot be sourced if initially absent [36].The velocity dispersion vector and tensor degrees of freedom are neglected here due to computational limitations.Their impact is expected to be subdominant compared to the scalar velocity dispersion degrees of freedom, as the vector and tensor modes can only be sourced non-linearly (if initially absent) from scalar modes due to the structure of the non-linear terms in the velocity dispersion equation (2.8).
fastest growing linear density contrast mode in the single-stream approximation. 7Further, it is sensible to introduce  + = d ln( + ) / d ln(), quantifying the deviation of  + from the EdS linear growth function. 8Finally, it is convenient to introduce the rescaled mean field σ = σ /( + H) 2 , such that ⟨  (, )⟩ =  3 δ()  2 σ (). 9 In terms of the field multiplet (3.1), the bare equations of motion (2.6), (2.7) and (2.8) can be written as where the linear dynamics is specified by the matrix after eliminating the gravitational potential using Poisson's equation (2.9).Here and in the following, the vector mode sector of objects with a matrix structure in field space is implicitly understood to carry an appropriate transverse projector.The non-linear dynamics are characterised by   , which can be formally written in a vertex expansion as with the (symmetric) vertices   ( 1 ,  2 ) =   ( 2 ,  1 ), some of which are explicitly given in Appendix B. 10  In the context of cosmology, one is interested in expectation values of composite fields of the field content (3.1) which obey the equations of motion (3.2) for Gaussian random initial conditions.As such, the initial state is completely characterised by the initial velocity dispersion mean field σ in and the power spectra  in  ().Due to statistical isotropy and homogeneity, the mean field is a zero mode and the power spectra depend on one wave number only, up to possible projectors for vector and tensor modes.Since the bare equations of motion (3.2) are non-linear, the fields naturally evolve away from their initial Gaussian shape and are completely characterised by an infinite set of correlation functions.
where 2  1 is the Gaussian hypergeometric function and  = 1/ m ( 0 ) − 1, with  0 the conformal time corresponding to today.Here, the growth function is normalised to  + ( 0 ) = 1 and in the Einstein-de Sitter (EdS) limit  m → 1 one obtains  + = . 8The approximation  m / 2 + = 1 allows to map a CDM to an EdS cosmology, where the linear equations of motion are time-translation invariant and often admit analytical solutions. 9Since the physical part of the mean field is σ , zero modes need to be treated cautiously and expressions of the form δ()  2  (, ) are only to be evaluated under an integral. 10The higher-order terms of the vertex expansion are due to the term     ln(1 + ) appearing in equation (2.7) and are formally only justified for  < 1.In a field-theoretic treatment, the higher-order vertices involve higher loop orders which is beyond the approximations investigated in Section 3.2.
In the functional formulation the system is characterised by the bare action where the first line characterise the (bare) dynamics, while the second line entails the Gaussian statistics of the initial state.
From the generating functional of connected correlation functions , one obtains the mean fields at vanishing source currents from while connected two-point functions are given by Here,   (,  ′ , ) is the power spectrum of two fields and  R  (,  ′ , ) is the (retarded) mean linear response function, in cosmology often called propagator, to which the advanced counterpart is related by In the following, it is convenient to also work with the one-particle irreducible (1PI) effective action , which is the generating functional of 1PI correlation functions.At vanishing source currents, the 1PI one-point functions are while the 1PI two-point functions take the form ) . (3.9) Here, () is the effective equation of motion for the mean field σ (),  A  (,  ′ , ) and  R  (,  ′ , ) are the inverse advanced and retarded propagators, respectively, and   (,  ′ , ) is the statistical 1PI two-point correlation function.
At vanishing source currents, the effective equation of motion reads the propagator equation is and the power spectrum equation is Equations (3.10), (3.11) and (3.12) are exact and although they cannot be straightforwardly solved since the 1PI one-and two-point functions are not explicitly known, they provide a generically non-perturbative way to compute the mean field and connected two-point correlation functions.
To investigate approximation schemes, it is sensible to split the 1PI one-and two-point functions into a linear contribution and a non-linear correction, where the former can be directly read off the bare action (3.5) and the latter parametrises the ignorance of the 1PI functions.For the effective equation of motion of the velocity dispersion mean field this amounts to where the second and third terms on the right-hand side are the Hubble drag term and the initial condition, respectively.The non-linear correction is parametrised in terms of the source (), which quantifies the back-reaction of fluctuations onto the mean field.Similarly, the 1PI two-point functions (3.9) can be decomposed as where the non-linear corrections are parametrised in terms of the self-energies  R  (,  ′ , ) and   (,  ′ , ).Further, it is convenient to split the inverse propagator self-energy into parts that are local and non-local in time, which in the following are referred to as Hartree and Fock self-energies, respectively.Since the statistical self-energy   (,  ′ , ) is naturally non-local in time, it is also regarded as a Fock-type self-energy.Eliminating the advanced in favour of the retarded propagator and plugging the expressions (3.13), (3.14) and (3.15) into the evolution equations (3.10), (3.11) and (3.12), one obtains the effective equation of motion the propagator equation and the power spectrum equation Up to here, no approximations have been made and equations (3.16), (3.17) and (3.18) are exact.The ignorance of the 1PI one-and two-point functions has simply been absorbed into the source  and the self-energies  H  ,  F  and   .Given explicit expressions for these, one can solve the evolution equations to obtain the mean field and connected twopoint correlation functions.However, in general it is not possible to provide exact and closed expressions for the source and self-energies since these typically involve higher-order correlation functions.Nonetheless, the expansion in terms of connected and 1PI correlation functions naturally allows for non-perturbative investigations with functional methods, such as the Dyson-Schwinger equations [19,42,43] and the renormalisation group [31][32][33][34][35], the former of which are studied in the next section.

Dyson-Schwinger equations
To obtain explicit expressions for the source and self-energies, the Dyson-Schwinger equations [45][46][47] of the statistical field theory are investigated.In the current context it reads where the dot product runs over all internal structures, namely the physical-response field content, time and space arguments.The Dyson-Schwinger equation (3.19) is an exact equation, which relates the 1PI onepoint correlation function to its bare counterpart, i.e. the first functional derivative of the bare action, evaluated on the field configuration given in the argument on the right-hand side.Equation (3.19) can be derived from the invariance of the functional integral measure with respect to additive shifts of the fields.For the type of bare action (3.5) the Dyson-Schwinger equation (3.19) in general involves all functional derivatives of the bare and effective action. 12his is due to the fact that the term   ln(1 + ) involves vertices of all orders in a vertex expansion.In the following, only the bare three-point vertex is taken into account and the terms originating from higher-order vertices are neglected.In this case the Dyson-Schwinger equation (3.19) depends on the third-order functional derivative of the bare action  (3) , that is up to normalisation the bare three-point vertex   , as well as the second-order functional derivative of the effective action  (2) .
By applying functional derivatives to the Dyson-Schwinger equation (3.19) one obtains equations for higher-order 1PI correlation functions, which in turn depend on 1PI correlation functions up to the next higher order, creating an infinite hierarchy of coupled equations.More explicitly, the Dyson-Schwinger equation for the 1PI one-point function reads while the equation for the 1PI two-point function is given by where the dot product and the trace run over all internal structures.The inverse of the second-order functional derivatives of the effective action is expressed in terms of the connected two-point correlation functions by virtue of equation (A.7) and should be understood as being field-dependent through the Legendre transform.Similarly, one can proceed for higher-order 1PI correlation functions to obtain an infinite tower of coupled equations.The Dyson-Schwinger equations (3.20) and (3.21) can now be used to assign explicit expressions to the source and self-energies.At vanishing source currents, the mean field source reads or equivalently using the diagrammatic representations introduced in Appendix A, From this expression it is evident that the mean field is sourced by the equal-time two-point cross-correlation of velocity and velocity dispersion fluctuations. 13he Hartree self-energy arises from evaluating  (2) [ ] at non-vanishing mean field and is given by or diagrammatically as and parametrises the presence of a non-vanishing mean field (velocity dispersion) onto the response of fluctuations.In contrast, the Fock self-energies contain loop expressions as parametrised by the trace term and read diagrammatically and The explicit formulae for the Fock self-energies are more involved and are not explicitly displayed in full form here, but are instead given later within the approximations used to close the Dyson-Schwinger hierarchy.

Hartree approximation
A simple truncation of the Dyson-Schwinger hierarchy is the Hartree approximation, where the system of equations is closed by neglecting the Fock self-energies,  F  = 0 and   = 0.
The equations of motion are then modified by the presence of the time-dependent velocity dispersion mean field and the propagator equation (3.17) reads Here, the velocity dispersion mean field acts very similar to the pressure term in Jeans instability through the Hartree self-energy (3.24), such that growing solutions are only allowed for wave numbers smaller than the free-streaming wave number [7]  fs () = √ 3 2 The power spectrum equation (3.18) is formally solved by due to the absence of the statistical self-energy.In the following, it is convenient to write the initial power spectrum as where  in  () is the initial density contrast power spectrum and  1 () = 1.To investigate the response of fluctuations, it is convenient to define the reduced propagator quantifying the mean linear response from initial conditions   ().

Linear mean field
To gain a better understanding of the physics captured in the Hartree approximation, it is sensible to first study a situation that can be tackled analytically.For a vanishing source term, the velocity dispersion mean field decays linearly with the Hubble expansion and one diagrammatically obtains From a physical point of view, the fluctuations evolve in the presence of a mean field, which does not feel the back-reaction of the fluctuations.Employing the approximation  m / 2 + = 1 and assuming initial conditions that are in the growing mode of the single-stream approximation,   = (1, 1, 0, 0, 0), the propagator equation (4.1) can be rewritten using the reduced propagator (4.5) to arrive at a single equation for the reduced density propagator, with the linearly decaying mean field σ () = σ in e −(− in ) .The general solution of this thirdorder differential equation can be constructed from the three independent solutions where 1  2 (; , ; ) are generalised hypergeometric functions.
The solutions  g and  d smoothly connect to the standard growing and decaying modes of the single-stream approximation in the limit  2 σ → 0, where the hypergeometric functions are unity, while the solution  * d is a new decaying mode due to the velocity dispersion degrees of freedom.Figure 1 displays that the growing mode is enhanced, while the decaying modes are suppressed and oscillate with increasing frequency at larger  2 σ .
It should be emphasised that even in this rather simple approximation, the propagator equation (3.11) is no longer time-translation invariant due to the time-dependent velocity dispersion mean field and therefore  R  (,  ′ , ) ≠  R  ( −  ′ , 0, ).

Sourced mean field
To capture the back-reaction of fluctuations onto the mean field, a non-vanishing source for the velocity dispersion mean field needs to be included.Diagrammatically this approximation reads That is, the velocity dispersion mean field is sourced by the back-reaction of fluctuations, which in turn evolve in the presence of the mean field.In this case, the propagator equation (3.17) has to be solved together with the mean field equation (3.16).The solution heavily depends on the initial conditions, especially on the initial velocity dispersion mean field [7,30].In the following, numerical solutions for different initial conditions σ in are investigated, while fluctuations are initialised in the single-stream approximation growing mode,   () = (1, 1, 0, 0, 0).In Figure 2, numerical solutions of the velocity dispersion mean field for different choices of σ in are shown.One finds that the smaller the initial velocity dispersion mean field is, the stronger it grows at later times.This is due to the fact that the power spectra, which source the mean field through the term (3.22), are less suppressed in the ultraviolet the smaller the value of the velocity dispersion mean field.Indeed, going to very small initial mean fields causes a rather involved numerical treatment, due to strongly suppressed and oscillating power spectra.
This can be understood from the wave number dependence of the reduced propagators shown in Figure 3, which directly determine the shape of the power spectra through equation (4.4).At larger wave numbers, the propagators oscillate with increasing frequency and decaying amplitude, where the scale is set by the free-streaming wave number (4.2), very similar to the decaying solutions shown in Figure 1.The velocity dispersion-velocity cross-spectra, sourcing the mean field through the term (3.22), increasingly extend into the ultraviolet for colder dark matter models with a larger free-streaming wave number.This leads to the observed strong growth for the colder candidates displayed in Figure 2.
The drawback of this approximation is that the source of the mean field is sensitive to the ultraviolet of power spectra, a regime that is not well described by linear theory.To obtain a realistic source and mean field, non-linearities for the evolution of power spectra need to be taken into account with an intrinsically non-perturbative approximation.

Hartree-Fock approximation
A more sophisticated approximation is obtained by closing the Dyson-Schwinger hierarchy by setting the 1PI three-point functions to its bare form, while keeping the full 1PI twopoint functions.This is also referred to as the self-consistent one-loop or Hartree-Fock approximation, in analogy to approximations of this type employed for quantum systems.
Because the bare vertex   couples one response to two physical fields, only the first two diagrams of equation ( 3 Reduced propagators for the density contrast (blue curves), velocity-divergence (red curves), isotropic velocity dispersion (yellow curves) and anisotropic velocity dispersion (violet curves) at redshift  = 0, normalised to the standard growing mode as a function of  2 σ in .Results are shown for two choices of initial velocity dispersion mean fields, 1 km 2 /s 2 (solid curves) and 10 −5 km 2 /s 2 (dashed curves).The damping scale responsible for the oscillations is set by the time-dependent mean field, shown in Figure 2, and is at much smaller  2 σ in for the colder dark matter model due to the strong late-time growth of the velocity dispersion mean field.
Hartree-Fock approximation.Diagrammatically one then has and More explicitly, the retarded Fock self-energy is given by With these expressions, equations (3.16), (3.17) and (3.18) form a closed system.

Large external wave number limit
Before turning to full numerical results, it is instructive to consider the limit of large external wave numbers in order to get an analytical understanding of the Hartree-Fock approximation.
In the limit of large external wave numbers, the propagator equation can be simplified to an extent where it can be solved in the absence of a mean field.
The inclusion of the retarded Fock self-energy  F  (,  ′ , ) introduces a new scale dependence in the propagator due to the non-linear coupling of modes.As long as the velocity dispersion degrees of freedom are comparably small, one expects the dominant non-linear contribution to be due to the random advection of density structures by a large-scale velocity field, also known as the sweeping effect [35,48,49].To isolate this contribution and investigate how it appears in the Hartree-Fock approximation, the large external wave number limit is studied in the absence of velocity dispersion degrees of freedom.
In the limit  → ∞, the leading-order contribution to the retarded Fock self-energy where (5.6) Notice that in the equal-time limit which is proportional to the mean-square velocity      (, , 0)/3 up to the factor ( + H) 2 , where      (,  ′ , | − |) is the covariance function of two velocity fields. 14Since  (,  ′ ) depends on the full velocity power spectrum, the general solution involves solving the power spectrum equation (3.18).To avoid this complication, the full velocity power spectrum may be approximated by the linear one.In this case one obtains where is often interpreted as a one-dimensional velocity dispersion in the sense of a mean-square velocity.
Consider the approximation where the propagator in the Fock self-energy (5.5) is evaluated at linear level.This is the 1PI one-loop approximation, which can be solved analytically [43], where  R  (,  ′ ) is the linear retarded propagator.In the case where the full propagator is kept, one obtains [19,43] where  1 is the first-order Bessel function of first kind.
Let us also compare this to results first obtained in the framework of renormalised perturbation theory [15,16] and thereafter also using renormalisation group methods [31,35].In the large external wave number limit one obtains a resummed propagator of the form with a very similar approximation using a linear instead of the full power spectrum. 15he large external wave number limit propagators (5.10), (5.11) and (5.12) can all be written as the linear retarded propagator multiplied by a correction factor that only depends on the dimensionless combination  v [e − in − e  ′ − in ] and are shown in Figure 4.
The propagator obtained from the functional renormalisation group calculation (5.12) is suppressed by a Gaussian factor while the 1PI one-loop and Hartree-Fock approximation propagators feature oscillations which in the latter case are suppressed at larger  v [e − in − e  ′ − in ].The characteristic scale of suppression respectively the frequency of the oscillations is determined by the non-linear wave number or more precisely by the difference of the inverse of two non-linear wave numbers at the two times of the propagator.The difference of the propagators is due to the different resummation schemes underlying the approximations in which they were derived.While the three methods are all non-perturbative in the sense that the perturbative series is resummed to infinite order, they correspond to different infinite partial resummations.This is best understood in the language of renormalised perturbation theory [15], where the propagator (5.12) can be obtained as a systematic resummation of the perturbative series.
The Gaussian suppression factor entering the propagator (5.12) can be written as the series where  =  v [e − in − e  ′ − in ] and !! is the double factorial.Here, (2ℓ − 1)!! are the number of diagrams contributing to the propagator at perturbative loop order ℓ while the factor 1/(2ℓ)! is due to all time integrations within a diagram [16].It is emphasised that these are not all diagrams, but rather only a subclass of diagrams that are assumed to be dominant in the large external wave number limit.This is also clear from functional renormalisation group calculations, where the propagator (5.12) was obtained in an approximation [35].As such, it corresponds to an infinite partial resummation of the perturbative series.
Applying the same expansion to the propagators (5.10) and (5.11), one can find the number of diagrams contributing at each loop order.For the 1PI one-loop approximation the expansion of the cosine reads while for the Hartree-Fock approximation one finds The number of diagrams contributing at each loop order to the different resummation schemes is summarised in Table 1.It is noteworthy that in the direct-interaction approximation similar results for the propagator hold in turbulence [50,51].
The different partial resummations can also be obtained from topological arguments.Consider the contributions to the propagator contained in the resummation of renormalised perturbation theory up to two-loop order.At zero-and one-loop order there is only a single diagram which is captured by all resummation schemes.At two-loop order though there are the three diagrams .
Based on the large external wave number limit, two approximation schemes interpolating between the perturbative small wave number and non-perturbative large wave number sector have been developed.In renormalised perturbation theory the propagator (5.12) is used for the interpolation [16], while in the direct-interaction approximation the propagator (5.11) is utilised [19].In the next section the numerical solution of the full propagator equation (3.17) is compared to these interpolation schemes, which in the following are referred to as 'CS' and 'TH' approximations, respectively.

Numerical solutions
In this section the Hartree-Fock approximation is solved with numerical methods and compared to -body simulations.In the following, it is convenient to work with dimensionless power spectra, ∆  (,  ′ , ) =  3 2 2   (,  ′ , ) . (5.18)

N -body simulations
Since -body simulations can only simulate a representative finite region of the Universe, one typically has to decide whether one wants to accurately reproduce large-scale structures, requiring a fairly large simulation volume, or is interested in resolving small-scale physics, which requires a rather high number density of particles.To test the performance of the Hartree-Fock approximation, numerical solutions are compared to the Horizon Run 2 (HR2) -body simulation [56], featuring a rather large simulation volume, as well as to a -body simulation of Buehlmann & Hahn (BH19) [59], featuring a comparably high number density of particles.The cosmological parameters and codes used for the simulations as well as the number of particles  p , the box side length  box and the initial redshift  in are summarised in Table 2. Three snapshots of the BH19 -body simulation were kindly provided by the authors M. Buehlmann and O. Hahn [59].The datasets contain estimates for the density, velocity and isotropic velocity dispersion, interpolated on a 1024 3 -mesh at redshifts  = 2.165,  = 0.994 and  = 0.The density field was estimated using the cloud-in-cell deposition algorithm [63], while the velocity and velocity dispersion fields were obtained from the Lagrangian tessellation method [64,65].
To derive an estimate for the power spectra, the fields are discretely Fourier transformed on a three-dimensional grid with fundamental wave number  f = 2/ box and  grid = 1024 grid points per dimension.Further, the density field is deconvolved using the cloud-in-cell window function [66] where   ∈ {1, … ,  grid }.The power spectrum estimate P  () is obtained from the spatial average of two fields over the simulation volume and an angular average of modes within spherical shells of thickness  f and radius  =  f for  ∈ {1, … ,  grid }.Finally, the statistical error is estimated as [67] where the first term on the right-hand side is the sample variance, while the second term is the Poisson shot noise.It is emphasised that the power spectrum and its error estimate should be taken with caution for small wave numbers due to systematic uncertainties related to the relatively small box size of the BH19 -body simulation.
Similarly, an estimate of the velocity dispersion mean field can be computed from the spatial average of the trace of the velocity dispersion tensor.

Numerical methods
The system of equations (3.16), (3.17) and (3.18) in the Hartree-Fock approximation (5.3) and (5.4) is numerically solved in time using an adaptive Runge-Kutta-Cash-Karp method [52][53][54] for Volterra integro-differential equations [55] and in wave number using a finite element method with B-splines of order unity as basis functions [20,43].To test the Hartree-Fock approximation against the HR2 and BH19 -body simulations, two sets of numerical solutions are computed.These account for the different cosmological parameters and initial redshifts listed in Table 2.
The initial power spectrum is generated on scales 10 −5 ℎ/Mpc ≲  ≲ 20 ℎ/Mpc using the class code [62] and extrapolated into the ultraviolet with the fitting formula  log() 2   s −4 , where  s is the scalar spectral index.To ensure comparability to the HR2 and BR19 -body simulations, wave numbers are interpolated on a grid with an infrared cut-off  min =  f and an ultraviolet cut-off  max =  Ny equal to a multiple of the Nyquist wave number  Ny =  grid / box .For the HR2 simulation these are given by  f ≈ 8.7 ⋅ 10 −4 ℎ/Mpc and  Ny ≈ 2.6 ℎ/Mpc, while for the BH19 simulation these are  f ≈ 2.1 ⋅ 10 −2 ℎ/Mpc and  Ny ≈ 11 ℎ/Mpc.More specifically, when comparing numerical solutions to the HR2 simulation, the cut-off is chosen as  max ≈ 5.2 ℎ/Mpc, while when comparing to the BH19 simulation the cut-off is varied between  max ≈ 11 ℎ/Mpc and  max ≈ 43 ℎ/Mpc. 16or the single-stream approximation a wave number grid with   = 100 interpolation points is used, while for the full set of fields (3.1) a grid with   = 50 is chosen.While it would be desirable to increase the number of interpolation points, it considerably slows down the solving algorithm since the arrays storing the Fock self-energies (5.3) and (5.4) scale cubic in   .17

Single-stream approximation
To gain a better understanding of the properties of the Hartree-Fock approximation, the equations (3.17) and (3.18) are solved in the absence of a mean field in the single-stream approximation, where the field content (3.1) reduces to the density and velocity-divergence only.
The density contrast propagator and power spectrum have also been studied using the two-particle irreducible method [42,43,48,68] and closure theory [19][20][21], which yield the same evolution equations as the Hartree-Fock approximation in the single-stream approximation.The purpose of the following is to identify how well the Hartree-Fock approximation captures non-linearities and to find out where the single-stream approximation fails.
Propagator.In Figure 5, the time (upper panel) and wave number (lower panel) dependencies of the density contrast and velocity-divergence reduced propagators are shown at constant wave number  ≈ 2 ℎ/Mpc and at fixed redshift  = 0, respectively.The numerical solution of the Hartree-Fock approximation is near to identical with the TH approximation.A slight deviation for wave numbers  ≳  nl () is observed, which is not too surprising since the TH approximation interpolates between the small and large wave number sector.As was already discussed in Section 5.1, the observed oscillations are not associated with real physical phenomena but are rather due to the partial resummation of the full perturbative series.
Density contrast equal-time auto-spectrum.In Figure 6, the dimensionless density contrast equal-time auto-spectrum, normalised to the linear spectrum, is shown at redshifts  = 2 (upper panel) and  = 0 (lower panel).At large wave numbers, the Hartree-Fock approximation does not overestimate the -body power spectrum as bad as the standard perturbation theory one-loop prediction.While the deviation from the -body power spectrum is fairly small at redshift  = 2, it grows towards redshift  = 0. Especially for wave numbers  ≲ 0.2 ℎ/Mpc the Hartree-Fock approximation seems to follow the standard perturbation theory one-loop prediction.Since the full 1PI three-point functions are set to their bare form in the Hartree-Fock approximation, it is likely that vertex corrections are needed to accurately capture the mildly non-linear regime.
In Figure 7, the dimensionless density contrast equal-time auto-spectrum is shown for a larger range of wave numbers at redshift  = 0.The Hartree-Fock approximation shows a much better convergence towards the -body power spectrum compared to the standard Single-stream approximation reduced density contrast (blue curves) and velocitydivergence (red curves) propagator time dependence at  ≈ 2 ℎ/Mpc (upper panel) and wave number dependence at  = 0 (lower panel).The numerical solution of the Hartree-Fock (HF) approximation (solid curves) is compared to the TH (dashed curves) and CS (dotted curves) approximation.The Hartree-Fock and TH approximations are near to identical since the latter is based on the large wave number limit of the former.perturbation theory one-loop prediction, which badly overestimates the power spectrum.The drop of the -body power spectrum at wave numbers near the Nyquist wave number is due to finite box size effects and should therefore not be trusted.
In Figure 8, the dimensionless density contrast equal-time auto-spectrum (upper panel) is shown at redshift  = 0.However, here it is computed with different cosmological parameters and cut-offs, as it is being compared to the BH19 -body simulation.Here, the standard perturbation theory one-loop prediction is actually performing better than the Hartree-Fock approximation, which is due to the comparably large infrared cut-off  min ≈ 2.1 ⋅ 10 −2 ℎ/Mpc that is resulting in smaller non-linear corrections to the power spectrum compared to the overestimation observed in Figure 7, where the infrared cut-off is much smaller.(lower panel).The numerical solution of the Hartree-Fock (HF) approximation (blue solid curves) is compared to the standard perturbation theory (SPT) one-loop prediction (black dashed curves) and data from the HR2 -body simulation (red solid curves).While the Hartree-Fock approximation shows better convergence properties at larger wave numbers, it overestimates the power spectrum at smaller wave numbers.

Velocity-divergence equal-time auto-spectrum.
In Figure 8, the velocity-divergence equal-time auto-spectrum (lower panel) is shown at redshift  = 0. Neither the Hartree-Fock nor the standard perturbation theory one-loop prediction can capture the sharp drop of the velocity-divergence power spectrum at small scales.Even worse, the Hartree-Fock approximation even fails to reproduce the power spectrum at relatively small wave numbers, where the standard perturbation theory one-loop prediction can capture at least the onset of the suppression.Indeed, the drop is associated with the energy transfer into vorticity modes [69] and is thus beyond the single-stream approximation.This shortcoming can be resolved by including vorticity and velocity dispersion degrees of freedom and is discussed in greater detail in Section 5.2.4.Single-stream approximation dimensionless density contrast equal-time auto-spectrum at redshift  = 0.The Hartree-Fock approximation captures the general shape of the non-linear power spectrum obtained from -body simulations much better than the standard perturbation theory oneloop prediction.
Unequal-time auto-spectra and equal-time cross-spectra.In Figure 9, the dimensionless density contrast unequal-time auto-spectrum (upper panel) at redshifts  = 0 and  ′ = 2.165, as well as the dimensionless velocity-divergence-density contrast equal-time crossspectrum (lower panel) at redshift  = 0, are shown.The oscillatory suppression of the unequal-time power spectrum is qualitatively captured by the Hartree-Fock approximation, although this is likely due to the partial resummation scheme that is also responsible for the unphysical oscillations in the propagator seen in Figure 5.
More interesting and important for the investigations in the next section that include velocity dispersion degrees of freedom is the cross-spectrum.The -body power spectrum stays close to the linear prediction up to  ≈ 1.0 ℎ/Mpc above which a sign change is observed.This is interpreted as a turnover from matter inflow to outflow and associated with shellcrossing [59,70].Naturally, this cannot be captured in the single-stream approximation and explains why neither the standard perturbation theory result nor the Hartree-Fock approximation are able to capture the -body power spectrum.
Summary.In summary, the insights gained within the single-stream approximation are: • The propagators are to very good approximation described by the large wave number limit (5.11), at least if the initial and final times are sufficiently far apart.
• The density contrast equal-time auto-spectrum converges well at small scales but fails to accurately capture the physics at mildly non-linear scales.This is most likely due to neglecting corrections to the bare vertex.
• The drop of the velocity-divergence equal-time auto-spectrum cannot be captured since it is associated to the energy transfer into vorticity modes.
• The oscillatory suppression of unequal-time auto-spectra is qualitatively captured due to the non-perturbative resummation of the propagators.Single-stream approximation dimensionless density contrast (upper panel) and velocitydivergence (lower panel) equal-time auto-spectrum at redshift  = 0.The numerical solution of the Hartree-Fock (HF) approximation (blue solid curves) is compared to the standard perturbation theory (SPT) one-loop prediction (black dashed curves) and data from the BH19 -body simulation (red solid curve).The Hartree-Fock approximation badly overestimates the -body velocity-divergence power spectrum, even in the small wave number regime.The drop of the velocity-divergence power spectrum is associated with the energy transfer into vorticity modes and is thus beyond single-stream approximation.
• The sign change in the velocity-divergence-density contrast equal-time cross-spectrum cannot be captured since it is associated with shell-crossing.

Including velocity dispersion degrees of freedom
Having seen the limitations of the single-stream approximation in the last section, the Hartree-Fock approximation is now studied for the field content (3.1), including the vorticity and velocity dispersion fields.In this case, the free-streaming of dark matter seen in the Hartree approximation as well as the non-linear coupling of modes encoded in the Fock self-energies contribute.To take both effects correctly into account, the two scales associated Single-stream approximation dimensionless density contrast unequal-time auto-spectrum (upper panel) at redshifts  = 0 and  ′ = 2.165 as well as velocity-divergence-density contrast equaltime cross-power spectrum (lower panel) at redshift  = 0.The numerical solution of the Hartree-Fock (HF) approximation (blue solid curves) is compared to the standard perturbation theory (SPT) oneloop prediction (black dashed curves) and data from the BH19 -body simulation (red solid curve).The oscillatory suppression of the unequal-time power spectrum is qualitatively captured by the Hartree-Fock approximation, while the sign change observed in the equal-time cross-spectrum cannot be captured since it is associated with shell-crossing.
with them, namely the free-streaming wave number  fs and the non-linear wave number  nl , need to be resolved.Numerically this is rather challenging for cold dark matter candidates because these scales are separated by several orders of magnitudes.Because the propagators are suppressed on small scales, the adaptive Runge-Kutta method takes smaller times steps as one solves at later times.The larger the ultraviolet cut-off wave number is, the more time steps are necessary to keep the errors under control.To obtain results in a sensible amount of time, but at the same time capture the effects of non-vanishing velocity dispersion, rather warm dark matter candidates are studied in the following.More specifically the initial velocity dispersion mean field is taken to be between σ in = 10 2 km 2 /s 2 and σ in = 10 5 km 2 /s 2 .Figure 10.Evolution of the velocity dispersion mean field in the Hartree-Fock approximation for different ultraviolet cut-offs that are multiples of the Nyquist wave number  Ny ≈ 11 ℎ/Mpc.Shown are colour-coded numerical solutions for differential initial conditions (violet, blue, green and yellow from coldest to warmest), the linearly decaying mean field (dotted curves) and three snapshots from the BH19 -body simulation (red diamonds).For warmer dark matter, where the free-streaming wave number is below the cut-off, the results are well converged, while for colder dark matter with a larger free-streaming wave number this is no longer the case.
The ultraviolet cut-off is varied between  max =  Ny and  max = 4 Ny to study the effect of including smaller scales.
Velocity dispersion mean field.In Figure 10, the velocity dispersion mean field for different ultraviolet cut-offs is shown.The different initial conditions are colour-coded in a scheme that is also used in all following power spectrum figures.The numerical solutions have been obtained for the ultraviolet cut-offs  max =  Ny and  max = 2 Ny as well as for  max = 4 Ny for the coldest candidate with an initial velocity dispersion mean field of σ in = 10 2 km 2 /s 2 .As long as the free-streaming wave number is below or of the order of the ultraviolet cut-off, the numerical solutions are well converged and insensitive to a larger ultraviolet cutoff.This is best seen for the warmest candidate with an initial velocity dispersion mean field of σ in = 10 5 km 2 /s 2 .The corresponding free-streaming wave number is  fs ≈ 2.1 ℎ/Mpc at  = 99 and grows to  fs ≈ 6.2 ℎ/Mpc at  = 0, well below the ultraviolet cut-off.For colder initial conditions, one naturally needs a larger ultraviolet cut-off to capture the smallscale power that sources the growth of the mean field.This can be seen best for the coldest candidate, where the free-streaming wave number at the point of turnover from decay to growth is roughly of the order  fs ≈ 210 ℎ/Mpc, which is way above the largest ultraviolet cut-off  max = 4 Ny ≈ 43 ℎ/Mpc used here.Naturally, this extends even more drastically to colder candidates, which increases the need to include ever larger wave numbers.
Comparing to the BH19 -body simulation, it is evident that the velocity dispersion mean field is even larger than the warmest dark matter candidates studied here.The similar velocity dispersion mean field values at late times for dark matter candidates with very different initial conditions suggest some universal fixed point at late times.The inability to quantitatively predict the simulation data is likely related to the Hartree-Fock approximation .Density contrast and isotropic velocity dispersion reduced propagators for a dark matter model with σ in = 10 5 km 2 /s 2 at redshift  = 0. Also shown is the ansatz (5.21), superposing the effect of free-streaming found in the Hartree approximation (4.8) and the ultraviolet suppression due to the sweeping effect in the large wave number limit (5.11).Qualitatively, the ansatz captures the numerical solution of the Hartree-Fock approximation which naturally captures both effects.
or the restricted set of degrees of freedom.On the other hand, testing this with very cold dark matter candidates, such as a weakly interacting particle with σ in ∼ 10 −10 km 2 /s 2 , would be necessary in order to rule out the indicated very strong growth behaviour that is observed in the Hartree approximation in Figure 2, which suggests that candidates that start out colder might become warmer at later times than candidates that initially start out warmer.Due to the above mentioned computational limitations it was not possible to verify this here.
Propagator.In Figure 11, the density contrast and isotropic velocity dispersion reduced propagators for the dark matter model with σ in = 10 5 km 2 /s 2 are shown at redshift  = 0. To qualitatively understand the wave number dependence of the propagator, the ansatz is also shown.It is based on the large external wave number limit (5.11),where G R  ( −  ′ ,  2 σ ) is the retarded propagator that is obtained in the Hartree approximation for a linearly decaying mean field and can be reconstructed from the three independent solutions (4.8).This naïve ansatz qualitatively captures the Hartree-Fock approximation and reproduces the superposition of oscillations from free-streaming encoded in the Hartree self-energy (3.24) and the decaying oscillations due to the coupling of modes encoded in the retarded Fock selfenergy (5.3). 18Although not shown here, the performance of the ansatz (5.21) worsens for colder dark matter models.Since the propagator G R  ( −  ′ ,  2 σ ) is strictly speaking only valid for a linearly decaying mean field, it fails to capture the effect of a strongly growing mean field, as is the case for colder dark matter models.10.They are compared to the linear power spectrum (black dotted curves), the standard perturbation theory (SPT) one-loop prediction (black dashed curves) and data from the BH19 -body simulation (red solid curve).Both spectra are suppressed for warmer dark matter models, although the velocity-divergence power spectrum shows a much stronger suppression.
Density contrast and velocity-divergence equal-time auto-spectra.Although only initially rather warm dark matter candidates are studied, the question arises whether the approximation including velocity dispersion degrees of freedom can overcome the shortcomings of the single-stream approximation.To see the impact on non-vanishing velocity dispersion, the dimensionless density contrast and velocity-divergence equal-time auto-spectra are shown in Figure 12.The solutions correspond to different velocity dispersion mean field initial conditions and are colour-coded in the same way as in Figure 10.Additionally, the free-streaming scales related to the corresponding velocity dispersion mean field are displayed.Since the numerical solutions all have a mean field of the same order at redshift  = 0, and therefore a similar free-streaming wave number, only a single free-streaming wave number is indicated here.It is clearly visible that the dark matter models starting out with a warmer initial condition show a large suppression in the density contrast as well as in the velocity-divergence power spectrum.Interestingly, one can observe that although the velocity dispersion is of similar order at redshift  = 0, the dark matter models starting out comparably cold show near to no suppression in the density contrast power spectrum and match the -body data rather well.Comparing with the evolution of the mean field shown in Figure 10, it seems that if the mean field obtains a large value only at late times, one does not see a suppression in the density contrast power spectrum.In contrast, consider the warmest dark matter model shown (orange lines) that is warm for a large part of its time evolution.Although the final mean field is even below the final mean field of the initially colder model, a significant suppression is observed in the density contrast power spectrum.This suggest that indeed the dark matter model should start out with a small velocity dispersion mean field in order to match the observed power spectrum.On the other hand, this also implies that a latetime large velocity dispersion mean field does not necessarily imply that the density contrast power spectrum needs to be suppressed.This is rather different for the velocity-divergence power spectrum, which shows a dip in the equal-time power spectrum for all dark matter models.Although the power spectrum does not match the -body simulations quantitatively, the qualitative behaviour of the dip can be understood from the numerical solutions obtained here.As the suppression of the velocity-divergence power spectrum in the mildly non-linear regime is related to the energy transfer into vorticity modes [69], the position of the suppression dip is sensitive to the vorticity power spectrum.As is shown in Figure 15 and argued below, the amplitude of vorticity spectrum is heavily dependent on the magnitude of the velocity dispersion mean field, which is underestimated in the numerical solutions as compared to the BH19 -body simulation.As a consequence, the vorticity power spectrum is of similar magnitude as the velocity-divergence power spectrum at larger wave numbers and thus the dip indicating the scale of energy transfer is shifted towards higher wave numbers.An accurate reconstruction of the velocity-divergence power spectrum and the suppression in the mildly non-linear regime therefore necessitates a correct description of the vorticity power spectrum, which in the current framework requires a realistic velocity dispersion mean field at late times.Isotropic velocity dispersion equal-time auto-spectra.In Figure 13, the time (upper panel) and wave number (lower panel) dependencies of the dimensionless isotropic velocity dispersion equal-time auto-spectrum are shown at constant wave number  ≈ 0.5 ℎ/Mpc and at fixed redshift  = 0, respectively.
At constant wave number (upper panel) the power spectrum grows in time and does so more strongly for initially colder dark matter models.Also shown are three snapshots of the BH19 -body simulation which exhibit a stronger growth at later times, roughly scaling ∝  8  + .The numerical solutions cannot reproduce the time dependence observed in the BH19 -body simulation since the time dependence is directly linked to the velocity dispersion mean field shown in Figure 10.The mean field fails to correctly reproduce the observed time dependence as initially rather warm dark matter models are studied here as discussed above.
At fixed redshift (lower panel) the power spectrum shows a power-law scaling in the infrared up to the free-streaming wave number.-body simulations suggest that this scale is in fact related to shell-crossing and the largest collapsed structures [59,70].As indicated, the free-streaming wave number is related to the velocity dispersion mean field in the current framework.It is clearly visible that for wave numbers larger than the free-streaming scale the  10.They are compared to the linear power spectrum (black dotted curves), the standard perturbation theory (SPT) one-loop prediction (black dashed curves) and data from the BH19 -body simulation (red solid curve).The sign change of the cross-spectra is associated with shell-crossing and captured by the Hartree-Fock approximation in both power spectra.contrast (upper panel) and velocity dispersion-density contrast (lower panel) equal-time crossspectra are shown at redshift  = 0.
As seen in Figure 9, the single-stream approximation badly fails to describe the velocitydivergence-density contrast equal-time cross-spectrum.This is related to the fact that the sign change observed in -body data is due to shell-crossing which cannot be captured by the single-stream approximation.In contrast, when including velocity dispersion degrees of freedom, the numerical solutions shown in Figure 14 exhibit the expected sign change.It is evident that the scale of this sign change is related to the free-streaming wave number (the sign change occurs roughly at 2 fs ) and thus to the velocity dispersion mean field.The numerical solutions shown here show a sign change at a higher wave number compared to pling of fluctuations and is therefore naturally generated in the Hartree-Fock approximation through the coupling of modes in the statistical self-energy.
One can observe the late-time power-law-like scaling in time ∝  7 + [14,70] in the BH19 -body data and the numerical solutions obtained here grow similarly, although depending on the specific dark matter model.This is related to the fact that vorticity is initially sourced by the coupling of density contrast and velocity dispersion fluctuations, the latter of which are sourced by the velocity dispersion mean field.Therefore, similar to the velocity dispersion power spectra shown in Figure 13, the time dependence of the vorticity power spectrum is sensitive to the evolution of the velocity dispersion mean field.
For the wave number dependence one finds a spectral index   ≈ 2.5 in the infrared from various -body simulations [14,69,70], corresponding to the scaling ∆ 55 ∝  5.5 of the dimensionless vorticity power spectrum, before dropping on scales smaller than the freestreaming scale.The numerical solutions obtained here show a spectral index of almost exactly   ≈ 2. 20 Indeed, there are arguments that if vorticity is absent initially, one expects the power spectrum to actually scale with a spectral index   = 2. Since the position space covariance function only has finite support due to causality, the Fourier transform needs to be analytic [71].Since the vorticity mode is transverse, the power spectrum carries an additional transverse projector P  () such that the simplest scaling required for analyticity is given by   = 2. Precisely such a scaling is observed in the infrared of the Hartree-Fock approximation.It has been argued [69] that a deviation from this scaling observed in body simulations could be due to missing contributions of large scale modes due to the finite simulation box size, rendering the -body power spectrum too steep.Along similar lines one could argue that the deviation of the isotropic velocity dispersion equal-time auto-spectrum from -body data shown in Figure 13 is due to this finite size effect.
It is remarked, that although the infrared is perturbative in the sense that fluctuations are small and the density contrast power spectrum is well described by perturbation theory, the generation of vorticity and velocity dispersion is highly non-perturbative.This is also the reason why perturbative methods do not correctly capture the scaling of the vorticity power spectrum [6,72].More recently, there have been investigations that reproduce a spectral index   = 2 at two-loop order in an approximation that involves also higher-order cumulants [9].It is interesting that the vorticity power spectrum shown in the lower panel of Figure 15 is enhanced in the ultraviolet, while the velocity-divergence power spectrum shown in the lower panel of Figure 12 significantly drops below the linear prediction with a pronounced dip near the free-stream wave number.Indeed, when looking at the dimensionful power spectra   and   , one actually finds that the peak of the vorticity power spectrum aligns with the position of the velocity-divergence power spectrum dip [70].This is interpreted as a consequence of angular momentum conservation converting power from velocity-divergence to vorticity modes, preventing the further infall of matter [69].

Conclusions
In summary, the Dyson-Schwinger equations have been investigated in the Hartree and Hartree-Fock approximations for dynamical cosmological structure formation caused by dark matter.To this end, an extension of the perfect pressureless fluid dark matter model was studied by including velocity dispersion degrees of freedom and in particular a non-vanishing velocity dispersion mean field.
The inclusion of velocity dispersion modifies the equation of motion of the velocity field.In particular, a time-local 'Hartree-type' contribution to the self-energy causes the propagator to depend on the solution of the velocity dispersion mean field and introduces a wave number dependence associated to the free-streaming of dark matter particles.Further, a non-local 'Fock-type' contribution to the self-energy captures the non-linear coupling of modes through loop effects.The resulting coupled equations for the mean field, propagators and power spectra were solved numerically in the Hartree-Fock approximation.While this approximation has a self-consistent one-loop structure, loop corrections to the three-and higher-point 1PI vertices are neglected.
A crucial ingredient is the inclusion of a time-dependent velocity dispersion mean field.At early times the mean field is decaying due to the Hubble expansion, however, it typically reaches a turnover point at later times, after which it starts to grow as it becomes non-linearly sourced by fluctuations.The magnitude of this growth depends on the initial value, with smaller initial values resulting in larger growth, such that the final order of magnitude of the velocity dispersion mean field seems to be largely independent of the initial value.It is worth noting that the numerical solutions obtained here lead to values that are still smaller in magnitude than those observed in -body simulations.
Overall, the numerical results agree well with expectations.More specifically, it was demonstrated that the equal-time density contrast power spectrum exhibits better convergence to -body simulations compared to standard perturbation theory.Further, as long as dark matter is initially sufficiently cold, the strong growth of velocity dispersion at late times does not seem to significantly affect the equal-time density contrast power spectrum.This is drastically different for the velocity-divergence power spectrum.Although the equal-time velocity-divergence power spectrum could not be reconstructed accurately, the qualitative behaviour of the dip in the non-linear regime could be captured by including velocity dispersion degrees of freedom and is inherently tied to the failure of the single-stream approximation.At shell-crossing, vectorial velocity modes are generated and energy is transferred into these modes from scalar velocity modes.Therefore, the scale of the dip in the velocity-divergence power spectrum is related to the free-streaming wave number and in turn related to the shape of the vorticity power spectrum.The amplitude of the vorticity power spectrum is directly linked to the velocity dispersion mean field, making an accurate description of the mean field crucial for correctly reproducing the dip.However, the numerical solutions for the mean field obtained here are still smaller than those observed in -body simulations.As a result, the vorticity spectrum has a smaller amplitude and consequently the dip in the velocity-divergence power spectrum appears at larger wave numbers compared to -body simulations.
By including a non-vanishing velocity dispersion mean field, vorticity vector modes are naturally generated through the non-linear coupling of scalar modes, even if initially absent.The numerical solutions grow in time and exhibit a power-law wave number dependence of the equal-time power spectrum in the infrared with a spectral index   ≈ 2, as expected from analytical arguments.This result is close to the spectra index   ≈ 2.5 observed in -Body simulations and agrees better than other (semi-)analytical methods with -body data, although it should be stressed that more recently a spectral index   = 2 was reproduced in a perturbative higher-order cumulant approach [9].
Finally, the equal-time power spectrum for the trace of the velocity dispersion tensor was computed, starting from vanishing initial conditions.The velocity dispersion fluctuations are naturally sourced through the velocity dispersion mean field and the overall amplitude is therefore closely tied to the evolution of the mean field, very similar to the vorticity power spectrum.A power-law wave number dependence was found in the infrared with a spectral index of   ≈ −1.8, which is slightly below the spectral index   ≈ −1 observed in -body simulations.
Due to computational limitation, only the scalar velocity dispersion degrees of freedom have been included in this work.Since the velocity dispersion vector and tensor modes can only be sourced non-linearly (if initially absent) through cross-correlations of velocity and velocity dispersion scalar modes, their impact is expected to be subdominant compared to the investigated scalar degrees of freedom.Further, the degree to which cold dark matter could be studied is also heavily limited by computational resources in the presented work.The dark matter candidates that were investigated are comparably warm as the colder candidates would necessitate an accurate resolution also far into the ultraviolet of the wave number regime.Here, the equations of motion are solved with a finite element method in wave numbers, which necessarily involves a rather large wave number grid for accurate results.For a precise treatment extending further into the ultraviolet, or even spanning the hole wave number domain, more advanced spectral methods might be more sensible to accurately capture the ultraviolet correlations sourcing the velocity dispersion mean field.
The findings of this work indicate that the velocity dispersion mean field is largely independent of the initial value, as it grows more strongly the colder the initial condition.The velocity and velocity dispersion power spectra are sensitive to the mean field, as it determines the free-streaming scale.Consequently, the velocity dispersion mean field is closely tied to the ultraviolet suppression of the power spectra, which in turn are responsible for the growth of the mean field.In other words, as the mean field becomes larger, the power spectra become more suppressed in the ultraviolet, resulting in a lesser contribution from fluctuations to the mean field.This chain of reasoning leads to the conjecture of a universal fixed point at late times for the velocity dispersion mean field.
It was shown, that for a viable description of dark matter, it is necessary to go beyond the single-stream approximation and perturbative methods.While the inclusion of velocity dispersion degrees of freedom allows for an effective multi-stream description that qualitatively captures effects related to shell-crossing, it is unclear whether the inclusion of higher-order velocity cumulants is needed for an accurate description of structure formation after shell-crossing.It is evident that non-perturbative methods are necessary to describe the non-linear regime, which can be understood from the perspective of an infinite but partial resummation of loop diagrams.Results of the self-consistent one-loop approximation suggest that corrections to higher-order 1PI vertices are needed to accurately describe even the mildly non-linear regime.To this end, the functional renormalisation group provides a method to systematically explore approximations and it is planned to investigate this direction further in the future.

A Functional approach to cosmological field theory
Here, the Martin-Siggia-Rose/Janssen-de Dominicis response field formalism for cosmology is reviewed.For a detailed review of the functional formalism for the dynamics of stochastic systems, the reader is referred to reference [73].

2 Figure 1 .
Figure1.Hypergeometric functions appear in the solutions (4.8) of the reduced density propagator as a function of  2 σ .The growing mode (solid curve) is enhanced at non-vanishing  2 σ , while the decaying modes (dashed and dotted curves) are suppressed and oscillate with increasing frequency and decaying amplitude.

5 Figure 2 .
Figure 2. Time evolution of the velocity dispersion mean field as a function of the scale factor for different choices of initial values.Solutions which are sourced by correlations (solid curves) grow in time compared to those decaying linearly with the Hubble expansion (dotted curves).Smaller initial values of the velocity dispersion mean field correspond to larger free-streaming wave numbers and lead to a stronger late-time growth by non-linear effects.

2 Figure 3 .
Figure 3.Reduced propagators for the density contrast (blue curves), velocity-divergence (red curves), isotropic velocity dispersion (yellow curves) and anisotropic velocity dispersion (violet curves) at redshift  = 0, normalised to the standard growing mode as a function of  2 σ in .Results are shown for two choices of initial velocity dispersion mean fields, 1 km 2 /s 2 (solid curves) and 10 −5 km 2 /s 2 (dashed curves).The damping scale responsible for the oscillations is set by the time-dependent mean field, shown in Figure2, and is at much smaller  2 σ in for the colder dark matter model due to the strong late-time growth of the velocity dispersion mean field.

Figure 4 .
Figure 4.The propagators (5.12), (5.10) and (5.11) (componentwise) normalised to the linear propagator as a function of  v [e − in − e  ′ − in ].While the functional renormalisation group (FRG) propagator (solid curve) shows a Gaussian suppression, the Hartree-Fock (HF) approximation (dashed curve) and the 1PI one-loop approximation (dotted curve) feature oscillations related to the different partial resummations of the perturbative series.

Figure 5 .
Figure5.Single-stream approximation reduced density contrast (blue curves) and velocitydivergence (red curves) propagator time dependence at  ≈ 2 ℎ/Mpc (upper panel) and wave number dependence at  = 0 (lower panel).The numerical solution of the Hartree-Fock (HF) approximation (solid curves) is compared to the TH (dashed curves) and CS (dotted curves) approximation.The Hartree-Fock and TH approximations are near to identical since the latter is based on the large wave number limit of the former.

Figure 6 .
Figure 6.Single-stream approximation dimensionless density contrast equal-time auto-spectrum, normalised to the linear spectrum (black dotted curves), at redshifts  = 2 (upper panel) and  = 0 (lower panel).The numerical solution of the Hartree-Fock (HF) approximation (blue solid curves) is compared to the standard perturbation theory (SPT) one-loop prediction (black dashed curves) and data from the HR2 -body simulation (red solid curves).While the Hartree-Fock approximation shows better convergence properties at larger wave numbers, it overestimates the power spectrum at smaller wave numbers.

2 Figure 7 .
Figure 7.Single-stream approximation dimensionless density contrast equal-time auto-spectrum at redshift  = 0.The Hartree-Fock approximation captures the general shape of the non-linear power spectrum obtained from -body simulations much better than the standard perturbation theory oneloop prediction.

3 Figure 8 .
Figure 8. Single-stream approximation dimensionless density contrast (upper panel) and velocitydivergence (lower panel) equal-time auto-spectrum at redshift  = 0.The numerical solution of the Hartree-Fock (HF) approximation (blue solid curves) is compared to the standard perturbation theory (SPT) one-loop prediction (black dashed curves) and data from the BH19 -body simulation (red solid curve).The Hartree-Fock approximation badly overestimates the -body velocity-divergence power spectrum, even in the small wave number regime.The drop of the velocity-divergence power spectrum is associated with the energy transfer into vorticity modes and is thus beyond single-stream approximation.

3 Figure 9 .
Figure 9.Single-stream approximation dimensionless density contrast unequal-time auto-spectrum (upper panel) at redshifts  = 0 and  ′ = 2.165 as well as velocity-divergence-density contrast equaltime cross-power spectrum (lower panel) at redshift  = 0.The numerical solution of the Hartree-Fock (HF) approximation (blue solid curves) is compared to the standard perturbation theory (SPT) oneloop prediction (black dashed curves) and data from the BH19 -body simulation (red solid curve).The oscillatory suppression of the unequal-time power spectrum is qualitatively captured by the Hartree-Fock approximation, while the sign change observed in the equal-time cross-spectrum cannot be captured since it is associated with shell-crossing.

Figure 11
Figure11.Density contrast and isotropic velocity dispersion reduced propagators for a dark matter model with σ in = 10 5 km 2 /s 2 at redshift  = 0. Also shown is the ansatz (5.21), superposing the effect of free-streaming found in the Hartree approximation (4.8) and the ultraviolet suppression due to the sweeping effect in the large wave number limit(5.11).Qualitatively, the ansatz captures the numerical solution of the Hartree-Fock approximation which naturally captures both effects.

3 Figure 12 .
Figure12.Dimensionless density contrast (upper panel) and velocity-divergence (lower panel) equaltime auto-spectrum at redshift  = 0.The numerical solutions of the Hartree-Fock (HF) approximation are colour coded according to the respective mean fields shown in Figure10.They are compared to the linear power spectrum (black dotted curves), the standard perturbation theory (SPT) one-loop prediction (black dashed curves) and data from the BH19 -body simulation (red solid curve).Both spectra are suppressed for warmer dark matter models, although the velocity-divergence power spectrum shows a much stronger suppression.

2 Figure 14 .
Figure14.Dimensionless velocity-divergence-(upper panel) and velocity dispersion-density contrast (lower panel) equal-time cross-power spectrum  = 0.The numerical solutions of the Hartree-Fock (HF) approximation are colour coded according to the respective mean fields shown in Figure10.They are compared to the linear power spectrum (black dotted curves), the standard perturbation theory (SPT) one-loop prediction (black dashed curves) and data from the BH19 -body simulation (red solid curve).The sign change of the cross-spectra is associated with shell-crossing and captured by the Hartree-Fock approximation in both power spectra.

Table 1 .
Number of diagrams included at loop order ℓ in the different resummation schemes.While the 1PI one-loop approximation captures only a single diagram at each loop order, the Hartree-Fock approximation improves on this, but captures only a subset those diagrams resummed in the functional renormalisation group (FRG).

Table 2 .
Details of the -body simulations used as a comparison to the numerical solutions of the Hartree-Fock approximation.