General relativistic non-ideal fluid equations for dark matter from a truncated cumulant expansion

A new truncation scheme based on the cumulant expansion of the one-particle phase-space distribution function for dark matter particles is developed. Extending the method of moments in relativistic kinetic theory, we derive evolution equations which supplement the covariant conservation of the energy-momentum tensor and particle number current. Truncating the cumulant expansion we obtain a closed, covariant and hyperbolic system of equations which can be used to model the evolution of a general relativistic non-ideal fluid. As a working example we consider a Friedmann-Lema\^itre-Robertson-Walker cosmology with dynamic pressure and solve for the time evolution of the effective equation of state parameter.


I. INTRODUCTION
Observations indicate that a substantial part of the energy budget of the Universe is made up of dark matter. To a good approximation, it seems to be non-interacting except through gravitational interactions [1,2]. While there are various ideas as to the nature of dark matter, there is no verified microscopic theory so far. In the absence of a direct detection, its properties can only be constrained through cosmological and astrophysical observations. These are necessarily observations in the macroscopic domain, while one would ultimately be interested in a microscopic understanding in terms of a quantum field theory.
In order to learn more about the properties of dark matter, one therefore turns to effective descriptions, such as fluid approximations or variants of kinetic theory. While isotropy and homogeneity constrain the cosmological background dynamics to be of perfect fluid type [3], gravitational collapse is caused by fluctuations which are expected to generate non-vanishing shear stress at smaller scales and therefore require a non-ideal fluid description. For a consistent coupling to gravity through the Einstein field equations and in order to have a large scale description that obeys relativistic causality, one should aim for a relativistic fluid description. This is also important to understand the possible interplay with modified theories of gravity.
Another striking example for the applicability of relativistic fluid dynamics are heavy-ion collisions, such as those studied at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC). At macroscopic scales these exhibit collective fluid-like behaviour, such that the bulk dynamics can be described in terms of relativistic non-ideal fluid dynamics. This should be understood as a low-energy effective description of quantum chromodynamics in a out-of-equilibrium situation of * erschfeld@thphys.uni-heidelberg.de † floerchinger@thphys.uni-heidelberg.de ‡ rupprecht@thphys.uni-heidelberg.de high energy density [4][5][6][7]. First attempts to construct a theory of relativistic fluid dynamics were made by Eckart [8] and later by Landau and Lifschitz [9]. However, these are known to be unstable and acausal [10][11][12] and can even have modes which propagate faster than light [13]. In general, many socalled first-order theories suffer from the fact that the equations of motion are parabolic in structure making them acausal [14]. Only recently, a class of first-order theories has been proposed that does not suffer from these problems [15][16][17][18]. While the energy-momentum tensor is expanded to first order in gradients around the ideal fluid form, this is done in such a way that the resulting equations of motion are second-order hyperbolic differential equations and can therefore obey relativistic causality.
As another possibility to remedy these issues, extended theories have been proposed. In general, hydrodynamical descriptions are possible due to conservation laws. In the relativistic context, this is the covariant conservation of the energy-momentum tensor, sometimes supplemented by a covariantly conserved particle number current. In such a framework one has ten degrees of freedom in the components of the energy-momentum tensor and four in the particle number current, but only five equations of motion from energy-momentum and particle number conservation. Assuming that all the degrees of freedom are independent and dynamical, the system of equations is not closed and additional information is needed.
The standard theory of hydrodynamics [9] provides such additional relations from an expansion of the conserved currents in terms of gradients of the fields that govern an ideal fluid, namely the fluid velocity and thermodynamic variables such an temperature and chemical potential. At a given finite order, this leads to additional constraint equations that related the components of the conserved currents to the hydrodynamic fields. Such a gradient expansion is particularly well motivated in the vicinity of thermal equilibrium or when interaction effects are so strong that they quickly drive the system back towards local equilibrium when the latter is violated as a consequence of fluid motion.
In the context of dark matter, this close-to-equilibrium gradient expansion is not particularly well motivated, especially at late times and at small scales. Indeed, when interactions besides gravity are absent, one cannot assume that local thermal equilibrium is restored quickly. However, additional information to close the system of equations can also be provided in the form of additional evolution equations. This leads to a formalism with more dynamical variables or fluid fields, for example shear stress, bulk viscous pressure, heat current or particle diffusion current.
Historically, such additional evolution equations have been first obtained from a kinetic theory approach where the system is characterised by a phase-space distribution function which obeys the relativistic Boltzmann equation. With the method of moments [19] one can derive additional equations of motion in order to evolve the dynamical degrees of freedom and close the system by a systematic truncation, such as a derivative expansion. This has been done, most notably in the non-relativistic case by Müller [20] and in the relativistic case by Israel and Stewart [21,22], keeping terms up to second order in gradients, so-called second-order theories. More modern approaches include derivative expansions, including all terms of second order in gradients [24] and formulations derived from the relativistic Boltzmann equation, keeping all terms in the moment expansion [23]. Recent applications of such approaches include references [25][26][27][28][29][30][31].
In the work presented here, we essentially adapt the method of moments and study a truncation of the cumulant expansion of the phase-space distribution function for dark matter particles without collisions. In particular, we concentrate on a description of dark matter as classical particles. A more direct relation of the fluid picture to a quantum field theoretic description is left for future work (see also references [46,47]).
In section II we present the relativistic Vlasov equation, which governs the evolution of the one-particle phase-space distribution function. We define moments and cumulants and study their behaviour under the Vlasov equation. Here, the first and second moment are the particle number current and energy-momentum tensor, respectively. We propose a truncation scheme in terms of the cumulant expansion of the distribution function and explicitly perform this truncation after the first and second cumulant. For these one can reconstruct a modified version of the phase-space distribution function. We discuss the non-consistency of these truncations, in the sense that they are not preserved by the Vlasov equa-tion.
In section III we derive closed equations of motion for the dynamical fluid fields parametrising the 14 independent degrees of freedom. To this end we use the covariant conservation of the first three moments and close the system of equations by neglecting the third cumulant. We discuss the hyperbolic structure of the system of equations, allowing for a causal description provided the characteristic propagation speeds are finite. By this we construct a closed, covariant and hyperbolic system of partial differential equations which can be used to model a general relativistic non-ideal fluid. Finally, we present a simple application of this truncation in the form of a Friedmann-Lemaître-Robertson-Walker (FLRW) cosmology.
Throughout this paper, except when otherwise stated, we work in natural units where c = = k B = 1.

