Onset of hydrodynamics for a quark-gluon plasma from the evolution of moments of distribution functions

The pre-equilibrium evolution of a quark-gluon plasma produced in a heavy-ion collision is studied in the framework of kinetic theory. We discuss the approach to local thermal equilibrium, and the onset of hydrodynamics, in terms of a particular set of moments of the distribution function. These moments quantify the momentum anisotropies to a finer degree than the commonly used ratio of longitudinal to transverse pressures. They are found to be in direct correspondence with viscous corrections of hydrodynamics, and provide therefore an alternative measure of these corrections in terms of the distortion of the momentum distribution. As an application, we study the evolution of these moments by solving the Boltzmann equation for a boost invariant expanding system, first analytically in the relaxation time approximation, and then numerically for a quark-gluon plasma with a collision kernel given by leading order 2 ↔ 2 QCD matrix elements in the small angle approximation.


Introduction
The evolution of the quark-gluon plasma (QGP) produced in high energy heavy-ion collisions is well described by the equations of relativistic hydrodynamics including viscous corrections (see [1] for a recent review). The success of such hydrodynamic descriptions suggests that the system of quarks and gluons which emerge shortly after the collisions is brought close to local equilibrium on a relatively short time scale, a process commonly referred to as thermalization. Understanding the detailed mechanisms by which such a thermalization occurs remains a theoretical challenge. At very early times, one may argue that the dynamics is dominated by classical color fields and the system evolution is governed by the classical Yang-Mills equations [2,3]. At a time scale of order 1/Q s , where Q s is the saturation momentum, the system may be sufficiently dilute for kinetic theory to become applicable [4]. Kinetic theory naturally fills the gap between the dynamics of classical field and that of dissipative fluids, and it offers the possibility to follow in details how the pre-equilibrium system evolves into a state of quasi local equilibrium well accounted for by viscous hydrodynamics. This is the framework that we shall consider in this work (see e.g. [5][6][7][8][9] for some recent representative works, and more specifically [10,11] for the analysis of the onset of hydrodynamics in a weakly coupled system using kinietic theory).
The transition from kinetic theory to hydrodynamics is commonly achieved by taking suitable moments of the kinetic equations, with lower moments encompassing conservation laws and higher moments various dissipative effects. For some particular geometries, it is JHEP11(2017)161 possible to recast the solution of the Boltzmann equation in terms of an infinite hierarchy of equations for a particular set of moments of the distribution function (see e.g. [12]). Provided this infinite hierarchy can be limited to the first few moments, this technique could represent a convenient alternative to the direct solution of the kinetic equation. Aside from this aspect, there is another interest in using moments. Defined in terms of integrals of the phase-space distribution function with suitable weights, the moments allow us to focus on the relevant (typically long wavelength) information, and wash out the irrelevant (short wavelength) one from the distribution function, thereby automatically implementing a strategy akin to that of effective field theories.
In this paper, we introduce a specific set of moments, defined as weighted integrals of the momentum distribution function f (p), L n ∝ p p 2 P 2n (cos θ)f (p), where P 2n is a Legendre polynomial, and cos θ = p z /p with p z the longitudinal momentum of a particle with total momentum p. These moments capture the deviation of the momentum distribution from an isotropic distribution and are therefore damped as the system approaches local equilibrium. They are tailored to take into account the (strong) effect of the longitudinal expansion. In addition, these moments are found to be in correspondence with viscous corrections in hydrodynamics, order by order in a gradient expansion. They can therefore provide an alternative view of the viscous corrections, in terms of the damping of high multipoles of the momentum distribution as the system approaches local equilibrium. Their relative magnitudes can be used as an indicator of the onset of hydrodynamics. In this work, we shall provide a description of the evolution of these moments, obtained by solving the Boltzmann equation for a longitudinally expanding system. The usefulness of these moments as a practical tool in solving the kinetic equation will be addressed in a separate publication.
This paper is organized as follows. The moments L n , and their usefulness in a boost invariant setting, are introduced in section 2. In this section, we also present some of the major results obtained in this work, concerning the correspondence between the L n 's and the transport coefficients used in second order viscous conformal hydrodynamics [13]. Relations to higher order viscous corrections are also derived. Then, in section 3, we discuss the evolution of the moments in the pre-equilibrium stage towards local equilibrium using kinetic theory. To do so we solve the Boltzmann equation in a boost invariant setting, using two different approximations for the collision kernel. In the first case, we use the relaxation time approximation, with a constant relaxation time, and obtain analytical solutions that allow us to verify the general relations to viscous hydrodynamics established in section 2. Then, in section 3.2, we consider a more realistic setting: we use leading order 2 ↔ 2 QCD matrix elements for the collision kernel, and solve the corresponding Boltzmann equation in the small scattering angle approximation in order to study the evolution of the moments. Summary and conclusions are given in section 4.

