Causal Relativistic Hydrodynamics of Conformal Fermi-Dirac Gases

TL;DR

This paper derives a causal relativistic hydrodynamics framework for Fermi-Dirac gases from kinetic theory, addressing divergence issues and analyzing stability of anisotropic configurations.

Contribution

It introduces a divergence type theory approach for Fermi-Dirac gases, providing a stable hydrodynamic description contrasting with Grad approximation predictions.

Findings

The derived hydrodynamics are stable for anisotropic configurations.

Divergence type theories regularize formal divergences in kinetic expressions.

Hydrodynamics from DTTs differ from Grad approximation in stability predictions.

Abstract

In this paper we address the derivation of causal relativistic hydrodynamics, formulated within the framework of Divergence Type Theories (DTTs), from kinetic theory for spinless particles obeying Fermi-Dirac statistics. The approach leads to expressions for the particle current and energy momentum tensor that are formally divergent, but may be given meaning through a process of regularization and renormalization. We demonstrate the procedure through an analysis of the stability of an homogeneous anisotropic configuration. In the DTT framework, as in kinetic theory, these configurations are stable. By contrast, hydrodynamics as derived from the Grad approximation would predict that highly anisotropic configurations are unstable.

Equations

$$T^{\mu \nu }=\int Dp{p}^{\mu }{p}^{\nu }f$$

$$Dp=\frac{2d{p}_0{d}^{3}p}{{\left(2\pi \right)}^3}\delta \left(p^2\right)\theta \left(p^0\right)$$

$$T^{\mu \nu }{u}_{\nu }=-\rho u^{\mu }$$

$$S^{\mu }=-\int Dp{p}^{\mu }\left\{\left(1-f\right)\ln \left(1-f\right)+f\ln f\right\}$$

$$f_{eq}=\frac{1}{e^{-{\beta }_{\mu }{p}^{\mu }}+1}$$

$$T_{eq}^{\mu \nu }$$

$$S_{eq}^{\mu }$$

$$p^{\nu }\frac{\partial }{\partial x^{\nu }}f=I_{col}\left[f\right]$$

$$\int Dp{p}^{\mu }{I}_{col}\left[f\right]=0$$

$$T_{,\nu }^{\mu \nu }=0$$

$$\int Dp\ln \left[\frac{1}{f}-1\right]I_{col}\left[f\right]\ge 0$$

$$S_{,\mu }^{\mu }\equiv \sigma \ge 0.$$

$$I_{col}\left[f\right]=\frac{1}{\tau }{u}_{\mu }{p}^{\mu }\left(f-f_{eq}\right).$$

$$T^{\mu \nu }={\sigma }_{SB}{T}^4\left[u^{\mu }{u}^{\nu }+\frac{1}{3}{\Delta }^{\mu \nu }\right]+{\Pi }^{\mu \nu }$$

$$f_{Ch-E}=f_{eq}\left[1+\tau \left(1-f_{eq}\right)\delta f_{Ch-E}\right]$$

$$\delta f_{Ch-E}=\frac{-1}{2T\left|u_{\rho }{p}^{\rho }\right|}{\sigma }_{\mu \nu }{p}^{\mu }{p}^{\nu }$$

$${\sigma }_{\mu \nu }={\Delta }_{\mu }^{\rho }{\Delta }_{\nu }^{\lambda }\left[u_{\rho ,\lambda }+u_{\lambda ,\rho }-\frac{2}{3}{\Delta }_{\rho \lambda }{u}_{,\tau }^{\tau }\right]$$

$${\Pi }^{\mu \nu }=-\eta {\sigma }^{\mu \nu }$$

$$f=f_{eq}\left[1+Z\right]$$

$$Z=\frac{1}{2T\left(-u_{\mu }{p}^{\mu }\right)}\left(1-f_{eq}\right){\xi }_{\mu \nu }{p}^{\mu }{p}^{\nu }$$

$${\Pi }_{\mu \nu }=\left(\frac{7{\pi }^2}{900}\right)T^4{\xi }_{\mu \nu }$$

$$\left[{\Delta }_{\mu }^{\rho }{\Delta }_{\nu }^{\lambda }-\frac{1}{3}{\Delta }^{\rho \lambda }{\Delta }_{\mu \nu }\right]\left\{A_{,\sigma }^{\sigma \mu \nu }-I^{\mu \nu }\right\}=0$$

$$A^{\sigma \mu \nu }$$

$$I^{\mu \nu }$$