A. The relativistic Vlasov equation
Kinetic theory for classical point particles can be formulated in general coordinates x and in curved spacetime. 1 In the absence of any force except for gravity, particles follow geodesics, where s is an affine parameter and Γ µ ρσ (x) are the Christoffel symbols of second kind of the space-time metric tensor g µν (x). 2 We denote the components of fourvectors by Greek letters and summation over the same co-and contravariant indices is understood. The state of the theory is described by a one-particle phase-space distribution function f (x, p). Here, the four-momentum for particles of mass m is p µ = m dx µ /ds.
The number of particles at position x, in the momentum range √ g d 4 p /(2π) 4 and in a mass range dm flowing through the hypersurface element is given by Here g(x) = − det(g µν (x)) is the determinant of the metric tensor, µνρσ is the total antisymmetric Levi-Civita symbol and we abbreviate the space-time inner product as p 2 = g µν p µ p ν . Further δ and θ denote the (onedimensional) Dirac delta and Heaviside unit step function, respectively. For fixed mass m, the momenta fulfil the on-shell constraint p 2 + m 2 = 0, but we allow more generally a distribution of masses ξ(m) which we take to be independent of time and space. As an example for a dark matter model with some distribution of masses one may think about primordial black holes. For a single species of particles with unique mass m * one has ξ(m) = δ(m − m * ).
The particle number current is given by and the energy-momentum tensor is Here and throughout we abbreviate the mass and momentum integrals as In the absence of scattering, the distribution function is conserved along a geodesic, df /ds = 0, which implies the Vlasov equation [48,49], Note however that this equation is somewhat more general than usual because f (x, p) is here not restricted to p 2 + m 2 = 0 being satisfied for a unique value of m. We allow a more general distribution of masses such that also the distribution function f (x, p) is more general, but equation (7) still determines its time evolution.
In the following it is convenient to introduce the modified distribution functioñ such that for example T µν = p p µ p νf (x, p). Equation (7) immediately implies that the modified distribution function also obeys the Vlasov equation,

