Coupled kinetic equations for quarks and gluons in the relaxation time approximation

Kinetic equations for quarks and gluons are solved numerically in the relaxation time approximation for the case of one-dimensional boost-invariant geometry. Quarks are massive and described by the Fermi-Dirac statistics, while gluons are massless and obey Bose-Einstein statistics. The conservation laws for the baryon number, energy, and momentum lead to two Landau matching conditions which specify the coupling between the quark and gluon sectors and determine the proper-time dependence of the effective temperature and baryon chemical potential of the system. The numerical results illustrate how a non-equlibrium mixture of quarks and gluons approaches hydrodynamic regime described by the Navier-Stokes equations with appropriate forms of the kinetic coefficients. The shear viscosity of a mixture is the sum of the shear viscosities of quark and gluon components, while the bulk viscosity is given by the formula known for a gas of quarks, however, with the thermodynamic variables characterising the mixture. Thus, we find that massless gluons contribute in a non-trivial way to the bulk viscosity of a mixture, provided quarks are massive. We further observe the hydrodynamization effect which takes place earlier in the shear sector than in the bulk one. The numerical studies of the ratio of the longitudinal and transverse pressures show, to a good approximation, that it depends on the ratio of the relaxation and proper times only. This behaviour is connected with the existence of an attractor solution for conformal systems.


I. INTRODUCTION
Comparisons between predictions of hydrodynamic models and exact kinetic-theory results have become an important method to verify the validity of hydrodynamic frameworks [1][2][3][4][5][6][7][8][9][10][11][12] which are now our basic tools to interpret the processes of heavy-ion collisions studied experimentally at RHIC and the LHC [13][14][15][16][17]. Such comparisons allow us also for deeper analyses of mutual relations between effective hydrodynamic models and microscopic, underlying theories [18][19][20][21], for a recent review see [22]. In this work we continue earlier studies on this topic and generalise previous results by studying a mixture of massive quarks and massless gluons forming a highly non-equilibrium system. Similarly to earlier works we restrict ourselves to boost-invariant systems [23].
Previous studies of mixtures [24][25][26][27] were restricted to the massless case and done mostly in the context of anisotropic hydrodynamics [28,29]. In this paper we restrict ourselves to the kinetic-theory study, leaving an anisotropic hydrodynamics context for a separate investigation. Nevertheless, we use here the results of the first-order Navier-Stokes hydrodynamics to demonstrate the process of hydrodynamization of the system [30]. We find that the hydrodynamization in the shear sector (equalisation of the longitudinal, P L , and transverse, P T , pressures) takes place earlier than the hydrodynamization in the bulk sector (equalisation of the average and equilibrium pressures).
In order to study the system behavior close to equilibrium we determine the shear and bulk viscosities of a mixture and find that the shear viscosity η is simply a sum of the quark and gluon shear viscosities, η = η Q + η G . On the other hand, the bulk viscosity of a mixture is given by the formula known for a massive quark gas, ζ. Nevertheles, we find that ζ depends on thermodynamic coefficients characterising the whole mixture rather than quarks alone, which means that massless gluons contribute in a non-trivial way to the bulk viscosity (provided the quarks are massive).
Interestingly, our studies of the time evolution of the ratio of the longitudinal and transverse pressures indicate that, to a very good approximation, it depends on the ratio of the relaxation and proper times only. This behaviour is related to the presence of an attractor which was found and discussed earlier for conformal systems [31][32][33][34][35][36], and, quite recently, also for non-conformal ones [37].
The paper is organised as follows: In Secs. II and III we introduce the system of kinetic equations for the quark-gluon mixture and study their momentum moments. This leads to two Landau matching conditions related to the baryon number and energy-momentum conservation. In Sec. IV we discuss an algebraic method useful for dealing with tensors describing our main observables. This method is used in Secs. V, VI, and VII to calculate various thermodynamic variables for systems exhibiting isotropic, anisotropic, and exact distribution functions, respectively. In Sec. VIII, which is central analytic part of this work, we discuss the conservation laws and present two integral equations used to determine the proper-time dependence of the effective temperature and baryon chemical potential. In Sec. IX we present our results describing proper-time dependence of various quantities, hydrodynamization process, and scaling properties of the P L /P T ratio. We summarise and conclude in Sec. X. Appendices A, B and C contain: details of the calculations of the generalised thermodynamic functions, discussion of the Navier-Stokes equations, and the explicit calculation of the shear and bulk viscosities for a quark-gluon mixture, respectively.
In this work we use x µ = (t, x, y, z) and p µ = (p 0 = E p , p x , p y , p z = p L ) to denote the particle space-time position and four-momentum. The longitudinal (z) direction corresponds to the beam axis. The transverse momentum is p T = p 2 x + p 2 y and particles are assumed to be always on the mass shell, E p = m 2 + p 2 T + p 2 L . The scalar product of two four-vectors is a µ b µ = a µ g µν b ν ≡ a · b where g µν = diag (1, −1, −1, −1) is the metric tensor. For the partial derivative we use the notation ∂ µ ≡ ∂/∂x µ . Throughout the paper we use natural units with c = k B =h = 1.
In numerical calculations we assume that τ eq is constant, which explicitly breaks conformal symmetry of the system. The other source of breaking of the conformal symmetry is a finite quark mass. The form of U µ (x) in (1) is defined by choosing the Landau hydrodynamic frame. We note, however, that for one-dimensional boost-invariant systems the structure of U µ (x) follows directly from the symmetry arguments, see Sec. VII A.
In Eq. (2) the functions f s,eq (x, p) are standard equilibrium distribution functions, which (unless specified otherwise) take the Fermi-Dirac and Bose-Einstein forms for (anti)quarks and gluons, respectively, Here T (x) is the effective temperature, µ(x) is the effective chemical potential of quarks, and The same value of T (x) appearing in Eqs. (3) and (4), as well as the same value of µ(x) appearing in the quark and antiquark distributions in Eq. (3) introduces interaction between quarks, antiquarks and gluons -all particles evolve toward the same local equilibrium defined by T (x) and µ(x). Since the baryon number of quarks is 1/3, we can use the relation with µ B being the baryon chemical potential.
All particles are assumed to be on the mass shell, p 2 = p · p = m 2 , so that the invariant momentum measure is where Θ is the Heaviside step function and t µ is an arbitrary time-like four-vector. Hereafter, the gluons are treated as massless, while quarks have a finite constant mass m.