$$\{\begin{array}{l}{\partial }_{\mu }{N}^{\mu }=0\\ {\partial }_{\mu }{T}^{\mu \nu }=0\end{array}$$

$${\partial }_{\mu }{A}^{\mu \nu \rho }=I^{\nu \rho },$$

$$N^{\mu }=\frac{\partial {\chi }^{\mu }}{\partial \alpha },T^{\mu \nu }=\frac{\partial {\chi }^{\mu }}{\partial {\beta }_{\nu }}\text{and}{A}^{\mu \nu \rho }=\frac{\partial {\chi }^{\mu }}{\partial {\zeta }_{\nu \rho }}.$$

$${\chi }_{eq}^{\mu }=-\int Dp{p}^{\mu }\text{ln}\left(1-f_{eq}\right),$$

$${\chi }^{\mu }=-\int Dp{p}^{\mu }\text{ln}\left(1-f\right)$$

$$f\left[p^{\mu },u^{\mu },{\zeta }^{\mu \nu },T\right]=\frac{1}{e^{-\frac{1}{T}{u}_{\mu }{p}^{\mu }-{\zeta }_{\mu \nu }{p}^{\mu }{p}^{\nu }}+1}$$

$$f_{eq}\left[p^{\mu },u^{\mu },T_{eq}\right]=\frac{1}{e^{-\frac{1}{T_{eq}}{u}_{\mu }{p}^{\mu }}+1}.$$

Full text

I Introduction

The successful application of relativistic hydrodynamics Eck40 ; LL6 ; ISR76 ; IsraelStewart ; ISR88 ; ReZa13 ; Sch14 to the description of high energy heavy ion collisions Rom10 ; Flor10 ; CalzettaVdL ; Monnai14 has led not only to a revival of this theory, but also to the demand of enlarging its domain of applicability to regimes where the system of interest is still away from local thermal equilibrium Str14b ; JeHe15 ; Rom16 , and so the usual strategy of deriving hydrodynamics as an expansion in deviations from ideal behavior is not available. Moreover, the best known implementations of this strategy, namely the Chapman-Enskog Chapman ; Enskog and Grad Grad1 ; Grad2 approximations, face severe problems, such as spurious instabilities, as it will be shown below.

The so-called divergence type theories (DTTs)liu72 ; LiuMullerRuggeri ; GerochLindblom1 ; GerochLindblom2 are an appealing alternative for the derivation of relativistic causal hydrodynamics because in this framework both the conservation laws for particle number and the energy-momentum tensor as well as the Second Law of Thermodynamics are rigorous properties of the theory, no matter how far from ideal behavior. For this reason the solutions of the theory may be trusted to be at least qualitatively faithful to the underlying kinetic theory. By contrast, a formalism that only enforces the Second Law in an approximate way could lead to unphysical results if the system makes a large excursion away from local thermal equilibrium, even if it is a transient one, and then the whole further evolution would be compromised. Furthermore, in more complete theories including gauge fieldsMaMr06 ; PRD12a ; AIP14 , these spurious instabilities could mask or get entangled with legitimate plasma instabilities Sch06 ; mama07 ; PRC13a ; MSS16 ; KC16 .

In this paper we shall analyze the derivation of DTT relativistic hydrodynamics from kinetic theory taking as test case a gas of spinless and massless particles obeying Fermi-Dirac statistics LewisRomatschke . In the Grad approach, this derivation consists on formulating an ansatz for the one particle distribution function (1pdf), parameterized by the hydrodynamic variables. Later on the hydrodynamic currents such as the particle number current and the energy-momentum tensor are derived as moments of the parameterized 1pdf, and the corresponding equations as moments of the Boltzmann equation DKR10 ; DMNR12a ; DMNR12b ; Cristian . This procedure may be replicated in the DTT framework, but it leads to formally divergent expressions. Therefore it is necessary to interpolate a process of regularization and renormalization by which these expressions become meaningful. The conclusion is that a DTT can be derived from kinetic theory, but not uniquely.

As a demonstration of the formalism we shall carry on the procedure, adopting a regularization and renormalization scheme that does not introduce new dimensionful parameters in the theory, preserves positive expressions, and gives the right results in equilibrium, where all relevant expressions are finite to begin with. We shall use the resulting DTT to investigate the stability of an anisotropic (though axisymmetric) homogeneous configuration (a precise characterization will be given below). These are always stable in kinetic theory, but we will show that the Grad approximation predicts an instability if the anisotropy exceeds a certain threshold. DTT agrees with kinetic theory predicting again stability. As it ought to be expected, the quantitative agreement worsens for larger deviations form equilibrium.

Let us be more specific about the contents of this paper. We consider a gas of massless, spinless particles obeying Fermi-Dirac statistics. In relativistic kinetic theory S71 ; I72 ; IsrSt79 ; dGvLW ; CalHu08 , the state of the gas is described by a 1pdf $f=f\left(x^{\mu},p_{\nu}\right)$, where the momentum variable is restricted to the positive mass shell $p^{2}=0=\vec{p}^{2}-p^{02}$, $p^{0}\geq 0$ (we adopt the $\left(-+++\right)$ signature for Minkowsky metric $\eta_{\mu\nu}$). From $f$ we derive the energy-momentum tensor

[TABLE]

where $Dp$ is the invariant measure

[TABLE]

Observe that $T^{\mu\nu}$ is traceless. $T^{\mu\nu}$ admits one (and only one) timelike eigenvector

[TABLE]

$u^{2}=-1$. We say $u^{\mu}$ is the (Landau-Lifshitz) fluid velocity LL6 , and $\rho$ the energy density. The other relevant current is the entropy flux

[TABLE]

and $s=-u_{\mu}S^{\mu}$ is the entropy density. In equilibrium, the 1pdf must maximize the entropy density for a given energy density. This obtains when $f$ is the Fermi-Dirac distribution

[TABLE]

where $\beta_{\mu}=u_{\mu}/T$, $T$ being the temperature. Thus in equilibrium

[TABLE]

where $\Delta^{\mu\nu}=\eta^{\mu\nu}+u^{\mu}u^{\nu}$ and $\sigma_{SB}=7\pi^{2}/240$ is Stefan-Boltzmann’s constant.

Out of equilibrium $f$ evolves according to the Boltzmann equation S71 ; I72 ; IsrSt79 ; dGvLW ; CalHu08

[TABLE]

The collision integral $I_{col}$ vanishes in equilibrium, and obeys

[TABLE]

which enforces energy-momentum conservation

[TABLE]

and the $H$ theorem

[TABLE]

which enforces the Second Law

[TABLE]

For concreteness we shall adopt the Anderson-Witting collision termAndersonWitting1 ; AndersonWitting2 ; TI10

[TABLE]

where $f_{eq}$ is the equilibrium distribution with the same velocity and energy density as the nonequilibrium 1pdf $f$. Another frequently used prescription is the Marle or BGK one Marle1 ; Marle2 , where $-u_{\mu}p^{\mu}$ in the right hand side is replaced by a power of temperature.

For a general 1pdf, $T^{\mu\nu}$ acquires a new term, the viscous energy-momentum tensor $\Pi^{\mu\nu}$,

[TABLE]

Since $\Pi^{\mu\nu}$ is traceless and transverse $\Pi^{\mu}_{\mu}=u_{\nu}\Pi^{\mu\nu}=0$ it has $5$ independent components, elevating the total number of degrees of freedom in $T^{\mu\nu}$ to $9$. The four conservation equations Eq. (9) are therefore not enough to predict the evolution of the energy-momentum tensor. The problem of relativistic hydrodynamics is to provide the missing equations.

The Chapman-Enskog approachChapman ; Enskog assumes that at every point $f$ is close to an equilibrium distribution, although with position dependent temperature and velocity. Then a solution of Eq. (7) is sought as a formal expansion in powers of the relaxation time $\tau$ introduced in Eq. (12)

[TABLE]

Inserting this into Eq. (7) with collision term (12), and using the conservation equations Eq. (9) to order $\tau^{0}$ to simplify the result, we obtain to lowest order

[TABLE]

where we introduced the shear tensor

[TABLE]

A straightforward computation yields

[TABLE]

where $\eta=\left(7\pi^{2}/900\right)\tau T^{4}$ is the shear viscosity. Thus in this approach the viscous energy-momentum tensor is slaved to the degrees of freedom that describe the ideal fluid at $\tau=0$. This eventually leads to a parabolic system of equations of motion, incompatible with relativistic causality HisLin83 ; HiscockLindblom ; HisLin88 .

To overcome this difficulty, the Grad approachGrad1 ; Grad2 proposes instead a 1pdf

[TABLE]

[TABLE]

$\xi_{\mu\nu}$ is traceless and transverse, and it is regarded as an independent tensorial degree of freedom. It is directly related to $\Pi_{\mu\nu}$, since

[TABLE]

Since $\xi_{\mu\nu}$ is not positive definite, the Grad approximation will lead to negative pressures if $\xi_{\mu\nu}$ is large enough, which underscores the unapplicability of the theory far from equilibrium. Moreover, we shall show below the theory has spurious instabilities even before that limit is reached.

To obtain a dynamics for these new $5$ degrees of freedom in the viscous energy-momentum tensor, one further moment of the Boltzmann equation is computed. The first moments yield energy-momentum conservation Eq. (9). Instead of the ten second moments, we only consider the traceless, transverse ones

[TABLE]

where

[TABLE]

The nonequilibrium current $A^{\sigma\mu\nu}$ is totally symmetric and traceless on any two indexes. This approach leads to a Maxwell-Cattaneo HeatWaves equation for $\xi^{\mu\nu}$ and enforces causality. However, it cannot be applied arbitrarily far from equilibrium, because, as we shall show below, it predicts instabilities that do not exist in the kinetic theory. For further discussion of the Chapman-Enskog and Grad approaches see Jais13

The DTT framework keeps Eqs. (9) and (21) as the fundamental equations, but now seeks a 1pdf which maximizes entropy density for given energy density and $A^{0ij}$ components in the rest frame. This leads to the introduction of a new tensor Lagrange multiplier $\zeta_{\mu\nu}$ besides $T$ and $u_{\mu}$. Assuming $\zeta_{\mu\nu}$ is symmetric, traceless and transverse, it is equivalent to $5$ new degrees of freedom. Thus the theory has the same number of degrees of freedom as the energy-momentum tensor, with the nonequilibrium current slaved to it (we will return to this point below). This means that in the DTT, two evolutions starting with the same energy-momentum tensor will remain identical, though it is known that they may diverge in kinetic theory Ollitrault . Even so, we will show that the DTT outperforms the Grad approximation, in the sense that it is free from the spurious instabilities that appear in the latter.

The problem is that, though the variational problem leading to the DTT 1pdf is easily solved, the formal expressions one obtains for the energy-momentum tensor and the nonequilibrium current diverge ReulaNagy . Thus it is necessary to regularize and renormalize them to make sense of the theory. This adds a new, non unique stage in the derivation of hydrodynamics from kinetic theory. Our goal is to show a concrete procedure to obtain finite quantities for the relevant currents, and then to use this procedure to demonstrate the stability of anisotropic, axisymmetric configurations, in agreement with kinetic theory.

The rest of the paper is organized as follows. In next section we provide some further background on the DTT framework, starting from the purely macroscopic point of view whereby it was first introduced, and then linking it to kinetic theory. Then we proceed to regularize and renormalize the formal expressions for the energy-momentum tensor and nonequilibrium current. Section 3 provides a first comparison of DTT and Grad hydrodynamics through the analysis of the pressure anisotropy; we show that while in Grad hydrodynamics the pressure anisotropy becomes negative when far from equilibrium, in DTT it is bounded below. Section 4 is the main part of this paper; here we discuss the stability of anisotropic homogeneous configurations, comparing the analysis made within DTT and Grad hydrodynamics to the one in kinetic theory. We conclude that while both kinetic theory and DTT predict anisotropic axisymmetric configurations are always stable, Grad hydrodynamics shows an instability if the anisotropy is large enough. We conclude with some brief final remarks.

The two appendices discuss important conceptual issues. Appendix A presents a general framework to analyze stability of conformal hydrodynamic theories of the type discussed in this paper. We show that the instability of the Grad approximation is a consequence of the linearization of the one-particle distribution function with respect to $\xi^{\mu\nu}$. In Appendix B we show how the regularization and renormalization scheme presented here may be applied to fluids obeying Maxwell-Jüttner or Bose-Einstein statistics.

II Divergence Type Theories

DTTs are theories in which all the dynamical equations can be written as divergences of tensor fields. They were originally developed by Liu, Müller and Ruggeriliu72 ; LiuMullerRuggeri as a response to the perceived flaws of the so-called “first order” relativistic hydrodynamics of Eckart Eck40 and Landau-Lifshitz LL6 . They were later extended by the works of Geroch and LindblomGerochLindblom1 ; GerochLindblom2 and Reula and NagyReulaNagy . They were applied to study relativistic hydrodynamic fluctuations in Calzetta_CQG , and free streaming flows in PRD15 . They were applied to the study of relativistic heavy ion collisions in PRC09 ; PRC10b ; IJMPD10 .

The main fields on a DTT are the particle-number current $N^{\mu}$ and the energy-momentum tensor $T^{\mu\nu}$ and their dynamics are governed by the conservation equations

[TABLE]

Closure of the system is achieved by the addition of the balance law of fluxes

[TABLE]

where $A^{\mu\nu\rho}$ and $I^{\nu\rho}$ are algebraic functions of $N^{\mu}$ and $T^{\mu\nu}$. This means we are not adding extra degrees of freedom.

To relate these currents among themselves, it is assumed not only that there is an entropy flux vector $S^{\mu}$ (cfr. Eq. (4)) whose divergence $\sigma$ (cfr. Eq. (11)) is positive, but moreover that both $S^{\mu}$ and $\sigma$ are algebraic functions of $N^{\mu}$ and $T^{\mu\nu}$, such that the positivity of $\sigma$ follows from Eqs. (23) and (24) alone. It can be shown liu72 that this implies the existence of a vector $\chi^{\mu}=\chi^{\mu}\left(\alpha,\beta_{\mu},\zeta_{\mu\nu}\right)$ and a source $I^{\mu\nu}=I^{\mu\nu}\left(\alpha,\beta_{\mu},\zeta_{\mu\nu}\right)$ in such a way that the fields $N^{\mu}$, $T^{\mu\nu}$ and $A^{\mu\nu\rho}$ can be computed as the following partial derivatives

[TABLE]

The variables $\alpha$ and $\beta_{\mu}$ are related to the chemical potential, the hydrodynamic velocity and the temperature and the symmetric tensor $\zeta_{\mu\nu}$ provides the necessary degrees of freedom to match any given energy-momentum tensor.

II.1 Dissipative Type Theories from Kinetic Theory

Since a DTT is totally defined by the generating function $\chi^{\mu}$, to establish a link with kinetic theory it is necessary to relate the generating function to the 1pdf. Let us consider the massless case from now on, so we shall drop $\alpha$ and the particle number current from the discussion.

In equilibrium, the energy-momentum tensor Eq. (6) may be recovered from the generating function

[TABLE]

where $f_{eq}$ is the Fermi-Dirac 1pdf Eq. (5). This suggests to generalize this to the dissipative case by writing

[TABLE]

and a source chosen to match the Anderson-Witting collision term, as in Eq. (22). In order that we may recover the nonequilibrium current as a derivative of $\chi^{\mu}$ we must write $f$ as a deformation of the Fermi-Dirac distribution ReulaNagy

[TABLE]

This distribution function maximizes the entropy density for given values of $T^{00}$ and $A^{0ij}-\left(1/3\right)\delta^{ij}A^{0k}_{k}$, as measured in the fluid rest frame. $f_{eq}$ in the collision term Eq. (12) reads

[TABLE]

with $T$, $\zeta^{\mu\nu}$ and $T_{eq}$ related by

[TABLE]

Although the theory can be formally defined as is and the energy-momentum tensor $T^{\mu\nu}$ and nonequilibrium current $A^{\mu\nu\rho}$ can be expressed as partial derivatives of the generating function $\chi^{\mu}$, there’s an obvious problem in Eq. (27). Since the quadratic form $\zeta_{\mu\nu}p^{\mu}p^{\nu}$ is not negative definite, there are values of $\zeta_{\mu\nu}$ such that

[TABLE]

making the generating function divergent. Since these are all states with occupation number one, this singularity can be interpreted as the Dirac Sea. There’s also a not-so-obvious singularity because of the behavior of the integrals near the manifold where the quadratic part of the argument of the exponential becomes zero. That is, the manifold defined by the equation

[TABLE]

This will be further clarified below.

II.2 Regularization

In order to take care of the singularities Eq. (31) and Eq. (32), let’s fix the dissipative tensor $\zeta_{\mu\nu}$. Motivated by the family of solutions first introduced by Romatschke and Strickland RomatschkeStrickland (see also Florkowski ; Tinti ), we choose the transverse, traceless and axisymmetric case

[TABLE]

($\zeta_{0}>0$) so the type of integrals to regularize are

[TABLE]

and

[TABLE]

where $g$ is a polynomial function in the variables $p$, $\text{cos}\,\theta$ and $\text{sin}\,\theta$ and $f$ is the dissipative Fermi-Dirac distribution Eq. (28),

[TABLE]

Note that the first singularity, Eq. (31), happens when $\text{cos}^{2}\theta<1/3$ and the second one, Eq. (32), when $\text{cos}^{2}\theta\rightarrow 1/3^{-}$ so the $pp_{z}$-plane gets divided into two sections delimited by $\text{cos}^{2}\theta_{0}=1/3$. On one of those sections, when $\text{cos}^{2}\theta>1/3$, the integral is regular and on the other one is where the singularities are located. Figure (1) shows this.

The Dirac Sea singularity, namely Eq. (31), can be eliminated by a simple integration by parts. The surface term, which is infinite, is discarded and the $f\left(p,\theta\right)$ that used to be in the integrand gets replaced by $f\left(p,\theta\right)\,\left[1-f\left(p,\theta\right)\right]$ in the remaining term, which goes to zero as $f$ goes to one. The other singularity, Eq. (32), requieres further analysis. The goal of this section is to introduce a regularization procedure which eliminates both singularities but without introducing new parameters to the theory, as to preserve the conformal invariance.

II.2.1 The $I_{1}$ Case

Let’s start by dividing the $I_{1}$ integral in $\theta_{0}=\text{cos}^{-1}\left(1/\sqrt{3}\right)$,

[TABLE]

$I_{1}^{<}$ is finite and is left as is. To study the $I_{1}^{>}$ term, define a function $G$ by

[TABLE]

so that $G\left(p,\pi/2\right)=0$ and

[TABLE]

Now integrate by parts to obtain

[TABLE]

it is evident that the first term is finite. Let’s call the second term $K_{Z}$. By using the sum of arguments relation of the hyperbolic cosine and realizing it could be written as a partial derivative, it is possible to rewrite $K_{Z}$ as

[TABLE]

Performing another integration by parts and defining the auxiliary function

[TABLE]

we arrive at the final expression for $K_{Z}$,

[TABLE]

The key here is to identify the ratio in Eq. (43) as the mean value

[TABLE]

where $F_{1}$ is a probability density function defined as

[TABLE]

if $u\in\left[0,1\right)$ and $F_{1}\left(u\right)=0$ otherwise. Although $F_{1}(1^{+})=0$, $F_{1}(1^{-})$ goes to infinity logarithmically and, therefore, $U$ is divergent. The first step to improve this behavior is to replace $U$ by its Cauchy principal value,

[TABLE]

which, as a consequence of the Sokhotski-Plemelj-Fox theorem SPFT , equals

[TABLE]

Now let $e^{W\left(t\right)}$ be the characteristic function of $F_{1}$,

[TABLE]

Then $U_{PV}$ can be written as

[TABLE]

Since $e^{W\left(t\right)}$ is the characteristic function of $F_{1}$, then $W$ is the cumulant-generating function. This means that $W$ has the formal power series expansion

[TABLE]

where $\kappa_{n}$ is the $n$th cumulant. That is, $\kappa_{1}=\left\langle u\right\rangle$, $\kappa_{2}=\sigma^{2}=\left\langle u^{2}\right\rangle-\left\langle u\right\rangle^{2}$ and so on.

II.2.2 The $I_{2}$ Case

This case is fairly similar to the first one the only difference being, due to possible divergent terms introduced by $G$, one less integration by parts is performed. We start by dividing the integral the same way as before,

[TABLE]

$I_{2}^{<}$ is finite so we focus our attention on $I_{2}^{>}$. Using the sum of arguments relation of the hyperbolic cosine we write

[TABLE]

Just as before we define the auxiliary function

[TABLE]

and the probability density function

[TABLE]

if $u\in\left[0,1\right)$ and $F_{2}\left(u\right)=0$ otherwise. Now $I_{2}^{>}$ can be written as the following mean value

[TABLE]

$\mathcal{K}$ is obviously finite but

[TABLE]

is not. In order to improve this behavior we replace $U$ by its Cauchy principal value,

[TABLE]

which could be written as

[TABLE]

where $W$ is the cumulant-generating function of $F_{2}$.

II.3 Renormalization

So far we have been able to rewrite our integrals in a way that singles out the divergent factors. We must now renormalize them in such a way as to obtain finite expressions.

The idea is to generate a series of expressions for $U_{PV}$ by replacing $F_{1}$ by another distribution function, better behaved than $F_{1}$ as $u\to 1$, but whose irreducible moments agree with those of $F_{1}$ up to a certain order. It is important that we replace $F_{1}$ by another positive function, since this preserves positivity, and it is important that the replacement pdf includes no dimensionful parameters not present in $F_{1}$, since otherwise conformal invariance would be spoiled. The simplest such replacement would be a $\delta$ function with support at $\left\langle u\right\rangle$. In this paper, we shall restrict ourselves to the next approximation, where both $\left\langle u\right\rangle$ and $\sigma^{2}$ are retained, and $F_{1}$ is replaced by a Gaussian pdf. This is equivalent to considering only the leading term in a Gram - Charlier approximation to $F_{1}$ Cramer57 ; it must be recalled that Gram - Charlier series, when truncated at higher orders, may not be a true distribution function because it may not be nonnegative definite. We shall discuss the accuracy of this lowest order approximation below.

Therefore, we keep up to the quadratic term in the expansion of $W$,

[TABLE]

This approximation leads to

[TABLE]

where $\mathcal{D}$ is the Dawson function, defined as

[TABLE]

The mean values $\left\langle u\right\rangle$ and $\left\langle u^{2}\right\rangle$ are

[TABLE]

and

[TABLE]

We finally arrive at the final result for the $I_{1}$-type integrals,

[TABLE]

We now turn to $I_{2}$-type integrals. The same way as before, we keep up to the quadratic term in the expansion of $W$,

[TABLE]

where $\mathcal{D}$ is the Dawson function defined in Eq. (61) and the relevant mean values are

[TABLE]

and

[TABLE]

We finally arrive at the final result for the $I_{2}$-type integrals,

[TABLE]

To conclude, let us discuss whether approximating $F_{1}$ and $F_{2}$ by a Gaussian distribution is quantitatively correct.

Let us call $\mathcal{F}_{i}$, $i=1,2$ the Gaussian approximation of $F_{i}$, given by

[TABLE]

where $\langle u\rangle_{i}$ and $\sigma^{2}_{i}=\langle u^{2}\rangle_{i}-\langle u\rangle^{2}_{i}$ are computed with $F_{i}$.

In the limit $\zeta_{0}\rightarrow 0$, both $F_{1}$ and $F_{2}$ converge to the Dirac Delta distribution,

[TABLE]

and, since $\langle u\rangle_{i}\rightarrow 0$ and $\sigma^{2}_{i}\rightarrow 0$, the Gaussian approximations converge to the same limit,

[TABLE]

Therefore, the limit $\zeta_{0}\rightarrow 0$ is exact.

Arbitrary values of $\zeta_{0}$ require numerical methods to analyze. Figure (2) shows a comparison between $F_{1}$ and its Gaussian approximation $\mathcal{F}_{1}$ as functions of $\zeta_{0}$, using a function $g$

[TABLE]

It shows that for small $\zeta_{0}$ values both $F_{1}$ and $\mathcal{F}_{1}$ tend to a Dirac Delta distribution centered at $u=0$ but for large $\zeta_{0}$ values $F_{1}$ diverges logarithmically at $u=1$ while $\mathcal{F}_{1}$ tends to a Gaussian distribution with constant mean and variance. Although approximating $F_{1}$ by a Gaussian allows the support of the new distribution to be different than the original interval $\left[0,1\right]$, figure (3), which is the integral of $\mathcal{F}_{1}$ over $\left[0,1\right]$, shows the area under the curve is mostly (at least $75\%$ of it) located in that interval. Similarly, figure (4) shows a comparison between $F_{2}$ and its Gaussian approximation $\mathcal{F}_{2}$ as functions of $\zeta_{0}$, using a function $g$

[TABLE]

For small $\zeta_{0}$ values both $F_{2}$ and $\mathcal{F}_{2}$ tend to a Dirac Delta distribution centered at $u=0$ but, unlike $F_{1}$, $F_{2}$ is finite for large $\zeta_{0}$ values. Figure (5) shows the area under the curve of $\mathcal{F}_{2}$ in $\left[0,1\right]$ is $80\%$ for small $\zeta_{0}$ values and reaches a constant value of $93\%$ for large $\zeta_{0}$ values.

We can conclude that the replacement of $F_{1}$ and $F_{2}$ by their corresponding Gaussian approximation effectively cuts off the integrals in a neighborhood of $u=1$ without the need to include an explicit cut-off, which would add a new dimensionfull parameter to the theory. Moreover, the procedure yields a quantitatively accurate approximation for small $\zeta_{0}$ (it is exact at $\zeta_{0}=0$). While there is a loss of accuracy for large values of $\zeta_{0}$, it must be observed that also the energy-momentum tensor is less sensitive to the exact value of $\zeta_{0}$ in that range, as we will show in next section.

III Pressure Anisotropy in DTTs and in the Grad Approximation

One interesting way to visualize the relationship of $\zeta^{\mu\nu}$ to the energy-momentum tensor in the fully nonlinear DTT is by considering the pressure anisotropy.

We consider an axisymmetric configuration where, in the rest frame, $T^{2}\zeta_{0}^{ij}=\mathrm{diag}\left(\zeta_{0},\zeta_{0},-2\zeta_{0}\right)$ with $\zeta_{0}\geq 0$. We define the pressure anisotropy as

[TABLE]

For comparison, under Grad approximation $\Pi^{\mu\nu}$ is given by Eq. (20). Assuming for $\xi^{\mu\nu}$ the same form as $\zeta^{\mu\nu}$, we get

[TABLE]

It becomes negative for $\xi_{0}>5/8$ and remains negative thereafter, approaching $\delta_{p}\to-2$ as $\xi_{0}\to\infty$.

By contrast, DTTs have a built-in lower limit for the pressure anisotropy, because the deformed Fermi-Dirac distribution Eq. (28) always has a finite dispersion in $p_{z}$. For example, and leaving out renormalization issues for the moment, when $\zeta_{0}\to\infty$ Eq. (28) becomes $f=\Theta\left(1-3\cos^{2}\theta\right)$, where $\Theta$ denotes the step function. Factoring out and cancelling a divergent radial integral this leads to $\delta_{p}=1/4$.

Of course, a correct evaluation of $\delta_{p}$ requires that $T^{\mu\nu}$ is computed by carrying out a proper renormalization procedure. When the integrals are regularized by the Gaussian approximation we have presented above, it is seen that $1/4$ is indeed the asymptotic value of $\delta_{p}$ as $\zeta_{0}\to\infty$, but that higher anisotropy is possible at finite values. A numerical evaluation shows that $\delta_{p}\geq 0.1466\approx 1/7$, which in the Grad approximation corresponds to $\xi_{0}=1/2$, as can be seen in figure (6).

Now let us study the changes in the slope of the DTT anisotropy shown in figure (6) left. If $\delta_{p}=p_{z}/p_{x}$ then the derivative with respect to $\zeta_{0}$ (denoted by primes from now on) is

[TABLE]

where $p_{z}$ and $p_{x}$ are given by

[TABLE]

[TABLE]

its derivatives $p_{z}^{\prime}$ and $p_{x}^{\prime}$ by

[TABLE]

[TABLE]

and the distribution function is as in eq. (28), with $T=1$ for simplicity,

[TABLE]

In the limit $\zeta_{0}\longrightarrow 0$, $\delta_{p}\longrightarrow 1$ and

[TABLE]

Since $p_{z}^{\prime}<0$ it follows that for sufficiently small $\zeta_{0}$, $\delta_{p}^{\prime}<0$. Conversely, in the limit $\zeta_{0}\longrightarrow\infty$, $\delta_{p}\longrightarrow 1/4$ and $f(1-f)$ goes to zero except when $x^{2}\approx 1/3$. Therefore, for big enough $\zeta_{0}$ values, the integrand in Eqs. (79) and (80) is concentrated in a neighborhood of $x^{2}=1/3$. On the one hand we have

[TABLE]

and on the other $p_{z}^{\prime}>0$ because $f\longrightarrow\Theta\left(1-3x^{2}\right)$, so $\delta_{p}^{\prime}>0$. We have proved that $\delta_{p}$ approaches its asymptotic value from below, in accordance with figure (6) left, which thus captures the general behavior of $\delta_{p}$. In particular, the anisotropy parameter reaches a minimum at some value of $\zeta_{0}$, regardless of the regularization and renormalization procedure.

This result also shows that the theory can not describe a configuration with anisotropy parameter less than $1/7$; since such configurations would be very extreme, we do not believe this is a significant drawback. At the same time, these are the cases where one would not expect the distribution could be accurately described by a few of its moments.

For this reason, it seems more important to us that, while an approximation such as Grad’s is prone to unphysical behavior in extreme ranges of parameters, the DTT has built in safety measures against such behavior; in this case, no matter how large $\zeta_{0}$ could become along the evolution, pressures will never become negative. Moreover, the pressure anisotropy is mostly insensitive to the value of $\zeta_{0}$ when it becomes large, underlying that inaccuracies in the approximations made not necessarily propagate to the physical predictions of the theory.

IV Perturbative Stability of Homogeneous Configurations

After outlining the procedure to obtain a well defined DTT out of kinetic theory, we are going to test the resulting theory by considering a problem we can solve both in kinetic theory and in the DTT, and also in a hydrodynamic formalism derived from Grad’s approximation. Concretely, we shall discuss whether non equilibrium, homogeneous, anisotropic (but axisymmetric) configurations are perturbatively stable. We shall show this is the case in kinetic theory, meaning that a non-homogeneous perturbation of such a background always decays in time (observe that the background itself is not a solution of the Boltzmann equation). Then we shall obtain a similar result in the DTT by considering the dynamics of linear perturbations to the $u^{\mu}$ and $\zeta_{\mu\nu}$ degrees of freedom. Finally, we shall show that the dynamics of the variables $u^{\mu}$ and $\xi_{\mu\nu}$ from Grad’s approximation is unstable if the background is anisotropic enough, even before the lowest pressure actually becomes negative.

We are going to assume an homogeneous temperature $T$ and energy density $\rho$. Following Romatschke and Strickland RomatschkeStrickland (see also Florkowski ; Tinti ), $\zeta_{\mu\nu}$ and $\xi_{\mu\nu}$ will be chosen symmetric, transverse and traceless and we will use a coordinate system such that in the unperturbed rest frame, in which $u_{0}^{\mu}=\delta^{\mu 0}$, the background part of $\zeta_{\mu\nu}$ can be written as $T^{2}\zeta^{\mu\nu}_{0}=\text{diag}\left(0,\zeta_{0},\zeta_{0},-2\zeta_{0}\right)$, with $\zeta_{0}\geq 0$, and similarly for $\xi_{\mu\nu}$. Since we are perturbing an homogeneous background the normal modes shall be plane waves $e^{st+ikz}$. Our goal will be to find the dispersion relation $s=s\left(k\right)$ by the three formalisms and compare them. An instability appears if for any $k$, $\mathrm{Re}(s)>0$.

IV.1 Kinetic Theory

We want to solve the kinetic equation for the 1pdf $f$ with an Anderson-Witting collision term Eq. (12) Sch06 ; TI10 . We shall investigate linearized fluctuations around an homogeneous background. To do so, we shall look for the solution of the Boltzmann equation with an initial condition given by a 1pdf of the DTT type Eq. (28), where moreover the parameters may be decomposed into an homogeneous plus a small, position dependence perturbation. We shall assume this dependence is of the form $e^{ikz}$, in the unperturbed fluid rest frame. At late times, when all transients have decayed, the solution will correspond to a normal mode of the Boltzmann equation.

We are going to assume a solution of the form

[TABLE]

where $f_{0}$ is the background,

[TABLE]

and $\delta f$ the perturbation. Likewise, $f_{eq}$ is given by

[TABLE]

where ${f_{eq}}_{0}$ is

[TABLE]

If we write the perturbations in the parameters at $t=0$ as

[TABLE]

and

[TABLE]

then the perturbation in the initial condition is

[TABLE]

For the perturbation in $f_{eq}$ we may write

[TABLE]

where $T$ and $T_{eq}$ are related by Eq. (30). Replacing Eqs. (84) and (86) in Eq. (12) we arrive at the solution

[TABLE]

where $\sigma$ and $F$ are given by

[TABLE]

and

[TABLE]

We will assume the only nonzero perturbations are $v_{1}$ and $z_{13}$. In order to preserve the transversality condition $v_{\mu}\zeta^{\mu\nu}=0$, $z_{01}$ should be nonzero as well and equal to $z_{01}=-v^{1}\zeta_{0}$. Under this assumptions the solution is

[TABLE]

In order to find the dispersion relation we are going to study the long-time behavior of the velocity perturbation. First, recall that the hydrodynamic velocity can be defined as the timelike eigenvector of the energy-momentum tensor, $u_{\mu}T^{\mu\nu}=-\rho u^{\nu}$. For the perturbations we are considering, there are no first order corrections to the energy density. Therefore, up to first order we have

[TABLE]

where $p_{0}$ is the equilibrium pressure and $\delta T^{\mu\nu}$ the non-equilibrium part of the energy-momentum tensor. Replacing Eq. (95) in Eq. (96),

[TABLE]

In the limit $t\rightarrow\infty$ we obtain the asymptotic behavior $v_{1}(t)\propto e^{st}$, provided

[TABLE]

where

[TABLE]

We introduce a new parameter $\gamma$ defined by

[TABLE]

With the aid of $\gamma$, the implicit relation Eq. (98) can be written in parametric form as

[TABLE]

where

[TABLE]

If $\gamma=0$, then the first term of Eq. (102) is equal to $p_{0}/\left(\rho_{0}+p_{0}\right)$, while the second is just $\rho_{eq}/\left(\rho_{0}+p_{0}\right)$. Due to the fact that $\rho=\rho_{eq}$, $\mathcal{H}\left[\gamma=0,\zeta_{0}\right]=1$. Therefore, for $\gamma=0$,

[TABLE]

Since $\mathcal{H}$ is a decreasing function of $\gamma$, we have $\tau s\leq 0$ and, as a consequence, $s$ is a decreasing function of $k$ with initial value $s(k=0)=0$. Finally, $s(k)\leq 0$ for any $\gamma$ and we arrive at the conclusion that the linear theory is stable.

IV.2 DTT

According to the DTT, the dynamics of the gas are governed by the set of equations (1), (9), (21) and (22), where $f$ and $f_{eq}$ are defined as in Eqs. (28) and (29), respectively, not just at $t=0$, but at all times. This means that we have the decompositions Eqs. (88), (89), (90) and (91) not just at $t=0$, but at all times. As before, we assume the only non zero perturbations are $v_{1}$, $z_{13}$ and $z_{01}=-v^{1}\zeta_{0}$. Linearization leads to the system of equations

[TABLE]

where the coefficients $C$, $D$, $E$, $F$, $G$ and $H$ are defined by the integrals

[TABLE]

By keeping the only relevant equations for the $v_{1}$ and $z_{13}$ perturbations (and $z_{01}$ because of the transversality condition), we arrive at the two-by-two linear system

[TABLE]

Since $T$ is the only dimensionful parameter, there is no loss of generality in setting $T=1$. The dispersion relation is given by the secular equation whose solutions are

[TABLE]

where $\Omega_{DTT}$ is defined by

[TABLE]

It is evident that the only interesting case is

[TABLE]

If $\Omega_{DTT}\left(\zeta_{0}\right)<0$ for some $\zeta_{0}$, then an instability arise. To see if such $\zeta_{0}$ exist, we use the following properties of $f_{0}$,

[TABLE]

and

[TABLE]

to write

[TABLE]

Using this identity we can integrate by parts to obtain

[TABLE]

and

[TABLE]

Since $F^{113}/H^{0113}>0$, we conclude that $\Omega_{DTT}>0$ and, just like kinetic theory, perturbations in DTT are stable.

IV.3 Grad’s Approximation

In order to use Grad’s approximation first we need to find Grad’s probability density function. As discussed in the Introduction, it takes the form Eqs. (18) and (19). As before, we assume that $T$ is unperturbed, Eq. (88) for the velocity, and

[TABLE]

where $\xi_{0}^{\mu\nu}=\mathrm{diag}\left(0,\xi_{0},\xi_{0},-2\xi_{0}\right)$, $\xi_{0}\geq 0$. As before, the only nonzero components of the perturbed variables are $v^{1}$, $x^{13}=x^{31}$ and $x^{01}=\xi_{0}v_{1}$.

The system’s dynamics are governed by Eqs. (1), (9), (21) and (22). Replacing Grad’s probability density function Eq. (18) in the previous equations and solving them for $v_{1}$ and $x_{13}$ up to first order we obtain

[TABLE]

The dispersion relation is given by the secular equation whose solutions are

[TABLE]

where $\Omega_{Grad}$ is defined by

[TABLE]

If we chose the negative sign in Eq. (116) and $\xi_{0}>1/2$, then perturbations show an exponential growth with coefficient

[TABLE]

There is not a minimum (or maximum) $k$ value for instabilities. They occur at every value of $k$ as long as $\xi_{0}>1/2$. As we have shown in the previous section, there is a range $1/2<\xi_{0}\leq 5/8$ where all three pressures are positive, but nevertheless this spurious instability appears. For larger $\xi_{0}$ the lowest pressure becomes negative, bringing the breakdown of the theory to the fore.

IV.4 Quantitative Comparison

The main difference between the three methods utilized before is the presence or absence of instabilities. While kinetic theory and DTT show no signs of them, Grad’s approximation has no stable solutions for $\xi_{0}>1/2$. This is a drawback for Grad since it shows its applicability is fairly limited.

All three theories predict an $s_{max}<0$ with its corresponding $k_{max}$ (also in Grad’s approximation $\xi_{0}$ must be less than $1/2$) such that if $k>k_{max}$ then the dispersion relation can be written as

[TABLE]

This means we have propagation in the form of damped waves. In kinetic theory this set of $s_{max}$ and $k_{max}$ are given by

[TABLE]

and

[TABLE]

Both DTT and Grad have the same $s_{max}$,

[TABLE]

but they have a different $k_{max}$. In DTT it’s given by

[TABLE]

while in Grad its value is

[TABLE]

Figure (7) shows a comparison between the three different values of $k_{max}$ as a function of the anisotropy $\zeta_{0},\xi_{0}$. It can be seen that DTT shows a qualitatively similar behavior to kinetic theory, better than Grad’s.

The isotropic case, that is $\zeta_{0}=\xi_{0}=0$, can be solved analytically. Up to second order in $k$ we have, in kinetic theory

[TABLE]

in DTT,

[TABLE]

and in Grad’s approximation,

[TABLE]

Figure (8) shows the full dispersion relations for $\zeta_{0}=\xi_{0}=0$, showing that even up to $\tau k\approx 0.9$ DTT looks very similar to kinetic theory.

Arbitrary values of $\zeta_{0}$ and $\xi_{0}$ require numerical methods to solve, always utilizing the regularization procedure defined in previous sections. Figure (9) shows the dispersion relations for all three theories as functions of $\zeta_{0}=\xi_{0}$. For both small and big values of $\zeta_{0}$, DTT is a good approximation to kinetic theory. Figure (10) shows the dispersion relation for $\zeta_{0}=\xi_{0}=0.45$ and $\zeta_{0}=\xi_{0}=0.55$, that is, before and after Grad’s instability.

V Final Remarks

Formulating a fully nonlinear hydrodynamics of dissipative relativistic fluids is a daunting challenge that we must nevertheless confront if we wish to make sense of the very early stages of relativistic heavy ion collisions Str14b ; JeHe15 ; Rom16 and also of the cosmic evolution in the period that goes, roughly, from reheating after inflation to the electroweak and QCD transitions BTW06 ; AHKK15 ; BVS06 ; CK10 ; CK14 . We believe this paper contributes to the ongoing effort to meet this challenge in two main ways. On the one hand, it delineates the boundary of applicability of a representative “second order” theory. These theories were introduced to solve the instability problems of the so-called “first order” theories HisLin83 ; HiscockLindblom ; HisLin88 ; nevertheless, as we have shown, they display spurious instabilities of their own. On the other hand, we show a definite way whereby a fully nonlinear hydrodynamics may be derived from kinetic theory in a systematic way.

Of course this is not the only strategy that is being tried out MNDR14 . The best known is simply to go to higher orders within the Chapman-Enskog or Grad approaches, as in the so-called Burnett’s equations Burnett ; the second order approximation is discussed in Jais13b ; Jais14 ; Jais15 ; Jais15b . These models very soon become extremely complex, which may become an issue if we consider that the kind of problems we have discussed in this paper are already an oversimplification of the problems we really want to study, and which include gauge and possibly the gravitational field as well.

A promising strategy is to duplicate the Grad approach but taking as zeroth order an already nonequilibrium state, as in the so-called anisotropic hydrodynamics Str14b ; BHS14 ; BNR16 . To the best of our knowledge this approach has been tried so far only in highly symmetric configurationsTRFS15 ; FRST16 ; FMR15 , so it is unclear what hurdles it could encounter in realistic scenarios.

In our view, the DTT framework we are advocating has two distinctive advantages. First, it enforces energy-momentum conservation and the Second Law in a rigorous way, contrary to “second order” theories in which it is enforced only to second order. In typical “second order” theories (see however Jaiswal and MartinezStrickland ), the entropy production, as computed from the entropy flux and the hydrodynamic equations, is nonnegative only if terms of order higher than second are neglected Cristian . In a DTT, the same procedure yields an strictly nonnegative expression to all orders in deviation from ideal behavior. Second, it can describe situations far from equilibrium without the addition of other degrees of freedom than those already present in $T^{\mu\nu}$; this puts a limit on how complex the theory may become, although of course it will never be as simple and compelling as the hydrodynamics of ideal fluids. Moreover, having a fully consistent theory to begin with gives one a solid framework whereby one can discuss simplifications in a systematic way (by contrast, observe that the Eckart expansion is known not to be convergent ChapmanEnskogDivergence ; Resurgence ). It could well be that the main value of the theory we have developed in this paper is that it exists, rather than its actual applications.

In last analysis, to be able to compare several alternative formalisms will be a definite asset for the community as we enter in this largely uncharted territory.

Appendix A: Stability in conformal hydrodynamics

We wish to provide a general template for the discussion of stability against incompressible perturbations in conformal hydrodynamical theories, where the fundamental equations are energy-momentum conservation and a new conservation law of type eq. (24) for some totally symmetric (and traceless on any pair of indexes) tensor $A^{\mu\nu\rho}$. Following Israel and Stewart IsrSt79 , we call an hydrodynamical theory one where the fundamental degrees of freedom are in one-to-one correspondence with the components of the energy-momentum tensor, at least in a neighborhood of the equilibrium states. Therefore, since we are restricting ourselves to conformal theories, the fundamental degrees of freedom can be chosen as a single dimensionful scalar $T$ (which becomes the temperature in equilibrium states), the Landau-Lifshitz velocity $u^{\mu}$ and a dimensionless, symmetric, traceless and transverse tensor $Z^{\mu\nu}$. Since $Z^{\mu\nu}$ is transverse, we cannot build new tensors by contracting it with $u^{\mu}$; the only other linearly independent transverse traceless tensor in the theory is $\hat{Z}^{\mu\nu}=Z^{2\mu\nu}-\left(1/3\right)\mathrm{tr}Z^{2}\;\Delta^{\mu\nu}$. It follows that we have the decomposition

[TABLE]

The scalars $A_{T},A_{A},B,C,D,B^{\prime},C^{\prime}$ and $D^{\prime}$ are functions of $Z^{\mu\nu}$ through invariants such as $\mathrm{tr}Z^{2}$ and $\mathrm{tr}Z^{3}$.

We consider linear perturbations to an homogeneous anisotropic background. This means quantity $X$ becomes $X=X_{background}+\delta X\;e^{st+ikz}$. In the background $u^{\mu}=\delta^{0\mu}$ and $Z^{\mu\nu}=\mathrm{diag}\left(0,Z_{0},Z_{0},-2Z_{0}\right)$ ($\mathrm{tr}Z^{2}=6Z_{0}^{2}$). The only perturbed component of the velocity is $\delta u^{1}=v$. The perturbed components of $Z^{\mu\nu}$ are $\delta Z^{13}=\delta Z^{31}=z$ and $\delta Z^{01}=\delta Z^{10}=Z_{0}v$, as demanded by transversality. It follows that $\delta\hat{Z}^{13}=\delta\hat{Z}^{31}=-Z_{0}z$ and $\delta\hat{Z}^{01}=\delta\hat{Z}^{10}=-Z_{0}^{2}v$, all other zero. $T$ and all the invariants constructed from $Z^{\mu\nu}$ are unchanged. The relevant equations of motion

[TABLE]

become

[TABLE]

leading to the dispersion relation (compare to eqs. (106) and (116))

[TABLE]

where

[TABLE]

If $\tau<0$ the theory is always unstable for long wavelenghts, which is clearly unphysical. If $\tau>0$, the theory becomes unstable for short wavelenghts if $\Omega<0$. So stability requires both $\tau>0$ and $\Omega>0$.

However this condition cannot be met if we force $T^{\mu\nu}$ and $A^{\mu\nu\rho}$ to be linear functions of $Z^{\mu\nu}$. This amounts to defining $B^{\prime}=C^{\prime}=D^{\prime}=0$ and $A_{T},A_{A},B,C$ and $D$ to be constants. Moreover stability at equilibrium implies that $\Omega\left(Z_{0}=0\right)=BA_{A}/4CA_{T}>0$, so $C/A_{A}$ and $B/A_{T}$ must have the same sign. Therefore $\Omega$ becomes a rational function which cannot be nonnegative everywhere. This is what happens in the Grad theory, where, after identifying $Z_{0}=\xi_{0}$, we obtain $A_{T}=\sigma_{SB}$, $B/A_{T}=4/15$ and $C/A_{A}=1/3$. $\Omega=\Omega_{Grad}$ is given by eq. (117), which is clearly negative for $\xi_{0}>1/2$.

In the DTT, on the other hand, we have $Z_{0}=\zeta_{0}$ and $\Omega=\Omega_{DTT}$ given by eq. (107), which as we have seen is indeed positive for $\zeta_{0}>0$.

We show both $\Omega_{DTT}$ and $\Omega_{Grad}$ as functions of $\zeta_{0}=\xi_{0}$ in figure (11). Both start from a positive value at $\zeta_{0}=0$ with a negative slope, but while $\Omega_{Grad}$ is monotonous and eventually reaches the asymptotic value of $-2$, $\Omega_{DTT}$ changes its tendency and remains positive for all $\zeta_{0}>0$. This underlies that the instability of the equations from the Grad approximation is an artifact of the linearization of $T^{\mu\nu}$ and $A^{\mu\nu\rho}$ with respect to $\xi^{\mu\nu}$.

Appendix B: Extension to Maxwell-Jüttner and Bose-Einstein statistics

In this appendix we will show that the same regularization and renormalization procedure we have applied to Fermi-Dirac (FD) particles may be used for particles obeying Maxwell-Jüttner (MJ) or Bose-Einstein (BE) statistics. Our starting point is the observation that all three statistics may be obtained as particular cases of theories described by the family of one-particle distribution functions

[TABLE]

where $a=0$ for FD, $a\to\infty$ yields MJ and $a=\pm i\pi$ gives BE. The idea is to obtain the expectation value of some function $g(p,\theta)$ (as in the main text) as a function of $a$ for $a$ real and positive, and then try and extend the result to the MJ and BE cases. We only consider case I above (cfr. eq. (34))

[TABLE]

As in the main text, we start by dividing the $I_{1a}$ integral in $\theta_{0}=\text{cos}^{-1}\left(1/\sqrt{3}\right)$,

[TABLE]

$I_{1a}^{<}$ is well defined in all three cases and will be left as is. To study the $I_{1a}^{>}$ term, define a function $G$ as in the main text (eq. (38)) and integrate by parts to obtain

[TABLE]

it is evident that the first term is finite in all three cases. Let’s call the second term $K_{Za}$. By using the sum of arguments relation of the hyperbolic cosine and realizing it could be written as a partial derivative, it is possible to rewrite $K_{Za}$ as

[TABLE]

Performing another integration by parts and defining the auxiliary function

[TABLE]

we arrive at the final expression for $K_{Za}$,

[TABLE]

The key here is to identify the ratio in Eq. (139) as the mean value

[TABLE]

where $F_{1a}$ is a probability density function defined as

[TABLE]

We leave open the range of $u$. As in the main text we replace $U_{a}$ by its Cauchy principal value,

[TABLE]

The idea is to approximate $F_{1a}$ by a Gaussian, for which we need the mean values $\left\langle u\right\rangle_{a}$ and $\left\langle u^{2}\right\rangle_{a}$

[TABLE]

and

[TABLE]

The limit $a\to\infty$ is not problematic. If $a=\pm i\pi$ we have the identities

[TABLE]

Again, the relevant expectation values are well defined.

Acknowledgments

Work supported in part by CONICET and University of Buenos Aires. It is a pleasure to thank R. Ferraro, A. Jaiswal, A. Kandus, F. Lombardo, P. Mininni, N. Mirón y C. Vega for discussions.