Chemical equilibration in weakly coupled QCD

We study thermalization, hydrodynamization, and chemical equilibration in out-of-equilibrium Quark-Gluon Plasma starting from various initial conditions using QCD effective kinetic theory, valid at weak coupling. In non-expanding systems gauge bosons rapidly lose information of initial state and achieve kinetic equilibrium among themselves, while fermions approach the equilibrium distribution only at later time. In systems undergoing rapid longitudinal expansion, both gluons and quarks are kept away from equilibrium by the expansion, but the evolution is well described by fluid dynamics even before local thermal equilibrium is reached. For realistic couplings we determine the ordering between the separate hydrodynamization, chemical equilibration and thermalization time scales to be $\tau_\text{hydro}<\tau_\text{chem}<\tau_\text{therm}$.


I. INTRODUCTION
How gauge theories pushed far from equilibrium thermalize is a central topic in the study of heavy-ion collisions [1]. To what extent the post-collisional debris created in the collision of two nuclei reach local thermal equilibrium before the system cools down, determines how well a fluid dynamical description of the system is applicable.
Our ability to perform first principles non-perturbative real-time calculations in QCD is limited by the infamous sign problem [2], and considerable efforts have been invested to understand thermalization and hydrodynamization in various approximations of the QCD. A prominent example is the N = 4 Super Yang-Mills theory in the limit of large number of colors and strong coupling, which has been studied extensively with holographic methods [3][4][5][6]. The holographic methodology can be applied only to a very limited set of gauge theories, and for generic theories-such as QCD-only weak coupling methods are available. So far the weak-coupling studies of thermalization of far-from-equilibrium systems have been limited to either pure gauge or scalar theories, and studies in QCD have been restricted only to near equilibrium systems [7,8]. Here, we extend the weak coupling treatment of [9,10] by including dynamical fermions, and study how far-from-equilibrium systems approach equilibrium in full leading order QCD description.
Introducing new degrees of freedom to the system adds new structures. It has been argued [11] that the off-equilibrium dynamics of quarks may be significantly slower than that of the gluons, owing partly to smaller group theoretic color factors, and partly to different spin * a.k@cern.ch † a.mazeliauskas@thphys.uni-heidelberg.de statistics and Pauli blocking. It may be then that the equilibration of quarks could be a bottleneck of thermalization as chemical equilibration may take place in a significantly longer time scale. In particular, in the weak coupling picture of heavy-ion collisions, the initial state in midrapidity is dominated by large number of gluons with only few fermions. If the production of fermions is delayed, this could have an impact on the fluid dynamical modelling of heavy-ion collisions, since the equation of state of Quark-Gluon Plasma is different than that of plasma consisting of gluons only. Furthermore, chemically equilibrated QGP is a standard explanation of strangeness enhancement in nucleus-nucleus collisions [12,13], so understanding fermion production from first principles provides an important theoretical validation of this picture. It has been been observed in several theories -both weakly and strongly coupled -that the hydrodynamical constitutive relations are approximately fulfilled in systems that have sizeable anisotropies, that is, hydrodynamization without thermalization [5,10,14,15]. Upon including quark degrees of freedom to the system we may ask the question when does the chemical equilibration happen with respect to the hydrodynamization and thermalization times τ hydro and τ therm .
The basic tool here is the set of Boltzmann equations that are applicable in isotropic systems any time the typical occupancies are smaller than 1/α s . Of course, many authors have already considered the evolution of quark-gluon systems under Boltzmann equations [16][17][18][19][20][21][22][23][24]. What sets our study apart from these works is that we use the effective Boltzmann equations derived by Arnold, Moore, and Yaffe [25], which account for all processes needed for a description that is accurate to leading order in α s . These processes include in-medium screening effects [26], and Landau-Pomeranchuk-Migdal [27][28][29][30] corrected splitting processes in far-from-equilibrium but isotropic systems. To do so we extend the previously developed setup of Refs. [6,9,10,31] to include quark arXiv:1811.03068v1 [hep-ph] 7 Nov 2018 degrees of freedom and we briefly summarize the key elements of the description in Sec. II.
In the following we will consider N c = 3 QCD plasma with N f = 3 flavours of massless quarks in different outof-equilibrium conditions, with and without expansion. While there is only one thermal equilibrium, there are many ways how a system can be out of equilibrium. In Sec. III A we will first study in detail a particularly simple non-expanding system where quarks are absent and gluons are in kinetic equilibrium among themselves. We discuss the processes that produce fermions and subsequently lead to chemical equilibration. In Sec. III B we then move on to discuss non-expanding systems which are initialized with an overoccupied gauge boson distribution. While we are unaware of a physical system where QCD would be found in these conditions, cosmological reheating may result in system of overoccupied gauge bosons (see e.g. [32]). Following the time evolution of the overoccupied system, we see that due to the slower dynamics of fermions, gauge bosons reach kinetic equilibrium before the chemical equilibration. Once the gauge bosons have reached kinetic equilibrium the evolution proceeds as in our first example. Finally, in Sec. IV we turn to a system of overoccupied gluons undergoing boost-invariant longitudinal expansion. This is the expected initial conditions in heavy-ion collisions in the asymptotic weak coupling limit [33,34]. For moderate values of the coupling constant α s ∼ 0.3 we observe a rapid memory loss of initial conditions and the chemical and hydrodynamical equilibrium is approached along a universal curve. Finally, we conclude with the discussion of the separate equilibration time scales in Sec. V.

