On the relative importance of second-order terms in relativistic dissipative fluid dynamics

In Denicol et al., Phys. Rev. D 85, 114047 (2012), the equations of motion of relativistic dissipative fluid dynamics were derived from the relativistic Boltzmann equation. These equations contain a multitude of terms of second order in Knudsen number, in inverse Reynolds number, or their product. Terms of second order in Knudsen number give rise to non-hyperbolic (and thus acausal) behavior and must be neglected in (numerical) solutions of relativistic dissipative fluid dynamics. The coefficients of the terms which are of the order of the product of Knudsen and inverse Reynolds numbers have been explicitly computed in the above reference, in the limit of a massless Boltzmann gas. Terms of second order in inverse Reynolds number arise from the collision term in the Boltzmann equation, upon expansion to second order in deviations from the single-particle distribution function in local thermodynamical equilibrium. In this work, we compute these second-order terms for a massless Boltzmann gas with constant scattering cross section. Consequently, we assess their relative importance in comparison to the terms which are of the order of the product of Knudsen and inverse Reynolds numbers.


I. INTRODUCTION AND CONCLUSIONS
Relativistic fluid dynamics has found widespread applications in heavy-ion physics, in modeling nuclear collisions at ultrarelativistic bombarding energies [1,2], in astrophysics, for instance in modeling binary mergers of compact stellar objects (see e.g. Ref. [3]), as well as in cosmology [4][5][6][7][8]. In the past, in order to solve the equations of motion of relativistic fluid dynamics, one has often made the assumption that the fluid is ideal, i.e., one demands instantaneous local thermodynamical equilibrium, which in turn allows to neglect all dissipative effects. However, there are no ideal fluids in nature, as can be seen for instance from the fact that the shear viscosity coefficient may attain a lower limit, but never vanishes [9][10][11].
A more realistic modeling of the dynamics of relativistic fluids thus demands that one uses the equations of relativistic dissipative fluid dynamics. The first attempts to formulate such equations were made by Eckart [12] and Landau and Lifshitz [13] based on a relativistic generalization of the nonrelativistic Navier-Stokes equations. However, their equations suffer from instabilities and acausal signal propagation [14]. The reason for this is the (erroneous) assumption that the dissipative quantities, like bulk viscous pressure Π, particle diffusion current n μ , and shear-stress tensor π μν , react instantaneously to the thermodynamic forces, like gradients of the fluid velocity or temperature and chemical potential. If one relaxes this assumption by introducing certain time scales τ Π , τ n , and τ π , on which the dissipative quantities are allowed to approach the values determined by the corresponding thermodynamic forces, these problems can be cured (provided the relaxation times fulfill certain conditions [15]). Similar to earlier works by Grad [16] and Müller [17][18][19] in the nonrelativistic context, Israel and Stewart (IS) were among the first to suggest equations of motion for relativistic dissipative fluids that were stable and causal [20][21][22].
In Ref. [24] a derivation of the equations of motion of relativistic dissipative fluid dynamics was presented, which is based on a systematic power-counting scheme in the Knudsen and inverse Reynolds number. The Knudsen number, Kn ¼ λ=L, is the ratio between a characteristic microscopic time/length scale, λ (e.g., the mean-free path between collisions), and a characteristic macroscopic scale of the fluid L. In this context, the inverse Reynolds numbers are the ratios of dissipative quantities and (local) equilibrium values of macroscopic fields, e.g. R −1 Π ¼ jΠj=P 0 , R −1 n ¼ jn μ j=n 0 , or R −1 π ¼ jπ μν j=P 0 , where P 0 is the thermodynamic pressure and n 0 is the particle density in equilibrium. The time scales τ Π , τ n , and τ π are identified with the slowest microscopic time scales of the Boltzmann equation.
The physical picture that emerges is that microscopic processes (i.e., in the case of the Boltzmann equation, binary collisions) occur on time scales smaller than (or, at most, as large as) τ Π , τ n , and τ π . These processes affect that the dissipative quantities Π, n μ , and π μν approach the values given by the (relativistic generalization of the) Navier-Stokes equations on the time scales τ Π , τ n , and τ π . Since microscopic physics influences the motion of the fluid only on short time scales, the term "transient fluid dynamics" was coined for such theories of relativistic dissipative fluid dynamics. The fact that the microscopic dynamics of the Boltzmann equation gives rise to relaxation-type equations of motion for the dissipative quantities, i.e., where these quantities exponentially decay towards the values given by the Navier-Stokes equations, was confirmed in Ref. [52]. It was also shown in that paper that approaches based on the AdS/CFT correspondence lead to equations of motion which are of the type encountered for an underdamped harmonic oscillator. The relaxation towards the Navier-Stokes values is then not exponential but oscillatory.
Let us recall the relaxation-type equations for the dissipative quantities derived in Ref. [24], where the overdot denotes the proper time derivative, Here, ζ is the coefficient of the bulk viscosity, κ the coefficient of particle diffusion (which is related to that of heat conduction) and η the coefficient of shear viscosity. Furthermore, with the fluid four-velocity u μ (chosen in the Landau frame), where u μ u μ ¼ 1, with Δ μν ¼ g μν − u μ u ν being the three-projector onto the subspace orthogonal to u μ , and with ∇ μ ¼ Δ μν ∂ ν being the threegradient, θ ¼ ∇ μ u μ is the expansion scalar, σ μν ≡ ∇ hμ u νi ¼ 1 2 ð∇ μ u ν þ ∇ ν u μ Þ − 1 3 θΔ μν is the shear tensor, while I μ ¼ ∇ μ α 0 is the gradient of α 0 ¼ μ=T, the ratio of chemical potential μ and temperature T.
For massless particles, the bulk viscous pressure vanishes identically, Π ¼ 0, and we do not have to solve Eq. (1). Also, terms in Eqs. (4)- (8) proportional to Π can be neglected, such that we do not need to compute the corresponding coefficients. Moreover, as mentioned above, in the 14-moment approximation all coefficients in Eqs. (5) vanish. Dividing Eq. (2) by n 0 , and Eq. (3) by P 0 , these two equations can be written in the following form, Here, we have made use of the equation of state of the massless Boltzmann gas, P 0 ¼ n 0 =β 0 , with β 0 ¼ 1=T. By dividing the dissipative quantities by n 0 or P 0 , respectively, we immediately identify terms which are proportional to the inverse Reynolds number. Furthermore, the coefficients of terms involving gradients (or time derivatives) all have dimension of time (or mean-free path) and are thus proportional to the Knudsen number. In this form, it is easy to apply power-counting arguments to estimate the order of magnitude of the various terms. The Navier-Stokes terms appearing in the first lines are of first order in the Knudsen number. The second lines contain terms proportional to the dissipative quantity that is evolved in the respective equation (n μ in the first and π μν in the second equation), while the third lines contain cross terms proportional to the other dissipative quantity that is not evolved (π μν in the first and n μ in the second equation). Terms in the second and third lines are of first order in the product of the Knudsen and inverse Reynolds number as well as of second order in the inverse Reynolds number. At this point, one cannot draw any further conclusion without making assumptions about the relative magnitude of the Knudsen and inverse Reynolds numbers. In the following, let us assume that all of them are of the same order of magnitude, Kn ∼ R −1 n ∼ R −1 π (situations where this is no longer the case were studied, e.g., in Ref. [53]). At least for asymptotically long times, when the values of the dissipative quantities approach their respective Navier-Stokes limits, this assumption is fulfilled.
Then we can simply assess the order of magnitude (and thus the importance) of the various terms by comparing the values of the coefficients accompanying them. These are listed in Tables I [for Eq. (9)] and II [for Eq. (10)]. Note that the coefficients of terms containing a power of the Knudsen number are given in units of the relaxation times τ n and τ π , respectively.
One immediately observes that the coefficients in the third line of Eq. (9) are (at least) an order of magnitude smaller than the ones in the second line. If the Knudsen and inverse Reynolds numbers are of the same order of magnitude, it is then a reasonable assumption to drop these terms altogether, including the last term in Eq. (9) that arises from nonlinear terms in the collision integral.
In the 14-moment approximation, all terms in the third line of Eq. (10), except for the last one, vanish identically. We remark, however, that this is accidental; going beyond the 14-moment approximation [24], all terms are of the same order of magnitude. From the terms in the second line, the last one (arising from the nonlinear terms in the collision TABLE I. The coefficients in the particle diffusion equation for a massless Boltzmann gas with constant cross section σ T in the 14-moment approximation (from Ref. [24]). λ mfp ¼ 1=ðn 0 σ T Þ is the mean-free path. κ n 0 ½λ mfp τ n ½λ mfp λ nn ½τ n δ nn ½τ n l nπ β 0 ½τ n τ nπ P 0 β 0 ½τ n λ nπ β 0 ½τ n φ 4 P 0 integral) is an order of magnitude smaller than the other second-order terms and hence can safely be neglected.
To summarize our results, we conclude that, in the 14moment approximation, for a massless Boltzmann gas, and for situations where the Knudsen and inverse Reynolds numbers are of the same order of magnitude, the equations of motion for the dissipative quantities read, to good approximation, as follows: In the application of relativistic dissipative fluid dynamics to describe the dynamics of nuclear collisions, this form of the equations should be applicable (far) above the QCD transition, i.e., (deep) in the quark-gluon plasma phase, where the velocity of sound is close to that of an ultrarelativistic gas, c 2 s ¼ 1=3, and the mass of the quasiparticles can be neglected. Nevertheless, the values of the transport coefficients may change from the values given in Tables I  and II when accounting for proper quantum statistics, the correct number of degrees of freedom, as well as a more realistic (angular-dependent) scattering cross section. Closer to the QCD phase transition, however, the breaking of conformality and the effect from nonvanishing bulk viscous pressure Π can no longer be neglected.
The remainder of this paper is organized as follows. In Secs. II A and II B we briefly recapitulate the derivation of the equations of relativistic dissipative fluid dynamics from the Boltzmann equation using the method of moments [24]. In Sec. II C we give a general derivation of the second-order corrections to the collision integral. Explicit results will be derived in Sec. III by reducing the number of moments to 14. Details of our calculations are relegated to several Appendices.