The method of moments and cumulants
Often one is not interested in the full distribution function f (x, p) but rather in its moments and cumulants. In analogy to probability distributions one can define fully symmetric four-momentum moments, the set of which completely characterise the distribution function. These can be conveniently derived from a momentgenerating function, by n-fold differentiation with respect to the source four-vector l µ . The normalisation z(x) = z(x; 0) is the zeroth moment and it is evident from equation (4) and (5) that the particle number current N µ (x) and the energy-momentum tensor T µν (x) are the first and second moment, respectively. Due to the on-shell Dirac delta function in equation (10) one obtains the relation which in the limit of a single particle species with unique mass m reads and relates moments which differ by two orders. For a spectrum of masses it is in general not possible to express the relation (12) uniquely in terms of finite moments.
Similarly the distribution function is completely characterised by the connected parts of the moments, the so-called cumulants C µ1...µn (x). These can be derived from the cumulant-generating function ln(z(x; l)) in the same manner as moments are derived from z(x; l).
The first few cumulants and moments are related by the expressions and vice versa Here we denote the symmetrisation with respect to a set of indices by a pair of parentheses around them. Similar relations hold for higher order moments and cumulants and can be straight forwardly derived from the corresponding generating functions. Combining the relations (14) and (15) one can express the n-th moment in terms of the lower order moments and the n-th cumulant, The Vlasov equation (7) implies the covariant conservation of all moments, where ∇ µ denotes the covariant derivative with respect to the coordinates x µ . It is immediately clear that the evolution equations are independent of each other at each order. 3 In contrast, the cumulants follow the non-linear evolution equation where the sums runs over all combinations of picking indices {α 1 , ..., α n } out of {µ 1 , ..., µ n }. 4 The non-linear terms in the evolution equation of the n-th cumulant involve all lower order cumulants and the sum of orders of the two cumulants in the quadratic terms always adds to n. This structure sets strong restrictions for a consistent truncation of the cumulant expansion as we discuss at the end of the next section. We remark that even though the evolution equation (18) involves all cumulants of lower order it does not couple to higher orders. This is to be seen in contrast to the non-relativistic limit, in which the evolution equations of the moments and cumulants depend on the next higher order moment or cumulant, creating the so-called Vlasov hierarchy [50][51][52].
This can be made explicit by taking the non-relativistic limit of the moments (10), which for simplicity we do for a single particle species of unique mass m. To do so we restore the speed of light c and evaluate the p 0 integral to arrive at where and Latin indices only range over spatial components.
In the limit c → ∞, where the metric is Minkowskian, the n-th moment involving k temporal indices and n − k spatial indices is Here M i1...i n−k nr (t, x) is the (n − k)-th non-relativistic moment and t denotes time and x position in configuration space. From the non-relativistic limit of evolution equation (17) it is then evident, that the evolution of the n-th (non-relativistic) moment also depends on the next higher (non-relativistic) moment. A similar argument holds true for the evolution of the non-relativistic cumulants.