II. EFFECTIVE KINETIC THEORY
The effective kinetic theory that we use to describe thermalization is the Effective Kinetic Theory (EKT) of Arnold, Moore and Yaffe [25], which is leading order accurate in the QCD coupling constant λ = g 2 N c = 4πα s N c in the combined limit of weak coupling (λ → 0) and nonperturbative occupancies (λf → 0) for modes whose momenta are larger than the thermal screening scale in the nonequilibrium system p 2 At leading order in the coupling constant, the EKT describes the time evolution of color/spin averaged distribution function f s with an effective 2 ↔ 2 scattering and a 1 ↔ 2 effective splitting terms. The resulting Boltzmann equation for homogeneous non-expanding system is with massless dispersion relation, p 0 = |p| = p. The index s refers to different particle species in the theory. Expanding upon previous implementations of pure gauge theories in Refs. [6,9,10,31] to QCD, s now stands for gluons and 2N f massless fermions (with quarks and antiquarks counted separately) 1 . The symmetrized 2 ↔ 2 collision terms in the right hand side of Eq. (1) reads where |M ab cd | 2 is a 2 ↔ 2 scattering amplitude-squared summed over all degrees of freedom of the external legs (ν q = 2N c for quarks and ν g = 2(N 2 c − 1) for gluons), abcd is a sum over all particle and antiparticle species, and p = d 3 p 2p(2π) 3 is a shorthand notation for Lorentz invariant momentum integral. The second line is the usual phase-space loss and gain terms, while the Kroneker and Dirac delta functions in the last line accounts for the possibility of particle s to be on any of the four external lines. Finally the numerical prefactors in front of the integral corrects the double counting of identical processes. 1 In this work we consider plasma at zero chemical potential with quark and anti-quark distributions being equal.
The effective matrix elements |M ab cd | 2 in Eq. (2) are for most kinematics the normal tree-level vacuum matrix element (see Table II in Ref. [25]). For soft small angle scatterings with energy transfer ω p, k, the treelevel Coulomb and Compton scatterings are infrared divergent elevating a set of diagrams with arbitrary number of loops to the same magnitude as the tree-level diagrams. These effects become important at the in-medium screening scale p ∼ m g , m q . Here the in-medium effective masses of gluon and quarks are given by For momentum transfer of this order the dispersion of the internal line in the computation of |M ab cd | 2 gets an O(1) correction from the in-medium physics. Therefore, for the problematic soft scattering we replace the matrix element with that computed in the Hard Thermal Loop (HTL) approximation that self-consistently treats the medium interaction correctly to leading order. We perform this substitution by removing the infrared divergent small angle approximation from the full matrix element and replace it with the small angle approximation of the full HTL rate [35]. Specifically for a soft gluon or fermion exchange with the momentum transfer q = |p − p| in tchannel, the divergent term (u − s)/t ∼ 1/q 2 is replaced by IR regulated term where ξ g = e 5/6 /2 and ξ q = e/2 are fixed such that the matrix element reproduces the full HTL results for drag and momentum diffusion properties of soft gluon scattering [35] and gluon to quark conversion gg → qq [36,37] at leading order for isotropic distributions. While the soft ω ∼ m g scatterings do not appreciably change the momentum state of the particle, they may bring the particle slightly off shell and make it kinematically possible for the particle to decay through nearly collinear splitting. This makes the effective 1 ↔ 2 matrix element a leading order effect. It is included as C 1↔2 [f ](p) on the right hand side of the Boltzmann equation Eq. (1) and explicitly where the unit vectorn =p/|p| defines the splitting direction and γ a bc (p; p , k ) is the effective collinear splitting rate including Landau-Pomeranchuk-Migdal [27][28][29][30] suppression of collinear radiation. Factoring out the kinematic splitting functions the rates γ g gg (p; p , k ) = γ q qg (p; p , k ) = are given by an effective vertex resuming an infinite number of possible soft interaction with the medium [25]. It is found by solving the following integral equation and F s (p; p , k ) is defined as In this work the strength of soft momentum background fluctuations A(q ⊥ ) are treated using an isotropic screening approximation [38] A the energy difference δE is defined as and the effective temperature T * is given by Instead of solving Eq. (10) directly, the required two dimensional integral Eq. (11) is expressed as the value at the origin of the Fourier transformed functionF s , which solves a Fourier transformed Eq. (10). We solve it using the basis function method [39] and parametrize the solution for the Monte-Carlo sampling of the collision kernel in Eq. (6). For the distribution functions we use previously developed discretization scheme, which preserves energy and number density [6,10,35]. The distribution functions are discretized in spherical polar coordinates on logarithmic momentum grid with typical momentum values p min = 0.01Q S and p max = 15Q s and N p = 100. For anisotropically expanding systems the longitudinal momentum fraction cos θ = p z /p is discretized on a uniform grid with typical value of N θ = 200. We considered azimuthally symmetric distributions. The collision integrals were calculated by Monte Carlo sampling of the phase space with importance sampling [6]. At each time step the 2 ↔ 2 collision integral Eq. (2) was calculated using N 2↔2 N p N θ randomly generated vector quadruplets p, k, p , k satisfying the momentum and energy conservation, where N 2↔2 = 50, 100. For 1 ↔ 2 collision integral Eq. (6) we used N 1↔2 N p samples of momentum p, p , k = p − p combinations, which were reused for each angular direction (here N 1↔2 = 50, 100).

