Evolution of dark matter velocity dispersion

TL;DR

This paper extends cosmological perturbation theory by including dark matter velocity dispersion, revealing late-time growth effects and implications for different dark matter candidates, which impacts the understanding of structure formation beyond shell-crossing.

Contribution

It introduces a tensor decomposition of dark matter velocity dispersion into scalar, vector, and tensor fields, and analyzes its late-time growth and effects on structure formation.

Findings

Velocity dispersion grows strongly at late times due to back-reaction.

Warm dark matter shows weaker late-time growth compared to cold dark matter.

Non-zero velocity dispersion breaks single-stream approximation, enabling beyond shell-crossing analysis.

Abstract

Cosmological perturbation theory for the late Universe dominated by dark matter is extended beyond the perfect fluid approximation by taking the dark matter velocity dispersion tensor as an additional field into account. A proper tensor decomposition of the latter leads to two additional scalar fields, as well as a vector and a tensor field. Most importantly, the trace of the velocity dispersion tensor can have a spatially homogeneous and isotropic expectation value. While it decays at early times, we show that a back-reaction effect quadratic in perturbations makes it grow strongly at late times. We compare sterile neutrinos as a candidate for comparatively warm dark matter to weakly interacting massive particles as a rather cold dark matter candidate and show that the late time growth of velocity dispersion is stronger for the latter. Another feature of a non-vanishing velocity dispersion expectation value is that it destroys the apparent self-consistency of the single-stream approximation and allows thereby to treat times and scales beyond shell-crossing.

Tables

Table 1. Table 1 : Cosmological parameters found by the Planck Collaboration [ 106 ] where h ℎ h is the reduced Hubble constant, Ω c , 0 subscript Ω 𝑐 0 \Omega_{c,0} and Ω b , 0 subscript Ω 𝑏 0 \Omega_{b,0} are the current cold dark matter and baryonic matter density parameters respectively, z re subscript 𝑧 re z_{\text{re}} is the redshift at reionisation, A s subscript 𝐴 𝑠 A_{s} is the primordial scalar amplitude at the pivot scale k ∗ = 0.05 h / k_{*}=0.05\,h/ Mpc and n s subscript 𝑛 𝑠 n_{s} is the scalar spectral index.

$h$	${\Omega }_{c,0}{h}^2$	${\Omega }_{b,0}{h}^2$	$z_{\text{re}}$	$A_s$	$n_s$
$0.6774$	$0.1188$	$0.0223$	$8.8$	$2.142\cdot {10}^{-9}$	$0.9667$

Equations

$$ds^{2}=a^{2}(\tau)\,\Big{[}-\big{(}1+2\,\phi(\tau,\mathbf{x})\big{)}\,d\tau^{2}+\big{(}1-2\,\phi(\tau,\mathbf{x})\big{)}\,d\mathbf{x}^{2}\Big{]}\;,$$

$$\frac{\dot{\mathcal{H}}(\tau )}{{\mathcal{H}}^2(\tau )}=1-\frac{3}{2}{\Omega }_m(\tau ),1={\Omega }_m(\tau )+{\Omega }_{\Lambda }(\tau ),$$

$${\Delta }_{\mathbf{x}}\varphi =\frac{3}{2}{\mathcal{H}}^2{\Omega }_m\delta .$$

$$\frac{d\mathbf{x}}{d\tau }=\frac{\mathbf{p}}{a(\tau )m},\frac{d\mathbf{p}}{d\tau }=-a(\tau )m{\nabla }_{\mathbf{x}}\varphi (\tau ,\mathbf{x}),$$

$${\partial }_{\tau }f+\frac{\mathbf{p}}{am}\cdot {\nabla }_{\mathbf{x}}f-am{\nabla }_{\mathbf{x}}\varphi \cdot {\nabla }_{\mathbf{p}}f=0.$$

$$\rho (\tau ,\mathbf{x})=\frac{m}{a^3(\tau )}{\int }_{{\mathbb{R}}^3}{d}^{3}pf(\tau ,\mathbf{x},\mathbf{p}).$$

$$\delta (\tau ,\mathbf{x})\equiv \frac{\rho (\tau ,\mathbf{x})-\bar{\rho }(\tau )}{\bar{\rho }(\tau )}.$$

$$\rho (\tau ,\mathbf{x})u_i(\tau ,\mathbf{x})=\frac{m}{a^3(\tau )}{\int }_{{\mathbb{R}}^3}{d}^{3}p\frac{p_i}{a(\tau )m}f(\tau ,\mathbf{x},\mathbf{p}),$$

$$\rho (\tau ,\mathbf{x})u_i(\tau ,\mathbf{x})u_j(\tau ,\mathbf{x})+T_{ij}(\tau ,\mathbf{x})=\frac{m}{a^3(\tau )}{\int }_{{\mathbb{R}}^3}{d}^{3}p\frac{p_i{p}_j}{a^2(\tau )m^2}f(\tau ,\mathbf{x},\mathbf{p}),$$

$$\tilde{f}(\tau ,\mathbf{x},\mathbf{p})\equiv \frac{1}{\bar{\rho }(\tau )}\frac{m}{a^3(\tau )}f(\tau ,\mathbf{x},\mathbf{p}),$$

$$M(\tau ,\mathbf{x};\mathbf{l})\equiv {\int }_{{\mathbb{R}}^3}{d}^{3}p\exp \left\{\frac{\mathbf{l}\cdot \mathbf{p}}{a(\tau )m}\right\}\tilde{f}(\tau ,\mathbf{x},\mathbf{p}),$$

$$m_{i_{1}...i_{n}}^{(n)}(\tau,\mathbf{x})=\frac{\partial^{n}M(\tau,\mathbf{x};\mathbf{l})}{\partial l_{i_{1}}...\partial l_{i_{n}}}\bigg{|}_{\mathbf{l}=\mathbf{0}}\;,$$

$$c^{(0)}(\tau,\mathbf{x})=\ln\big{(}1+\delta(\tau,\mathbf{x})\big{)}\;,\qquad c_{i}^{(1)}(\tau,\mathbf{x})=u_{i}(\tau,\mathbf{x})\;,\qquad c_{ij}^{(2)}(\tau,\mathbf{x})=\sigma_{ij}(\tau,\mathbf{x})\;.$$

$${\partial }_{\tau }C+\mathcal{H}\mathbf{l}\cdot {\nabla }_{\mathbf{l}}C+{\nabla }_{\mathbf{x}}\cdot {\nabla }_{\mathbf{l}}C+{\nabla }_{\mathbf{x}}C\cdot {\nabla }_{\mathbf{l}}C+\mathbf{l}\cdot {\nabla }_{\mathbf{x}}\varphi =0.$$

$$C(\tau,\mathbf{x};\mathbf{l})=c^{(0)}(\tau,\mathbf{x})+l_{i}\,c_{i}^{(1)}(\tau,\mathbf{x})\;,\qquad M(\tau,\mathbf{x};\mathbf{l})=\big{[}1+\delta(\tau,\mathbf{x})\big{]}\,e^{\mathbf{l}\cdot\mathbf{u}(\tau,\mathbf{x})}\;,$$

$$\tilde{f}(\tau,\mathbf{x},\mathbf{p})=\big{[}1+\delta(\tau,\mathbf{x})\big{]}\,\delta_{\text{D}}^{(3)}\big{(}\mathbf{p}-a(\tau)\,m\,\mathbf{u}(\tau,\mathbf{x})\big{)}\;,$$

$$\tilde{f}(\tau,\mathbf{x},\mathbf{p})=\sum_{s}\big{[}1+\delta_{s}(\tau,\mathbf{x})\big{]}\,\delta_{\text{D}}^{(3)}\big{(}\mathbf{p}-a(\tau)\,m\,\mathbf{u}_{s}(\tau,\mathbf{x})\big{)}\;,$$

$$C(\tau ,\mathbf{x};\mathbf{l})=c^{(0)}(\tau ,\mathbf{x})+l_i{c}_i^{(1)}(\tau ,\mathbf{x})+\frac{1}{2}{l}_i{l}_j{c}_{ij}^{(2)}(\tau ,\mathbf{x}),$$

$$M(\tau,\mathbf{x};\mathbf{l})=\big{[}1+\delta(\tau,\mathbf{x})\big{]}\,e^{l_{i}\,u_{i}(\tau,\mathbf{x})+\frac{1}{2}\,l_{i}\,l_{j}\,\sigma_{ij}(\tau,\mathbf{x})}\;.$$

$$\tilde{f}(\tau,\mathbf{x},\mathbf{p})=\frac{1+\delta}{\sqrt{(2\pi)^{3}\,\det(a^{2}\,m^{2}\,\boldsymbol{\sigma})}}\,\exp\left\{-\frac{1}{2}\,\big{(}p_{i}-a\,m\,u_{i}\big{)}\,\frac{\sigma_{ij}^{-1}}{a^{2}\,m^{2}}\,\big{(}p_{j}-a\,m\,u_{j}\big{)}\right\}\;,$$