III. MOMENTS OF THE KINETIC EQUATIONS
We introduce the n-th moment operator in the momentum spacê with the zeroth moment operator defined aŝ Acting withÎ µ 1 ···µn on the distribution functions f s (x, p) and multiplying them by the degeneracy factors k s , one obtains the n-th moments of the distribution functions Here k s ≡ g s /(2π) 3 , with g Q ± = 3 × 2 × N f and g G = 8 × 2 being the internal degeneracy factors for (anti)quarks and gluons, respectively. In our calculations we assume that we deal with two (up and down) quark flavors with equal mass, which reflects the SU(2) isospin symmetry.
With the above definitions, the first and second moments of the distribution functions which are identified with the particle number current and the energy-momentum tensor of the species "s", respectively. In addition, we define the baryon number current where q s = {1/3, −1/3, 0} is the baryon number for quarks, antiquarks, and gluons, respectively. The total particle number current and total energy-momentum tensor read We now consider the n-th moments of the kinetic equations (1), which are obtained by acting with the operatorÎ µ 1 ···µn given by (8) on their left-and right-hand sides and multiplying them by the degeneracy factors k s . The zeroth and first moments have the form which, using Eqs. (10)- (12), may be rewritten as Taking difference between s = Q + and s = Q − components of Eqs. (18) we obtain the baryon current evolution equation On the other hand, when taking the sum over "s" components of Eqs. (19) one gets the total energy and momentum conservation equation In order to have the baryon number conserved it is required that the left-hand side of Eq. (20) vanishes, ∂ µ B µ = 0. The latter implies vanishing of the right-hand side of Eq. (20), which leads to the Landau matching condition for baryon current Analogously, the energy and momentum conservation means that the left-hand side of Eq. (21) vanishes, ∂ µ T µν = 0. This condition results in vanishing of the right-hand side of Eq. (21), which leads to the Landau matching condition for energy and momentum