Formulation of moments in expanding systems
In this section, after a brief review of the main features of expanding boost invariant systems, we define a set of moments of the distribution function that are suited to the JHEP11(2017)161 description of such systems. We show that, near the hydrodynamic regime, these moments are in correspondence to the viscous corrections that emerge from a gradient expansion.

Kinetic theory in a boost invariant expanding system
The quark-gluon plasma produced in the very early stages of a heavy ion collision experiences a strong expansion along the collision axis (referred to as the longitudinal direction). A simple description of the system (Bjorken flow) is obtained when one assumes boost invariance in the longitudinal direction and translational invariance in the transverse directions [14]. It becomes then convenient to use in place of the usual space-time coordinates, the proper time τ = √ t 2 − z 2 and the space-time rapidity tan −1 (z/t) , where z is the longitudinal coordinate and t the time. In fact, boost invariance makes it possible to focus on a slice of the fluid located around the plane z = 0, where τ = t. There, the distribution function depends only on time and the three components of the momentum, f (t, p T , p z ), and it obeys the kinetic equation with C the collision term. 1 Averages of various physical quantities with the phase space distribution function play an important role in this paper, and we shall denote them with double brackets 2 For instance, the energy density is given by e(t) = (p 0 ) 2 = p 2 , where we have used p 0 = p, assuming massless particles. By multiplying eq. (2.1) by p 2 , and integrating over momentum, we obtain, where we have used the fact that the collisions conserve energy, so that the contribution of the collision term vanishes. The quantity P L is the longitudinal pressure, We define similarly the transverse pressure where in the last equality, we have used e(t) = P L + 2P T and eq. (2.3).
1 For simplicity, we consider in this section a single-species system. In section 3.2 we shall introduce distinct distributions for quarks and gluons and take proper care of degeneracy factors. 2 Throughout this paper, we use bold lower case letters to denote a three vector, such as p, with the magnitude of the vector denoted by the corresponding normal lower case letter, e.g., |p| = p.

JHEP11(2017)161
In the hydrodynamic regime, i.e. when local equilibrium is achieved, the pressure is isotropic and simply related to the energy density, P = P L = P T = e/3, and eq. (2.3) becomes a closed equation for the energy density. It yields e(t) ∼ 1/t 4/3 ∼ T 4 , with T (t) the local temperature. Before reaching this regime, viscous corrections need to be taken into account. The first order correction involves the shear viscosity η, and yields a relaxation equation for the difference of the pressures Viscous corrections are accompanied by an entropy increase, given in leading order by the equation (with T s = e + P, where s is the entropy density and T the temperature) Note that ts represents the total entropy in a (expanding) covolume, or equivalently the entropy density in the transverse plane. In the absence of viscosity, this is contant. Our purpose in this paper is to explore specific features of the approach to the hydrodynamical regime, studying in particular how the deviations of the momentum distribution from an isotropic distribution are damped.

Thermalization, isotropization and moments
Quite generally, the effect of collisions is to wash out the anisotropy of the momentum distribution, leading eventually to a spherically symmetric distribution. In the case of expanding boost invariant systems, this effect is counterbalanced by the strong longitudinal expansion. The competition between the two effects is commonly investigated in terms of the ratio between the longitudinal pressure P L and the transverse pressure P T , local equilibrium being established when P L /P T = 1. As an illustration, and anticipating on the discussion of section 3.2, we show in figure 1 the typical evolution of P L /P T obtained from the numerical solution of the Boltzmann equation, with leading order 2 ↔ 2 scatterings of QCD matrix elements, in the small scattering angle approximation. The trend towards isotropization is clearly visible. However, this is a slow process, and complete local thermal equilibrium is not reached in the time span of the simulation, with quark production delaying the approach to equilibrium even further as compared to the case of a purely gluonic plasma. 3 Such a slow approach to isotropy of the momentum distribution seem to be quite generic for a longitudinally expanding system undergoing Bjorken flow [11]. However, it has also been realized that the complete isotropy of the pressures, characterized by P L /P T = 1, may not be necessary for hydrodynamics to be applicable, since viscous corrections can accommodate rather large differences between P L and P T .
3 Some readers may find it surprising to see this trend towards isotropization arise from elastic scattering alone, as it seems to go against a prediction from the bottom-up scenario [15] emphasizing the role of inelastic collisions in isotropization. The parametric analysis of the bottom-up scenario relies however on values of the coupling constant much smaller than those used in the present simulation (αs ≈ 0.3). In fact, Tanji and Venugopalan [16] have recently completed a systematic study of the same equations as ours, in which they show explicitly how deviations form the bottom-up scenario develop as the coupling increases. This was first revealed by strong coupling calculations [17], where it was shown that viscous hydrodynamics can handle pressure ratios P L /P T 0.5. This has been later verified in a number of calculations (see ref. [11] and references therein).
In this paper, in order to describe more details of the isotropization of an expanding quark-gluon plasma, we introduce the following moments of the distribution function L n ≡ p 2 P 2n (cos θ) ≡ p p 2 P 2n (cos θ)f (p), (2.8) where P 2n is a Legendre polynomial of order 2n, and cos θ = p z /p. For an expanding system with Bjorken geometry, odd order moments vanish as a consequence of the invariance of the distribution function under parity. There are two distinct advantages of these moments. First, except for the n = 0 moment which corresponds to the energy density, as already mentioned, all higher order moments defined in eq. (2.8) naturally quantify the details of the longitudinal momentum anisotropy. For instance, the information concerning P L /P T is contained in the n = 1 moment, More precisely, L 1 → 0 is equivalent to P L /P T → 1. Similarly, the moments L n of higher order are associated to finer structures of the momentum anisotropy of the distribution function. Second, the moments defined in eq. (2.8), with the specified weight p 2 , are closely related to hydrodynamics through Landau's matching condition