$$\tilde{f}(\tau ,\mathbf{x},\mathbf{p})=\frac{1+\delta }{a^3{m}^3}{\left(\frac{m}{2\pi k_{\text{B}}T}\right)}^{\frac{3}{2}}\exp \left\{-\frac{m}{2k_{\text{B}}T}\frac{{\mathbf{p}}^2}{a^2{m}^2}\right\}.$$

$$\dot{\delta }+u_i{\delta }_{,i}+(1+\delta )u_{i,i}=0,$$

$${\dot{u}}_i+\mathcal{H}{u}_i+u_j{u}_{i,j}+{\sigma }_{ij,j}+{\sigma }_{ij}\ln (1+\delta {)}_{,j}+{\varphi }_{,i}=0,$$

$${\dot{\sigma }}_{ij}+2\mathcal{H}{\sigma }_{ij}+u_k{\sigma }_{ij,k}+{\sigma }_{jk}{u}_{i,k}+{\sigma }_{ik}{u}_{j,k}=0.$$

$$⟨{\sigma }_{ij}(\tau ,\mathbf{x})⟩={\delta }_{ij}\bar{\sigma }(\tau ).$$

$${\varsigma }_{ij}(\tau ,\mathbf{x})\equiv {\sigma }_{ij}(\tau ,\mathbf{x})-{\delta }_{ij}\bar{\sigma }(\tau ).$$

$$\dot{\bar{\sigma }}(\tau )+2\mathcal{H}(\tau )\bar{\sigma }(\tau )+\frac{1}{3}⟨u_i(\tau ,\mathbf{x}){\varsigma }_{jj,i}(\tau ,\mathbf{x})⟩+\frac{2}{3}⟨{\varsigma }_{ij}(\tau ,\mathbf{x})u_{i,j}(\tau ,\mathbf{x})⟩=0.$$

$$⟨T_{ij}(\tau ,\mathbf{x})⟩={\delta }_{ij}\bar{P}(\tau ).$$

$$\bar{P}(\tau )=\bar{\rho }(\tau )\bar{\sigma }(\tau )+\frac{\bar{\rho }(\tau )}{3}⟨\delta (\tau ,\mathbf{x}){\varsigma }_{ii}(\tau ,\mathbf{x})⟩.$$

$${\dot{u}}_i+\mathcal{H}{u}_i+u_j{u}_{i,j}+\bar{\sigma }{\varsigma }_{ij,j}+\bar{\sigma }\ln (1+\delta {)}_{,i}+\bar{\sigma }{\varsigma }_{ij}\ln (1+\delta {)}_{,j}+{\varphi }_{,i}=0,$$

Full text

1 Introduction

The matter distribution throughout the Universe shows structures on a wide range of scales. Surveys indicate that matter is distributed homogeneous and isotropic on the larges scales, commonly referred to as the cosmological principle. But observations also show a rich variety of dense cosmological objects, from stars to galaxies, galaxy groups to galaxy clusters, forming filaments separated by immense voids making up the so-called cosmic web. It is assumed that these structures are the outcome of gravitational instabilities of near to scale-invariant and Gaussian initial density fluctuations in the matter distribution in an early epoch where the Universe is dominated by non-relativistic matter, dubbed cold dark matter. The collapse of matter is a competition between the expansion of the Universe and gravity, and current observations favour a scenario in which small scale structures collapse and virialise earlier than larger scales structures, also known as the bottom-up scenario [1, 2].

From a theoretical point of view, structure formation can be described by following the evolution of an ensemble of self-gravitating dark matter particles. Commonly this is done by considering Newtonian gravity because the scales under consideration are much smaller than the Hubble horizon, where general relativistic corrections need to be taken into account. One often uses a description in terms of classical kinetic theory for the phase-space distribution of dark matter particles. This uses the collisionless limit of the Boltzmann equation with gravitational fields, generically known as Vlasov equation. The further description is done in either of two frames. Eulerian coordinates are fixed to the expanding Friedmann-Lemaître-Robertson-Walker background cosmology, while Lagrangian coordinates follow the trajectories of individual dark matter particles or the so-called displacement field [3]. Working in the former, one can consider moments or cumulants of the dark matter velocity distribution to obtain an infinite tower of coupled equations. Often they are truncated by setting the velocity dispersion tensor and higher cumulants to zero. The resulting equations describe dark matter as a perfect pressureless fluid. This is the so-called single-stream approximation since it assigns a unique velocity value to each point in space. Although the evolution equations are non-linear, these non-linearities are assumed to be small at early times and large scales and can therefore be treated perturbatively. For a review of standard cosmological perturbation theory we refer to reference [4].

Linear perturbation theory has had a tremendous success in describing the evolution of cosmic fields at early times and large scales. However, at later times and smaller scales, gravitational collapse becomes non-linear. A common way to deal with non-linearities in structure formation are numerical $N$-body simulations [5, 6]. While these allow to study structure formation in detail over a large amount of scales, they are limited by computational power. For this reason and to gain further insight into the mechanisms at play, different theoretical schemes have been put forward. Higher orders in perturbation theory have been calculated [7, 8, 9, 10, 11, 12, 13, 14, 15], various resummation schemes have been devised [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43], effective theories for a description of a larger range of scales have been developed [44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73] and the renormalisation group for a (statistical) field theoretic description has been studied [21, 74, 75].

However, to the extend that these methods are based on the single-stream approximation, they are ultimately limited to early times or very large scales. The approximation breaks down when multiple streams of matter coexist at the same region in space, a phenomenon known as shell-crossing. By definition, the single-stream approximation can not account for multiple velocities at the same point in space. In Lagrangian coordinates, the gravitational force becomes intrinsically non-local after stream crossing and thus the density field can no longer be approximated as an expansion of the determinant of the local deformation tensor.

There have been a number of works addressing this fundamental flaw in the description of dark matter. In an approach to estimate the effects of shell-crossing in Eulerian perturbation theory, reference [76] computed the back-reaction of vorticity and velocity dispersion generation on the density and velocity-divergence fields using numerical simulations. In references [77] and [78], velocity dispersion has been accounted for explicitly in the derivation from the Vlasov equation as source terms which are provided by numerical $N$-body simulations. Similarly, reference [49] formulated an effective field theory by extending Lagrangian perturbation theory with sources. These encode small scale effects such as shell-crossing and the long wavelength expansion is fitted using numerical simulations. All these approaches are in some way ‘effective’ in the sense that they encode the small scale effects in parameters which need to be fitted using numerical simulations. In order to extent the description of dark matter to times and scales beyond shell-crossing, a Vlasov-Einstein system [79, 80] as well as Schrödinger-Poisson methods have been studied [81, 82, 83, 84, 85]. More recently, advances have been made in the Lagrangian picture of the Vlasov-Poisson system [86, 87], especially in terms of 1+1 dimensional gravity [88, 89, 90, 91, 92, 93], where Zel’dovich dynamics is exact up to shell-crossing, and recently even beyond the 1+1 dimensional setting [91]. Rather than pursuing down this road we extent the Eulerian framework of the Vlasov-Poisson system, similarly as has been started in reference [94].

Let us recall that the insufficiency of the single-stream approximation is due to the truncation of the infinite Vlasov hierarchy. The single-stream approximation seems to be self-consistent if velocity dispersion is absent, but is unstable under perturbations. Including the second cumulant of the velocity distribution, namely velocity dispersion, introduces new degrees of freedom which allow to account for multiple streams within the fluid description. While velocity dispersion effects should be relatively small for cold dark matter, their corrections are expected to kick in at small enough scales. As we show, even a small seed of velocity dispersion can lead to accountable corrections at late times and small scales.

The other aspect we want to start investigating in this paper are cosmological vector and tensor fields. The single-stream approximation is described by two scalar fields, namely the density and velocity-divergence, while the only vector field present is the curl of the velocity field, the vorticity field. It is well known that at linear level the vorticity field decays due to the expansion of the Universe. In the single-stream approximation there are no non-linear terms which source vector fields from scalar initial conditions only. However, after including a non-vanishing velocity dispersion tensor we find that vector and tensor fields are non-linearly sourced by scalar fields. In the present paper we treat the full set of four scalar, two vector and one tensor matter degrees of freedom on linear level. In future follow-up work we plan to address non-linear dynamics in this setup.

2 The Vlasov-Poisson system for dark matter

2.1 Dark matter distribution function, moments and cumulants

We assume that dark matter consists of non-relativistic, classical point particles of mass $m$ which interact only gravitationally, at least to leading approximation. We work in an approximation where the metric is in the conformal Newtonian gauge [95]

[TABLE]

with the scale factor $a(\tau)$ and peculiar Newtonian gravitational potential $\phi(\tau,\mathbf{x})$. The conformal time $\tau$ is related to standard cosmic time $t$ by $dt=a(\tau)\,d\tau$. The gravitational potential is assumed to be relatively small on the scales we are interested in and for $\phi(\tau,\mathbf{x})=0$ the space-time has the symmetries of a flat Friedmann-Lemaître-Robertson-Walker cosmology.111These symmetries actually become statistical symmetries for the full theory where the expectation values are either taken as spatial averages over large enough volumes or for conveniently defined ensembles of cosmologies. For a cosmology accounting for a cosmological constant and one dark matter species only, the evolution of the scale factor is determined by the Friedmann equations [96]