IV. TENSOR DECOMPOSITION
It is convenient to introduce the four-vector basis [42,43] which in the local rest frame (LRF) reads Using Eqs. (25) one may express the metric tensor as follows [44] The projector on the space orthogonal to the four-velocity, and satisfies the conditions U µ ∆ µν = 0, ∆ µ α ∆ αν = ∆ µν and ∆ µ µ = 3. The basis (25) is a unit one in the sense that and complete so that any four-vector may be decomposed in the basis A (α) . In particular, one may express the particle number flux as follows where the coefficients n s A , due to Eqs. (28), are given by the projections with A 2 = A · A (note that A 2 = −1 for space-like four-vectors of the basis (25)). The tensorial basis for the rank-two tensors is constructed using tensor products of the basis four-vectors A µ (α) . Thus the decomposition of the energy-momentum tensor takes the form with the components of T µν s (x) defined in the following way Using Eqs. (11) and (12) in Eqs. (30) and (32) in one gets

V. ISOTROPIC DISTRIBUTIONS
In the case of momentum-isotropic distribution functions (in particular, in the case of equilibrium distribution functions f s,eq (x, p) = f s,eq (p·U (x)), as defined by Eqs. (3) and (4)), which are invariant with respect to SO(3) rotations in the three-momentum space, by the symmetry of the integrands in Eqs. (33) and (34) one has n s,eq so that for the momentum-isotropic state Eqs. (29) and (31) have the following structure T µν s,eq (x) = E s,eq U µ U ν − P s,eq ∆ µν , with N s,eq = n s,eq U , E s,eq = t s,eq U U , P s,eq = t s,eq being the particle density, energy density, and pressure in equilibrium. Explicit forms of these expressions are given in App. A 2.

VI. ANISOTROPIC ROMATSCHKE-STRICKLAND DISTRIBUTIONS
It is also useful to consider anisotropic phase-space distributions introduced by Romatschke and Strickland in [45]. In the covariant form they read [24] f where is the quark transverse-momentum scale, and λ(x) is the non-equilibrium baryon chemical potential of quarks. Similarly, ξ G (x) is the gluon anisotropy parameter and Λ G (x) is the gluon transverse-momentum scale. The anisotropy parameters ξ s vary in the range −1 < ξ s < ∞, with the cases −1 < ξ s < 0, 0 < ξ s < ∞ and ξ s = 0 corresponding to the prolate, oblate and isotropic momentum distribution, respectively.
The distributions defined by Eqs. (40) and (41) are invariant only with respect to SO (2) rotations around the z direction in the three-momentum space. In this case one still has and Eqs. (11) and (12) have the following structure [46] N µ s,a (x) = N s,a U µ , with Here ∆ µν is the projection operator orthogonal to U and Z. Explicit forms of Eqs. (46) are given in App. A 1.

VII. EXACT SOLUTIONS OF THE KINETIC EQUATIONS
In order to solve Eqs. (22) and (23) we need to know the form of the distribution functions f s (x, p) being solutions of the kinetic equations (1). In general, such solutions are difficult to find and Eqs. (1) may be at best solved numerically. However, it is possible to find formal analytic solutions of Eqs. (1) in the case where the system is boost invariant and transversally homogeneous. Below, we discuss this case in more detail.

A. Boost-invariance and transversal homogeneity
Hereafter, we assume that the considered system is boost-invariant in the longitudinal (beam) direction and homogeneous in the transverse direction. In such a case we may choose [42] U µ = (t/τ, 0, 0, z/τ ), where τ is the (longitudinal) proper time As a result the system becomes effectively one-dimensional. Since its evolution is governed completely by the proper time that mixes t and z, one usually refers to such a system as (0+1)-dimensional.