A. General variables
In relativistic kinetic theory of dilute single-component gases, an ensemble of particles with mass m and four-momenta , at a given spacetime point x μ ¼ ðt; xÞ, is characterized by the invariant singleparticle distribution function f k ðt; xÞ ≡ f k . In the absence of external forces or fields, the spacetime evolution of f k is described by the relativistic Boltzmann equation [54,55], where C½f is the collision term. In the case of binary elastic collisions it is given by Heref k ¼ 1 − af k , with a ¼ 1ð−1Þ for fermions (bosons) and a ¼ 0 for classical particles, dK ¼ gd 3 k=½ð2πÞ 3 k 0 is the Lorentz-invariant momentum-space volume, where g denotes the number of internal degrees of freedom, and W kk0→pp0 is the Lorentz-invariant transition rate. The factor ν ¼ 2 takes into account that the particles are indistinguishable, while the transition rate satisfies detailed balance, i.e., it is symmetric for time-reversed states W kk0→pp0 ¼ W pp0→kk0 . The conservation of particle number and of energy and momentum in individual collisions leads to the following continuity equations for the particle four-current, N μ , and energy-momentum tensor, T μν [54,55]: The currents N μ ¼ N μ ðt; xÞ and T μν ¼ T μν ðt; xÞ are identified as the first and second moments of the single-particle distribution function, respectively. Without loss of generality, they can be tensor decomposed in terms of the fluid four-velocity as The coefficients in the shear-stress tensor equation for a massless Boltzmann gas with constant cross section σ T in the 14-moment approximation (from Ref. [24]). λ mfp ¼ 1=ðn 0 σ T Þ is the mean-free path.
Using Eqs. (17)- (18) we are able to identify the fundamental fluid-dynamical quantities as where n is the particle number density, ε is the energy density, P is the local isotropic pressure, n μ is the particle diffusion current, and π μν is the shear-stress tensor. It is customary to separate the isotropic pressure into two components, P ¼ P 0 þ Π, with P 0 being the thermodynamic pressure and Π the bulk viscous pressure. The thermodynamic and bulk viscous pressures are defined with respect to the local equilibrium distribution function, and The temperature and chemical potential, μ ¼ Tα 0 , introduced in f 0k are defined by the so-called matching conditions which impose that the particle number density and energy density are given by their respective values in a fictitious local thermodynamic equilibrium state, i.e., n 0 ≡ hE k i 0 ¼ n and ε 0 ≡ hE 2 k i 0 ¼ ε. In this state, there is an equation of state of the form P 0 ðT; μÞ, such that n 0 ¼ ∂P 0 =∂μ, s 0 ¼ ∂P 0 =∂T, and the fundamental thermodynamical relation ε 0 ¼ Ts 0 þ μn 0 − P 0 is fulfilled.
It is also convenient to introduce the irreducible moments of δf k , Such irreducible moments are constructed to be symmetric, traceless, and orthogonal to the four-velocity, with the symmetrized,traceless,andorthogonalprojectionsbeingdefinedas k hμ 1 :::::k μ l i ¼ Δ μ 1 …μ l ν 1 …ν l k ν 1 …k ν l : The details of construction and the properties of such tensors can be found in Appendix B. Note that the bulk viscous pressure, particle diffusion current, and shear-stress tensor are also irreducible moments of δf k , Furthermore, the matching conditions and the definition of the local rest frame can also be expressed using irreducible moments. The matching conditions correspond to while the Landau definition of the fluid four-velocity leads to B. Moment expansion of f k and the equations of motion for the irreducible moments Following Ref. [24], δf k is parametrized as The function ϕ k is then expanded in terms of a series in the irreducible tensors given in Eq. (28), By expanding the tensor λ hμ 1 …μ l i k using a set of orthogonal polynomials, it is straightforward to prove that H ðlÞ kn ρ μ 1 …μ l n ; (35) where N l denotes the order at which the expansion is truncated (for the coefficient of rank l) and ρ μ 1 :::::μ l n is the irreducible moment defined in Eq. (27) using Gram-Schmidt orthogonalization. The measure ω ðlÞ depends on the rank l of the tensor being expanded and reads where W ðlÞ is a normalization constant defined as For more details, see Refs. [24,25].
Using the Boltzmann equation, one can derive the general equations of motion satisfied by ρ μ 1 :::::μ l r . This is accomplished by explicitly taking the comoving derivative of the corresponding irreducible moment, i.e., _ ρ dKE r k k hν 1 …k ν l i δf k , and using the Boltzmann equation to express the comoving derivative of δf k in terms of the collision term, f 0k and its derivatives, and spatial derivatives of δf k . The details of this derivation as well as the general form of the resulting equations of motion are contained in Refs. [24,25]. For the three lowest-rank moments, these equations of motion read Here we define the following thermodynamic quantities: where h 0 ¼ ðε 0 þ P 0 Þ=n 0 denotes the enthalpy per particle and The variables I nþr;q ðα 0 ; β 0 Þ and J nþr;q ðα 0 ; β 0 Þ correspond to thermodynamic integrals defined as More details can be found in Appendix A.
We also introduced the generalized irreducible collision integral C μ 1 :::::μ l r and its symmetric, traceless, and orthogonal projection, As was shown in Ref. [24], the moment equations (41)-(43) reduce to the fluid-dynamical equations for the dissipative variables when the fast-varying modes of the Boltzmann equation can be neglected and, simultaneously, the Knudsen number(s) and inverse Reynolds number(s) are small. In this case, the linear parts of the collision integrals introduced above determine the relaxation times for the dissipative variables, while their nonlinear parts give rise to the terms that are of second order in the inverse Reynolds number(s), i.e., the tensors R, R μ , and R μν that appear in Eqs. (6)- (8). In Ref. [24] the existence of such nonlinear terms was pointed out, but the explicit calculation of the corresponding transport coefficients was left for future work. In the next sections we shall complete this task.