[TABLE]

assuming dark matter to be pressureless at large enough scales. Here $\mathcal{H}(\tau)\equiv\dot{a}(\tau)/a(\tau)$ is the conformal Hubble parameter where we denote partial derivatives with respect to $\tau$ by a dot. Moreover, $\Omega_{m}(\tau)$ and $\Omega_{\Lambda}(\tau)$ are the dark matter and cosmological constant density parameters at time $\tau$, respectively. On scales much smaller than the Hubble horizon $1/\mathcal{H}(\tau)$, the gravitational potential is in first order approximation subject to Poisson’s equation

[TABLE]

To simplify notation we suppress time and space arguments when unambiguous. The dark matter mass density fluctuation field $\delta(\tau,\mathbf{x})$ is defined below.

Dark matter is described in the framework of a kinetic theory approximation by the one-particle phase-space distribution function $f(\tau,\mathbf{x},\mathbf{p})$ such that the number of particles in the comoving phase-space volume $d^{3}x\,d^{3}p$ at $(\mathbf{x},\mathbf{p})$ and time $\tau$ is $f(\tau,\mathbf{x},\mathbf{p})\,d^{3}x\,d^{3}p$. Single dark matter particles move along the characteristics [97]

[TABLE]

such that the dark matter distribution function is conserved along them which leads to the Vlasov equation

[TABLE]

The Vlasov-Poisson system of equations (2.3) and (2.5) is closed by the dark matter mass density field

[TABLE]

In the following it is convenient to split the density field into a spatially homogeneous expectation value $\bar{\rho}(\tau)\equiv\left\langle\rho(\tau,\mathbf{x})\right\rangle\propto a^{-3}(\tau)$ which determines the dark matter density parameter $\Omega_{m}(\tau)$ and the locally varying mass density fluctuation field

[TABLE]

The latter enters Poisson’s equation (2.3) such that the Vlasov-Poisson system of equations is non-linear.

The Vlasov-Poisson system of equations specifies the time evolution of the distribution function $f(\tau,\mathbf{x},\mathbf{p})$ and the metric parametrised by the gravitational potential $\phi(\tau,\mathbf{x})$. A theoretical drawback when considering the distribution function is that it is a function of seven variables which makes it naturally difficult to solve the Vlasov-Poisson system of equations. Further, in a statistical description where one considers expectation values $\left\langle f(\tau,\mathbf{x},\mathbf{p})\right\rangle$ or correlation functions $\left\langle f(\tau,\mathbf{x},\mathbf{p})\,f(\tau^{\prime},\mathbf{x}^{\prime},\mathbf{p}^{\prime})\right\rangle$, the symmetries reduce the complexity only little. On the other side, from an observational point of view, one is often not interested in the full phase-space distribution function but rather in moments and cumulants with respect to the momentum argument. The zeroth moment of the distribution function is the mass density field (2.6). The first moment is given by the momentum density field

[TABLE]

where $u_{i}(\tau,\mathbf{x})$ is the peculiar velocity field. The second moment is the momentum flux density field

[TABLE]

where $T_{ij}(\tau,\mathbf{x})$ is the stress tensor field. Higher moments of the distribution function can be defined very similar to (2.6), (2.8) and (2.9) and in close analogy to the usual definitions for probability distributions. It is convenient to split the stress tensor field $T_{ij}(\tau,\mathbf{x})=\rho(\tau,\mathbf{x})\,\sigma_{ij}(\tau,\mathbf{x})$ where $\sigma_{ij}(\tau,\mathbf{x})$ is the peculiar velocity dispersion tensor field which quantifies the deviation of particle velocities from the velocity field $u_{i}(\tau,\mathbf{x})$. Since only the density fluctuation field $\delta(\tau,\mathbf{x})$ enters Poisson’s equation (2.3), it is convenient to discuss moments and cumulants with respect to the momentum argument of the redefined distribution function

[TABLE]

which obeys the Vlasov equation (2.5) because it is linear in $f(\tau,\mathbf{x},\mathbf{p})$ and the density background decays $\bar{\rho}(\tau)\propto a^{-3}(\tau)$. Moments $m_{i_{1}...i_{n}}^{(n)}(\tau,\mathbf{x})$ of any order $n\in\mathbb{N}_{0}$ of the distribution function with respect to the momentum argument parametrise spatial matter degrees of freedom of the fluid. They can be obtained from the moment-generating function

[TABLE]

by applying derivatives with respect to $\mathbf{l}$ and evaluating at $\mathbf{l}=\mathbf{0}$,

[TABLE]

where the indices $i_{j}\in\{1,2,3\}$ are spatial and totally symmetric for all $j\in\{1,...,n\}$. Similarly, cumulants $c_{i_{1}...i_{n}}^{(n)}(\tau,\mathbf{x})$ of the distribution function are statistically independent quantities which represent the connected part of the moments. They can be obtained from the cumulant-generating function C(\tau,\mathbf{x};\mathbf{l})\equiv\ln\big{(}M(\tau,\mathbf{x};\mathbf{l})\big{)} in the same manner moments are derived from $M(\tau,\mathbf{x};\mathbf{l})$ in equation (2.12). The first few cumulants are given by

[TABLE]

With the definition of the cumulant-generating function and the Vlasov equation (2.5) one can derive the equation of motion for the cumulant-generating function [76, 85],

[TABLE]

From the structure of equation (2.14) one can already infer the most important properties of the time evolution of cumulants of any order $n$:

• •

The second term on the left-hand side generates the Hubble drag term accounting for the expansion of the Universe, which scales with the number of spatial indices $n$.

• •

The third term on the left-hand side involves the next higher order cumulant, creating an infinity hierarchy of coupled equations.

• •

The fourth term on the left-hand side is quadratic in the cumulants and generates a sum over the next higher and all lower order cumulants.

• •

The last term on the left-hand side only contributes at $n=1$ and couples Poisson’s equation (2.3) to the evolution of the first cumulant.

To have a finite set of closed equations one needs to truncate the Vlasov hierarchy in an appropriate way. From the non-linear term in the equations of motion (2.14) it is clear, that the only (apparent) consistent truncation is neglecting all cumulants of order $n\geq 2$, because there are no source terms in the corresponding equations of motion which solely depend on lower order cumulants. This corresponds to assuming dark matter to behave like a perfect pressureless fluid and amounts to the generating functions

[TABLE]

which correspond to the distribution function

[TABLE]

where $\delta_{\text{D}}^{(n)}(x)$ is the Dirac delta function in $n$ dimensions and we enforce the Einstein summation convention. This so-called single-stream approximation assigns a single velocity field $\mathbf{u}(\tau,\mathbf{x})$ to each point in configuration space. While the single-stream approximation is self-consistent, it is unstable under perturbations. Any deviation from the distribution function (2.16) naturally generates all higher order moments [76]. Further, the single-stream approximation can not account for multiple streams within the fluid and is therefore unable to describe the phenomenon of shell-crossing, when trajectories of particles meet in configuration space. At this point the velocity field in equation (2.16) is multivalued and all roots of the Dirac delta function contribute. This amounts to a distribution function of the form

[TABLE]

where $s$ runs over all streams with density and velocity field, $\delta_{s}(\tau,\mathbf{x})$ and $\mathbf{u}_{s}(\tau,\mathbf{x})$, respectively. From the corresponding cumulant-generating function it is clear that this form of distribution function induces cumulants of any order [76].

In order to overcome the perfect pressureless fluid approximation we truncate the Vlasov hierarchy at the next higher cumulant, thereby accounting for velocity dispersion within the fluid. This amounts to the cumulant-generating function

[TABLE]

and thus to the moment-generating function

[TABLE]

The corresponding distribution function is the three-dimensional normal distribution

[TABLE]

where this form only holds if the isotropic velocity dispersion is sufficiently larger than anisotropic velocity dispersion, such that $\sigma_{ij}(\tau,\mathbf{x})$ is positive definite. A few comments are in order:

• •

The momenta of single dark matter particles are normal distributed around the mean velocity field $u_{i}(\tau,\mathbf{x})$ with covariance $\sigma_{ij}(\tau,\mathbf{x})$.

• •

The distribution function (2.20) can capture the phenomenon of shell-crossing.

• •

The single-stream approximation (2.16) is recovered in the limit $\sigma_{ij}(\tau,\mathbf{x})\to 0$.

• •

In thermal equilibrium we expect $u_{i}(\tau,\mathbf{x})=0$ and $\sigma_{ij}(\tau,\mathbf{x})=\delta_{ij}\,k_{\text{B}}T/m$, such that momenta are Maxwell-Boltzmann distributed

[TABLE]

There are indications that a Maxwell-Boltzmann distribution also describes the dark matter velocity distribution within bond structures such as halos [98, 99].

From the discussion of the equation of motion (2.14) we know that the Vlasov equation (2.5) does not conserve this approximation, because higher order cumulants are dynamically generated. However, cum grano salis we assume this truncation to be appropriate for times and scales where the dark matter distribution function is perturbed around thermal equilibrium sufficiently close to (2.20) and henceforth set all higher order cumulants zero.