JHEP11(2017)161
of which the moment L 1 provides the simplest illustration (see eq. (2.6)). These moments are therefore expected to acquire a physically transparent meaning at late times when the system approaches the hydrodynamic regime. Note that one may generalize the definition in eq. (2.8) to space-time symmetries other than the boost invariant setup considered in this paper. For a system with spatial SO(3) rotational symmetry, for instance, one could replace the Legendre polynomials with associated Legendre polynomials to further characterize momentum anisotropies in the azimuthal direction.

The moments L n in the hydrodynamic regime
We now analyze more precisely the correspondence of the moments (2.8) to viscous hydrodynamics. At this point covariant notation is helpful. Four vectors are denoted by normal upper case letters. Thus, for instance, U µ is the fluid four-velocity, normalized so that U 2 = −1. The operator ∆ µν = g µν + U µ U ν , with g µν the metric tensor, 4 projects on the space orthogonal to U µ . The energy momentum tensor is written as where π µν represents the viscous corrections. In leading order, π µν = −ησ µν , where η is the shear viscosity, and Here, and in the following, single brackets . . . around a tensor implies that the tensor has been made symmetric, traceless and transverse to U µ . For simplicity, we assume throughout this work that the fluid is made of massless constituents, so that conformal symmetry implies that the energy-momentum tensor is traceless, T µ µ = 0. Hydrodynamics effectively applies to systems whose evolution is dominated by long wavelength modes. Corrections to ideal hydrodynamics are then naturally searched for in a gradient expansion. The successive terms in such an expansion can be obtained by applying the Chapman-Enskog technique.
As an illustration of the method, let us repeat the derivation of eq. (2.6) within the relaxation time approximation for the collision term. Assuming that the viscous correction corresponds to a small deviation δf of the phase-space distribution function from its local equilibrium value f eq , we linearize the Boltzmann equation and obtain, in covariant notation, where τ rel is the relaxation time. Note that P · U reduces to −P 0 in the local rest frame of a fluid cell. We shall allow the relaxation time to depend on the energy of the particle, i.e, we assume τ rel = τ rel (P · U/T ), where T is the local temperature (which enters the local equilibrium distribution). Various ansatz have been considered in the literature, in particular a so-called "linear" ansatz, corresponding to a constant relaxation time, and a "quadratic" ansatz, corresponding to a linear dependence of τ rel on the energy [20]. The

JHEP11(2017)161
origin of this terminology is that, in the first case, δf /f eq ∼ p/T , while δf /f eq ∼ (p/T ) 2 when τ rel is linear in p. Note that the present analysis does not rely on the specific dependence of the relaxation time upon energy. To proceed, we find it convenient to definẽ (2.14) Since the local equilibrium distribution is a function of P · U/T , the effect of the operator P µ ∂ µ when acting on f eq is to generate the structure P µ P ν σ µν , i.e., where the prime on f eq indicates a derivative with respect to P ·U/T . Ignoring the derivative of δf in the left hand side of eq. (2.14), one then identifies the first order viscous correction to the phase-space distribution function as well as the first viscous correction to the energy momentum tensor, A simple calculation then allows us to determine the shear viscosity [18] In the case of Bjorken flow, the contraction of the irreducible tensors in eq. (2.16) is easily calculated in terms of the Legendre polynomials. One gets By using this expression in eq. (2.16), and the expression (2.18) of the shear viscosity, one can then calculate the n = 1 moment, eq. (2.9), and obtain eq. (2.6). The higher order corrections are obtained iteratively along the same line. In the second order correction δf (2) needed to the calculation of the transport coefficients of conformal viscous hydrodynamics [13], the following new tensor structures appear [18,19],