C. Expansion of the collision integral
In this section we show how to express the collision integrals in terms of irreducible moments of δf k . Substituting the expression of the collision term for binary elastic collisions (14) into the expression for the irreducible collision integral (51), one obtains Substituting the distribution function f k ¼ f 0k þ f 0kf0k ϕ k into the above formula and using together with the equality f 0k f 0k 0f 0pf0p 0 ¼ f 0p f 0p 0f 0kf0k 0 , the part that is linear in ϕ k reads L hμ 1 ::::: where we abbreviated f 0k f 0k 0f 0pf0p 0 and used the fact that the collision term vanishes for the local equilibrium distribution function C μ 1 …μ l r ½f 0 ¼ 0. Inserting the expression for ϕ k from the moment expansion (34) into the previous equation, we obtain L hμ 1 ::::: It was shown in Ref. [24] that the linear part of the collision integral simplifies to where While using the properties of the irreducible projection tensors, one can show that The coefficient A ðlÞ rn is the (rn) element of a ðN l þ 1Þ × ðN l þ 1Þ matrix, A ðlÞ , and, in the linearized case, contains all the information of the underlying microscopic theory. We remark that, for l ¼ 0, the second and third rows and columns (r, n ¼ 1, 2) and, for l ¼ 1, the second row and column (r, n ¼ 1) are zero, because the moments ρ 1 , ρ 2 , and ρ μ 1 vanish due to the matching conditions and our choice of frame.
The computation of the nonlinear part of the collision integral is analogous. Inspecting the previous formulas we observe that the collision integral is a quartic function of ϕ k . However, in this paper we shall restrict our calculations to the case of Boltzmann statistics (a ¼ 0), in which case the dependence on ϕ k becomes quadratic. The collision integral can be written as where the quadratic contribution to the collision integral reads This nonlinear contribution can be further simplified to where we defined the following tensor of rank l þ m þ m 0 : In comparison with Eq. (61), we have split the double sum P m m 0 ¼0 , and subtracted the superfluous term m ¼ m 0 in the last sum with the help of a Kronecker delta. Then we interchanged indices m↔m 0 , n↔n 0 in the second sum.
The tensor ðN rnn 0 Þ μ 1 …μ l α 1 …α m β 1 …β m 0 is symmetric under permutations of μ-type, α-type, and β-type indices, and depends solely on equilibrium distribution functions and the corresponding cross section(s). The equilibrium distribution function contains only one four-vector, i.e., the fluid four-velocity u μ . Therefore, ðN rnn 0 Þ μ 1 …μ l α 1 …α m β 1 …β m 0 must be constructed from tensor structures made of u μ and the metric tensor g μν , or, equivalently, u μ and Δ μν . Furthermore, ðN rnn 0 Þ μ 1 …μ l α 1 …α m β 1 …β m 0 must be orthogonal to u μ , which implies that it can only be constructed from combinations of elementary projection operators, Δ μν . This already constrains the rank of the tensor, l þ m þ m 0 , to be an even number. Finally, it must satisfy the following property: For our purposes it is sufficient to calculate terms that are of second order in the inverse Reynolds number, i.e., the terms R, R μ , and R μν . Therefore, we only need to consider the cases l ¼ 0, l ¼ 1, and l ¼ 2. Since the actual deduction of the nonlinear collision integrals is complicated, this task is relegated to Appendix C and here we shall only give the final results.
The scalar nonlinear collision integral from Eq. (62) is given by where C 0ðm;mÞ rnn 0 is the special case l ¼ 0 of a more general coefficient Similarly, the nonlinear collision term for l ¼ 1 becomes, where the coefficient C Finally, the rank-2 tensor terms are obtained taking l ¼ 2, where C 2ðm;mþ2Þ rnn 0 can be calculated from Eq. (66); we introduce another coefficient, The normalization d ðmÞ is complicated and is discussed in Appendix C together with other details of the derivation of the nonlinear collision term.