Thus, the set of equations we are taking into account are the continuity equation

[TABLE]

the Cauchy momentum equation

[TABLE]

and the velocity dispersion equation

[TABLE]

Here we use a notation where a comma indicates partial derivates with respect to the components of $\mathbf{x}$.

2.2 Statistical description and background-fluctuation splitting

We are interested in a statistical description of dark matter formulated in terms of expectation values and correlation functions of the cosmological fields. The latter are given by the gravitational potential $\phi(\tau,\mathbf{x})$ as well as the matter fields defined by the zeroth moment (2.6) and the first two cumulants (2.13) of the dark matter distribution function $f(\tau,\mathbf{x},\mathbf{p})$.

The statistical field theory we employ can be seen as describing an ensemble of cosmological histories with stochastic initial conditions or, equivalently for this purpose, an ensemble of different spatial subvolumes of a single cosmological history. We assume that the symmetries of the idealised flat Friedmann-Lemaître-Robertson-Walker cosmology are realised as statistical symmetries for this field theory which concerns in particular spatial translations and rotations.

In the following it is convenient to split the cosmological fields into a background and a fluctuation part. The background is taken to be spatially homogeneous and isotropic and corresponds to the expectation value of the fields. For the density field we have already done such a splitting in definition (2.7). For the velocity field $u_{i}(\tau,\mathbf{x})$ the statistical rotation symmetry forbids an expectation value such that it is a pure fluctuation field. For the velocity dispersion tensor field defined in (2.9) one may have a non-vanishing isotropic expectation value,

[TABLE]

Note that $\bar{\sigma}(\tau)$ is positive semi-definite and in the following we assume $\bar{\sigma}(\tau)>0$, although the value could be very small. Similar to the density fluctuation field (2.7) we define the peculiar velocity dispersion tensor fluctuation field

[TABLE]

The velocity dispersion background $\bar{\sigma}(\tau)$ obeys the expectation value of the trace of the velocity dispersion equation (2.24), given by

[TABLE]

The evolution of the velocity dispersion background depends on equal-time two-point correlation functions of the velocity and the velocity dispersion fluctuation fields.

Note that the expectation value (2.25) differs from a possible expectation value of the pressure of the dark matter fluid which would instead be given by

[TABLE]

While one might wonder if directly working with the pressure background is more sensible, we do not see an a priori reason why this should be the case. In terms of the background-fluctuation splitting (2.7) and (2.26) the dark matter pressure and velocity dispersion background can be related

[TABLE]

This pressure background enters the Friedmann equations (2.2), but is suppressed by powers of the speed of light.

In terms of the background-fluctuation splitting (2.25) and (2.26), the Cauchy momentum equation (2.23) and velocity dispersion equation (2.24) can be written as

[TABLE]

and

[TABLE]

The last term on the left side of equation (2.31) is a background term which, together with the evolution equation (2.27), ensures that $\langle\dot{\varsigma}_{ij}(\tau,\mathbf{x})\rangle=0$. Note that $\bar{\sigma}(\tau)$ enters the equations of motion (2.30) and (2.31) and that the evolution equation (2.27) depends on the two-point correlation function of the fluctuation fields. In this sense, the background and fluctuation fields are coupled and need to be solved simultaneously.

2.3 Scalar, vector and tensor decomposition

For further analysis we decompose the matter fields according to their transformation properties under spatial rotations into scalar, vector and tensor fields. This is most conveniently done by moving to Fourier space using the convention

[TABLE]

for any integrable field $f(\mathbf{x})$ where we abbreviate configuration and Fourier space integrals by

[TABLE]

The velocity fluctuation field $u_{i}(\tau,\mathbf{k})$ can be decomposed into an irrotational and a solenoidal vector field. The irrotational part of the decomposition can be parametrised in terms of the scalar peculiar velocity-divergence field

[TABLE]

and describes the potential flow of the dark matter fluid. The remaining solenoidal part can be parametrised by the peculiar vorticity field

[TABLE]

which quantifies the local rotation of fluid elements. The vorticity field $\omega_{i}(\tau,\mathbf{k})$ is a pseudovector field which is transverse to the flow direction $\mathbf{k}$. In terms of the velocity-divergence and vorticity fields, the decomposition takes the form

[TABLE]

Similarly, the symmetric second-rank velocity dispersion tensor fluctuation field $\varsigma_{ij}(\tau,\mathbf{k})$ can be decomposed into two scalar fields, a solenoidal vector field and a symmetric, transverse and traceless second-rank tensor field [100]. We parametrise the scalar fields in terms of the trace field

[TABLE]

and the off-trace field

[TABLE]

which quantify isotropic and anisotropic velocity dispersion, respectively. The solenoidal vector field we parametrise by

[TABLE]

and the symmetric, transverse and traceless second-rank tensor field by

[TABLE]

both of which describe anisotropic velocity dispersion degrees of freedom. Here we introduced the symmetric, transverse and traceless projector

[TABLE]

defined in terms of the transverse projector

[TABLE]

The decomposition of the velocity dispersion tensor fluctuation field $\varsigma_{ij}(\tau,\mathbf{k})$ can then be written as

[TABLE]

All ten matter degrees of freedom are described by the four scalar fields \big{(}\delta,\theta,\varsigma,\vartheta\big{)}, the two solenoidal vector fields \big{(}\omega_{i},\vartheta_{i}\big{)} and the symmetric, transverse and traceless second-rank tensor field $\vartheta_{ij}$.

In order to rewrite the evolution equation (2.27) in terms of expectation values of the decomposed fluctuation fields it is convenient to define the equal-time power spectra

[TABLE]

The evolution equation (2.27) can then be rewritten to

[TABLE]

where we defined the dimensionless function

[TABLE]

From equation (2.45) it is evident that the velocity dispersion background is sourced by the covariance of the velocity and velocity dispersion fields, or, physically speaking, by inhomogeneities in the velocity and velocity dispersion distribution. In this sense we henceforth refer to the function $Q(\tau)$ as source function.

To cast the Cauchy momentum equation (2.30) into a form which is a polynomial of finite degree in the fluctuation fields, we expand the logarithmic density field $\ln(1+\delta(\tau,\mathbf{x}))$ around the background configuration $\delta(\tau,\mathbf{x})=0$ by a Taylor series up to terms cubic in $\delta(\tau,\mathbf{x})$ and neglect all higher orders.222Higher order terms in the fluctuation fields do not enter one-loop calculations within the framework of renormalised cosmological perturbation theory [17]. By splitting the equations of motion (2.22), (2.30) and (2.31) for $\mathbf{k}\neq 0$ according to the decompositions (2.36) and (2.43) and using Poisson’s equation (2.3) as well as the evolution equation (2.45) one is left with the continuity equation

[TABLE]

the velocity-divergence equation

[TABLE]

the vorticity equation

[TABLE]

and the four velocity dispersion equations

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here we abbreviated terms that are quadratic and cubic in the fluctuation fields which are listed in appendix A.

The assumption of a non-vanishing $\bar{\sigma}(\tau)$ breaks the apparent self-consistency of the single-stream approximation since already at linear level the velocity dispersion fields $\varsigma(\tau,\mathbf{x})$, $\vartheta(\tau,\mathbf{x})$ and $\vartheta_{i}(\tau,\mathbf{x})$ are sourced by the velocity fields $\theta(\tau,\mathbf{x})$ and $\omega_{i}(\tau,\mathbf{x})$, respectively. Further, the single-stream approximation does not allow to generate the vorticity field from scalar initial conditions because $C_{\theta\omega}^{ij}(\mathbf{k}_{1},\mathbf{k}_{2})$ and $C_{\omega\omega}^{ijk}(\mathbf{k}_{1},\mathbf{k}_{2})$ are the only vertices present. In contrast, by including velocity dispersion, it is possible to source vector and tensor fields non-linearly from scalar initial conditions through the vertices $C_{\varsigma\delta}^{i}(\mathbf{k}_{1},\mathbf{k}_{2})$, $C_{\vartheta\delta}^{i}(\mathbf{k}_{1},\mathbf{k}_{2})$, $C_{\varsigma\delta\delta}^{i}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})$, $C_{\vartheta\delta\delta}^{i}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})$, $F_{\varsigma\theta}^{i}(\mathbf{k}_{1},\mathbf{k}_{2})$, $F_{\vartheta\theta}^{i}(\mathbf{k}_{1},\mathbf{k}_{2})$, $G_{\varsigma\theta}^{ij}(\mathbf{k}_{1},\mathbf{k}_{2})$, $G_{\vartheta\theta}^{ij}(\mathbf{k}_{1},\mathbf{k}_{2})$.

Note that $\bar{\sigma}(\tau)$ enters the velocity-divergence and vorticity equations (2.48) and (2.49) in the combination $k^{2}\bar{\sigma}(\tau)$, as is required by dimensional arguments. This introduces a scale dependence even at linear order in the fluctuation fields which is crucial for the convergence of the source function $Q(\tau)$ as we explain in section 4.2.

2.4 Compact notation