JHEP11(2017)161
For Bjorken flow, the vorticity tensor Ω µν = 1 2 [∇ µ U ν − ∇ ν U µ ] vanishes. One can also prove that contractions of irreducible tensors of odd ranks do not contribute due to parity. Finally, the relevant structures in eq. (2.20) are those arising from the contraction of irreducible tensors of even ranks, P µ P ν P α P β σ µν σ αβ = 32 35 One can thus rewrite the viscous corrections to the phase space distribution up to second order in the gradient expansion in terms of Legendre polynomials, where, for convenience, we have defined the dimensionless variablesχ p = f eqC p andp = p/T . Ellipses in the brackets of eq. (2.22) stand for terms of higher power in 1/tT . One recognizes in eq. (2.22) a gradient expansion, with δf expressed in powers of 1/tT . The latter quantity may be viewed as a measure of the Knudsen number, that is, as the ratio between a microscopic and a macroscopic length scale characterizing the fluid. Here the typical microscopic length scale is the inverse of the temperature, while the macroscopic length scale can be taken as the inverse of the local expansion rate, which, for a medium system expanding according to Bjorken flow, is simply the time t. The expansion (2.22) is also an expansion in Legendre Polynomials, the term multiplying P 2n (cos θ) having a gradient expansion starting with a leading contribution of order n, that is, ∼ 1/(tT ) n . 5 More explicitly, the evolution of the moments of order n = 1 and n = 2 can be expressed in terms of the transport coefficients that enter (conformal) viscous hydrodynamics: and where τ π , and λ 1 are second order transport coefficients [13]. In obtaining the above results, we have used the following expressions of these transport coefficients in kinetic JHEP11(2017)161 theory (generalizing eq. (2.18) for the shear viscosity) Again, we emphasize that the momentum dependence of the relaxation time only affects the values of these transport coefficients, as given by eqs. (2.25), but it does not alter the form of relations such as those given in eq. (2.23) and eq. (2.24). Viscous corrections in hydrodynamics get more complicated in higher orders, which makes it more involved to derive an explicit correspondance between transport coefficients and the p 2 -moments associated with higher order Legendre polynomials. Nevertheless, it is possible to generalize the analyses in eqs. (2.23) and (2.24) to higher orders, at least for the leading terms. Indeed, it can be shown that, after linearizing the Boltzmann equation in order to construct the gradient expansion, there is a term with highest power in momentum which appears on the left hand side of the Boltzmann equation as a result of applying P µ ∂ µ to f eq (P · U/T ) iteratively. This highest power appears also as the leading contribution to the coefficient of the Legendre polynomial P 2n (cos θ) in δf , through contraction of irreducible tensors, Therefore, one has The coefficient c n is expected to be identified to some combinations of transport coefficients of order n, since 1/t n represents an n-th order contribution in the gradient expansion with a definite angular structure. For a conformal and classical gas, c n can be analytically determined when the relaxation time is linear in p/T (quadratic ansatz), One may regard eq. (2.28) as an analytical prediction of the n-th order transport coefficient of a conformal system, with η/s a constant input. In particular, one easily verifies that c 1 = −2η, as expected. More details regarding eq. (2.26) and eq. (2.28) are given in appendix A.

Evolution of the L n moments in expanding systems
We now demonstrate how the L n moments evolve in the pre-equilibrium stage of heavy ion collisions by calculating them using kinetic theory. To do so, we solve the Boltzamnn JHEP11(2017)161 equation (2.1) for a boost invariant system. In line with the color-glass picture (CGC) [21], we take an initial momentum distribution function of the form where f 0 is a free parameter characterizing the typical initial gluon occupation number, and Q s is the saturation momentum. We assume that the kinetic description applies at an initial time t 0 ∼ 1/Q s . In this work, we choose a value f 0 = 0.1, for which the approach to equilibrium occurs smoothly, without encountering Bose-Einstein condensation [7,8]. 6 Especially, with this small initial occupancy, variation of particle number would have negligible effect regarding the 2 ↔ 2 QCD matrix elements within a small angle apprximation. The parameter ξ controls the initial momentum anisotropy. For ξ > 1, one has initially P L /P T < 1. Throughout this work, we take for definiteness a fixed value, ξ = 1.5, corresponding to an initial momentum anisotropy P L /P T ≈ 0.5. One may check that, initially, the moments calculated form the momentum distribution (3.1) are ordered such that |L n+1 | < |L n |. Furthermore, all odd order moments are negative, while even order moments are positive. We shall present the results of two calculations. We start with the simple case where the collision term is written in the relaxation time approximation. For a constant value of the relaxation time (linear ansatz), an analytical solution to the Boltzmann equation can be obtained [24]. For a somewhat more realistic analysis, we then proceed to the numerical solution of the Boltzmann equations for quarks and gluons, considering 2-to-2 scatterings among gluons and quarks within a small scattering angle approximation [8].