B. Boost-invariant Bialas-Czyz variables
In the case of (0+1)-dimensional system exhibiting symmetries discussed in the previous section it is convenient to use the variables w and v which are defined as follows [47,48] Due to the fact that particles are on the mass shell w and v are related by the formula v(τ, w, p T ) = w 2 + (m 2 + p 2 T ) τ 2 .
Equations (52) and (53) can be inverted to express the energy and longitudinal momentum of a particle in terms of w and v, namely The Lorentz invariant momentum-integration measure can be written now as For boost-invariant systems, all scalar functions of space and time, such as the effective temperature T and quark chemical potential µ, may depend only on τ . In addition, one can check that the phase-space distribution functions (which are Lorentz scalars) may depend only on the variables w, τ and p T . We use these properties in the next section.

C. Formal solutions of the kinetic equations
With the help of the variables w, v and p T we can rewrite (1) in a simple form [1,2,49,50] where the boost-invariant versions of the equilibrium distribution functions are straightforward to find using (52) and (53) Below we assume that distribution functions f s (τ, w, p T ) are even functions of w, and depend only on the magnitude of p T 1 , The formal solutions of Eqs. (57) have the form [1,2,49, 50] where f 0 s (w, p T ) ≡ f s (τ 0 , w, p T ) is the initial distribution function (we have introduced here the notation τ eq = τ eq (τ ) for the general case where the equilibration time may depend on the proper time).

D. Damping function
In Eq. (61) we have introduced the damping function The function D(τ 2 , τ 1 ) satisfies the two differential relations and converges to unity if the two arguments are the same, D(τ, τ ) = 1. These properties imply the identity [26] For a constant relaxation time used in this work Eq. (62) reduces to

E. Initial distributions
In what follows we assume that the initial distributions f 0 s (w, p T ) are given by the anisotropic Romatschke-Strickland (RS) forms f 0 s,a (w, p T ) which follow from Eqs. (40) and (41), Here ξ 0 s ≡ ξ s (τ 0 ), Λ 0 s ≡ Λ s (τ 0 ), and λ 0 ≡ λ(τ 0 ) are initial parameters. In view of the form (61), the use of Eqs. (40) and (41) implies that the decomposition of the particle current and the energy-momentum tensor for (61) has the form of Eqs. (44) and (45), namely with Hereafter, we refer to results obtained with the solution (61) as the kinetic or exact ones.

A. Baryon number conservation
Using the expression for the baryon number current (13) and the decompositions (37) and (67) one may rewrite Eq. (22) as where we define the equilibrium and exact baryon number densities as The explicit formula for B eq (τ ) is derived in App. A 2, see Eq. (A33), becomes an integral equation Equation (73) is a single equation for two functions, T (τ ) and µ(τ ). The second necessary equation required for their determination is obtained from the Landau matching condition for the energy, which we discuss in the next section.
Meanwhile, it is interesting to notice that Eq. (73) can be rewritten as an integral equation for the function B(τ ), namely By differentiating (74) with respect to τ we get which is nothing else but the form of baryon number conservation law valid for the Bjorken geometry (in the original Bjorken paper [23] the same equation was obtained for the conserved entropy current). Equation (75) has scaling solution Combining (70) and (72) with (76) we find the equation which allows to determine µ in terms of T and τ for a given initial baryon number density.
Unfortunately, in the general case we study (Fermi-Dirac statistics for quarks) Eq. (77) is an implicit equation for µ. The situation simplifies in the case of classical statistics, where the function H B becomes independent of µ.