To cast the equations of motion (2.47) – (2.53) into a more symmetric and compact form we introduce the time evolution parameter

[TABLE]

corresponding to the number of $e$-folds of the linear growth function in the single-stream approximation [101]

[TABLE]

normalised to unity at some initial time $\tau_{\text{in}}$.333We assume a cosmology with a cosmological constant and one sufficiently cold dark matter species without radiative component. Here the growth function is given in terms of the Gaussian hypergeometric function ${}_{2}F_{1}(a,b;c;x)$ and the dark matter density parameter is written as $\Omega_{m}(\tau)=\Omega_{m,0}/a^{3}(\tau)$ where $\Omega_{m,0}$ is the value at $a(0)=1$, corresponding to today. Further it is useful to define the dimensionless function

[TABLE]

which parametrises the logarithmic deviation of $D_{+}(\tau)$ from the Einstein-de Sitter linear growth function $a(\tau)$ and obeys the evolution equation [77]

[TABLE]

In the following we use $\eta(\tau)$ rather than $\tau$ which are related through definition (2.54).

Grouping scalar, vector and tensor fields in a condensed notation, we define the field vector

[TABLE]

The component fields are referred to by $\Psi_{a}(\eta,\mathbf{k})$ where the index $a$ runs over all fields and implicitly carries possible spatial indices which are only displayed if necessary. We use the Einstein summation convention which implies summation over possible carried vector or tensor indices. Scalar, vector and tensor fields are referred to by the indices $a_{s}$, $a_{v}$ and $a_{t}$ which run over the scalar, vector and tensor field components only. When only the first two components of $\Psi_{a}(\eta,\mathbf{k})$ are considered, that is the density and velocity-divergence fields, the analysis reduces to that of perfect pressureless fluid dark matter [103, 17, 18].

Similar to before we define the equal-time power spectra

[TABLE]

which depend only on the wave number $k$ due to the statistical homogeneity and isotropy symmetry. The power spectrum (2.59) is naturally block diagonal for scalar, vector and tensor fields and is understood to implicitly carry a possible projector for the vector and tensor fields, i.e.

[TABLE]

It is convenient to also define the rescaled background

[TABLE]

and introduce the rescaled source

[TABLE]

such that the evolution equation (2.45) can be written as

[TABLE]

The equations of motion (2.47) – (2.53) can be cast into the form

[TABLE]

The coefficients of the terms linear in the fluctuation fields are given by the matrix $\Omega_{ab}(\eta,k)$, which is block diagonal for scalar, vector and tensor fields. The scalar submatrix is given by

[TABLE]

while the vector and tensor submatrices are

[TABLE]

Similar to the power spectra (2.60) the submatrices $\Omega_{a_{v}b_{v}}(\eta,k)$ and $\Omega_{a_{t}b_{t}}(\eta,k)$ are understood to implicitly carry a vector or tensor projector, respectively. The non-linear terms quadratic in the fluctuation fields are given by

[TABLE]

whereas the terms cubic in the fluctuation fields are

[TABLE]

Here $Y_{abc}(\eta,\mathbf{k},\mathbf{k}_{1},\mathbf{k}_{2})$ and $X_{abcd}(\eta,\mathbf{k},\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})$ are symmetrised versions of the vertices appearing in equations (A.1) – (A.4).

3 Evolution of the background and linear fluctuations

3.1 Linear response

At early times and large scales we assume the fields $\Psi_{a}(\eta,\mathbf{k})$ to be small compared to the homogeneous and isotropic background. Therefore, it seems natural to study how the fields evolve in the linear regime where we neglect the coupling of different Fourier modes. In the limit where we turn off interactions, the equations of motion (2.64) reduce to

[TABLE]

While the equations of motion (3.1) describe linear fields, in the sense that no interaction between Fourier modes is taken into account, they nevertheless form a set of non-linear differential equations. Because the velocity dispersion background $\hat{\sigma}(\eta)$ enters $\Omega_{ab}(\eta,k)$ the evolution equation (2.63) is needed to close the system of equations, which in turn depends on the correlation of the fluctuation fields through the source $\hat{Q}(\eta)$. In this situation one cannot expect to find solutions in a closed form in terms of standard mathematical functions.

To fully determine the system we need to specify initial conditions at some initial time evolution parameter ${\eta_{\text{in}}}\equiv\eta(\tau_{\text{in}})$.444By definition (2.54) the initial time evolution parameter is ${\eta_{\text{in}}}\equiv 0$. In the following we use the approximation $\Omega_{m}(\tau)/f^{2}(\tau)\approx 1$, which is sufficiently accurate for the times we are interested in [102]. For an initial velocity dispersion background $\hat{\sigma}_{\text{in}}$, the evolution equation (2.63) can be solved with

[TABLE]

The fields $\Psi_{a}(\eta,\mathbf{k})$ can formally be written as a linear response to the initial conditions $\Phi_{a}(\mathbf{k})$,

[TABLE]

The response function $g_{ab}^{R}(\eta,{\eta_{\text{in}}},k)$ is subject to the defining equation

[TABLE]

together with the causal boundary conditions

[TABLE]

and is henceforth called the retarded linear propagator. With the initial power spectrum

[TABLE]

the source (2.62) can be written as

[TABLE]

We emphasise that equation (3.4) is non-linear and the retarded linear propagator is not a Green’s function in the standard sense of linear differential equations, because $g_{ab}^{R}(\eta,{\eta_{\text{in}}},k)$ depends on the initial velocity dispersion background $\hat{\sigma}_{\text{in}}$ and is a functional of the initial power spectrum $P_{ab}^{\text{in}}(k)$. As a consequence, the linear fields do not grow independently for different Fourier modes in contrast to the single-stream approximation, because $\hat{Q}(\eta)$ is sensitive to the evolution of the fluctuation fields for all wave numbers $k$. Nevertheless, equation (3.3) can be seen as a linear response to a small perturbation in the initial random field.

3.2 Initial conditions

In the following we assume that the initial vector and tensor fields can be neglected and that the initial scalar fields can be written as proportional random fields,

[TABLE]

where $\delta^{\text{in}}(\mathbf{k})$ is the initial dark matter mass density fluctuation field. The vector $w_{a_{s}}(k)$ quantifies the proportionality of the initial scalar fields to $\delta^{\text{in}}(\mathbf{k})$ and we chose it to be the eigenvector of the growing solution given in section 3.3. Further, we assume $\delta^{\text{in}}(\mathbf{k})$ to be a Gaussian random field which is completely characterised by the initial scalar power spectrum

[TABLE]

It is particularly convenient to discuss solutions in terms of the linear growth functions

[TABLE]

such that the equal-time linear evolved scalar power spectra can be written as

[TABLE]

For a vanishing initial vector power spectrum the source for the velocity dispersion background takes the form

[TABLE]

which depends on the scalar linear growth functions only.

Since $\Omega_{ab}(\eta,k)$ is block diagonal for scalar, vector and tensor fields, and because $\hat{Q}(\eta)$ depends on the scalar linear growth functions only, we can solve the scalar part of equation (3.4) together with the velocity dispersion background (3.2) independent from the vector and tensor parts. Subsequently the vector and tensor part can be solved.

3.3 Scalar growth factors and free-streaming scale

Because the scalar retarded linear propagator is a functional of the initial power spectrum, we need to specify $w_{a_{s}}(k)$ in order to solve the scalar part of equation (3.4). To make a sensible choice we consider the eigenvalues of the submatrix $-\Omega_{a_{s}b_{s}}(\eta,k)$ which determine the growth of the solutions at time $\eta$ and wave number $k$ and are henceforth called growth factors. Three growth factors are given by the roots of the polynomial

[TABLE]

whereas the fourth is $s_{s}(\eta,k)=-1$. In the limit $k\to 0$, one can easily find the roots of the polynomial (3.13) and obtains

[TABLE]

To fully characterise the solutions the corresponding eigenvectors are needed. These are given by

[TABLE]

for $n\in\{1,2,3\}$.

In the limit $k\to 0$ the density and velocity-divergence fields decouple from the velocity dispersion fields. In turn the former evolve according to the dynamics of the single-stream approximation and obey a growing and decaying solution proportional to $e^{\eta-{\eta_{\text{in}}}}$ and $e^{-\frac{3}{2}(\eta-{\eta_{\text{in}}})}$ [103] while the latter decay proportional to $e^{-(\eta-{\eta_{\text{in}}})}$. The growth factors $s_{s}^{(1)}(\eta,0)$ and $s_{s}^{(2)}(\eta,0)$ can be identified with the single-stream approximation solutions, whereas $s_{s}^{(3)}(\eta,0)$ and $s_{s}^{(4)}(\eta,0)$ correspond to the two new decaying solutions.