Relaxation time approximation
The Boltzmann equation (2.1), with the relaxation time approximation for the collision term, reads We only require here energy conservation, 7 so that the local equilibrium distribution function f eq is a Bose-Einstein distribution with a vanishing chemical potential, For a constant τ rel , the solution to eq. (3.2) can be written formally as, Including the inelastic processes would change the behavior of the distribution function at very small momenta, as discussed in [22], and analyzed in detail recently in [23]. However, because of the weighting factor p 2 in their definition, the moments Ln are not sensitive to this particular region, and we do not expect that including inelastic, number changing, processes would alter the main conclusions of the present paper. 7 If one would also require conservation of the number of constituents, the local equilibrium distribution would also depend on a chemical potential [25]. This situation will be considered in the next subsection. We recognize in the first term in the right hand side of eq. (3.4) a contribution that represents free-streaming from the initial condition f (t 0 , p T , p z ), eq. (3.1). The time dependence of the temperature T in f eq in the second term is fixed by the condition of energy conservation This relation, together with eq. (3.4), completely determines the solution. The resulting energy density exhibits the expected transition from the early free streaming regime, where e(t) ∼ 1/t (see below), to the hydrodynamic regime at late times where e(t) ∼ 1/t 4/3 . The evolutions with time of the energy density and the pressures obtained from eq. (3.4) are illustrated in figure 2, and compared to the solution of first order viscous hydrodynamics (which we also refer to as Navier-Stokes hydrodynamics). As shown by this figure, the energy density is well accounted for by viscous hydrodynamics for times t 15τ rel . This is also the case for the pressure, although in this case, the existence of significant viscous corrections at the latest times is attested by the fact that the longitudinal pressure is not yet equal to the transverse pressure, P L /P T 0.9 at t = 50 Q −1 s . The moments can be calculated from the distribution function given in eq. (3.4). Consider first the free-streaming regime (τ rel → ∞). In this case, one has