B. Four-momentum conservation
Using the expression for the energy-momentum tensor (15) and the decompositions (38) and (68) one may rewrite Eq. (23) as where E eq and E contain contributions from quarks, antiquarks, and gluons Using Eqs. (A26), (A30), (A36), and (A40) we obtain where the functionsH ± are defined by Eqs. (A9) and r is the ratio of the degeneracy factors Equations (73) and (81) are two integral equations that are sufficient to determine the proper-time dependence of the functions T (τ ) and µ(τ ). This is done usually by the iterative method [51]. The two initial, to large extent arbitrary, input functions T in (τ ) and µ in (τ ) are used on the right-hand sides of (73) and (81)  Our use of the two coupled integral equations is similar to the case studied previously in [52]. We find that it is more straightforward than using (81) together with (77). However, the situation is different in the case of classical statistics, where (77) can be used to determine analytically µ/T . In this case, the expression for µ/T obtained from (77) may be substituted into (81) and we are left with a single integral equation for the function T (τ ).
One may check, using (A18) and (A19), that Eq. (81) is consistent with the formula where P L = P Q + L + P Q − L + P G L is the total longitudinal momentum of the system. Equation (83) holds in general for the Bjorken expansion. It follows directly from the conservation law in the form ∂ µ T µν = 0.

IX. RESULTS
In this section we present the results of our numerical calculations. In all studied cases we use a constant equilibration time τ eq = 0.25 fm, which is the same for quark and gluon components. 2 The starting proper time is τ 0 = 0.1 fm and the evolution continues till τ f = 5.0 fm (or τ f = 10 fm in several cases). The initial transverse momentum scales of quarks and gluons are taken identical and always fixed to Λ 0 Q = Λ 0 G = 1 GeV. The initial non-equilibrium chemical potential λ 0 is chosen in such a way that the initial baryon number density is either B 0 = 0.001 fm −3 or B 0 = 1 fm −3 , see Eq. (A43).
Other initial conditions correspond to different values of the anisotropy parameters. We use three sets of the values for ξ 0 Q and ξ 0 G : i) ξ 0 Q = 1 and ξ 0 G = 10, ii) ξ 0 Q = −0.5 and ξ 0 G = 10, and iii) ξ 0 Q = −0.5 and ξ 0 G = −0.25. They correspond to oblate-oblate, prolate-oblate, and prolate-prolate initial momentum distributions of quarks and gluons, respectively. Such initial values for ξ 0 Q and ξ 0 G were used previously in Ref. [27]. We note that different values of ξ 0 Q , ξ 0 G , and λ 0 imply different initial energy and baryon number densities, hence, due to matching conditions, also different initial values of T 0 and µ 0 . We also note that the oblatedue to additional integral in Eq. (62). oblate initial configuration is supported by the microscopic calculations which suggest that the initial transverse pressure is much higher than the longitudinal one [30,53].
We perform our calculations for three different choices of the particle statistics and the quark mass: in the first case both quarks and gluons are described by the classical, Boltzmann statistics 3 and the quark mass is equal to 1 MeV 4 , in the second case we use again the classical statistics but the quark mass is 300 MeV, finally, in the third case the quarks are described by the Fermi-Dirac statistics and have the mass of 300 MeV, while the gluons are described by the Bose-Einstein statistics. The gluon mass is always set equal to zero. The case with classical statistics, B 0 = 0.001 fm −3 , and negligibly small quark mass of 1 MeV agrees well with the exact massless case studied in Ref. [27]. This agreement is used as one    Fig. 4.