For finite values of $k^{2}\,\hat{\sigma}(\eta)$ the picture is more complicated. In figure 1 we show the real and imaginary parts of the scalar field growth factors as a function of the combination $k^{2}\,\hat{\sigma}(\eta)$. The natural extension of the single-stream approximation growth factor $s_{s}^{(1)}(\eta,k)$ is the only positive growth factor in the range $0\leq k^{2}\,\hat{\sigma}(\eta)<3/2$ and turns negative for larger values of $k^{2}\,\hat{\sigma}(\eta)$. While $s_{s}^{(1)}(\eta,k)$ and $s_{s}^{(4)}(\eta,k)$ are real for all $k^{2}\,\hat{\sigma}(\eta)$ the growth factors $s_{s}^{(2)}(\eta,k)$ and $s_{s}^{(3)}(\eta,k)$ develop an imaginary part at some $k^{2}\,\hat{\sigma}(\eta)>0$ from where on their real parts stay negative and coincide.

We emphasis that due to the time dependence of $\Omega_{ab}(\eta,k)$ the retarded linear propagator is not time translation invariant and thus the eigenvectors generally differ from those of $\Omega_{ab}(\eta,k)$. Therefore exciting a solution proportional to an eigenvector (3.15) at initial time ${\eta_{\text{in}}}$ naturally evolves to a superposition of different eigenvectors at time $\eta$. Further, the full time dependence of a solution is determined by the growth factors as well as their corresponding eigenvectors.

In the following we concentrate on the growing solution characterised by the growth factor $s_{s}^{(1)}(\eta,k)$. We assume that initially only the growing solution is present and the corresponding eigenvector determines the initial scalar field configuration, i.e.

[TABLE]

Before solving the system of equations (3.2) and (3.4) by numerical means, it is sensible to infer some properties of the solution. The growth factor $s_{s}^{(1)}(\eta,k)$ is a monotonically decreasing function of $k$ and has a zero-crossing at the comoving free-streaming wave number

[TABLE]

and is negative for all $k>k_{\text{fs}}(\eta)$. The wave number $k_{\text{fs}}(\eta)$ determines the scale below which the fluctuation fields do not grow due to the free-streaming of dark matter particles. The corresponding solution is therefore damped at smaller scales and decays for wave numbers larger than $k_{\text{fs}}(\eta)$. In order to qualitatively understand how $k_{\text{fs}}(\eta)$ evolves in time, we consider the limit of early and late times. At early times the fluctuation fields are small in amplitude and one can neglect the enhancement of the velocity dispersion background (3.2) due to $\hat{Q}(\eta)$ as long as the variance of the velocity-divergence and scalar velocity dispersion fields is smaller than the velocity dispersion background,

[TABLE]

In this limit the velocity dispersion background decays in time and in turn the free-streaming wave number is raised which allows structures to form at smaller scales. At later times, however, the fluctuation fields grow in amplitude and the validity of the limit (3.18) breaks down. Eventually, at times where

[TABLE]

the velocity dispersion background starts to grow and in turn the free-streaming wave number is lowered which stalls the formation of structure at small scales.

3.4 Vector and tensor growth factors

Similar to the case of scalar fields it is sensible to study properties of the vector retarded linear propagator in terms of the eigenvalues of the submatrix $-\Omega_{a_{v}b_{v}}(\eta,k)$. The growth factors for the vector fields are given by

[TABLE]

and the corresponding eigenvectors are

[TABLE]

for $n\in\{1,2\}$. In the limit $k\to 0$ one obtains

[TABLE]

Since the limit $k\to 0$ restores the equations of motion of perfect fluid dark matter by decoupling the vorticity field from the vector velocity dispersion field, one can identify $s_{v}^{(1)}(\eta,0)$ with the decaying solution of the single-stream approximation vorticity field proportional to $e^{-\frac{1}{2}(\eta-{\eta_{\text{in}}})}$. The growth factor $s_{v}^{(2)}(\eta,0)$ corresponds to a new decaying solution proportional to $e^{-(\eta-{\eta_{\text{in}}})}$, similar to the scalar case.

For finite values of $k^{2}\,\hat{\sigma}(\eta)$, the real and imaginary parts of the vector field growth factors are shown in figure 2. Both growth factors $s_{v}^{(1)}(\eta,k)$ and $s_{v}^{(2)}(\eta,k)$ are real and negative for sufficient small $k^{2}\,\hat{\sigma}(\eta)$. For larger $k^{2}\,\hat{\sigma}(\eta)$ they develop an imaginary part from where on their real parts coincide and take a constant negative value. The growth factor $s_{v}^{(1)}(\eta,k)$ is monotonically decreasing while $s_{v}^{(2)}(\eta,k)$ is monotonically increasing for larger $k^{2}\,\hat{\sigma}(\eta)$ until they develop an imaginary part.

Because the vector part of equation (3.4) is a first-order linear ordinary differential equation, any general solution can be constructed by superposition of the solutions characterised by the growth factors $s_{v}^{(n)}(\eta,k)$ for $n\in\{1,2\}$.555This is only true as long as we neglect initial vector fields which would enter the source of the velocity dispersion background by definition (2.62) and thus couple the scalar and vector part of equation (3.4). The corresponding initial field configurations are determined by the eigenvectors, i.e.

[TABLE]

Finally, the tensor part of equation (3.4) together with the submatrix (2.66) is independent of $k$ and the tensor retarded linear propagator can therefore be given in the closed form

[TABLE]

The corresponding tensor growth factor is $s_{t}(\eta,k)=-1$ for which we simply chose $w_{a_{t}}(k)=\hat{\sigma}_{\text{in}}$ as the initial tensor field configuration.

3.5 Approximation for small velocity dispersion and soft power spectrum

In the following we discuss an approximate analytical solution that becomes available in the limit where the velocity dispersion is small and where the initial scalar power spectrum $P_{s}^{\text{in}}(k)$ does not extend too far into the ultraviolet, in a sense that becomes clear below.

Assuming that the velocity dispersion background and fluctuation fields are small, we consider the equations (3.4) at wave numbers $k^{2}\,\hat{\sigma}(\eta)\ll 1$, $k^{2}\,\varsigma(\eta)\ll 1$, $k^{2}\,\vartheta(\eta)\ll 1$ and $k^{2}\,\vartheta_{i}(\eta)\ll 1$. Since the evolution of the velocity dispersion background depends on $\hat{Q}(\eta)$, which involves an integral over all wave numbers, this assumption does on first sight not seem sensible. However, as discussed in section 4, the initial scalar power spectrum $P_{s}^{\text{in}}(k)$ is heavily suppressed for wave numbers larger than the comoving free-streaming wave number at matter-radiation equality $k_{\text{fs,eq}}$. The above proposed approximation can therefore be used as long as $k_{\text{fs,eq}}^{2}\,\hat{\sigma}(\eta)\ll 1$, $k_{\text{fs,eq}}^{2}\,\varsigma(\eta)\ll 1$, $k_{\text{fs,eq}}^{2}\,\vartheta(\eta)\ll 1$ and $k_{\text{fs,eq}}^{2}\,\vartheta_{i}(\eta)\ll 1$.

By dropping the $k$ dependent terms in equations (3.4) and denoting objects in this approximation by a tilde, the equations of motion for the linear growth functions can be written as

[TABLE]

where $\tilde{\Omega}_{ab}(\eta)$ is given by the submatrices

[TABLE]

The source function (3.12) takes the form

[TABLE]

where

[TABLE]

is the dimensionless variance of the initial density fluctuation field and is the only quantity appearing in the approximation which is sensitive to the initial scalar power spectrum.666Notice the difference to the single-stream approximation which depends even at one-loop level only on the dimensionful parameter $\sigma_{v}^{2}\equiv 4\pi/3\int_{0}^{\infty}dq\,P_{s}^{\text{in}}(q)/(2\pi)^{3}$ [104]. For a perfectly cold dark matter power spectrum as in equation (4.1), the integral (3.28) diverges approximately logarithmically.

The initial conditions are taken to be the $k\to 0$ limit of the initial field configuration $w_{a}(k)$ specified in section 3.3 and 3.4. For scalar fields excited in the growing solution with $s_{s}^{(1)}(\eta,0)=1$ these are given by

[TABLE]

The vector fields corresponding to $s_{v}^{(1)}(\eta,0)=-1/2$ and $s_{v}^{(2)}(\eta,0)=-1$ are excited by

[TABLE]

respectively and for the tensor field corresponding to $s_{t}(\eta,0)=-1$ we use $w_{a_{t}}(0)=\hat{\sigma}_{\text{in}}$.

To close the system of equations (3.25) the velocity dispersion background

[TABLE]

is needed.

The linear growth functions of the density and velocity-divergence fields obey the growing solution

[TABLE]

while the scalar velocity dispersion linear growth functions are

[TABLE]

Here we defined the function

[TABLE]

where $C\equiv\sqrt{2\,\sigma_{d}^{2}/3}$ is a real positive constant. The vorticity field linear growth function is given by

[TABLE]

and the vector velocity dispersion linear growth function is

[TABLE]

for $n\in\{1,2\}$. Here we introduced the functions

[TABLE]

and

[TABLE]

where $\text{erf}(x)$ and $\text{erfi}(x)$ denote the real and imaginary error function, respectively. The tensor velocity dispersion linear growth function has the solution

[TABLE]

which coincides with $\tilde{g}_{6}^{R}(\eta)$ for $n=2$. Finally, the velocity dispersion background is given by