JHEP11(2017)161
where F n (x) is a function defined from the following integral (0 ≤ x ≤ 1) This function has the following limits: F n (x) → πP 2n (0)/2 as x → 0, and F n =0 (x) → 0 as x → 1. Thus, for asymptotically large t, t t 0 /ξ, F n (x) reduces to a constant, and the moments decay as 1/t. When t = t 0 /ξ, the moments with n = 0 vanish, which implies in particular that they vanish at t = t 0 if there is no initial momentum anisotropy (ξ = 1). The energy density is given by the zeroth moment, with F 0 (0) = π/2 and F 0 (1) = 2. To within the slowly varying function F 0 (t 0 /ξt), the energy density exhibits the expected behavior in 1/t. It can also be verified that the longitudinal pressure drops rapidly, as ∼ 1/t 2 , so that at times t t 0 /ξ, the distribution function is peaked around p z = 0, and the energy density is dominated by transverse degrees of freedom.
For a finite τ rel , eq. (3.4) leads to where ζ(n) is the Riemann-zeta function. In this equation, the first term represents the contribution of the free-streaming of the initial distribution. This is suppressed in a time scale τ rel , i.e., when collisions start to play a significant role. One thus expects the evolution of the moments to exhibit a transition between the free-streaming regime at short time, t τ rel , and the late time regime, dominated by collisions and represented by the second term in eq. (3.8). Figure 3 displays the evolution of the absolute values of the normalized moments |L n /L 0 | up to n = 4 (recall that L 1 and L 3 are negative). Also shown are the moments of the pure free-streaming solution, eq. (3.6), which saturate at late times to their corresponding asymptotic values determined by F n (0), namely, Note that this limit is independent of the initial pressure anisotropy paramater ξ. Clearly, the isotropization of the momentum distribution, signaled by the vanishing of the moments, can only be achieved by the collisions. Indeed the effect of the collisions starts to be visible around the time scale (t − t 0 ) ∼ τ rel where they begin to compete with the free streaming and later drive the momentum distribution to isotropy. Another effect of the collisions is to reduce the overall magnitude of the higher moments. As can be seen on figure 3, while the free streaming moments evolve towards comparable (within a factor 2) values at late times, the hierarchy of moments present in the initial condition is preserved but after a time t 15 τ rel only the moment L 1 remains significant. At this time the evolution of the system is well accounted for by viscous hydrodynamics.  This is confirmed by a more detailed study of the late-time evolution of the moments in eq. (3.8). First we consider the moment L 1 . According to eq. (2.23), the ratio between this moment and the entropy density, more precisely −tL 1 /2s, approaches η/s in the hydrodynamic regime. This ratio is plotted in figure 4, where we purposely rescaled the time in figure 4 by 1/τ rel so that they evolve on the same time scale. As time increases, −tL 1 /2s indeed approaches the kinetic theory expectation: η/s = 1 5 T τ rel , and reaches it for t 15 τ rel . Figure 5 completes the discussion and illustrates the system evolution in terms of higher order moments. Again, we consider specified dimensionless combinations of moments which are supposed to reduce to dimensionless ratios of higher order transport coefficients at late times. The following combination plotted in figure 5(a) is related to second order transport coefficients λ 1 and τ π . In a similar way, the combinations  drodynamics 9 it is nevertheless interesting to estimate their ratios from the asymptotic behaviors: these are the black arrows in figure 5. Note that the asymptotic value of (λ 1 + ητ π )/η 2 /(e + P) ≈ 8.57 in figure 5(a) is consistent with what one expects from kinetic theory with a linear ansatz for the relaxation time [26]. The asymptotic values of the ratios in figure 5(b) and 5(c) are found to be approximately 1.3 and 1.0, respectively.
Actually, the combinations plotted in figure 5 are nothing but double ratios among three consecutive orders of transport coefficients (to within simple numerical constants). In terms of the c n 's defined in eq. (2.27), these are Asymptotic values of these double ratios can be obtained analytically in the case of the quadratic ansatz for the relaxation time (eq. (2.28)), or numerically for the linear ansatz. Results for these two cases are shown in figure 6 for n ≤ 20 and n ≤ 5, respectively, 10 9 Note however that the third order coefficients could in principle be extracted from the works on third order viscous hydrodynamics in refs. [27,28]. We thank A. Jaiswal for alerting us about this. 10 For the linear ansatz we limit ourselves to n ≤ 5, since the evaluation of the double ratio is challenged by precision in numerical integrations at asymptotic large t in eq. (3.8), which become less stable for higher orders.  Figure 5. Evolution of higher order moments, plotted in specified combinations which correspond to dimensionless ratios of transport coefficients in hydrodynamics regime. Black arrows indicate asymptotic values. Note that for (λ 1 + ητ π )/(η 2 /(e + P)), the asymptotic value is consistent with what was expected from kinetic theory [26].
as filled points and open symbols. For large n, the double ratio is seen to approach unity from above, which implies a saturation of transport coefficients of asymptotically high orders. This behavior reflects the fact that the gradient expansion leading to viscous hydrodynamics is only asymptotic [29,30] (see also refs. [31,32] for a recent discussion of this issue within the context of Bjorken expansion).

Quark-gluon system with small angle approximation
We now apply the analysis of the previous section to the evolution of a quark-gluon plasma containing N f = N c = 3 flavors of massless quarks, and introduce separate distributions for quarks (f q ) and gluons (f g ). We consider QCD tree-level 2 to 2 scatterings, within the small angle approximation. Following the strategy taken in ref. [8], we reduce the coupled Boltzmann equations for the gluon distribution function f g and the quark distribution f q to coupled Fokker-Plank equations, where the currents J q,g and sources S q,g are given by While the currents are associated to conservation laws (the particle numbers), and hence to elastic processes, the sources arise from the inelastic processes that correspond to quark creations and annihilations. Note that all the processes considered conserve the total number of particles. The quantities I a , I b and I c are the following integrals, For a thermal system in equilibrium, I a /I b is identical to the temperature. Finally the quantity denotes a logarithmically divergent integral (the Coulomb logarithm) which arises from small angle scatterings. For a QCD plasma close to thermal equilibrium, ≈ ln g −1 (with g the strong coupling constant). This Coulomb logarithm is treated as a constant in the present work.
With all the quantities thus specified, one may check that the collision kernel in eq. (3.13) conserves energy-momentum and constituent number (number of gluons, quarks and anti-quarks), which implies that the local equilibrium distribution is, in the local rest frame, a Bose-Einstein distribution for the gluons and a Fermi-Dirac distribution for the quarks and the antiquarks, with given temperature and chemical potential (see ref. [8] for details).