C. Truncated cumulant expansion
In general, the distribution function f (x, p) has an infinite amount of moments, as well as cumulants, the complete set of which encode the same information. Therefore a characterisation in terms of moments or cumulants are two sides of the same coin and are a matter of taste or problem at hand. However, since cumulants of different orders are statistically independent of each other it seems natural to study the cumulant expansion of the distribution function in order to put forth a truncation scheme which allows to describe a finite set of independent degrees of freedom. Assuming one could truncate the expansion after the n-th cumulant one would be left with independent degrees of freedom. Here, the first sum represents the degrees of freedom from the (fully symmetric) cumulants up to order n. However, these are not all independent of each other due to the on-shell constraints (12) which reduce the independent degrees of freedom. This is represented by the second sum. From the full set of cumulants one can in principle reconstruct the distribution functionf (x, p). However, this is not so for a finite truncation, with the exceptions of a truncation after the first or second cumulant, corresponding to a degenerate or normal distribution, respectively.
To work out the consequences of such truncations in more detail it is useful to decompose the particle number current as and the energy-momentum tensor as Here n is the particle number density, u µ is the local fluid four-velocity normalised to u µ u µ = −1 and ν µ is the diffusion current orthogonal to the fluid velocity, u µ ν µ = 0. Further, is the energy density in the local rest frame, p the thermodynamic pressure which is related to n and by the equilibrium expression, π bulk is the bulk viscous pressure, ∆ µν = u µ u ν + g µν is a projector orthogonal to the fluid velocity, π µν is the shear stress tensor, which is symmetric, transverse to the fluid velocity u µ π µν = 0 and traceless π µ µ = 0 and q µ is the heat current which is orthogonal to the fluid velocity, u µ q µ = 0. In the following we abbreviate the sum of thermodynamic and bulk viscous pressure as an effective pressure, p eff = p + π bulk . These rather general relations can be specialised to a frame by fixing the fluid four-velocity. Common choices are the Landau frame [9], where the fluid four-velocity is a time-like eigenvector of the energy-momentum tensor and the heat current vanishes, or the Eckart frame [8], in which the fluid four-velocity is defined by the direction of the particle number current so that one has vanishing diffusion current.
In the most simple non-trivial case one truncates the cumulant expansion after the first order and obtains the so-called single-stream approximation. It is characterised by four independent degrees of freedom and for a single particle species with unique mass m the on-shell constraint (13) leads to z = n/m. The modified distribution function is a degenerate distribution, where δ (4) denotes the four-dimensional Dirac delta function. In this case one has = nm and p eff = ν µ = q µ = π µν = 0, or in other words, the particle number current and the energy-momentum tensor are the ones of an perfect pressureless fluid.
Including the second cumulant and truncating the expansion at third order one obtains a normal distribution characterised by ten independent degrees of freedom. 5 Again specialising to a single particle species with unique mass m, the second moment on-shell constraint g µν T µν = m 2 z leads to The third moment on-shell constraint projected along the fluid velocity u ρ g µν M µνρ = m 2 u ρ N ρ gives while the projection orthogonal to the fluid velocity The modified distribution function is given by, where u · p = u µ p µ and due to the constraint (28) the second cumulant is purely transverse to the fluid velocity, and therefore the parallel part collapses to a Dirac delta function. The constraints (26) and (27) can be combined to give reducing the independent degrees of freedom to ten. Interestingly enough this truncation implies vanishing diffusion and heat current and thus the Landau and Eckart frame are the same. Formally we have now a fluid with non-vanishing bulk viscous pressure and shear stress. This truncation scheme can be straight forwardly generalised to higher orders. However, there are no distributions which are characterised by a finite number of cumulants beside the degenerate and normal distributions and one is therefore not able to explicitly reconstruct f (x, p). Further one cannot find a corresponding distribution function f (x, p) because one is not able reconstruct the δ(p 2 + m 2 ). Nevertheless, as an approximate model for a more complex form, equation (29) may be quite reasonable. A nice feature is that the cumulants that govern it have rather transparent equations of motion that can be solved in the spirit of a non-ideal fluid approximation as we discuss in the next section. While we are in principle able to truncate the cumulant expansion of the distribution function, the questions arises whether these truncations are preserved by the Vlasov equation (18). Interestingly enough the single-stream approximation (25) is (apparently) preserved since no higher order cumulants are sourced by terms solely depending on the zeroth or first cumulant. However, this apparent self-consistency is not stable under perturbations, that is as soon as any other cumulant obtains a non-vanishing value, cumulants of all orders are generated. This can be checked from the combinations of cumulants appearing in the quadratic terms of equation (18). This also implies that any truncation beyond the single-stream approximation is not preserved by the Vlasov equation and higher order cumulants are generated throughout evolution in time. The self-consistency of the single-stream approximation is only apparent due to the phenomenon of shell-crossing, when multiple streams of matter coexist at the same region in space. At this point in configuration space the velocity field in equation (25) is multivalued and the corresponding distribution functionf (x, p) has non-vanishing second and higher order cumulants [50].