[TABLE]

With the analytical solutions (3.32), (3.33), (3.35), (3.36), (3.39) and (3.40) we can infer the early and late time behaviour of the linear growth functions and velocity dispersion background in the approximative scenario.

As for the single-stream approximation, the density and velocity-divergence fields grow at all times exponentially with the time evolution parameter $\eta$ while the vorticity field decays proportional to $e^{-\frac{1}{2}(\eta-{\eta_{\text{in}}})}$. The evolution of the velocity dispersion fields strongly depends on the value of $C$, but since we assume the power spectrum to not extent too far into the ultraviolet, we study the early and late time behaviour for $C\ll 1$. The scalar velocity dispersion fields grow at early times exponentially,

[TABLE]

and growth with double exponential power at late times

[TABLE]

According to the definitions (2.26), (2.37), (2.38), (2.39), (2.40) and (2.58) the velocity dispersion linear growth functions describe perturbations relative to the velocity dispersion background. The latter is therefore needed for a complete picture. According to the approximation (3.40) it decays exponentially at early times,

[TABLE]

while at late times one finds in the formal limit $\eta\to\infty$ the dominant behaviour to be a double exponential growth,

[TABLE]

This qualitatively agrees with the early and late times regimes (3.18) and (3.19) discussed in section 3.3.

4 Numerical results

4.1 Initial scalar power spectrum and velocity dispersion background

To solve the system of equations (3.2) and (3.4) numerically we need to specify the velocity dispersion background $\hat{\sigma}_{\text{in}}$ as well as the scalar power spectrum $P_{s}^{\text{in}}(k)$ at an initial time $\tau_{\text{in}}$. We set $\tau_{\text{in}}$ deep within the matter dominated era of the Universe at a redshift of $z(\tau_{\text{in}})\equiv 1/a(\tau_{\text{in}})-1=100$ and assume that the non-linear evolution of perturbations is negligible at earlier times, such that we can use a linearly evolved power spectrum as initial condition. Further we assume that the matter fluctuation fields are initially small enough in amplitude to be in the regime (3.18) where the velocity dispersion background simply decays in time.

We generate a linear evolved cold dark matter scalar power spectrum $P_{\text{CDM}}(k)$ using the CLASS code [105] for wave numbers $k\in\left[10^{-5},\,10\right]h$/Mpc assuming a $\Lambda$CDM cosmology using the parameters given in table 1. To extrapolate to smaller and larger wave numbers we use the prominent Eisenstein & Hu fitting formula [107]

[TABLE]

where

[TABLE]

and $\alpha,\beta,e$ are real positive constants.

For dark matter with non-vanishing velocity dispersion, structure formation is stalled at small scales due to free-streaming of dark matter particles. In the matter dominated era, structures are erased for wave numbers larger than the comoving free-streaming wave number at matter-radiation equality $k_{\text{fs,eq}}$ [108]. For wave numbers $k\gtrsim k_{\text{fs,eq}}$ the power spectrum is strongly suppressed and for simplicity we employ the ultraviolet cutoff $\Lambda\equiv k_{\text{fs,eq}}$ and use the initial scalar power spectrum

[TABLE]

where $\Theta_{\text{H}}^{(n)}(x)$ denotes the $n$-dimensional Heaviside unit step function.

For our numerical investigations we consider two dark matter candidates, namely sterile neutrinos of mass $m\sim 1$ keV and weakly interacting massive particles (WIMPs) of mass $m\sim 100$ GeV. The free-streaming wave number and velocity dispersion background at matter-radiation equality are given by

[TABLE]

for sterile neutrino dark matter and

[TABLE]

for WIMP dark matter [110]. We evolve the velocity dispersion background to $\tau_{\text{in}}$ using the evolution equation (2.63) in an Einstein-de Sitter cosmology ($\Omega_{m}(\tau)=1$, $f(\tau)=1$) and neglect the enhancement from $\hat{Q}(\eta)$. The initial velocity dispersion background is then $\hat{\sigma}_{\text{in}}=(a_{\text{eq}}/a_{\text{in}})\,\hat{\sigma}_{\text{eq}}$ where $a_{\text{eq}}\approx 1/3229$ is the scale factor at matter-radiation equality. This approximation is sensible as long as the initial time $\tau_{\text{in}}$ is well within the matter dominated era and $\hat{Q}({\eta_{\text{in}}})/\hat{\sigma}_{\text{in}}\ll 1$. Indeed, for $z(\tau_{\text{in}})=100$ we find $\hat{Q}({\eta_{\text{in}}})/\hat{\sigma}_{\text{in}}\approx 10^{-3}$ for sterile neutrino dark matter and $\hat{Q}({\eta_{\text{in}}})/\hat{\sigma}_{\text{in}}\approx 0.05$ for WIMP dark matter.

4.2 Results

In order to compare sterile neutrino and WIMP dark matter we show the former in the left panel and the latter in the right panel in all subsequent figures. Figures 3 and 4 display the time evolution of the velocity dispersion background and the zero mode of the fluctuation fields. Numerical solutions are shown as solid colored lines and the analytical approximations derived in section 3.5 as dashed black lines. Figures 5, 6 and 7 show the scale dependence of the density, velocity and velocity dispersion fluctuation fields, respectively, at three different redshifts.

Figure 3 shows that initially the velocity dispersion background decays in time for both dark matter candidates, as expected for the regime (3.18) where the contribution from fluctuations can be neglected for the evolution of $\hat{\sigma}(\eta)$. For later times however, the fluctuations grow in amplitude such that the regime (3.19) is entered and $\hat{\sigma}(\eta)$ receives an enhancing contribution from $\hat{Q}(\eta)$. Physically speaking, the formation of inhomogeneities and anisotropies of the velocity and velocity dispersion degrees of freedom source the growth of the velocity dispersion background. From figure 3 it is evident that the growth of $\hat{\sigma}(\eta)$ strongly depends on the dark matter model. For warm dark matter the free-streaming wave number $k_{\text{fs,eq}}$ is smaller than for cold dark matter and correspondingly $\hat{Q}(\eta)$ receives less contribution from larger wave numbers. For sterile neutrino dark matter the velocity dispersion background stays below its initial value while for WIMP dark matter it undergoes an immense growth at late times. For the latter we only show the evolution for $a(\tau)\leq 2/3$ because at times where $\hat{\sigma}(\eta)\gtrsim 1/k_{\text{fs,eq}}^{2}$ the validity of the analytical approximation breaks down and the source function $\hat{Q}(\eta)$ involves an integral over strongly oscillating power spectra which become difficult to control.

In figure 4 we show how fluctuations evolve in the large-scale limit $k\to 0$. Here the density and velocity-divergence growth functions (red) as well as the vorticity $n=1$ mode growth function (blue) exactly obey the dynamics of the perfect pressureless fluid approximation and thus obey a growing and decaying solution, respectively.777In the limit $k\to 0$ the $n=2$ mode of the vorticity field is not excited. The scalar velocity dispersion growth functions (yellow) are growing in time, but depend on the evolution of the velocity dispersion background. For sterile neutrino dark matter the growth of the scalar velocity dispersion is smaller than for the density and velocity-divergence fluctuations while for WIMP dark matter the former exceeds the growth of the latter by many orders of magnitude, originating from the immense growth of $\hat{\sigma}(\eta)$. Similarly, the vector $n=1$ mode velocity dispersion growth function (purple) is decaying in time for sterile neutrino dark matter while for WIMP dark matter the enhancement from $\hat{\sigma}(\eta)$ allows for a growing behaviour at later times. Finally, the vector $n=1$ mode and tensor velocity dispersion growth functions (red) are independent of the velocity dispersion background and simply decay in time for both dark matter candidates.

Knowing how the zero modes evolve in time we now turn to the scale dependence of the fluctuations. To this end figures 5, 6 and 7 display the growth functions normalised to their zero mode as a function of $k/k_{\text{fs,eq}}$. From the discussion of the growth factors in section 3.3 we expect fluctuations to be suppressed at wave numbers $k\gtrsim k_{\text{fs}}(\eta)$ due to free-streaming. For sterile neutrino dark matter $k_{\text{fs}}(\eta)>k_{\text{fs,eq}}$ because $\hat{\sigma}(\eta)<\hat{\sigma}_{\text{in}}$ at all times and therefore the growth functions depend weakly on redshift. For WIMP dark matter on the other hand $k_{\text{fs}}(\eta)\ll k_{\text{fs,eq}}$ at late times because $\hat{\sigma}(\eta)\gg\hat{\sigma}_{\text{in}}$ and concordantly the growth functions are sensible to small changes in redshift. For sterile neutrino dark matter we show results at redshift $z(\tau)=0$ (green line), $z(\tau)=40$ (yellow line) and $z(\tau)=80$ (red line) while for WIMP dark matter we choose $z(\tau)=0.5$ (green line), $z(\tau)=1$ (yellow line) and $z(\tau)=5$ (red line).