Shear sector
The results shown in Fig. 4 suggest that the non-equilibrium dynamics of the system enters rather fast the hydrodynamic regime described by the NS equations (at the stage where deviations from local equilibrium are still substantial). Such a phenomenon was identified first in the context of AdS/CFT calculations [30] and is known now as the hydrodynamization process. To illustrate this behaviour in our case, we show in Fig. 5 see Eqs. B3. The effective shear viscosity (solid red line in Fig. 5) is compared with the standard shear viscosity coefficient, η, valid for the system close to equilibrium. For the quark-gluon mixture the latter is defined as the sum of the quark and gluon coefficients, 6 where following [55], see also [56][57][58], we use The coefficient η is calculated as a function of T and µ obtained either from the perfect-fluid (green dot-dashed line in Fig. 5) or NS hydrodynamic calculation (navy blue dashed line 5 We use the notation where calligraphic symbols such as E, P T or P L refer to exact values obtained from the kinetic theory. In the situations where the system is close to equilibrium and described by the Navier-Stokes hydrodynamics we add the subscript N S. The standard kinetic coefficients describe the systems close to equilibrium, hence, the shear viscosity is defined by the formula η = τ 2 (P T − P L ) NS and the bulk viscosity by ζ = −τ Π NS , see App. B. If we use the exact kinetic-theory values on the right-hand sides of these definitions we deal with effective values, which should agree with the standard definitions for systems being close to local equilibrium. 6 For general collision kernels, the total shear viscosity (although written formally as a sum of the individual contributions) may not be a simple sum of independent terms, for example, see [54].
in Fig. 5). In the two cases we find that η eff agrees very well with η for τ > 0.5 fm which is about two times the relaxation time. Thus, in the shear sector we observe a very fast approach to the hydrodynamic NS regime. It is important to notice that the agreement with the NS description is reached when η is significantly different from zero, which supports the idea that the hydrodynamic description becomes appropriate before the system thermalises, i.e., before the state of local thermal equilibrium with P T ≈ P L is reached.
In panel (a) in Fig. 6 Fig. 4. The results shown in Fig. 6 show that the shear viscosity of the mixture is dominated by the shear viscosity of quarks thoughout the system evolution. The information complementary to panel (a) is provided in the panel (b) in Fig. 6 where we present the ratio η Q /η G as a function of m/T and µ/T (contour lines) together with the system trajectory (colored line).

Bulk sector
Similarly to the shear-viscosity effects we can analyse the bulk sector, where we define the effective bulk viscosity by the expression where Π is the exact bulk pressure The time dependence of the effective bulk viscosity is compared in Fig. 7 with the time dependence of the bulk viscosity coefficient given by the expression where κ 1 and κ 2 are defined by the thermodynamic derivatives For a simple fluid with zero baryon density, the coefficient κ 1 becomes equal to the sound velocity squared. The steps leading to Eq. (90) are described in more detail in Appendix C.
The form of (90) agrees with that given in [55] for fermions. There is, however, one important difference between our approach and that of [55]. In Ref. Similarly as in the shear sector, we can see in Fig. 7 that ζ eff (τ ) approaches ζ(τ ), however, the agreement is reached for significantly larger times, τ > 2 fm. This means that the hydrodynamization of the bulk sector is slower and follows the hydrodynamization of the shear sector. Observations that the hydrodynamization in the shear sector may happen before the hydrodynamization in the bulk sector have been done recently in Ref. [59] within the non-conformal models using the gauge/gravity correspondence, where the hydrodynamization in the bulk sector has been dubbed the EoSization process. In this scenario first P L and P T tend to a common valueP = P eq and, subsequently,P approaches P eq , which signals establishing equation of state of the system.
To visualize the importance of the gluon degrees of freedom in expressions (91) for the bulk viscosity of the mixture in Fig. 8 we show the bulk viscosity coefficient ζ and compare it with the coefficient ζ 0 that has been calculated in the same way as ζ except that the thermodynamic coefficients κ 1 and κ 2 of the former were calculated only for the quark component. We find that neglecting the gluon contribution in κ 1 and κ 2 changes substantially the values of ζ making it significantly smaller. This finding indicates that gluons, although, massless, contribute to the bulk viscosity of a quark-gluon mixture. The necessary requirement for this effect is, however, that quarks are massive. 3. P T /E and P L /P T ratios Figures 9, 10, and 11 correspond to Figs. 1, 2, and 3, respectively, and show the time dependence of the ratios P T /E (upper panels) and P L /P T (lower panels). In the case of quarks with a very small mass (green dot-dashed lines) the ratios P T /E tend to 1/3 as expected for massless systems approaching local equilibrium. The ratios P L /P T in all studied cases tend to unity which again reflects equilibration of the system. Interestingly, the ratios P L /P T very weakly depend on the quark mass and the choice of the statistics. the results for the oblate-oblate, prolate-oblate, and prolate-prolate initial conditions. The blue dotted line describes (P L /P T ) NS obtained from the Navier-Stokes hydrodynamics.