A. Evolution equations
The particle number density and energy-momentum tensor have 14 independent degrees of freedom which we parametrise in terms of the fields introduced in equation (23) and (24). Evolution equations for these can be obtained from the covariant conservation of the first three moments, by expressing the third moment in terms of lower order moments and the third cumulant, as done in equation (16). The conservation of the particle number current yields an evolution equation for the particle number density, Projecting the conservation of the energy-momentum tensor along the fluid four-velocity, u µ ∇ ν T µν = 0, gives an evolution equation for the energy density, while a projection orthogonal, ∆ µ ν ∇ ρ T νρ = 0, yields essentially an evolution equation of the fluid four-velocity, A projection onto the orthogonal parts of the covariant conservation of the third moment, ∆ µν ∇ ρ M µνρ = 0, gives an evolution equation for the effective pressure, The evolution equation for the three degrees of freedom parametrised by the diffusion current ν µ and heat current q µ is obtained from the projection ∆ α µ u ν ∇ ρ M µνρ = 0, Finally, the evolution for the shear stress tensor π µν follows from the projection P αβ µν ∇ ρ M µνρ = 0, where the symmetric, traceless and transverse projector reads It is evident that the system of equations (33) - (38) is not closed, since they couple to the normalisation z and the third cumulant C µνρ . In order to close the system of equations we impose a truncation of the cumulant expansion after the second cumulant, C µνρ = 0, and use the energy-momentum tensor on-shell constraint (26) to eliminate z. Strictly speaking not all of the equations are independent of each other due to equations (27) and (28). These reduce the independent degrees of freedom to ten in accordance with formula (22). Nevertheless, we assume that these constraints can be neglected with the reasoning as follows. We know from the discussion at the end of section II that the imposed truncation is not preserved by the Vlasov equation and the third cumulant is generated throughout the evolution in time. Therefore the constraints (27) and (28) are not strictly satisfied but rather proportional to the third cumulant as indicated on the right-hand side. That is, we do not assume the third cumulant to be exactly zero throughout time evolution, but rather to stay small enough to be neglected in the evolution equations. A similar truncation of the cumulant expansion after the second order has also been used in the non-relativistic limit to model dark matter with non-vanishing velocity dispersion [52,53].
Collecting the 14 degrees of freedom in a superfield Φ a (x) the evolution equations (33) - (38) can be written in the quasi-linear form Here A a and B ρ a are field dependent matrices and u ρ B ρ a = 0. The index a carries the appropriate amount of covariant indices, e.g. in the Eckart frame one has Φ a = (n, , p eff , u µ , q µ , π µν ). Since the expressions for the matrices A a and B ρ a are less transparent and rather cumbersome if written explicit, we restrain from displaying them. We checked that the matrix A a is diagonalisable and invertible, making the system hyperbolic [54]. Finally, we are left with a closed, covariant and hyperbolic system of equations which describe the 14 degrees of freedom introduced in equations (23) and (24).