In figure 5 we show the density linear growth function. These are stronger suppressed at lower redshift and, as expected from the above discussion, for larger wave numbers. Sterile neutrino dark matter is weakly suppressed while WIMP dark matter is heavily suppressed for wave numbers $k\leq k_{\text{fs,eq}}$. We observe oscillations on top of the suppression for WIMP dark matter, which are not accounted for in the growing solution growth factor nor eigenvector. As mentioned in section 3.3 the eigenvectors of $\Omega_{ab}(\eta,k)$ and $g_{ab}^{R}(\eta,k)$ are in general not identical and thus the solution dynamically evolves away from being characterised by a single growth factor. The superposition of solutions of different growth factors lead to the observed oscillations which are due to the imaginary parts of the growth factors shown in figure 1.

Figure 6 displays the velocity-divergence (upper panel), the vorticity $n=1$ mode (middle panel) and the vorticity $n=2$ mode (lower panel) growth functions. Since the vorticity $n=2$ mode growth function vanishes in the $k\to 0$ limit, we show it not normalised. The velocity-divergence and vorticity $n=1$ mode growth functions are stronger suppressed for lower redshift and larger wave numbers. The vorticity $n=2$ mode growth function is suppressed for lower redshift and smaller wave numbers. This behaviour qualitatively matches the behaviour of the growth factors $s_{s}^{(1)}(\eta,k)$ and $s_{v}^{(1)}(\eta,k)$ which are suppressed for larger wave numbers while $s_{v}^{(2)}(\eta,k)$ is enhanced. Interestingly, on top of the suppressive and oscillatory behaviour of the WIMP dark matter fluctuations discussed above we observe an enhancement of the velocity-divergence growth function in the wave number range $3\cdot 10^{-3}\lesssim k/k_{\text{fs,eq}}\lesssim 4\cdot 10^{-1}$. This enhancement is responsible for the strong late time growth of the velocity dispersion background and the deviation from the analytical approximation. After discussing the velocity dispersion growth functions we comment more on this, since these are needed for a complete picture of the source function $\hat{Q}(\eta)$ which depends on the velocity-velocity dispersion correlations.

Figure 7 displays the isotropic scalar (upper panel), the vector $n=1$ mode (middle panel) and the vector $n=2$ mode (lower panel) velocity dispersion growth functions. The anisotropic scalar velocity dispersion growth function is not separately displayed since it is near to identical with the isotropic one. The scalar and vector $n=1$ mode velocity dispersion growth functions are stronger suppressed for lower redshift and larger wave numbers. Contrary the vector $n=2$ velocity dispersion mode growth function is suppressed for higher redshift and smaller wave numbers. Similar to the velocity-divergence growth function the scalar velocity dispersion growth function is enhanced in the wave number range $4\cdot 10^{-3}\lesssim k/k_{\text{fs,eq}}\lesssim 2\cdot 10^{-1}$. Together with the velocity-divergence growth function the strong growth of the velocity dispersion background may be qualitatively understood. For identical isotropic and anisotropic scalar velocity dispersion growth functions (at linear level these differ below the percent level) the source function (3.12) takes the form

[TABLE]

At large wave numbers the integral runs over the heavily oscillating growth functions which are enhanced compared to the analytical approximation (3.40). That is, at smaller scales the matter fluctuation fields oscillate due to the competition of gravitational collapse and free-streaming of dark matter and enhance the growth of the velocity dispersion background. Further, the suppression of the growth functions at large wave numbers is necessary for the convergence of $\hat{Q}(\eta)$. This becomes evident when taking perfectly cold initial conditions, $\Lambda\to\infty$, because the dimensionless variance (3.28) diverges in this limit. The scale dependence at linear level which induces the free-streaming suppression is therefore needed to obtain finite results.

5 Conclusion

In summary, we have studied an extension of the standard pressureless perfect fluid approximation for dark matter by taking the velocity dispersion tensor as an additional field into account. Neglecting all interaction of dark matter besides gravity, and assuming a kinetic description in terms of classical particles, we have obtained an equation of motion for the velocity dispersion tensor by taking the second cumulant of the Vlasov equation.

The set of evolution equations is closed by neglecting the third cumulant of the local dark matter velocity distribution. In a subsequent step, we have decomposed the velocity dispersion tensor into scalar, vector and tensor fields. Together with the density and velocity fields one obtains in total four scalar fields, two solenoidal vector fields and one symmetric, transverse and traceless second-rank tensor field.

Maybe the most significant modification of the single-stream approximation based on a perfect pressureless fluid picture comes from a homogeneous expectation value of the trace of the velocity dispersion tensor. On a technical level, the presence of such an expectation value has the consequence that it destroys the apparent self-consistency of the single-stream approximation. Perturbations in the scalar fields that parameterise the velocity dispersion tensor can be directly generated by other scalar perturbations. The extension to a fluid with velocity dispersion is important because it allows to account for shell-crossing and to address late times and small scales.

Interestingly, the evolution equation for the velocity dispersion expectation value contains a back-reaction term at quadratic order in perturbations. Physically, this term describes how velocity dispersion grows again at late times (after a decrease due to the cosmological expansion at earlier times) due to large inhomogeneities in particular at small scales. Technically, this back-reaction effect is proportional to an integral over the power spectrum of scalar perturbations that is strongly dominated by the small scale or ultraviolet regime. Beyond our present approximation, there is also a similar contribution involving an integration over the power spectrum of vector perturbations. This ultraviolet sensitivity makes it particularly interesting to compare different dark matter candidates.

We concentrate here in particular on a scenario where dark mater consists of rather light sterile neutrinos of mass $m\sim 1$ keV and weakly interacting massive particles with mass in the range $m\sim 100$ GeV. The latter are much colder than the former in the sense that they have initially, at matter-radiation equality, a much smaller velocity dispersion. However, this implies that the scalar power spectrum extends much further into the ultraviolet. As we show, the non-linear back-reaction effect is then rather strong at late times, leading to a strong increase of the velocity dispersion expectation value.

A similar computation was undertaken by reference [94]. Here a non-vanishing velocity dispersion background was introduced and evolved in a similar way as we did, by considering one scalar velocity dispersion fluctuation field, corresponding to the combination $(\varsigma(\tau,\mathbf{k})+\vartheta(\tau,\mathbf{k}))/2$. Reference [94] found a similar suppression of the power spectrum as well as a strong late time growth of the velocity dispersion background. Further a first non-linear approximation was made by including the dominant terms which showed that the growth of the velocity dispersion background persists even for non-linear mode coupling. Another path in the Lagrangian picture was pursed by reference [86]. Here the velocity dispersion tensor was solved for perturbatively and it was found that at third order growing contributions are possible. While this approach also features a velocity dispersion background as the zeroth order solution, it differs from our approach in the sense that we solve the background evolution equation fully non-linear, while solving the fluctuation fields linearly. The feature of a suppressed power spectrum due to a characteristic free-streaming scale was also found by reference [86].

In the present paper we have concentrated on a relatively simple approximation where the background is evolved non-linearly including the back-reaction effects from integrals over the power spectra of perturbations. The perturbations, or deviations from the homogeneous and isotropic background, are propagated linearly, however. While many qualitative effects can be well studied on this basis, one should keep in mind that non-linear effects in the evolution of perturbations can change the picture quantitatively, in particular at small scales and late times. Because of the ultraviolet dominance described above, this can also affect the evolution of the velocity dispersion expectation value. In future work we plan to take such effects into account. This could be done perturbatively but also using resummation schemes [16, 17, 18, 19, 20, 21, 22, 30, 23, 24, 25, 26, 27, 28, 29, 31, 32, 33, 34, 37, 38, 39, 40] and the non-perturbative renormalisation group [74, 75]. We are particularly motivated by the perspective that such a theoretical framework can go beyond the limitations of the conventional pressureless perfect fluid approximation when it comes to shell-crossing at late times and small scales. This could extend the range of applicability of semi-analytical techniques for cosmological structure formation substantially and thereby lead to an improved understanding of dark matter and late time cosmology.

Appendix A Non-linear terms of the equations of motion

The terms quadratic in the fluctuation fields in the continuity equation (2.47) are abbreviated by

[TABLE]

Similarly, for the velocity-divergence and vorticity equations (2.48) and (2.49) we abbreviate

[TABLE]

and in the velocity dispersion equations (2.50) – (2.53) we use

[TABLE]

Finally, the terms cubic in the fluctuation fields are abbreviated as

[TABLE]

All fields are evaluated at time $\tau$ and in order to simplify notation we denote a field with wave vector argument $\mathbf{k}_{1}$, $\mathbf{k}_{2}$ or $\mathbf{k}_{3}$ by a subscripted $1$, $2$ or $3$, respectively. The vertices $A,B,C,D,E,F$ and $G$ are the convolution kernels of the Fourier transformation of the non-linear terms of the equations of motion of the corresponding fluctuation fields. They carry an overall momentum conserving Dirac delta function and depend on the wave vectors which are integrated over. We provide the explicit form of these vertices in a future publication.

Acknowledgments

This work is part of the DFG Collaborative Research Centre “SFB 1225 (ISOQUANT)” as well as the DFG project BE 2795/4-1 and EXC-2181/1 – 390900948 (the Heidelberg STRUCTURES Excellence Cluster).