A. Kinetically equilibrated initial conditions
In order to gain intuition to the far-from-equilibrium dynamics, we start with a particularly simple system where gauge bosons and fermions are initialized at time but with different initial temperatures T g,init. = T q,init. . In such situation we call that quarks and gluons are in kinetic equilibrium among themselves, but not in thermal equilibrium with each other. The system will relax into state in which both the quarks and the gluons are in equilibrium with the same temperature T final . The energy conservation dictates that the final temperature will be given by If we start with pure gluon initial state, i.e., T q,init. = 0, the fermion number is initially zero. In order to reach chemical equilibrium, the fermion number has to be subsequently generated by pair production either through medium induced g → qq-splitting processes or alternatively through gg → qq conversions. Multiplying the Boltzmann equation Eq. (1) by 2N f ν q p 2 /λ 2 T 3 for quarks, we obtain the equation for the rate of change of fermion number (per momentum) and similarly for gluons. In Fig. 1 we show the rates C s 22 and C s 12 separately and as a sum for coupling constant λ = 0.1 and temperature T = T g,init . We see that the 2 ↔ 2 processes dominate fermion production around p ∼ T (blue dashed line), while for p ∼ m D = √ 2m g the splitting processes become roughly equally important (blue dotted line). The changes in the gluon distribution mirror the fermionic ones. We see that gg → qq conversion reduces the number of gluons at the same momentum scale p ∼ T where fermions are created (blue and red dashed lines), while soft collinear radiation from p ∼ T gluons (red dotted line) produces soft fermions. Importantly, the resulting fermion spectrum is non-thermal.
We now turn to study how the above system evolves toward equilibrium as a function time. For the same λ = 0.1 as above, Fig. 2 displays the time evolution of several different effective temperatures T eff g/q,α , defined through the αth moments of the fermion (F = 1) and boson (F = 0) distribution functions   The effective temperatures are normalised such that when the system is in thermal equilibrium all T eff g/q,α (for all α) are equal to the equilibrium temperature T eff g/q,α = T final , that is .
Lower α values are more sensitive to the infrared of the distribution, whereas larger α values describe the UV part of the distribution function.
We display the effective temperatures as a function of relaxation time of the final thermalized system where η/s is the specific shear viscosity, whose relation to λ is discussed in Appendix A. This relaxation time is parametrically of order of the transport mean free path of the thermalized system τ R ∼ (λ 2 T final ) −1 , but it has been observed in [6,10,40,41] that expressing the coupling constant λ in terms of the specific viscosity accounts for large numerical corrections that go beyond the parametric expression and leads to better scaling behaviour for different values of λ.
We first note that the effective temperatures of gluons (red lines in Fig. 2) decrease while those of the quarks (blue lines) increase. During this evolution the gluon effective temperatures approximately overlap signifying that gluons remain close to the kinetic equilibrium through the whole evolution. In contrast to gluonsand consistent with the non-equilibrium spectrum of the fermion production rate in Fig. 1-the fermion effective temperatures differ until full chemical equilibration is achieved. During this evolution the fermion spectrum is harder than that in kinetic equilibrium as seen from the ordering of the effective temperatures T eff 4 > T eff 3 > . . . > T eff −1 . In chemical equilibrium N f = 3 fermions constitute e q,eq /e q,total ≈ 0.66 of the total equilibrium energy density (see Eq. (16)). It can be seen in Fig. 3 that by the time t chem ≈ 1.4τ R for λ = 0.1 the fermion energy density e q has reached 90% of its equilibrium value e q,final , which we take as our somewhat arbitrary definition of the chemical equilibration time, i.e.
To study the coupling dependence of the chemical equilibration, we repeat the above calculation with the same initial conditions, but with several different values of λ = 0.1, 1, 10, corresponding to η/s ≈ 1, 35, 1900, respectively. Figure 3 shows the time-evolution of fermion energy fraction for these different values of 't Hooft couplings, for which the chemical equilibration time varies by three orders of magnitude. However, as is seen from the figure, the functional forms of the time-evolutions of the energy densities follow closely a common form when described in terms of the relaxation time τ R , with the system reaching chemical equilibrium around for all studied values of λ.
We note however that for moderate values λ 10, we expect substantial NLO order corrections to transport properties of QGP [42]. Nevertheless, while the kinetic theory suffers from large systematic uncertainties for these values of λ, the fact that the kinetic theory calculation itself does not fail catastrophically may allow one to extrapolate the results for semi-quantitative order of magnitude estimates even at larger values of λ using τ R scaling, where the dependence on the coupling constant only enters through the specific shear viscosity η/s. Note that other transport coefficients in units of η/s, e.g. τ π /(η/(sT )), are much less sensitive to the changes of the coupling constant [43].
We note in passing that similar qualitative features of the equilibration are seen by varying the initial starting gluon and fermion temperatures. In particular, the gluons remain in approximate kinetic equilibrium throughout evolution even if one starts with fermion dominated initial conditions (data not shown). In this case the fermion energy fraction approaches equilibrium ratio from above (see Fig. 3) and the chemical equilibration takes place at t chem /τ R ∼ 0.5 − 0.9, where now e q (t chem )/e q,final = 1.1 in this case.