JHEP11(2017)161
Assuming Bjorken flow, and integrating over azimuthal angle, one transforms eq. (3.13) into where J g , J q , S g and S q are the currents and sources obtained from eqs. (3.14) after a simple rescaling that eliminates the common factor 4πα 2 s . That is, in deriving eq. (3.17) we have absorbed the factor α 2 s in a redefinition of the time, t → t(α 2 s ), and we have redefined the integrals (3.15) by multiplying them by a factor 4π. With these redefinitions, the time scale of the simulation is measured in units of τ s ≡ (α 2 s Q s ) −1 . For technical reasons related to the numerical precision of the calculated moments, our simulations end at t = 400τ s . We initialize the system at a time t 0 = 1/Q s , and assume that, at this time, the system contains only gluons, with a distribution function given by eq. (3.1). In order to be able to compare with the previous subsection, we shall make a rough estimate of the relaxation time. To do so, we focus on the gluon sector, and the diffusion piece of the current (the one proportional to the gradient of the distribution). By writing the linearized version of the gluon kinetic equation (eq. (3.13)) as δf g /τ rel ∼ δf g /p 2 , we get For p ≈ Q s , N c = 3, this yields We first plot in figure 7 the absolute values of the moments L n divided by the energy density, for n = 1 up to n = 4 for a pure gluon system (solid lines) and a QGP (dashed lines). The pattern seen here is analogous to that in figure 3 corresponding to the relaxation time approximation. The curves present a (broad) peak at the transition between freestreaming at early times and a collision dominated regime at late times. Clearly, the elastic collisions, treated in the small angle approximation, isotropize the system: all moments tend to vanish in the collision dominated regime. The higher moments are damped on a time scale that decreases with increasing order of the moment. Overall, the evolution of moments in figure 7 is comparable to that in figure 3, even semi-quantitatively (to within a factor ∼ 2) if one takes into account the crude estimate of the relaxation time presented above. Note that for a quark-gluon plasma, the isotropization process is significantly delayed and the transition from the free-streaming regime to the collision dominated regime is slower than for the pure gluon case.
In order to compare with hydrodynamics we need to take into account the fact that our Boltzmann equation conserves the number of constituents. The equation for the energy density needs therefore to be completed by a continuity equation for the number density, which, for Bjorken flow, has the simple solution n(t) = n(t 0 )t 0 /t. This also implies that the local distribution function depends on a temperature and a chemical potential. Details on the hydrodynamical calculations are given in appendix B. The fact that particle number is conserved can modify the shear viscosity from its value at vanishing chemical potential. When µ = 0, it is known in weak coupling QCD [33,34] that where κ is a constant depending on the number of colors and quark flavors. For instance, for a pure gluon system, κ ≈ 0.17. We do not have (yet) a precise determination of the viscosity appropriate to the present setting. In the present work, we shall therefore rely on a crude approximation whereby we use eq. (3.20) for η, with κ left as an adjustable parameter. We initialize the hydrodynamical calculation at t = 400τ s , which corresponds to the largest time for which we can solve accurately the Boltzmann equation. At that time, the energy density and the number density are identified with those values obtained from the solution of the Boltzmann equation. Since after t = 200τ s , the high order (n ≥ 2) moments are relatively less significant, as shown in figure 7, one expects the system evolution to follow first order viscous hydrodynamics at later times, t 200τ s . Note that the ratio P L /P T ≈ 0.85 and 0.62 at t = 400τ s for the pure gluon gas and the QGP, respectively. This value for the pure glue system (0.85) is to be compared with that obtained in the previous subsection in the relaxation time approximation (0.9), at times where we found the system to be in the (viscous) hydrodynamic regime. The evolution of the moment L 1 is shown in figure 8, and compared to the solution of the Navier-Stokes hydrodynamics. This comparison allows for the determination of the parameter κ in eq. (3.20). We obtain κ = 0.62 and 0.67 respectively for the pure gluon gas and the QGP. These values are about four times larger than the values expected from weak coupling estimates. These estimates, however, do not take into account the effect of a non vanishing chemical potential, and the formula (3.20) is from that point of view only an approximation in the present context, as already stated. This approximation is sufficient for our purpose: with the viscosity given by eq. (3.20) the pattern observed in figure 8 is similar to that observed in figure 4, showing a brief free-streaming at early times, followed by a slow approach to the viscous hydrodynamic regime at late times. The evolution of the energy density obtained from the solutions of eqs. (3.13) is presented in figure 9(a) for a pure gluon gas (red solid line) and a QGP (blue dashed line). The energy density is rescaled by a factor t 4/3 to be plotted on a visible scale, and also to make clearer the approach to ideal hydrodynamics where the quantity should be constant. One sees that this ideal regime is not reached at t = 400 τ s , which was also indicated by the non zero value of L 1 at that time (see figure 4). Figure 9(b) displays the evolution of the entropy density in the transverse plane. It goes to a constant when the ideal hydrodynamical regime is reached. As expected in kinetics, a significant amount of entropy is produced at early times, when the system starts to redistribute momenta from the initial step function toward a smoother distribution, with in particular a rapid growth of soft modes [7,8]. More entropy is generated in the QGP (blue line), as compared to the pure gluon gas (red line), due to the additional quark degrees of freedom which, for three flavors, contribute about twice as much as the gluons to the total entropy (chemical equilibrium is nearly reached at t = 400τ s ). One sees that Navier-Stokes hydrodynamics describes well the entropy production already at times t 100τ s , for both the gluon gas and the QGP.