B. FLRW cosmology
As a simple cosmological working example we consider a spatially flat Friedmann-Lemaître-Robertson-Walker metric with line element where a(t) is the dimensionless scale factor which parametrises the relative spatial expansion of the Universe. Due to the symmetries of the metric, namely homogeneity and isotropy, the particle number current and energy-momentum tensor have the form of a perfect fluid, Here the energy density (t) and effective pressure p eff (t) are functions of time only and the four-velocity is given by u µ = (1, 0, 0, 0). We define the effective equation of state parameter as the ratio of the effective pressure and energy density, ω eff = p eff / . The evolution equations (33) and (34) read while equation (36) gives and H = ∂ t a/a is the Hubble rate. These equations can be solved in terms of the scale factor, n = n 0 a −3 , where quantities subscripted with a 0 are the values at a(t) = 1, corresponding to today. The particle number density decays as expected with the expansion of space while the decay of the energy density also depends on the effective equation of state parameter. state parameter ω eff < 1/3 the solution evolves towards the attractive non-relativistic solution fixed point as is expected for matter that is non-interacting except for gravity. For an effective equation of state parameter ω eff > 1/3 the solution exhibits a strong growth, but the physical interpretation behind such values of the equation of state parameter is not clear. In a next step it would be interesting to treat nonhomogeneous solutions and cosmological structure formation. The full set of non-ideal fluid fields, including the peculiar fluid velocity, shear stress and diffusion or heat current, are then expected to become non-trivial. Interestingly, perturbations in these fields can also influence the overall cosmological expansion through a dissipative back-reaction [55].

IV. CONCLUSIONS
Starting from a general relativistic kinetic theory approach for a system of collisionless classical point particles we presented the method of moments and cumulants of the one-particle phase-space distribution function. We studied how the moments and cumulants evolve under the relativistic Vlasov equation and found in particular that all moments are covariantly conserved while the cumulants follow a more complex non-linear evolution equation. For the first and second moment these are the common covariant conservation laws related to the particle number current and energy-momentum tensor, but the evolution of higher order moments furnish additional evolution equations.
We showed how the cumulant expansion of the distribution function can be truncated at finite order and explicitly performed this truncation after the first and second cumulant, corresponding to the single-stream and a Gaussian approximation, respectively. In particular the Gaussian approximation is capable of describing a nonideal fluid with non-vanishing bulk viscous pressure and shear stress. We discussed that these kind of truncations are not preserved by the Vlasov equation since higher order cumulants are naturally generated by lower order cumulants.
From the covariant conservation of the first three moments we derived a closed, covariant and hyperbolic system of equations by neglecting the third cumulant. The equations give the evolution of the 14 degrees of freedom of the particle number current and energy-momentum tensor and can be used to describe a general relativistic non-ideal fluid. As a working example we considered a Friedmann-Lemaître-Robertson-Walker cosmology with non-vanishing dynamic pressure and solve its time evolution. We find that the solution has an attractive and repulsive fixed point, corresponding to non-relativistic matter and ultra-relativistic radiation, respectively.
The equations of motion for a non-ideal fluid approximation to dark matter that we have derived can be extended in different directions. One would be to include the effects of dark matter self interaction or interactions with baryonic matter. Another would be to include quantum effects. Moreover, it would of course be highly interesting to solve the evolution equations for the non-ideal fluid fields or cumulants directly, either numerically or with further analytical techniques such as perturbation theory or field theoretic methods similar to those developed in references [56][57][58]. In particular we are curious to see whether the truncated cumulant expansion developed here, or an extension of it, can agree with numerical solutions of the Vlasov equation through N -body simulations. This might lead to a rather useful framework to study extensions of the collisionless cold dark matter model and for a comparison to observational data.