C. Scaling properties
Each panel of Figs. 9, 10, and 11 shows our results obtained for different values of the quark mass and particle statistics but for the same initial anisotropies. In Figs. 12, 13, and 14 we rearrange this information showing in each panel our results obtained for different initial anisotropies, i.e., for oblate-oblate, prolate-oblate, and prolate-prolate initial quark and gluon distributions. Figures 12, 13, and 14 collect the results for different mass and statistics. The most striking feature of our results presented in these figures is that the P L /P T ratios (shown in lower panels) converge to the same values, although they describe the system evolutions starting from completely different initial conditions. The origin of this behaviour can be found if we analyse the NS formula for the P L /P T ratio. Let us first consider the massless case where we may neglect the bulk viscosity and write P L P T NS = P Q,eq − 4η Q /(3τ ) + P G,eq − 4η G /(3τ ) P Q,eq + 2η Q /(3τ ) + P G,eq + 2η G /(3τ ) .
Assuming in addition that the baryon number density is zero, we may use the following relations connecting the shear viscosity with equilibrium pressure 7 : It is interesting to note that the coefficient 4/5 is the same for quarks and gluons, hence  temperature, Eq. (94) indicates that (P L /P T ) NS depends on the product of τ and T , which is expected for conformal systems and related to the existence of a hydrodynamic attractor for such systems [31][32][33][34][35][36]. It turns out that the inclusion of the finite mass and baryon chemical potential (with the values studied in this work) affects very little Eqs. (93) connecting the shear viscosity with pressure. The main difference is that the coefficient 4/5 is slightly changed. It should be replaced by an effective value obtained for the studied range of T and µ.
To analyse the (P L /P T ) NS ratio in a general case in Fig. 15 we plot it as a function of two variables, τ eq /τ and m/T , for a fixed value of µ. The left panel of Fig. 15 shows the contour plot of (P L /P T ) NS in the case where quantum statistics are used and µ = 0. The fact that the contour (red dashed) lines have horizontal shapes indicates that (P L /P T ) NS depends effectively only on τ eq /τ (except for the region where τ ≈ τ eq and T ≈ m/5). The the studied, rather broad range of τ eq /τ and m/T . These observations explain similarities of the close-to-equilibrium behavior of (P L /P T ) NS in the left panels of Figs. 9, 10, and 11.
The right panel of Fig. 15 shows the contour plot of (P L /P T ) NS for µ/T = 2. In this case we find again a weak dependence on m/T as compared to the case of classical statistics and µ = 0 represented by the solid black lines. Again this helps to understand the similarities of the right and left panels of Figs. 9, 10, and 11.

D. Remarks on non-conformal attractor
In a very recent paper [37] it has been suggested by Romatschke to look for attractor behaviour by studying the quantities  and as functions of the variable (97) Note that in Eqs. (95)-(97) we used boost invariance to simplify our notation.
In Fig. 16 we show the function A 1 (Γ) obtained for three different initial anisotropies studied in this work. To get the connection with [37] we consider the case with negligible baryon number density. Otherwise, we include the finite mass of quarks and quantum statistics. Figure 16 shows that the lines corresponding to three different initial conditions converge and later approach the Navier-Stokes line. This observation supports the existence of a non-conformal attractor for A 1 in our system.