B. Overoccupied initial conditions
We now turn to a system starting with an overoccupied distribution of gluons and no fermions. This system has initially too many gluons, f g 1, with too soft momenta compared to a thermal ensemble with the same energy density, that is T eff 2 T eff 1 T eff 0 . As the fermions cannot be overoccupied such a system is necessarily dominated by the gluons in the early stages. As discussed in detail in [9,35,[44][45][46][47][48][49], overoccupied gluonic systems thermalize via a self-similar cascade which carries the energy and particle number from the infrared to the ultraviolet via elastic and inelastic scattering. As long as the system is parametrically overoccupied f g 1, the cascade is self-similar, that is, the gluon distribution function at a given time t is given by an approximately stationary scaling functionf which is insensitive to the initial conditions where Q 4 = 2π 2 λ p pf g (p). This scaling form is reached in a time that is proportional to the scattering rate of the initial condition [9], which is parametrically faster than the thermalization time for a parametrically overoccupied system. The approach to the scaling form is discussed in detail in [49]. Once the typical momentum of the cascade p max reaches the thermal scale T final , the system equilibrates. How the presence of fermions may affect the cascade has been studied using (semi-)classical Yang-Mills simulations in the regime where the gluons are still highly overoccupied and p max T final [50]. We now answer the question what happens when p max ∼ T .
The green lines in Fig. 4 show the time evolution of a pure N c = 3 gauge theory with over-occupied gluon initial conditions 2 . For the effective temperatures, this scaling form Eq. (22) corresponds to a power-law behaviour If the system is parametrically overoccupied, the effective temperatures are also parametrically separated T eff α T eff α+1 , but when the system thermalizes around t ≈ 0.3 τ R , the effective temperatures again collapse. (N c = 3, N f = 3) with the same initial condition of overoccupied gluons with no quarks present. The time evolution differs now from the pure glue system as the fermion number is dynamically generated. At early times the quark effective temperatures are small and gluonic temperatures approximately follow the classical powerlaws. The cascade ends and the gluonic sector reaches a kinetic equilibrium among themselves again around where we have defined the kinetic equilibration timein analogy with the chemical equilibration time-by demanding This timescale is significantly faster than the timescale for chemical equilibration t chem . Indeed at the time the kinetic equilibrium among gluons is reached there are only few quarks present and the state of the system is approximately the same as in the case of thermal initial conditions studied in the previous section. From this point on, the chemical equilibrium follows the same pattern which was described in the previous section. To emphasize this point, in Fig. 5 we show the fermion energy fraction for over-occupied initial conditions (colored lines) on top of the thermal initial conditions (gray lines), which were shown in Fig. 3.

IV. CHEMICAL EQUILIBRATION IN EXPANDING SYSTEMS
In the weak coupling description of heavy-ion collisions the state of the system right after the initial particle production is given by an overoccupied distribution of gluons [33,34]. The novel feature compared to the previous section is that the geometry of the collision system is such that the overoccupied matter is undergoing a rapid, approximately boost invariant longitudinal expansion. Such system has been studied in detail in pure gauge theory [10,31] and it forms the link between the initial condition and the hydrodynamic stage in phenomenological multi-stage simulations of heavy-ion collisions [40,41]. Here, we study the expanding plasma in full QCD.
Assuming boost invariant form of the distribution function, the Boltzmann equation can be written in the form [17] where τ is the Bjorken time τ = √ t 2 − z 2 . The expansion redshifts the distribution in p z direction making it more anisotropic along the longitudinal momentum, while 2 ↔ 2 scatterings act to isotropize the system [51]. While anisotropic systems could suffer from the presence of unstable plasma modes [52][53][54] affecting the kinetic dynamics [48,55], the detailed 3+1D classical-statistical Yang-Mills simulations [47,56] found that late time evolution of anisotropic systems is in agreement with kinetic theory expectations neglecting plasma instabilities [51]. In the absence of general non-equilibrium formulation of QCD kinetic theory [48], we use QCD kinetic theory with isotropic approximations, which remove the unstable modes, to study the equilibration in expanding systems.
While the expansion conserves total energy, the local energy density in a given rapidity slice decreases as a function of time. At late times when the system is close to local thermal equilibrium, the time evolution of the temperature is given by ideal hydrodynamics with constant T (τ )τ 1/3 . As the target temperature to which the out-of-equilibrium system aims to thermalize changes, so does the kinetic relaxation time. In the following we follow the practice of [40] and for each simulation, we determine the asymptotic value of T (τ )τ 1/3 | τ →∞ , and use the ideal hydrodynamics relation to extract what the temperature of the system would have been at earlier times if the full time evolution of the system were described by ideal fluid dynamics and use that in our definition of time dependent kinetic relaxation time We consider expanding system with over-occupied initial condition motivated by Color-Glass-Condensate where the values of A and Q 0 , and ξ are adjusted to reproduce the mean transverse momentum squared p 2 T and energy density e(τ 0 ) at the initial time τ 0 extracted from the classical lattice simulations of initial stages of the collision [10,57]. The anisotropy parameter ξ determining the ratio of longitudinal to transverse pressure is chosen such that P L P T and is set in the following ξ = 10. The same initial conditions have been studied in pure gauge theory in [10,31,40]. Here, as in [10,31,40], we use τ 0 = 1/Q s , Q 0 = 1.8Q s and we set A = 5.24. Here Q s is the saturation scale of Color-Glass-Condensate, and is of order Q −1 s ∼ 0.1fm. In Fig. 6 we show the time evolution of the effective temperatures T eff α for the initial conditions Eq. (29) and λ = 5, which corresponds to η/s ≈ 2.75. The rapid longitudinal expansion quickly inverts the hierarchy of temperatures from an overoccupied state f g 1 (T eff α > T eff α+1 ) to an underoccupied state f g 1, which is well understood as the first stage of the bottom-up termalization [51]. This transition takes place well before a substantial number of fermions are produced as can be seen from the good overlap of the time evolution of the effective temperatures in pure gauge theory with the full theory in Fig. 6. This suggests that during this early phase the presence of fermions does not significantly affect the evolution of bulk quantities. Note that the effective temperatures are sensitive only to the angular averaged distribution function, but because of the longitudinal expansion the system is highly anisotropic.
Next we look at the chemical equilibration of fermion energy fraction of the equilibrium energy density. In Fig. 7 we show the approach to chemical equilibrium for different values of the coupling constant λ and overoccupied initial state. Using the same ad-hoc criterion of chemical equilibration as in non-expanding system Eq. (20), we find that chemical equilibration happens for τ chem ∼ 1 − 2τ R for a wide range of coupling constants 0.5 ≤ λ ≤ 20 and only for very small coupling λ = 0.1, we get τ chem ≈ 2.8τ R . We would like to remind that parametrically the chemical equilibration scales as ∼ (λ 2 T ) −1 , so in physical units the equilibration times change by several orders of magnitude. Scaling with relaxation time Eq. (28), reduces this vast separation of scales and for couplings λ > 1 we observe the collapse of equilibration dynamics to the same curve and the chemical equilibration is reached at τ chem ≈ 1.2τ R . The difference in equilibration times for λ ≤ 1, which was not present in the non-expanding case, arise from the additional scale -the expansion rate 1/τ . Since the starting time τ 0 = 1/Q s is kept fixed, for smaller values of the coupling constant, the system experiences a longer phase where the anisotropy is still increasing, i.e. the first stage of bottom-up, which delays equilibration.
We note that the chemical equilibration takes place when the system is still highly anisotropic. This is seen in Fig. 8(a) depicting the time evolution of the ratio of the longitudinal and transverse pressures The system becomes isotropic only at very late times. of λ. However similarly to the pure gauge theory, the system's time evolution is well described by fluid dynamics well before pressure anisotropies become small, that is the system exhibits "hydrodynamization without thermalization".
We quantify the approach to thermal equilibrium and hydrodynamization by defining two additional timescales τ therm and τ hydro , in analogy to Eq. (20). We require the combined gluon and fermion energy density e = e g + e q to be within 10% of ideal and viscous hydrodynamic estimates of energy density Here T id is the ideal estimate of the local temperature Eq. (27) and T 1st is the 1st order viscous hydrodynamic solution of longitudinally expanding system with the same late time asymptotics as T id . The analytical solution for first order conformal hydrodynamical equations of motion can be written in units of τ R as [58,59] T 1st (τ ) In Fig. 8(b) we show the time evolution of the total energy density scaled by ideal estimate, i.e. e/e id. = (T /T id. ) 4 . We find that 90% of equilibrium energy is reached at τ therm ∼ 2τ R for the coupling constant values 0.5 ≤ λ ≤ 20. Only for λ = 0.1 the approach to equilibrium is slower and this criterion is satisfied at τ therm ∼ 2.5. Next, we use the 1st order hydrodynamic estimate for temperature Eq. (34) and compare it to the kinetic theory 3 . We achieve agreement with full kinetic theory evolution at very early times and by the time τ hydro 0.5τ R the criterion Eq. (33) is satisfied for 0.5 ≤ λ ≤ 20. For λ = 0.1 thus defined hydrodynamization takes place somewhat later at τ ∼ 1.3τ R .

V. CONCLUSIONS
In this paper we presented a complete simulation of chemical equilibration in leading order QCD kinetic theory in stationary and expanding systems with infinite transverse extent. By analysing how out-of-equilibrium plasma of N c = 3 gluons and N f = 3 quarks relax to the common thermal equilibrium for different values of the coupling constant λ, we determined the chemical equilibration time in non-expanding isotropic systems, which we define by requiring quark energy fraction to be within 10% of their equilibrium value. For initial conditions with no quarks present, thus defined chemical equilibrium is reached at time t chem ∼ 1.1−1.4τ R , where τ R = (4πη/s)/T final is the kinetic relaxation time and η/s(λ = 0.1, 1, 10) ≈ 1900, 35, 1. We also note faster gluon dynamics, which results in gauge bosons reaching kinetic equilibrium among themselves before thermalizing with fermions. Consequently, for the case of the over-occupied gluon initial state, gluons first equilibrate through a self-similar cascade as in pure glue theory at times t g kinetic ∼ 0.2τ R , and then the equilibrium quark densities are produced by the quasi-thermal gluon background. We would like to emphasise that the chemical equilibration time dependence on the coupling constant λ-the only microscopic parameter of the theory -is very well captured by the specific shear viscosity η/s, which is the macroscopic QGP property. As the relaxation time τ R changes by orders of magnitude, the equi-libration dynamics in rescaled units t/τ R remains largely unchanged.
Next we studied the QCD equilibration in homogeneous, but longitudinal expanding systems, which is a relevant situation for heavy ion phenomenology. There the expansion prevents the system from ever reaching static thermal equilibrium and one instead may define thermalization time τ therm by requiring the total energy density to be within 10% with the value given by ideal hydrodynamic evolution e id. (τ ) = (eτ 4/3 ) ∞ τ −4/3 . Here we note that in expanding case the effective temperature of the plasma is decreasing and so we defined time-dependent kinetic relaxation time τ R (τ ) ∼ τ 1/3 which grows in time, see Eq. (28). For Color-Glass-Condensate motivated, anisotropic and overoccupied initial conditions, Eq. (29), such thermalization happens at τ therm ∼ 2τ R (τ ) for a range of the coupling constants 0.5 ≤ λ ≤ 20, but is somewhat delayed for λ = 0.1. The chemical composition of the plasma changes rapidly in the first couple units of relaxation time. Similarly to the nonexpanding systems, gluons undergo a kinetic equilibration faster than fermions, in agreement with two-stage QGP equilibration argued in Ref. [60]. However only when the expansion rate slows down and viscous corrections to the particle distribution function are small enough, particle distributions are well approximated by the equilibrium Bose-Einstein or Fermi-Dirac distributions. Keeping in mind that the effective kinetic relaxation time τ R is growing in time as the temperature is decreasing, the chemical equilibration in longitudinally expanding systems for moderate values of the coupling constants λ = 5, 10, 20 (α s ∼ 0.1−0.5) proceeds very similarly in rescaled units to non-expanding case and chemical equilibration is reached at τ chem ∼ 1.2τ R (τ ). For smaller values of the coupling constant λ ≤ 1 we do not see the collapse to the same universal curve and the chemical equilibration (in units of τ R ) takes place later.
In summary, the chemical composition is an important property of the expanding QGP fireball, which is not captured by hydrodynamic modelling of heavy ion collisions, but is essential for the hadrochemistry, photon production, and determines which equation of state best describes the medium. For realistic values of the coupling constant α s ∼ 0.3, we find that even in expanding systems the coupling constant dependence can be factored out by rescaling time with kinetic relaxation time τ R = (4πη/s)/T id. , which results in the following ordering of hydrodynamization, chemical equilibration and thermalization timescales according to criteria given in Eqs. (20), (33) and (32). Such universality allows one to convert the dimensionless time τ /τ R to physical units by matching the late time constants (τ 1/3 T ) ∞ and η/s from hydrodynamical modelling of heavy ion collisions and which is explored in our companion paper [61]. Shear viscosity over entropy ratio as a function of the coupling constant λ = Ncg 2 in our leading order kinetic theory implementation with two regularization schemes of the elastic collision kernel (see text). The dashed line corresponds to the next-to-leading-log result from Ref. [7].