III. TRANSPORT COEFFICIENTS IN THE 14-MOMENT APPROXIMATION
In this section we calculate the previously introduced coefficients A in the 14-moment approximation. As shown in Refs. [24,25], this corresponds to the truncation N 0 ¼ 2, This implies that the following irreducible moments appear: As one can see, they are uniquely related to the dissipative quantities.
Before proceeding and for the sake of later convenience, we reexpress the coefficients H ðlÞ kn using Eqs. (36) and (37) as Note that, for n ¼ 0, the second sum in Eq. (70) identically vanishes, which greatly simplifies the calculation of the collision integral. Furthermore, from the definition of the irreducible moments and using Eqs. (34)-(37) together with the orthogonality condition (B8) we obtain the following general result: ρ μ 1 :::: (72) where we used Eq. (70) in the last step. Therefore, truncating the above general result in the 14-moment approximation we obtain 20 J rþ2;0 ÞΠ; 10 J rþ3;1 Þn μ ; 20 , are calculated from Eq. (71) and listed in Appendix D. These linear relations between the moments are the main result of the 14-moment approximation, which was also obtained in Ref. [25].
It is straightforward to show using Eqs. (58), (59), and (70) that the A ðlÞ rn coefficients of the linear collision term can be expressed in terms of A ðlÞ rn . For l ¼ 0 where, in the 14-moment approximation, where the integrals proportional to A ð0Þ 00 Þ ¼ 0 vanish due to particle number and energy conservation in binary collisions. Here, we introduce the rank-4 tensor which is symmetric upon the interchange of indices ðμ; νÞ and ðα; βÞ, i.e., X μναβ ðrÞ ¼ X
Comparing the above result to Eq. (6) and taking into account the corresponding relaxation time from Eq. (1), we obtain Similarly, the vector term (67) in the 14-moment approximation leads to the following formula: Now, recalling Eq. (7) we obtain Finally, the tensor term (68) is and a comparison with Eq. (8) yields In order to calculate these coefficients we introduce five new tensors, which are similar to X μναβ ðrÞ and are given as follows: Note that the Y iðrÞ;j terms in the previous equations are different contractions of these five tensors. Our notation is such that the i index specifies the tensor while the j index labels a particular contraction. More details are given in Appendix F. We have shown earlier that the coefficients in the equations of motion depend explicitly on the choice of the moment, i.e., the index r. Therefore, once the 14-moment approximation is enforced, any moment of the Boltzmann equation leads to a closed set of equations, but when calculating the coefficients of the nonlinear collision integrals one has to account for the exact form of the relaxation equations which follow from Eqs. (41)- (43) using Eqs. (73)-(75). As an example we quote the equations for the particle diffusion current and shear-stress tensor for arbitrary r, Note that in order to recover Eqs. (1)-(3) one must take and τ π ¼ τ 0 π . Using the above relations, it was already shown in Ref. [25] how to derive the equations of motion and calculate the transport coefficients for different choices of the moments corresponding to the traditional method by Israel and Stewart [22] and to the one proposed by Denicol, Koide, and Rischke (DKR) [23]. Here, we shall follow this recipe and calculate the coefficients of the nonlinear collision integral in both cases.
The equations of motion derived by Israel and Stewart [22] can be obtained by the choice (73)-(75) and substituting these values into the equations of motion (41)- (43). Therefore, in the IS theory all coefficients need to be calculated with r ¼ 3 for the scalar moments, r ¼ 2 for vector moments, and r ¼ 1 for second-rank tensor moments. In contrast, the choice of DKR is to use r ¼ 0 for all equations and coefficients.
We explicitly compute some of these coefficients in the ultrarelativistic limit, mβ 0 → 0, for a classical gas (a ¼ 0) with fixed cross section. Since in this limit Π ¼ 0, we do not need to compute the coefficients φ 1 , φ 2 , and φ 3 appearing in the term R, Eq. (6), which enters the equation of motion (1) for the bulk viscous pressure. Furthermore, φ 5 , φ 6 are coefficients in terms which are proportional to Π, and thus also need not be computed. The remaining, nonvanishing coefficients corresponding to the DKR choice r ¼ 0 are simply denoted as φ 4 , φ 7 , and φ 8 , while the ones corresponding to the IS choice are denoted by φ IS 4 for r ¼ 2, while φ IS 7 and φ IS 8 for r ¼ 1. They read We observe that only φ 8 (multiplying n hμ n νi ) differs in sign between the DKR and IS choices, with the absolute magnitude of the latter being half as large. The coefficients φ 7 (multiplying π λhμ π νi λ ) are approximately of the same magnitude for both choices, while φ 4 (multiplying n ν π μν ) is more than six times smaller in DKR than in IS theory. The implications of these results have already been discussed in the Introduction and Conclusions for the DKR choice.
Finally, for further reference, we also quote the coefficients of particle diffusion and shear viscosity, where λ mfp ¼ 1=ðn 0 σ T Þ is the mean-free path and σ T is the total cross section. Note that all remaining transport coefficients from Eq. (4) were already computed in Ref. [25] for both the DKR and IS choices and it was shown that the differences are of the order of 10%-30%. Therefore, the only coefficients that change considerably from one formalism to the other are φ 4 and φ 8 .
In closing we remark that the numerical solutions of both the IS and DKR theories were compared to the numerical solutions of the Boltzmann equation in various cases [23,53,[59][60][61][62]. These investigations showed the advantages of the DKR choice for the corresponding coefficients, which leads to a far better agreement with numerical solutions of the Boltzmann equation than the IS theory.
The b nq coefficient is equal to the number of permutations in the set ℘ μ , which lead to identical tensor products of the u μ and Δ μν projectors [22], see Eq. (A2) of Ref. [22]. The thermodynamic integrals I nq and J nq were defined in Eqs. (49), (50), Replacing ðΔ αβ k α k β Þ q ¼ ðm 2 − E 2 k Þ q we get the following recursion relations for 0 ≤ q ≤ n=2, I nþ2;q ¼ m 2 I n;q þ ð2q þ 3ÞI nþ2;qþ1 ; (A5) while an integration by parts of Eq. (A3) leads to the following relation: with a similar relation for dJ nq ðα 0 ; β 0 Þ.
for any given distribution of μ-type indices in the product Δ μ 2qþ1 ν 2qþ1 :::::Δ μ n ν n , there are ðn − 2qÞ! possible ways to distribute the ν-type indices, which lead to the same product of projectors.

APPENDIX C: REDUCTION OF COLLISION TENSORS
In this Appendix, we show how to derive the general structure of the collision integrals introduced in the main text. As already discussed, the tensor structure of ðN rnn 0 Þ μ 1 :::::μ l α 1 :::::α m β 1 :::::β m 0 can only be constructed from tensors formed using projection operators Δ μν . We start by collecting all possible combinations of projection operators that can appear in ðN rnn 0 Þ μ 1 :::::μ l α 1 :::::α m β 1 :::::β m 0 : (i) Terms where all μ-type indices pair up on projectors, all α-type indices pair up on projectors, and all β-type indices pair up on projectors, e.g. All possible permutations of the μ, α, β-type indices among themselves are allowed. In this case, l, m, and m 0 must all be even. (ii) Terms where at least one μ-type index pairs with an α-type index on a projector, or one μ-type index pairs with a β-type index on a projector, or one α-type index pairs with a β-type index on a projector, e.g.
In this case, only projectors of the type Δ μ i α j , Δ μ i β j or Δ α i β j exist, with no leftover projectors containing indices of the same type. Such terms have the form β q Δ α r β s :::::: Again, all permutations of the μ, α, β-type indices among themselves are allowed. It is important to emphasize that terms of the type (i) and (ii) by themselves do not satisfy the property (64), since they are not traceless. This can also be seen from the fact that any term which contains at least one projector of the type Δ μ i μ j , Δ α p α q , or Δ β r β s vanishes when contracted with . Thus, ðN rnn 0 Þ μ 1 :::::μ l α 1 :::::α m β 1 :::::β m 0 cannot be solely constructed from terms of type (i) and (ii) and there must be at least one term of type (iii).

APPENDIX E: THE COEFFICIENTS OF THE COLLISION TERMS
In this Appendix, we calculate the coefficients of the nonlinear collision integral in the 14-moment approximation. The nonlinear scalar term N r from Eq. (65) is expanded with the help of the following coefficients: where terms proportional to ðA ð0Þ 00 Þ 2 and A ð0Þ 00 A 10 vanish on account of energy conservation in binary collisions. The above result can be reexpressed using the X ðrÞ and Y iðrÞ tensors given in Eqs. (77) and (106)-(110). Thus, after some calculation we obtain C 0ð00Þ r00 ¼ ðA After this step we still need to evaluate the terms Y μναβ 2ðr−3Þ u μ u ν u α u β , Y μναβκλ 5ðr−3Þ u μ u ν u α u β u κ u λ , etc. This is relegated to Appendix F, for example Y μναβ 2ðr−3Þ u μ u ν u α u β ¼ Y 2ðrÞ;1 as shown in Eq. (F9). In a similar fashion we repeat the calculation for all components and later, in Appendix H, we calculate them in the massless limit.

APPENDIX F: TENSOR DECOMPOSITIONS
In this Appendix, we discuss the decompositions of all collision tensors required in the 14-moment approximation. The collision tensor X μναβ ðrÞ ¼ X ðμνÞðαβÞ ðrÞ from Eq. (77) is symmetric upon the interchange of indices ðμ; νÞ and ðα; βÞ, and it is also traceless in the latter indices, X μναβ ðrÞ g αβ ¼ 0, which follows from the massshell condition, k μ k μ ≡ p μ p μ ¼ m 2 . Using these properties one can show that X μναβ ðrÞ is a spatially isotropic tensor which can be constructed using the four-velocity u μ , the projector Δ μν , and different scalar coefficients x ij . The most general decomposition of X μναβ ðrÞ which is symmetric upon the interchange of indices ðμ; νÞ and ðα; βÞ is The indices of the scalar coefficients x ij are chosen such that the second index (j) denotes the number of projection tensors belonging to the respective coefficient while the first index (i) counts the number of such coefficients. For example x 10 is the coefficient without any projection tensor, while x 12 is the first coefficient which contains two elementary projection tensors.

APPENDIX G: COLLISION INTEGRALS IN THE BOLTZMANN LIMIT
In this Appendix we calculate the previously defined collision tensors Eqs. (77)-(110) in the massless Boltzmann limit. These tensors are based on the collision integral C½f defined in Eq. (14), where the Lorentzinvariant transition rate, W kk0→pp0 , may depend only on the following collision invariant: s ≡ ðk μ þ k 0μ Þ 2 ¼ ðp μ þ p 0μ Þ 2 : In the center of mass (CM) frame, where the sum of threemomenta vanishes k þ k 0 ¼ p þ p 0 ¼ 0 and so k 0 ¼ k 00 ¼ p 0 ¼ p 00 , the scattering angle θ s is given with the help of another collision invariant, t ≡ ðk μ − p μ Þ 2 ¼ ðk 0μ − p 0μ Þ 2 , such that The above defined collision invariants, together with u≡ ðk μ − p 0μ Þ 2 ¼ ðk 0μ − p μ Þ 2 , are the so-called Mandelstam variables, satisfying s þ t þ u ¼ 4m 2 , where m is the mass of the particles.