Anisotropic distributions
The forms of the generalised thermodynamic functions for anisotropic distributions are given by the following integrals: P s,a L ≡ t s,a ZZ = k s dP (p · Z) 2 f s,a p · U, p · Z .
With b = 0 the functions H 2 (a, b) reduce to the functions H(a), H L (a) and H T (a) used in [2].
The integrals in (A10) are analytic [4]: For gluons one has: where ζ is the Riemann zeta function (the coefficient ζ(3) is known as Apéry's constant).
The expressions on the right-hand sides of Eqs. It is useful to notice that the functions H 2 and H 2L are related by the expression hence, we also have ∂H ± (a, y, z) ∂a =H ± (a, y, z) +H ± (a, y, z) a . (A19) We can use (A19) to derive (83) from (81).
We close this section with the formula for the baryon number density valid for anisotropic RS systems where H B (y, z) ≡ 1 4 ∞ 0 r 2 dr 1 cosh r 2 + y 2 + cosh z . (A21)

Isotropic distributions
The forms of the thermodynamic functions for the equilibrium state are commonly known, nevertheless, we quote them here for completeness. They are given by the formulas N s,eq ≡ n s,eq Their explicit forms for quarks and anti-quarks may be obtained from Eqs. (A5)-(A8) as a special case of ξ s → 0, Λ s → T , and λ s → µ, Note that H 2L (1, b) = 2/ 3 √ 1 + b 2 and H 2T (1, b) = 4/ 3 √ 1 + b 2 which means that P Q ± ,eq T = P Q ± ,eq L ≡ P Q ± ,eq , as expected for the isotropic state. Analogous results may be obtained for gluons where to get the last expressions on the right-hand sides we again assumed the Bose-Einstein statistics. Here, similarly as for quarks P G,eq T = P G,eq L ≡ P G,eq . Similarly to anisotropic case the baryon number density is Using above definitions and repeating the calculation from Sec. A 1, gives for quarks and for gluons. We define the baryon number density for the exact solution of the kinetic equation as follows

Appendix B: Navier-Stokes hydrodynamics
The results of our kinetic-theory calculations are compared with the viscous hydrodynamic results obtained by solving the Navier-Stokes (NS) hydrodynamic equations. The latter have the form d dτ E Q,eq + E G,eq = − E Q,eq + E G,eq + P Q,eq + P G,eq + Π NS − π NS τ , dB eq dτ + B eq τ = 0.
Here E Q,eq = E Q + ,eq + E Q − ,eq is the equilibrium energy density of quarks and antiquarks, P Q,eq = P Q + ,eq + P Q − ,eq is the equilibrium pressure of quarks and antiquarks, Π NS is the bulk pressure, and π NS is the shear pressure (both used in the close-to-equilibrium limit).
All the functions appearing in (B1) and (B2) depend on T and µ, hence Eqs. (B1) and (B2) are two coupled equations that can be used to determine T (τ ) and µ(τ ). One can easily notice that Eq. (B1) may be written in the form of Eq. (83) once we identify P L = P eq − π NS + Π NS , (B3) where P eq = P Q,eq + P G,eq . Within NS approach, the shear and bulk pressures are expressed by the kinetic coefficients η and ζ, The expressions for η Q , η G , and ζ are given by Eqs. (86), (87), and (90).
For the moment let us denote T (τ ) and µ(τ ) obtained from the kinetic theory as T KT (τ ) and µ KT (τ ), while those obtained from the NS hydrodynamics as T NS (τ ) and µ NS (τ ). We expect that T KT (τ ) and µ KT (τ ) agree well with T NS (τ ) and µ NS (τ ) in the late stages of the evolution, when the system approaches local equilibrium. To check this behaviour we choose such initial conditions for hydrodynamic equations (B1) and (B2) that for the final time τ = τ f we match the temperature and chemical potential in the two approaches: T NS (τ f ) = T KT (τ f ), µ NS (τ f ) = µ KT (τ f ). Then, we check if the functions T NS (τ ) and µ NS (τ ) smoothly approach T KT (τ ) and µ KT (τ ) if τ → τ f . By neglecting the bulk and shear pressures in (B1) we can also make comparison with perfect fluid hydrodynamics and check if the system approaches local equilibrium.