Summary and conclusions
In this paper, we have proposed a set of moments, L n , that provide a simple characterization of the momentum anisotropy of the momentum distribution function and take into account Bjorken's symmetry of the longitudinally expanding quark-gluon plasmas created in relativistic heavy ion collisions. At late times, i.e., when the hydrodynamic regime is approached, these moments can be associated to viscous corrections to the energy-momentum tensor, order-by-order in a gradient expansion. The explicit correspondence with the transport coefficients of viscous conformal hydrodynamics has been established. The generalizations to cases without conformal symmetry can be carried out in a similar way. These moments offer a convenient tool to characterize the onset of hydrodynamics.
The time evolution of the moments was studied first with the help of a Boltzmann equation solved within a (constant) relaxation time approximation. This has allowed us to verify the correspondence of the L n 's with viscous hydrodynamics up to order n = 4. The onset of hydrodynamics was subsequently analyzed in terms of these L n 's up to n = 4. The approach to hydrodynamics manifests itself as the decrease of the magnitudes of the L n with n larger than 1. In fact the moments L n are damped more and more rapidly as n increases, which reflects the fact that long wavelength modes thermalize faster than short wavelength ones. Note that the preeminence of the moment L 1 is reminiscent of the arguments used in anisotropic hydrodynamics [35,36].
Similar patterns are observed when we use the more realistic 2 ↔ 2 QCD matrix elements in the Boltzmann equation, which we have solved numerically, within the small scattering angle approximation. Again, we have determined the evolution of L n up to order JHEP11(2017)161 n = 4. In the calculations that we have presented, quarks are entirely generated via gluonquark conversion according to tree-level QCD processes. Although quarks are produced on a relatively short time scale, this production is hindered by Pauli blocking effects which tend to delay significantly the thermalization and the onset of hydrodynamics. Finally, the present analysis ignores the inelastic collisions (soft and collinear splitting of gluons) whose role in the isotropization is emphasized in particular in the bottom-up scenario. We defer their study to future work.

A Contraction of irreducible tensors and Legendre polynomials
The equation (2.19) expresses in terms of a Legendre polynomial the tensorial structure that appears in the leading order gradient expansion of the distribution function. This can be generalized to higher orders. In particular, in the n-th order viscous corrections, the term with the highest power in momentum involves the contractions of irreducible tensors of rank 2n.
Generally, the contraction of a two irreducible tensor of rank 2n constructed from the product of 2n vectors A µ and 2n vectors B µ can be written in terms of Legendre polynomials as follows (cf. [37]), A µ 1 A µ 2 . . . A µ 2n B µ 1 B µ 2 . . . B µ 2n = |A| 2n |B| 2n N 2n P 2n (Â ·B) . (A.1) In this equation, |A| and |B| denote the modulus of the vectors, whileÂ andB are unit vectors. The constant coefficient is As a check, let us apply eq. (A.1) to P µ P ν σ µν . First, we note the following identity, We also know that the tensor ∇ µ U ν with respect to Bjorken flow has only non-zero component at ηη as ∇ η U η = τ , it is thus able to write equivalently

B Hydrodynamics with a with non-zero chemical potential
When comparing the solution of the Boltzmann equation in section 3.2 to viscous hydrodynamics, we need to pay attention to the fact that the QCD 2-to-2 scatterings conserve the number of constituents. As a result the hydrodynamic equations should be complemented with a continuity equation, i.e., we need to solve with energy-momentum tensor T µν given in eq. (2.11) and n µ = nU µ + I µ . 2) can be solved independently and yields τ n(τ ) = cste. The equation for the energy density requires the knowledge of η(T, µ). In this work we ignore the explicit dependence of η on the chemical potential, and use simply eq. (3.20) for η, with κ an adjustable parameter. The comparison with hydrodynamics is done by fixing "initial" conditions at t = 400τ s , and running the equations backwards in time. That is, the initial density and energy density allow us to determine the initial temperature and chemical potential, and hence the viscosity. This process can be repeated to evolve backwards in times steps. The entropy is determined from the thermodynamic relation T s = e + P − µn, which does not involve explicitly the viscosity. This provides a check of the overall consistency of the calculation.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.