Equilibration of weakly coupled QCD plasmas

We employ a non-equilibrium Quantum Chromodynamics (QCD) kinetic description to study the kinetic and chemical equilibration of the Quark-Gluon Plasma (QGP) at weak coupling. Based on our numerical framework, which explicitly includes all leading order processes involving light flavor degrees of freedom, we investigate the thermalization process of homogeneous and isotropic plasmas far-from equilibrium and determine the relevant time scales for kinetic and chemical equilibration. We further simulate the longitudinally expanding pre-equilibrium plasma created in ultrarelativistic heavy-ion collisions at zero and non-zero density of the conserved charges and study its microscopic and macroscopic evolution towards equilibrium.


I. Introduction
Non-equilibrium systems are ubiquitous in nature and of relevance to essentially all disciplines of modern physics. Despite the appearance of non-equilibrium phenomena in a variety of different contexts, there is a rather limited number of theoretical methods to study the realtime evolution of quantum systems, most of which rely on a set of approximations to study microscopic and macroscopic real-time properties of complex many-body systems.
Specifically, for fundamental theories of nature, the question of understanding and describing nonequilibrium processes in the strongly interacting Quantum Chromodynamics sector of the standard model, has gained considerable attention in light of highenergy heavy-ion collision experiments at the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC). Somewhat surprisingly, it turns out that the complex space-time dynamics of high-energy heavy-ion collisions on space-and timescales ∼ 10fm, can be rather well described by modern formulations of relativistic viscous hydrodynamics [1], which has become the primary tool of heavy-ion phenomenology [2,3]. Nevertheless, due to the limited availability of theoretical approaches, it remains to some extent an open question how the macroscopic hydrodynamic behavior emerges from the underlying non-equilibrium dynamics of QCD, albeit significant progress in this direction has been achieved in recent years [4][5][6][7][8][9][10][11][12][13][14][15].
Beyond high-energy heavy-ion collisions similar questions arise in Cosmology, where the non-equilibrium dynamics of QCD and QCD-like theories can certainly be expected to play a prominent role in producing a thermal abundance of standard model particles between the end of inflation and big bang nucleosynthesis. However, at the relevant energy scales, the field content of the early universe is not necessarily well constrained, and a detailed understanding of the thermalization of the early universe at least requires the knowledge of the coupling of the standard model degrees of freedom to the inflation sector, which makes this problem significantly more difficult. Nevertheless studies of the thermalization of the isolated QCD sector still bear relevance to this question, as some of the basic insights into the thermalization process of QCD or QCD-like plasmas can be adapted to Cosmological models, as recently discussed e.g. in Ref. [16][17][18].
Even though Quantum Chromodynamics (QCD) exhibits essentially non-perturbative phenomena such as confinement at low energies, strong interaction matter becomes weakly coupled at asymptotically high energies owing to the renowned property of asymptotic freedom. Specifically, for thermal QCD properties, it is established from first principles lattice QCD simulations, that above temperatures T pc ∼ 155MeV [19][20][21][22] hadronic bound states dissolve into a Quark-Gluon plasma (QGP) and the approximate chiral symmetry of light-flavor QCD is restored. While (resummed) perturbative approaches to QCD are able to describe the most important static thermal properties of high-temperature QCD down to approximately ∼ 2T pc [23], the perturbative description appears to be worse for dynamical properties, where e.g. next-to-leading order calculations of transport coefficients [24,25] yield large corrections to the leading order results [26,27], indicating a poor convergence of the perturbative expansion. Nevertheless, it is conceivable that at energy scales corresponding to ≳ 4T pc , achieved during the early stages of high-energy heavy-ion collisions [28], perturbative descriptions can provide useful insights into the early-time non-equilibrium dynamics of the system. Besides the potential relevance to earlyuniverse Cosmology and Heavy-Ion phenomenology, it is also of genuine theoretical interest to understand the unique microscopic dynamics of thermalization processes in QCD or QCD-like plasmas.
During the past few year, significant progress in understanding thermalization and "hydrodynamization", i.e. the onset of hydrodynamic behavior, in high-energy heavy-ion collisions has been achieved, within the limiting cases of weakly coupled QCD [4,[6][7][8][9] and stronglycoupled holographic descriptions [11,[29][30][31]. Despite clear microscopic differences, a common finding is that the evolution of macroscopic quantities, such as the energy momentum tensor, follows a hydrodynamic behavior well before the system reaches an approximate state of local thermal equilibrium.
Specifically, for weakly-coupled QCD plasmas, a detailed microscopic understanding of the thermalization process has also been established, as described e.g. in the recent reviews [32,33]. Different weak-coupling thermalization scenarios based on parametric estimates [34][35][36][37], distinguish between two broadly defined classes of non-equilibrium systems, commonly referred to as overoccupied or under-occupied [32], which undergo qualitatively different thermalization processes. While the thermalization of over-occupied QCD plasmas proceeds via a self-similar direct energy cascade [38][39][40][41], as is the case for many far-from equilibrium systems [42][43][44], underoccupied QCD plasmas undergo the so-called "bottomup" scenario [34] where thermalization proceeds via an inverse energy cascade, which is in many ways unique to QCD and QCD-like systems. Earlier parametric estimates have now been supplemented with detailed simulations of the non-equilibrium dynamics based on classicalstatistical lattice gauge theory [38][39][40][41]45] and effective kinetic theory [46][47][48][49][50]. However, with the exception of Ref. [47,50], all of the aforementioned studies have been performed for SU (N c ) Yang-Mills theory, i.e. only taking into account the bosonic degrees of freedom and neglecting the effect of dynamical fermions.
Central objective of this paper is to extend the study of thermalization processes of weakly coupled non-abelian plasmas, to include all relevant quark and gluon degrees of freedom. Based the leading order effective kinetic theory of QCD [51], we perform numerical simulations of the non-equilibrium dynamics of the QGP, to characterize the mechanisms and time scales for kinetic and chemical equilibration processes. By explicitly taking into account all light flavor degrees of freedom, i.e. gluons (g) as well as u,ū, d,d, s,s quarks/anti-quarks, we further investigate the non-equilibrium dynamics of QCD plasmas at zero and non-zero values of the conserved u, d, s charges.
We organize the discussion in this paper as follows. We begin with an brief explanation of the general setup in Sec. II, where we discuss the characterization of weakly coupled non-equilibrium QCD plasmas in Sec. II A, and outline their effective kinetic description in Sec. II B. Based on this framework, we study different thermalization mechanisms of the QGP, starting with the chemical equilibration of near-equilibrium systems in Sec. III. Subsequently, in Sec. IV we investigate kinetic and chemical equilibration processes in far-from equilibrium systems considering the two stereotypical examples of overoccupied systems in Sec. IV A and under-occupied systems in Sec. IV B. In Sec. V we continue with the study of longitudinally expanding QCD plasmas, which are relevant to describe the early time dynamics of high-energy heavy-ion collisions. Here, we mainly focus on the microscopic aspects underlying the isotropization of the pressure, and evolution of the QGP chemistry at zero and non-zero net-baryon density, and refer to our companion paper [52] for additional discussions on the implications of our findings in the context of relativistic heavy-ion collisions. We conclude in Sec. VI with a brief a summary of our most important findings and a discussion of possible future extensions. Several Appendices A, B, C contain additional details regarding the details of our numerical implementation of the QCD kinetic equations.

II. Non-equilibrium QCD
Generally the description of non-equilibrium processes in Quantum Chromo Dynamics (QCD) represents a challenging task, and at present can only be achieved in limiting cases, such as the weak coupling limit. We employ a leading order kinetic description of QCD [51], where the non-equilibrium evolution of the system is described in terms of the phase-space density f (⃗ x, ⃗ p, t) of on-shell quarks and gluons. We will focus on homogenous systems, for which the phase-space density f (⃗ x, ⃗ p, t) = f (⃗ p, t) only depends on momenta and time, and investigate the non-equilibrium dynamics of the QGP, based on numerical solutions of the QCD kinetic equations. Below we provide an overview of the relevant ingredients, with additional details on the numerical implementation provided in Appendices A, B, C.

A. Non-equilibrium properties of the Quark-Gluon Plasma
Before we address the details of the QCD kinetic description, we briefly introduce a few relevant quantities that will be used to characterize static properties and interactions in non-equilibrium systems. We first note that both equilibrium, as well as non-equilibrium systems can be characterized in terms of their conserved charges, which for the light flavor degrees of freedom of QCD correspond to the conserved energy density e, and the conserved net-charge densities ∆n u , ∆n d , ∆n s of up,down and strange quarks. Evidently in thermal equilibrium, the maximal entropy principle uniquely determines the phase-space distribution of gluons and quarks with well-defined temperature T eq and chemical potential µ f,eq determined by the values of the densities of the charges according to where we denote z f = f,eq Teq for the three light flavors f = u, d, s, which we will treat as massless throughout this work. Even though a non-equilibrium system can no longer be characterized uniquely in terms of its conserved charges, it is nevertheless useful to associate effective temperatures T ldm and chemical potentials µ f,ldm with the system, which can be determined via the so called Landau matching procedure of determining T ldm , µ f,ldm from the conserved charges according to the relations in Eq. (2). Specifically for systems with conserved energy and charge densities, T ldm and µ f,ldm will ultimately determine the equilibrium temperature T eq = T ldm and chemical potential µ f,eq = µ f,ldm once the system has thermalized.
Besides the densities of the conserved quantities, there is another set of important quantities relevant to describe the interactions in non-equilibrium QCD plasmas [51]. Specifically, this includes the in-medium screening masses of quarks and gluons, which in the case of the gluon can be expressed as in terms of the Debye mass with non-equilibrium gluon and quark distributions f g (⃗ p), f q (⃗ p), fq(⃗ p). Similarly, the thermal quark masses m 2 Q f for f = u, d, s quarks also enter in the kinetic description and can be expressed as While the screening masses m 2 D and m 2 Q f determine the elastic scattering matrix elements, the calculation of the effective rates for inelastic processes also requires the asymptotic masses of quarks and gluons, m 2 ∞,a which to leading order in perturbation theory can be related to the respective screening masses according to m 2 ∞,g = m 2 D 2 and m 2 ∞,Q f = 2m 2 Q f . Since inelastic interactions are induced by elastic collisions, their effective in-medium rates are also sensitive to the density of elastic interaction partners which receives the usual Bose enhancement f g (⃗ p)(1 + f g (⃗ p)) and Fermi blocking f q q (⃗ p)(1 − f q q (⃗ p)) factors. Since we will frequently characterize the non-equilibrium evolution of the QGP in terms of the above dynamical scales, we further note that the quantity g 2 T * characterizes the rate of small angle scatterings in the plasma, with T * defined such that in equilibrium T * (eq) = T eq corresponds to the equilibrium temperature.

B. Effective Kinetic Theory of Quark-Gluon Plasma
We adopt an effective kinetic description of the QGP, which at leading order includes both "2 ↔ 2" elastic processes as well as effective "1 ↔ 2" collinear inelastic processes. Specifically for a spatially homogeneous system, the time evolution of the phase-space density of quarks and gluons is then described by the Boltzmann equation for QCD light particles "a = g, u,ū, d,d, s,s" With regards to the numerical implementation, we follow previous works and solve the QCD Boltzmann equation directly as an integro-differential equation using pseudo spectral methods [31,53]. Our numerical implementation of the non-equilibrium dynamics rely on a discretized form of the Boltzmann equation based on a weight function algorithm [54], which is described in detail in Appendices A, B. Based on Eq. (7), the discretized form of the Boltzmann equation for species "a" can be written as where in accordance with Eq. (7), C 2↔2 a [n](i p , j θ , k φ , t) and C 1↔2 a [n](i p , j θ , k φ , t) correspond to discretized moments of the collision integral. Based on a suitable choice of the weight functions w , the discretization of the collision integrals is performed such that it ensures an exact conservation of the particle number for elastic collision, as well as an exact conservation of energy for both elastic and inelastic collisions.

Elastic Collisions
Within our effective kinetic description, we include all leading order elastic scattering processes between quarks and gluons, where following previous works [4,7,46,50] the relevant in-medium scattering matrix elements are determined based on an effective isotropic screening assumption.
a. Collision Integral We follow the notation of Arnold, Moore and Yaffe (AMY) [51], where the elastic collision integrals for particle a with momentum ⃗ p 1 participating in scattering process a, b → c, d with p 1 , p 2 ↔ p 3 , p 4 takes the form with dΠ 2↔2 denoting the measure and ν G = 2(N 2 c − 1) = 16, ν Q = 2N c = 6 denoting the number of gluon and quark degrees of freedom. By M ab cd (⃗ p 1 , ⃗ p 2 ⃗ p 3 , ⃗ p 4 ) 2 we denote the square matrix element for the process "a, b ↔ c, d" summed over spin and color for all particles, while F ab is the statistical factor for the "a, b ↔ c, d" scattering process where "±" provides a Bose enhancement (+) for gluons and Fermi blocking (−) for quarks, such that the first term in Eq.(11) represents a loss term, whereas the second term in Eq.(11) corresponds to a gain term associated with the inverse process.

b. Scattering matrix elements
Elastic scattering matrix elements for the various 2 ↔ 2 processes can be calculated in perturbative QCD (pQCD) [51], with the corresponding leading order matrix elements listed in Table (I), where g is the gauge coupling, s = (p 1 +p 2 ) 2 , t = (p 1 −p 3 ) 3 and u = (p 1 −p 4 ) 2 denote the usual Mandelstam variables and C F = c − 1 = 8 denote the group theoretical factors. However, due to the enhancement of soft t, u-channel gluon and quark exchanges, the vacuum matrix elements in Table I give rise to divergent scattering rates, which inside the medium are regulated by incorporating screening effects via the insertions of the Hard-Thermal loop (HTL) self-energies, as discussed in detail in [51]. Even though it should in principle be possible to include the full HTL self-energies in the calculation of the in-medium elastic scattering matrix elements (at least for homogenous and isotropic systems), this would represent yet another significant complication as the corresponding expressions would have to be re-evaluated numerically at each time step and we did not pursue this further. Instead, we follow previous works [4,7,46,50], and incorporate an effective isotropic screening, where soft t− and u− channel exchanges are regulated by screening masses m 2 D and m 2 Q f for different species of internal exchange particles, by replacing t and u in the singly and doubly underlined expressions in Table I with where is the spatial momentum of the exchanged particle, and the parameters ξ g = e 5 6 2 √ 2 , ξ q = e √ 2 have been determined in [50] by matching to leading order HTL results. Based on the above expressions for the collision integrals and scattering matrix elements, the corresponding integrals for the discretized moments C 2↔2 a [n](i p , j θ , k φ , t), is then calculated at each time step by performing a Monte-Carlo sampling described in detail in Appendix B 1.

Inelastic Collisions
Within our effective kinetic description, we include all leading order inelastic scattering processes between quarks and gluons, where following previous works [4,7,46,50] the relevant in-medium scattering matrix elements are determined within the formalism of Arnold, Moore and Yaffe [51], by solving an integro-differential equation for the effective collinear emission/absorption rates to take into account account coherence effects associated with the Landau-Pomeranchuk-Migdal (LPM) effect [55][56][57].

a. Collision Integral
Generally, the inelastic collision integral for particle "a" with momentum ⃗ p 1 participating in the splitting process a → b, c (p 1 ↔ p 2 , p 3 ) and the inverse joining process where dΠ a↔bc 1↔2 and dΠ ab↔c 1↔2 denote the measures where again "±" provides a Bose enhancement (+) for gluon and Fermi blocking (−) for quarks. Since for ultra-relativistic particles, the 1 ↔ 2 processes require collinearity to be kinematically allowed, the collision integral in Eq.(13) can be recast into an effectively one-dimensional collinear process Based on the formalism of AMY [51], the effective inelastic rate can be expressed in the following factorized form where the matrix element for the collinear splitting is expressed in terms of the leading-order QCD splitting functions (DGLAP) [59][60][61], Secondly, the factor ∫ encodes the relevant aspects of the current-current correlation function inside the medium, and satisfies the following integral equation for particles a, b, c carrying momentum where m 2 ∞,a , m 2 ∞,b , m 2 ∞,c denote the asymptotic masses of particles a, b, c, i.e m 2 ∞,g = m 2 D 2 for gluons and m 2 ∞,q f = 2m 2 Q f for quarks, C R a , C R b , C R c denote the Casimir of the representation of, i.e. C R q = C F for quarks and C R g = C A for gluons, and dΓ el d 2 ⃗ q denotes the differential elastic scattering rate stripped by its color factor, which is given by 1 We note that several different notations for the rate dΓ a bc dz exist in the literature, and we refer to the Appendix of Ref. [58] for a comparison of different notations. at leading order [51].
Self-consistent solutions to Eq. (19) can be efficiently constructed in impact parameter space, i.e. by performing a Fourier transformation w.r.t. to ⃗ p b (see e.g. [62]), and resum the effects of multiple elastic scatterings during a single emission/absorption. Since the effective inelastic rates depend on the kinematic variables p, z as well as the time dependent medium parameters , in practice we tabulate the rates as a function of p, z for different values of T * , m 2 D , m 2 such that for small variations T * , m 2 D , m 2 Q f which occur in every time step we interpolate between neighboring points, whereas for larger variations of T * , m 2 D , m 2 Q f which occur over the course of many time steps the entire database gets updated. Similar to the elastic scattering processes, the discretized versions C 1↔2 a [n](i p , j θ , k φ , t) of the inelastic collision integrals in Eq. (16) are then calculated using a Monte-Carlo sampling, as described in more detail in Appendix B 2.
Even though we will always employ Eq. (19) to calculate the effective inelastic rates in our numerical studies, it proves insightful to briefly consider the two limiting cases where the formation time t form ∼ 1 δE (p,z) (⃗ p b ) of the splitting is small or large compared to inverse of the (small angle) elastic scattering rate 1 Γ el ∼ 1 (g 2 T * ) and closed analytic expressions for the effective inelastic rates can be obtained. We first consider the limit of small formation times, commonly referred to as the Bethe-Heitler regime [63], where radiative emissions/adsorptions are induced by a single elastic scattering and Eq. (19) can be solved perturbatively (see e.g. [32]), yielding where ,a , such that the effective inelastic rate is essentially determined by the small angle elastic scattering rate (∼ g 2 T * ).
Since the typical transverse momentum acquired in a single scattering is ⃗ p 2 b ∼ m 2 D the validity of this approximations requires the formation time t form ∼ 2pz(1−z) m 2 D to be small compared to the mean free time between small angle scatterings t mfp ∼ Γ −1 el ∼ 1 g 2 T * , giving rise to a characteristic energy scale ω BH = m 2 D g 2 T * such that for 2pz(1 − z) ≲ ω BH radiative emissions/adsorptions typically occur due to a single elastic scattering. Conversely, for 2pz(1 − z) ≳ ω BH the radiative emission/adsorption occurs coherently over the course of many elastic scatterings, leading to the famous Landau-Pomeranchuk-Migdal suppression [55][56][57] of the effective inelastic interaction rate. Specifically, in the high-energy limit 2pz(1−z) ≫ ω BH , the effective rate can be approximated as [32,58,64] , where in contrast to Eq. (22) the effective rate is determined by the formation time t −1 form ∼ q 2pz(1−z) of the splitting/merging rather than the elastic scattering rate.

III. Chemical equilibration of near-equilibrium Systems
Before we address kinetic and chemical equilibration of non-abelian plasmas which are initially far-from equilibrium, we will address the conceptually simpler case of studying the chemical equilibration of systems, where initially there is only one species of particles present. While it is conceivable that such kind of states could be created in a cosmological environment, whenever the QCD sector is selectively populated via the coupling to e.g. the standard model Higgs or other BSM particles, our primary goal is to understand and characterize the dynamics underlying chemical equilibration of the QGP, and we do not claim relevance to any particular physics application. We will for simplicity assume that, e.g. due to the interaction with other non-QCD particles, the particle species that is present initially is already in thermal equilibrium at a given temperature T 0 and chemical potential µ 0 , such that over the course of the chemical equilibration process the energy of the dominant species needs to be redistributed among all QCD degrees of freedom, until eventually the final equilibrium state with a different temperature T eq and chemical potential µ eq is reached.
Since the leading order kinetic description of massless QCD degrees of freedom is manifestly scale invariant, we can express the relevant momentum and time scales in terms of an arbitrary chosen unit. Naturally, for this kind of investigation, we will express our results in terms of the final equilibrium temperature T eq and chemical potential µ eq , such that the corresponding estimates of the physical time scales can be obtained by evaluating the expressions for the relevant temperatures and densities. Even though we employ a leading order weak-coupling description, we will investigate the behavior for different values of the QCD coupling strength 2 , typically denoted by the t'Hooft coupling λ = g 2 N c , and frequently express the dependence on the coupling strength in terms of macroscopic quantities, such as the shear-viscosity to entropy density ratio η s ∼ 1 g 4 [26,27].

A. Chemical equilibration at zero density
We first consider the case of chemical equilibration at zero (net-) density of the conserved u, d, s charges, where the systems features equivalent amounts of quarks and antiquarks, resulting in zero chemical potentials for all quark flavors. We distinguish two cases, where in the first case the system initially features a thermal distribution of gluons, without any quarks or antiquarks present at initial time, whereas in the second case the system is initially described by the same distribution of quarks/antiquark for all flavors, without gluons present in the system. Specifically, for the first case with thermal gluons only, we have where due to energy conservation, the initial parameter T 0 can be related to thermal equilibrium temperature T eq by ν g π 2 30 T 4 0 = (4ν g + 7ν q N f ) π 2 120 T 4 eq according to Eq. (2). Similarly, for the second case where only quarks/antiquarks are initially present in the system, we have and the initial parameter T 0 has the following relation to final equilibrium temperature T eq by ν q N f 7π 2 120 T 4 0 = (4ν g + 7ν q N f ) π 2 120 T 4 eq according to Eq. (2). Since the final equilibrium temperature T eq is a constant scale, it is then natural to express other scales in terms of T eq , or alternatively in terms their corresponding equilibrium values, such as m 2 D (T eq ), m 2 Q (T eq ), ⋯. Besides providing a reference scale for static equilibrium quantities, the inverse of the equilibrium temperature ∼ 1 Teq also provides a natural time scale for the evolution of the system, and it is convenient to express the time evolution in units of the near-equilibrium kinetic relaxation time  p/T eq t=0 t=10 FIG. 1: Evolution of gluon fg(t, p) and quark fq(t, p) distribution for gluon dominated initial conditions (λ = 1) at different times 0 ≤ t ≤ 2τ R expressed in units of the equilibrium relaxation time τ R in Eq. (27). Spectra of anti-quarks fq(t, p) are identical to the spectra of quarks fq(p) at zero density and not depicted in the figure. where η s is the constant shear viscosity to entropy density ratios, with η s ≃ 1900, 35, 1 for t'Hooft couplings λ = g 2 N c = 4πα s N c = 0.1, 1, 10 [50]. 3

Spectra Evolution
We first investigate the evolution of the phase-space distribution of quarks and gluons over the course of the chemical equilibration of the QGP. We present our results in Figs. 1 and 2, where we depict the evolution of the spectra of quarks and gluons for initially gluon (Fig. 1) and quark (Fig. 2) dominated systems.
Starting with the evolution of the gluon dominated system in Fig. 1, one observes that the gluon spectrum only varies modestly over the course of the chemical equilibration of the system, such that throughout the evolution the spectrum can be rather well described by an p/T eq t=0 t=10 FIG. 2: Evolution of gluon fg(t, p) and quark fq(t, p) distribution for quark/anti-quark dominated initial conditions (λ = 1) at different times 0 ≤ t ≤ 2τ R expressed in units of the equilibrium relaxation time τ R in Eq. (27). Spectra of antiquarks fq(t, p) are identical to the spectra of quarks fq(p) at zero density and not depicted in the figure.
effectively thermal distribution f g (⃗ p, t) ≃ 1 exp(p Tg(t))−1 , with a time dependent temperature T g (t), decreasing monotonically from the initial value T g (t = 0) = T 0 to the final equilibrium temperature T g (t → ∞) = T eq . Due to soft gluon splittings g → qq and elastic quark/gluon conversion gg → qq, the quark/antiquark spectra quickly built up at soft scales p ≲ T eq , as can be seen from the spectra at early times (t ≪ τ R ) in the bottom panel of Fig. 1. The quark/antiquark follows a power-law spectrum f q q (⃗ p, t) ∝ 1 p 2 associated with Bethe-Heitler spectrum. While the production of quark/antiquark at low momentum continues throughout the early stages of the evolution, the momentum of previously produced quarks/antiquarks increases due to elastic interactions, primarily qg ↔ qg andqg ↔qg scattering, such that by the time t ≃ 0.5τ R the spectrum of produced quarks/antiquarks extends all the way to the temperature scale p ∼ T eq and eventually approaches equilibrium on a time scale on the order one to two times the kinetic relaxation time τ R . Similar behavior can be observed for the quark/antiquark dominated scenario, which is depicted in Fig. 2. While quarks/antiquarks feature approximately thermal spectra f q q (⃗ p, t) ≃ 1 exp(p Tq(t))+1 , gluons are initially produced at low momentum mainly due to the emission of soft gluon radiation q → gq, which at early times (t ≪ τ R ) gives rise to a power law spectrum f g (⃗ p, t) ∝ 1 p 3 associated with the Bethe-Heitler spectrum. Subsequently, elastic and inelastic processes lead to a production of gluons with momenta p ∼ T eq until the system approaches equilibrium on a time scale on the order of the kinetic relaxation time τ R .

Collision Rates
While the evolution of the spectra in Figs. 1 and 2 provides an overview over the chemical equilibration process, we will now investigate how the individual QCD processes contribute to the evolutions of the gluon and quark/antiquark spectra in Figs. 1 and 2. We provide a compact summary of our findings in Figs. 3 and 4, where we present result for the collision rates for initially gluon dominated (Fig. 3) and initially quark dominated scenarios (Fig. 4). Different columns in Figs. 3 and 4 show the collision rates of individual processes at the initial time t = 0, at an intermediate time t = 0.1τ R and near-equilibrium at time t = 0.5τ R . We note that due to the zero net-density of u, d, s quarks, the quark and antiquark collision rates in Figs. 3 and 4 are identical and briefly remind the reader, that according to our convention in Eq. (6), positive contributions to the collision rate represent a number loss and negative collision rates exhibit a number gain for the specific particle. a. Gluon dominated scenario Starting with the collision rates for the gluon dominated scenario in Fig. 3, one observes that at initial time t = 0, the gluon splitting process g → qq shown by the dark blue curve is dominating the production of quarks/antiquarks. By comparing the collision rates for quarks and gluons, one finds that gluons with momenta p ≃ 1 − 2 T eq copiously produce soft quarks/antiquarks at low momenta p ≪ T eq . Since the individual splittings are typically asymmetric with z(1 − z) ≪ 1, the energy of thermal gluons is re-distributed to soft quarks/antiquarks, and the splittings fall into the Bethe-Heitler regime as typically pz(1 − z) ≲ ω BH ∼ T eq . In addition to the inelastic splitting, elastic conversion processes gg → qq shown as a lime curve evenly redistribute the energy of gluons with momenta p ≃ T eq into quarks/antiquarks at an intermediate scale p ≃ T eq . Due to the absence of quarks and antiquarks at initial time, the contributions of all other processes involving quarks/antiquarks in the initial state vanish identically at initial time, as do the collision rates for processes involving only gluons due to the detaily balanced in the gluon sector.
Subsequently, as quarks/antiquarks are produced at low momenta, additional scattering processes involving quarks/antiquarks in the initial state become increasingly important, as can be seen from the second column of Fig. 3, where we present the collision rates at t = 0.1τ R . While the rate of the initial quark/antiquark production processes g → qq, gg → qq decrease, as the corresponding inverse processes qq → g, qq → gg start to become important, elastic scattering of quarks and gluons qg → qg (orange curve) and gluon absorption gq → q (light blue curve) become of comparable importance. Specifically, in each of these processes, the previously produced quarks/antiquarks at low momentum p ≪ T eq gain energy via elastic scattering or absorption of a gluons, resulting in an increase of the spectrum for p ≳ 1.5T eq . By inspecting the collision rates for gluons in the top panel of Fig. 3, one observes that the depletion of soft gluons (p ≪ T eq ) due to gluon absorption by quarks gq → q is primarily compensated by the emission of soft gluon radiation due to the g → gg process (black curve). Beside the aforementioned process, the elastic scattering of gluons gg → gg (red curve) also plays an equally important role in re-distributing energy among gluons, clearly indicating that over the course of the chemical equilibration process the gluon distribution also falls out of kinetic equilibrium.
Eventually, the chemical equilibration process proceeds in essentially the same way, until close to equilibrium the collision rates of all processes decrease as the corresponding inverse processes start to become of equal importance, as seen in the right column of Fig. 3 where we present the collision rates at t = 0.5τ R . By the time t ≃ 2τ R , which is no longer shown in Fig. 3, all the collision rates decrease by at least one order of magnitude as the the system gradually approaches the detailed balanced chemical and kinetic equilibrium state.
b. Quark/antiquark dominated scenario Next we will analyze the collision rates in the quark/antiquark dominated scenario shown in Fig. 4 and compare the underlying dynamics to the gluon dominated scenario in Fig. 3.
Starting from the collision rates at initial time shown again in the left panel, one finds that in addition to quark/antiquark annihilation via elastic qq → gg (lime) and inelastic qq → g (dark blue) processes, soft gluons are copiously produced by q → gq Bremsstrahlungs processes initiated by hard quarks/antiquarks with momenta p ≳ 3T eq . Noteably the q → gq process also leads to the re-distribution of the energy of quarks/antiquarks from momenta p ≳ 3T eq , to lower momenta p ≲ 3T eq ; however the negative collision rate for the q → gq process partially cancel against the positive contribution from qq → g processes, such that there is effectively no increase/decrease of the quark/antiquark distributions at very low momenta p ≪ T eq . Similar to the processes involving only gluons at t = 0 in Fig. 3, processes involving only quarks and antiquarks (green, pink) in Fig. 4 vanish identically at t = 0 due to cancellations of gain and loss terms in the statistical factor, while other processes gg → gg, qg → qg, g → gg are exactly zero due to the absence of gluons in the initial state. By comparing the collision integrals for quarks and gluons in Figs. 3 and 4, one also observes that inelastic processes are initially much more dominant for the quark/antiquark dominated  scenario in Fig. 3 as compared to the gluon dominated scenario in Fig. 4.
Similarly to the evolution in the gluon dominated scenario, the energy of the soft gluons produced in the previous stage increases through successive elastic and inelastic interactions, as can be seen from the middle column of Fig. 4, where we present the collision rates at the intermediate time t = 0.1τ R for the quark dominated case. By inspecting the collision rates for gluon in more detail, one finds that quark-gluon scattering qg → qg (orange) as well as gg → g (black) are the dominant processes that increase the number of hard p ≳ T eq gluons. Elastic scattering between gluons gg → gg (red) plays a less prominent role for the evolution of the gluons, as do elastic (green) qq → gg and inelastic (dark-blue) qq → g conversions. With regards to the collision rates for quarks and antiquarks, one finds that elastic qq → gg (lime) and inelastic qq → g (dark blue) annihilation processes as well as q → gq Bremsstrahlungs processes continue to deplete the number of hard quarks/antiquarks. However at this stage of the evolution, qg → qg scattering processes (orange) also lead to an efficient energy transfer from quarks to gluons, depleting the number of hard quarks p ≳ 2T eq in the system. While the non-vanishing quark/antiquark scattering rates (green, pink) reveal slight deviations of quark/antiquark spectra from kinetic equilibrium, these processes clearly have a subleading effect.
Subsequently, the evolution continues along the same lines as illustrated in the right column for t = 0.5τ R , with the collision rates of all processes decreasing as the system approaches kinetic and chemical equilibrium.

Scale Evolutions
Beyond the characterization of the microscopic processes in terms of spectra and collision rates, it is instructive to investigate the evolution of the characteristic scales m 2 D , m 2 Q and T * defined in Sec. II A, which further provides a compact way to compare the time scales of the chemical equilibration process at different coupling strength. Corresponding results are presented in Fig. 5, where we compare the evolution of the various scales for quark and gluon dominated initial condition at two different coupling strengths λ = 1, 10. By taking into account the corresponding change in the equilibrium relaxation rate τ R (c.f. Eq. (27)), one finds that the time evolutions of the various scales are quite similar, and rather insensitive to the coupling strength, such that by the time t = 1 ∼ 2 τ R all relevant dynamical scales are within a few percent of their equilibrium values.
During the earlier stages, t ≲ τ R , some interesting patterns emerge in the evolution of m 2 D , m 2 Q and T * , which can be readily understood from considering the evolution of the spectra in Figs  different effects that quarks and gluons have on each of these quantities. Since the occupancy of soft quarks is always limited to below unity, soft gluons contribute more significantly to in-medium screening, such that in the gluon dominated case the screening masses m 2 D and m 2 Q in Fig. 5 decrease monotonically, whereas in the quark dominated case on observes a monotonic increase of the same quantities. The effective temperature T * which characterizes the strength of elastic interactions inside the medium, drops throughout the chemical equilibration process for the gluon dominated case, whereas for quark dominated initial conditions, the evolution of T * shows a non-monotonic behavior featuring a rapid initial drop followed by a gradual increase of T * towards its equilibrium value. By careful inspection of the spectra in Fig. 2, one finds that this rather subtle effect should be attributed to the effects of Bose enhancement and Fermi suppression in the determination of T * . Besides the evolution of the characteristic scales m 2 D , m 2 Q and T * , it is also important to understand how the overall energy is shared and transferred between quark and gluon degrees of freedom over the course of the evolution. A compact overview of the energy transfer during the chemical equilibration process is provided in Fig. 6, where we show the evolution of the energy density of gluons as well as quarks and antiquark for the two scenarios. Starting from a rapid energy transfer at early times, the flattening of the individual energy densities towards later times eventually indicates the approach towards chemical equilibrium. Even though the evaluation of an exact chemical equilibration time depends on the quantitative criterion for how close the energy densities (or other scales) are compared their equilibrium values are, the figures still speak for themselves indicating the occurrence of chemical equilibration roughly at the same time scale as kinetic equilibration, with subject to mild variations for the two different coupling strengths.

B. Chemical equilibration of finite density systems
So far we have investigated the chemical equilibration of charge neutral QCD plasmas, and we will now proceed to study the chemical equilibration process of QCD plasmas at finite density of the conserved u, d, s charges, featuring an excess of quarks to antiquarks (or vice versa). Since a finite net charge density of the system can only be realized in the presence of quarks/anit-quarks, we will focus on quark dominated initial conditions and modify and T * (t) (black) during the chemical equilibration process for the quark dominated (fine dashed) and gluon dominated (long dashed) scenarios at two different coupling strengths λ = 1 (lighter colors) and λ = 10 (darker colors). Scales are normalized to their respective equilibrium values, while the evolution time t is normalized to the equilibrium relaxation time τ R in Eq. (27) in order to take into account the leading coupling dependence.
the corresponding initial conditions as where for simplicity, we consider equal densities of u, d and s quarks. Similar to Eqs. (25) and (26), the initial parameters T 0 , µ 0 can be related to corresponding equilibrium temperature T eq , and equilibrium chemical potential µ eq via the Landau matching procedure in Eq. (2). Due to energy and charge conservation, T eq and µ eq , then determine the final equilibrium state of the system, and we will characterize the different amounts of net charge in the system in terms of the ratio µ eq T eq , with µ eq T eq = 0 corresponding to the charge neutral plasma considered in the previous section. When comparing the evolutions at different coupling strengths, we follow the same procedure as discussed above and express the evolution time in units of the kinetic relaxation time which in accordance with the last equality reduces to the same expression for a charge neutral system (µ = 0) in Eq. (27). The effective temperature is evaluated as = T eq . Since we did not explicitly determine the dependence of the shear-viscosity η(T, µ) on the chemical potential µ for all coupling strengths λ, we will approximate η(T,µ)T by the corresponding value of η(T,µ=0) s at vanishing density of the conserved charges, which are quoted below Eq. (27).

Spectra Evolutions
We follow the same logic as in the charge neutral case and first investigate the evolution of the spectra of quarks, antiquarks and gluons, which is presented in Fig. 7 for the chemical equilibration of a system with quark chemical potentials µ eq T eq = 2.5. Similar to the quark dominated scenario at zero density, we find that the spectra for quarks and antiquarks are always close µ eq /T eq =2.50  to a thermal distribution with the expected moderate deviation at intermediate times. Specifically, the antiquark spectra in the low momentum sector p ≲ 0.3T eq are depleted at intermediate times t ≲ 0.5τ R , due to elastic qq → gg and inelastic qq → g conversions. Besides quark/antiquark annihilations, the radiative emission of gluons due to q → qg andq →qg processes leads to a rapid population of the soft gluon sector seen in the top panel of Fig. 7. By comparing the results in Figs. 2 and 7, one finds that the soft gluon sector builds up even more rapidly at finite density as compared to zero density, such that already by the time t = 10 −3 τ R , the gluon distribution at low momentum p ≲ 0.1T eq features a quasi-thermal spectrum f (p ≪ T eq ) ≃ T eq p, whereas the high momentum tail is yet to be populated. Eventually on a time scale of t ≃ 1.5τ R , a sufficiently large number of hard gluons has been produced and the spectra of all particle species relax towards equilibrium, such that significant deviations from the thermal distributions are no longer visible for t = 1.5τ R in Fig. 7.

Collision Rates
Beyond the evolution of the of the spectra, it again proves insightful to investigate the collision rates in Fig. 8 in order to identify the microscopic processes that drive chemical and kinetic equilibration of gluons, quarks and antiquarks at different stages.
Similar to the results for the charge neutral case in Fig. 4, the initial gluon production in Fig. 8 is still dominated by soft radiation q → gq +q →qg (light blue), with even more substantial contributions due to the larger abundancies of quarks. Conversely, the gluon production from elastic qq → gg (lime) and inelastic qq → g (dark blue) quark/antiquark annihilation processes is markedly suppressed due to the shortage of antiquarks. Similar differences between the evolution at zero and finite density can also be observed in the collision rates for quarks and antiquarks, where in the case of the quark the emission of gluon radiation leads to a depletion of the hard sector p ≳ 3T eq , along with an increase of the population of softer quarks with typical momenta p ∼ T eq . While elastic qq → gg (lime) and inelastic qq → g (dark blue) processes initially contribute at a much smaller rate, such that the inelastic q → qg process dominates the evolution of the quarks, a manifestly different picture emerges for the collision rates of antiquarks. Due to the large abundancies of quarks, elastic qq → gg (lime) and inelastic qq → g (dark blue) quark/antiquark annihilation initially occur at essentially the same rate as gluon radiation off antiquarksq →qg (light blue), resulting in a net-depletion of the antiquark sector across the entire range of momenta. Besides the aforementioned processes, the collision rates of all other processes vanish identically at initial time for all particle species due to cancellations of the statistical factors.
Subsequently, for t = 0.1τ R depicted in the central column of Fig. 8 a variety of different processes becomes relevant as soft gluons have been copiously produced during the previous evolution. Besides the processes involving quark-gluon interactions q → gq (light blue), qg → qg (orange), inelastic absorptions of soft gluons gg → g (black) also have an important effect for the thermalization of the gluon sector, whereas elastic scattering of gluons gg → gg (red) as well as elastic qq → gg (lime) and inelastic qq → g (dark blue) quark/antiquark annihilation processes are clearly subleading. By comparing the results at zero and finite density in Figs. 4 and 8, one further notices an increment of the gg → g collision rates, indicating a more rapid gluon production from quarks at finite density, consistent with the observations of the spectra in Figs. 2 and 7. Due to the fact that at finite density there are more quarks present in the system, the collision rates for quarks are generally larger compared to the zero density case. Nevertheless, the underlying dynamics remains essentially the same as compared to the zero density case, with gluon radiation q → gq (light blue) and quark-gluon scattering qg → gq providing the dominant mechanisms to transfer energy from hard quarks to softer gluons. Due to the larger abundance of quarks at finite density, elastic scattering processes qq → qq involving quarks of the same (green) and different flavors (pink), also play a more prominent  role in restoring kinetic equilibrium in the quark sector, while they were more or less negligible at zero density. Surprisingly small changes appear in the collision rates for antiquarks between the initial time t = 0 and t = 0.1τ R , where at later times the inelasticq →qg process becomes suppressed due to the fact the inverse process of absorbing a soft gluonqg →q becomes increasingly likely. Similarly, elastic scattering processesqg →q (orange) between antiquarks and gluons also contribute to the energy transfer from the antiquark to the gluon sector.
Eventually for t = 0.5τ R , the energy transfer from quarks to gluons due to elastic qg → qg (orange) and inelastic q → gg (light blue) becomes smaller and smaller, so do the collision rates for inelastic gluon absorptions g → gg (black) and elastic scatterings between quarks/antiquarks (pink and green), which are primarily responsible for restoring kinetic equilibrium in the gluon and quark sectors. Beyond the time scales shown in Fig. 8, the evolution of the system continues in essentially the same way, with continuously collision rates decreasing until eventually gluons, quarks and antiquarks all approach their respective equilibrium distribution.

Scale Evolutions
Now that we have established the microscopic processes underlying the chemical equilibration of finite density, we again turn to the evolution of the dynamical scales m 2 D , m 2 Q and T * , which serve as a reference to determine the progress of kinetic and chemical equilibration. We present our results in Fig. 9, where we compare the evolution of the dynamical scales in systems with a different amount of net-charge density, as characterized by the ratio µ eq T eq = 0, 0.14, 1.34, 2.5 of the equilibrium chemical potential over the equilibrium temperature. By comparing the evolution of the various quantities in Fig. 9, one observes that for larger chemical potentials m 2 D , m 2 Q as well as T * are generally closer to their final equilibrium values over the course of the entire evolution. While the smaller deviations of m 2 D , m 2 Q and T * can partly be attributed to the fact that in the finite density system the initial values for these quantities are already closer to the final equilibrium value, it also appears that the ultimate approach towards equilibrium occurs on a slightly shorter time scale. We attribute this to the fact that for larger values of µ eq T eq at a fixed temperature, the system features a larger energy density (c.f. Eq. (2)), which should effectively speed up the various collision processes.
Similar phenomena can also be observed in Fig. 10, where we present the evolution of the energy and num- µ eq /T eq =2.50 µ eq /T eq =1.34 µ eq /T eq =0.14 µ eq /T eq =0.00 ber density of gluons, quarks and antiquarks over the course of the chemical equilibration process at different densities µ eq T eq = 0, 0.14, 1.34, 2.5. While initially there is always a rapid production and energy transfer to the gluon sector, the flattening of the curve at later times show the relaxation towards chemical equilibrium, which occurs roughly on the same time scale as the kinetic equilibration of the dynamical scales m 2 D , m 2 Q and T * . By comparing the results for different µ eq T eq , one again observes that the chemical equilibration happens slightly earlier for larger chemical potential, consistent with the observations from spectra in Fig. 2, Fig. 7, collision rates in Fig. 4, Fig. 8 and from the scale evolutions in Fig. 9. Nevertheless, we believe that at least for the range of µ eq T eq considered in Fig. 9, our estimate of the kinetic and chemical equilibration time scales in Eq. (29), remains valid also at finite density.

IV. Equilibration of Far-From-Equilibrium Systems
We will now analyze the equilibration process of QCD systems which are initially far from equilibrium. By focusing on systems which are spatially homogeneous and isotropic in momentum space, we can distinguish two broad classes of far-from equilibrium initial states which following [32,36] can be conveniently characterized by considering the initial average energy per particle ⟨p⟩ 0 in relation to the equilibrium temperature T eq of the sys- 0.14 0.14 Energy density fraction: e a /e t/τ R µ eq /T eq =2.50 µ eq /T eq =1.34 µ eq /T eq =0.14 µ eq /T eq =0.00 0.14 0.14 Number density fraction: n a /n eq t/τ R µ eq /T eq =2.50 µ eq /T eq =1.34 µ eq /T eq =0.14 µ eq /T eq =0.00 Gluon Quark Antiquark tem. Specifically, for far-from equilibrium initial states, we can consider a situation where the average energy per particle is initially much smaller than the equilibrium temperature, i.e. ⟨p⟩ 0 ≪ T eq , such that the energy is initially carried by a large number f 0 ≫ 1 of low momentum gluons. Such over-occupied initial states typically appear as a consequence of plasma instabilities [36,41,65] and they also bear some resemblance with the saturated "Glasma" initial state created in high-energy collisions of heavy nuclei [40,[66][67][68][69][70][71], although the detailed properties of this state are quite different as the system is highly anisotropic and rapidly expanding in the longitudinal direction as discussed in more detail in Sec. V. While for ⟨p⟩ 0 ∼ T eq , the system is in some sense close to equilibrium and one would naturally expect kinetic and chemical equilibration to occur on the time scales of the equilibrium relaxation time ∼ τ R , there is a second important class of far-from equilibrium initial states corresponding to under-occupied states. In under-occupied systems the average energy per particle is initially much larger then the equilibrium temperature ⟨p⟩ 0 ≫ T eq , such that the energy is initially carried by a small number f 0 ≪ 1 of highly energetic particles, as is for instance the case for an ensemble of high-energy jets. While earlier works [40,46] have established the equilibration patterns of such systems for pure glue QCD, we provide an extension of these studies to full QCD with three light flavors, as previously done for over-occupied systems in [50].

A. Equilibration of Over-occupied Systems
We first consider over-occupied systems characterized by a large occupation number f 0 ≫ 1 of low-energy ⟨p⟩ 0 ≪ T eq gluons, 4 and we may estimate the energy density of the over-occupied system as e 0 ∼ f 0 ⟨p⟩ 4 0 . Since the total energy density is conserved, we have e eq = e 0 , such that with e eq ∼ T 4 eq the final equilibrium temperature T eq ∼ f 1 4 0 ⟨p⟩ 0 ≫ ⟨p⟩ 0 is much larger than the average initial momentum ⟨p⟩ 0 . Due to this separation of scales, energy needs to be re-distributed from low momentum to high momentum degrees of freedom, which as will be discussed shortly is achieved via a direct energy cascade from the infrared to ultraviolet in momentum space.

Theoretical Aspects
Due to the large population of low momentum gluons, interaction rates for elastic and inelastic processes are initially strongly enhanced, such that e.g. the large angle elastic scattering rate Γ el ∼ g 2 T * m 2 D ⟨p⟩ 2 is initially much larger than in equilibrium Γ 0 el ∼ g 4 f 2 0 ⟨p⟩ ≫ Γ eq el ∼ g 4 T eq . Even though the time scale for the actual equilibration process is eventually controlled by the equilibrium rate ∼ 1 Γ eq el , the system will therefore encounter a rapid memory loss of the details of the initial conditions on a time scale ∼ 1 Γ 0 el , and subsequently spend a significant amount of time in a transient non-equilibrium state, where the energy transfer from the infrared towards the ultraviolet is accomplished.
Since the dynamics remains gluon dominated with f g ≫ 1 ≥ f q,q all the way until the system eventually approaches equilibrium, one should expect that the evolution of the over-occupied Quark-Gluon plasma follows that of pure-glue QCD, where it has been established [38,39,41,45,53], that for intermediate times 1 Γ 0 el ≪ t ≪ 1 Γ eq el , the evolution of the gluon spectrum follows a self-similar scaling behavior of the form where t 0 ≃ 1 Γ 0 el , ⟨p⟩ 0 are the characteristic time and momentum scales, f 0 is the initial occupancy and f S (x) is a universal scaling function up to amplitude normalization and we adopt the normalization conditions f S (x = 1) = f ′ S (x = 1) = 1. We note that the emergence of self-similar behavior as in Eq. (32), is by no means unique to QCD, and in fact constitutes a rather generic pattern in the equilibration of far-from-equilibrium quantum systems, with similar observations reported in the context of relativistic and non-relativistic scalar field theories [42,72]. Specifically, the scaling exponents α = −4 7, β = −1 7 follow directly from a dimensional analysis of the underlying kinetic equations [36][37][38]41], and describe the energy transport from the infrared towards the ultra-violet due to a direct energy cascade [73].
Based on Eq. (32), we can further estimate the evolutions of some physical quantities knowing that gluon are dominant f g ≫ 1 ≥ f q,q in the self-similar scaling regime. In particular, the average momentum ⟨p⟩ increases as a function of time according to while the typical occupancies of hard gluons decrease as Similarly, one finds that the screening mass decreases, such that the system dynamically establishes a separation between the soft (∼ m D ) and hard (∼ ⟨p⟩) scales over the course of the self-similar evolution [36,45]. Since the effective temperature also decreases according to (f g (p, t) ≫ 1) the large-angle elastic scattering rate Γ el (t) ∼ g 2 T * m 2 D ⟨p⟩ 2 ∼ g 4 f 2 0 ⟨p⟩ 0 (t t 0 ) −1 decreases over the course of the selfsimilar evolution and eventually becomes on the order of the equilibrium rate Γ el (t) ∼ g 4 T eq at the same time t t 0 ∼ f 7 4 0 when the occupancies of hard gluons f t, ⟨p⟩(t) become of order unity, and the average momentum ⟨p⟩(t) becomes on the order of the equilibrium temperature T eq ∼ ⟨p⟩ 0 f 1 4 0 , indicating that the energy transfer towards the ultra-violet has been accomplished and gluons are no longer dominant for eq . Beyond this time scale, the system can be considered as close to equilibrium, and should be expected to relax towards equilibrium on a time scale on the order of the kinetic relaxation time τ R ∼ g −4 T −1 eq , which is parametrically of the same order as the time it takes to accomplish the energy transfer towards the ultra-violet.
We first consider an over-occupied system with a relatively large scale separation ⟨p⟩ 0 T eq =0.2 at weak coupling λ=0.1, and investigate the evolutions of the spectra of quarks and gluons depicted in the top panel of Fig. 11. Starting from a large phase-space occupancy of soft gluons, the initial spectra undergo a quick memory loss at very early times and then gradually evolve into harder spectra through a direct energy cascade, pushing the low momentum gluons towards higher momenta. In order to illustrate the self-similarity of this process, we follow previous works [32,38] and show re-scaled versions of the gluon spectra in the bottom panel of Fig. 11. By re-scaling the phase-space distribution as f S (x) ≃ (t t 0 ) −α f g (⟨p⟩ 0 (t t 0 ) −β x, t), and plotting it against the re-scaled momentum variable x = (t t 0 ) β p ⟨p⟩ 0 , one indeed finds that in the relevant scaling window, which corresponds approximately to times 10 −6 ≤ tλ 2 T eq ≤ 10 1 for this particular set of parameters, the spectra at different times overlap with each other to rather good accuracy, clearly indicating the selfsimilarity of the underlying process. Beside the gluons, all species of quarks/antiquarks are produced democratically over the course of the evolution from elastic gg → qq conversions and inelastic splitting g → qq processes. Generally, one finds that the quark/antiquark spectra follow the evolution of the gluon spectra, albeit due to their Fermi statistics the number of quarks/antiquarks in the system remains negligibly small compared to the overall abundance of gluons during the self-similar stage of the equilibration process.
Eventually for times t ≳ 10 2 λ 2 T eq the self-similar cascade in Fig. 11 stops as the occupancies of hard gluons fall below unity and the system subsequently approaches thermal equilibrium on time scales ∼ 10 3 λ 2 T eq for the parameters chosen in Fig. 11. It is interesting to point out, that due to the negligible abundance of quarks and antiquarks in the system, the evolution of the gluon spectra slightly overshoots the equilibrium temperature at times t ≳ 10 2 λ 2 T eq , and subsequently relaxes back towards equilibrium as the correct equilibrium abundance of quarks and antiquarks is produced along the lines of our previous discussion of gluon dominated systems in Sec. III.
Next we will discuss the evolution of the average momentum ⟨p⟩ (t),the screening mass square m 2 D (t) and the effective temperature T * (t) summarized in Fig. 12, where the upper panel shows the results for ⟨p⟩ 0 T eq = 0.2, λ = 0.1, i.e. the same parameters as in Fig. 11, while the middle and bottom panels show the results for a smaller scale separation ⟨p⟩ 0 T eq = 1, at larger values of coupling λ = 1, 10. By comparing the evolutions of the various scales with the theoretical predicted powerlaw scaling (dashed line) in the turbulent regime (c.f. Eqns. (33), (35), (36)), one finds that the scaling behavior ⟨p⟩ ∝ t 1 7 , T * ∝ t −3 7 and m 2 D ∝ t −2 7 associated with the turbulent energy transport towards the ultra-violet is indeed realized during intermediate times. Due to the large separation of scales for ⟨p⟩ 0 T eq = 0.2, λ = 0.1, the scaling window in the top panel of Fig. 12 extends over a significant period of time 10 −6 ≤ tλ 2 T eq ≤ 10 1 , consistent with the scaling of the gluon distribution observed in Fig. 11. Even though the scaling window shrinks significantly for the smaller scale separations ⟨p⟩ 0 T eq = 1 shown in the middle and bottom panels of Fig. 12, it is remarkable that the same turbulent mechanism appears to be responsible for the energy transfer even for such moderately strongly coupled systems.
Even though a significant amount of time is spent to accomplish the turbulent energy transfer, the logarithmic representation in Fig. 12 spoils the fact, that it is in fact the ultimate approach towards equilibrium which requires the largest amount of time. Beyond the investigation of the dynamical evolutions of various scales, it is therefore useful to consider the evolutions of the ratios of different scales compared to their equilibrium values, as indicators of the equilibration progress. We present our results in Fig. 13, where the upper panel shows the evolutions of the energy densities of gluons and quarks, approaching their equilibrium limits around t ≃ 1.5 − 2τ R , similar to near-equilibrium systems shown in Fig. 6. The next two panels of Fig. 13 show the screening mass square evolutions of m 2 D (t) and m 2 Q (t), which rapidly decrease at early times, an eventually approach their equilibrium values at t ≃ 1 − 1.5τ R . Similar observations also hold for the effective temperature T * (t) shown in the forth panel of Fig. 13. Due to the delayed chemical equilibration of the system, the average momentum ⟨p⟩ (t) shown in the bottom panel has a non-monotonic behavior, where the rapid increase at early times due to the direct energy cascade overshoots the equilibrium value, before ⟨p⟩ (t)'s gradual decrease at later times as energy is re-distributed between quarks and gluons, eventually approaching the equilibrium limit around t ≃ 1.5 − 2τ R .
Since in Fig. 13 the ultimate approach towards equilibrium is mostly insensitive to the initial scale separations ⟨p⟩ 0 T eq and coupling strength λ in Fig. 12 when expressed in units of the kinetic relaxation time τ R , we can estimate the equilibration time of an over-occupied system as where, as usual, the exact numerical value depends the detailed criteria chosen to define the equilibration time.

B. Equilibration of Under-occupied Systems
Next we consider the opposite limit of an underoccupied system, where the energy density is initially carried by a small number f 0 ≪ 1 of high energetic particles, with average momentum ⟨p⟩ 0 ≫ T eq . While there can be a large separation of scales, one finds that in contrast to over-occupied systems the final equilibrium temperate T eq ∼ f 1 4 0 ⟨p⟩ 0 ≪ ⟨p⟩ 0 is much smaller than the average initial momentum for under-occupied systems. Since the scale hierarchy is reversed, the thermalization process for an under-occupied system requires an energy transport from the ultra-violet to the infrared, which as we will discuss shortly will be accomplished via an inverse turbulent cascade of successive radiative emissions. While the qualitative features of this "bottom-up" thermalization mechanism have been established a long time ago [34], recent works in the context of thermalization and jet quenching studies [74][75][76] have provided a more quantitative description of the different stages and clarified the relation to turbulence. Based on our effective kinetic description of QCD, we will extend previous findings in pure glue QCD [46] to full QCD at zero and non-zero densities.

Theoretical Aspects
Before we turn to our numerical results, we briefly recall the basic features of the bottom up thermalization in QCD plasmas following the discussion in [32]. Starting from a dilute population of f 0 ≪ 1 highlyenergetic particles with ⟨p⟩ 0 ≫ T eq , elastic interactions between primary hard particles induce the emission of soft gluon radiation, which accumulates at low momenta. Due to the fact that elastic and inelastic interactions are more efficient at low momentum, the initially overpopulated soft sector eventually thermalizes on a time scale t ∼ g −4 f −1 3 0 ⟨p⟩ −1 0 , before the highly-energetic primary particles have had sufficient time to decay. Even though at this time most of the energy is still carried by the hard primaries, the soft thermal bath begins to dominate screening and scattering, such that in the final stages of bottom-up equilibration, the few remaining hard particles loose their energy to the soft thermal bath, much like a jet loosing energy to a thermal medium [32,34,46,76]. Based on recent studies [74][75][76], the energy loss of hard primaries is accomplished by a turbulent inverse energy cascade, where the hard primary quarks/gluons, undergo successive splittings until the momenta of the radiated quanta becomes on the order of the temperature T soft (t) of the soft thermal bath. Specifically, at intermediate scales T soft (t) ≪ p ≪ ⟨p⟩ 0 , the distributions of quarks/antiquarks and gluons can be expected to feature the Kolmogorov-Zakharov spectra of weakwave turbulence [75,76] which describe a scale-invariant energy flux from the ultra-violet ∼ ⟨p⟩ 0 to the infrared ∼ T soft (t), ensuring that the energy of the hard particles is deposited in the thermal medium without an accumulation of energy at intermediate scales.
Due to the energy loss of the hard primary particles, the temperature of the soft thermal bath increases until eventually the hard primaries have lost most of their energy to the thermal bath and the system approaches equilibrium. We note that due to the parametric sup- pression of inelastic rates for high-energy particles 5 Γ eq inel (⟨p⟩ 0 ) ∼ g 4 T eq Teq ⟨p⟩0 , the energy loss of the hard primaries is slow compared to the equilibration of the soft sector, such that for sufficiently large scale separations ⟨p⟩0 Teq ≫ 1 the thermalization of the system occurs on time scales t ∼ g −4 T −1 eq ⟨p⟩0 Teq , which can be significantly larger than the kinetic relaxation time τ R ∼ g −4 T −1 eq .

Bottom Up Thermalization of Quark-Gluon Plasma
When considering the dynamics of under-occupied QCD plasmas, we need to specify the initial conditions for the momentum distribution and we can further distinguish different chemical compositions of the plasma. 5 Since quasi-democratic z ∼ 1 2 splittings dominate the turbulent energy transfer [32,75], this can be seen by evaluating Eq. (24) for z ∼ 1 2 T eq =0 τλ 2 T eq =10 0 τλ 2 T eq =10 1 τλ 2 T eq =10 2 τλ 2 T eq =10 3 τλ We will limit our investigation to the following three cases, corresponding to (1) an initially under-occupied plasma of gluons, (2) an initially under-occupied plasma of quarks/antiquarks, and (3) an initially under-occupied plasma of quarks.
a. Under-occupied gluons We start by analyzing the evolutions of underoccupied gluon systems in order to provide a direct and intuitive understanding of the bottom up thermalization scenario. The evolution of the momentum spectra of quarks and gluons during the thermalization process is presented in Figs. 14, 15, 16 and 17 for weakly coupled plasmas λ = 1 with different average initial momenta ⟨p⟩ 0 T eq = 3 in Fig. 14, ⟨p⟩ 0 T eq = 10 in Fig. 15, ⟨p⟩ 0 T eq = 30 in Fig. 16 and ⟨p⟩ 0 T eq = 100 in Fig. 17. Different panels show the evolutions of the gluon distributions f g (p) and quark/antiquark distributions f q (p), while different curves in each panel correspond to different evolution times tλT eq with vertical arrows marking the characteristic Bethe-Heitler frequency at each stage of the evolution.

FIG. 19: Evolution of the energy densities (top) and average momenta (bottom) of quarks (solid) and gluons (dashed) in
an under-occupied gluon system at coupling strengths λ = 1 for different scale separation ⟨p⟩0 Teq = 3, 10, 30, 100. Energy densities and average momenta are normalized to their respective equilibrium values, while the evolution time t is normalized to τ R p 0 Teq in order to take into account the leading dependence on the initial energy ⟨p⟩0. and ⟨p⟩ 0 T eq = 100 in Fig. 17, one clearly observes that soft radiation processes g → gg and g → qq rapidly build up a large population of soft quarks and gluons with typical momenta p ≲ ω BH . Even though at early times, such as e.g. tλ 2 T eq ≪ 1 in Fig. 17, the soft sector is over-occupied and thus highly gluon dominated, one finds that for sufficiently large scale separations, the over-occupation is depleted and the soft sector thermalizes before the hard primaries loose most of their energy to the soft thermal bath. Since at intermediate scales ω BH ≪ p ≪ ⟨p⟩ 0 the emission is in the LPM regime, the spectra of gluons and quarks initially feature a characteristic power law behavior f g ∼ p −7 2 , f q ∼ p −5 2 for momenta ω BH ≪ p ≪ ⟨p⟩ 0 , associated with the single emission spectra of the g → gg and g → qq processes. Subsequently, the energy of the hard primaries is transferred to the soft thermal bath, via an inverse turbulent cascade due to multiple successive g → gg, g → qq and q → qg branchings, giving rise to the characteristic Kolmogorov-Zakharov spectrum f g q ∼ p −7 2 in both the gluon and quark sector. Since the energy injected into this cascade by the hard primaries at the scale ∼ ⟨p⟩ 0 , is transmitted all the way to the soft bath ∼ ω BH the temperature of the soft bath increases monotonically, as seen e.g. in Fig. 17, until eventually the hard primaries have lost nearly all of their energy and the system thermalizes. During the final stages of the approach towards equilibrium, a small number of hard primaries continues to loose energy giving rise to high momentum tails of the quark and gluon spectra seen for tλ 2 T eq = 10 3 in Figs. 15, 16, 17. Notably, the under-occupied system initially maintains a memory of the momentum distribution of hard primaries until the final stages of the thermalization process, which then closely resembles the mechanism of jet energy loss in a thermal medium [76].
Even for the smallest separation of scales ⟨p⟩ 0 T eq =3 shown in Fig. 14, some of the characteristic patterns of bottom up thermalization are still clearly visible, although in this case radiative emissions occur in the Bethe-Heitler regime. Nevertheless, hard gluons with momenta p ∼ ⟨p⟩ 0 still radiate soft gluons via g → gg, leading to the formation of a soft thermal spectrum of gluons at low momenta. Even though quarks/antiquarks are also produced via g → qq branching, one observes that the evolution in the quark sector is slightly slower than in the gluon sector, indicating once again that the energy transfer from gluons to quarks associated with the chemical equilibration of the system can cause a delay in the equilibration of the system. Now in order to compare the evolutions of the different systems, we again consider the evolutions of the characteristic dynamical scales m 2 D , m 2 Q , T * and ⟨p⟩. Since in accordance with the discussion in Sec. III we anticipate that, for sufficiently large scale separations, the equilibration time of the system will be delayed by a factor ⟨p⟩ 0 T eq , relative to the equilibrium relaxation time τ R , we will consider normalizing the evolution time to τ R ⟨p⟩ 0 T eq when comparing the results for different average initial momenta ⟨p⟩ 0 T eq in Figs. 18 and 19. Since the different scales m 2 D , m 2 Q , T * and ⟨p⟩, exhibit different sensitivities to the hard and soft components of the plasma, their time evolutions are actually quite different. While for scale separations ⟨p⟩ 0 ≳ 10T eq , screening masses m 2 D , m 2 Q are very quickly dominated by the soft thermal bath, and subsequently experience a strong rise as the soft bath heats up, the scale T * characterizing the strength of elastic interactions, receives significant contributions from the hard primaries at early times, before it is eventually dominated by the soft bath. Since the hard primaries carry most of the energy of the system until they eventually equilibrate, the average energy per particle ⟨p⟩ is always dominated by the hard sector, and decreases monotonically over the course of the evolution. Besides the equilibration of the various scales, it is also interesting to consider the chemical equilibration of the system in Fig. 19, where we present the energy fractions and average momenta separately for quarks and gluons. While for large scale separations, chemical equilibration in Fig. 19 occurs on the same time scales as kinetic equilibration in Fig. 18, one finds that for smaller scale separations the energy transfer from gluons to quarks requires additional time, delaying the equilibration of the system.
Generally, for scale separations ⟨p⟩ 0 T eq ≳ 10, one finds that the scaling of the evolution time with ⟨p⟩ 0 T eq , leads to comparable results for the equilibration time  the scale separations considered in our study.

FIG. 24: Evolution of the energy densities (top) and average momenta (bottom) of quarks (solid) and gluons (dashed)
in an under-occupied quark/anti-quark system at coupling strengths λ = 1 for different scale separations ⟨p⟩0 Teq = 3, 10, 30. Energy densities and average momenta are normalized to their respective equilibrium values, while the evolution time t is normalized to τ R p 0 Teq in order to take into account the leading dependence on the initial energy ⟨p⟩0.
ified, and there is no longer a significant difference between under-occupied gluon systems and under-occupied quark/antiquark systems.
By comparing the results for the evolutions of the dynamical scales m 2 D , m 2 Q , T * and ⟨p⟩ in Fig. 23 for the under-occupied quark/antiquark systems to the corresponding results for under-occupied gluons, one again observes essentially the same qualitative patterns. However, it is interesting to see, that for under-occupied systems of quarks and antiquarks, the approach towards equilibrium appears to occur on a somewhat larger time scale ≳ 0.5τ R ⟨p⟩ 0 T eq as compared to under-occupied gluon systems, where by 0.5τ R ⟨p⟩ 0 T eq all the scales m 2 D , m 2 Q , T * and ⟨p⟩ are already close to their respective equilibrium values. Based on our discussion in Sec. II B 2 b, we believe that this discrepancy at intermediate times can be attributed to the different color factors in the inelastic interactions rates for the hard primary quarks/antiquarks and gluons, as discussed in detail in the context of jet quenching in [76]. However, if one is concerned with the ultimate approach towards equilibrium, one should take into account the fact that at p/T eq tλ 2 T eq =0 tλ 2 T eq =10 0 tλ 2 T eq =10 1 tλ 2 T eq =10 2 tλ 2 T eq =10 3 tλ 2 T eq =10 4 thermal late times the quark/gluon composition is significantly modified, such that under-occupied systems of quarks and gluons can ultimately be expected to equilibrate at the same rate. Next, in order to investigate the chemical equilibrations of the under-occupied quark/antiquarks systems, we present our simulation results for the evolutions of the energy fraction of quarks and gluons, and their average momenta in Fig. 24. Interestingly, one finds that in contrast to the behavior for under-occupied gluon systems in Fig. 19, the energy fractions of quarks and gluons in the system exhibit a non-monotonic behavior. Even though initially all the energy is carried by the hard primary quarks and antiquarks, it turns out that for larger scale separations ⟨p⟩ 0 T eq =10, 30, gluons dominate the energy budget before the chemical equilibration of the system. By inspecting also the behaviors of the average momenta in the lower panel of Fig. 24, one finds that these gluons are typically soft, with the average momenta ⟨p⟩ close to the equilibrium value. We believe that this behavior can be attributed to the fact that gluon radiation dominates the energy transfer from the hard to the soft sector, such that the soft thermal bath absorbs the energy pre-dominantly in form of gluons, before the energy is eventually re-distributed among quarks and gluons.
c. Under-occupied quarks So far we have investigated the equilibrations of charge neutral systems of under-occupied gluons, quarks/antiquarks. Next we consider the equilibrations of under-occupied systems of quarks, which in accordance with Eq. (42) carry non-zero densities of the conserved u, d, s charges. Since in the presence of finite charge densities, the evolutions of quarks and antiquarks will be different, we first study the evolutions of spectra of gluons, quarks and antiquarks, which are depicted in Fig. 25 for ⟨p⟩ 0 T eff = 3 and in Fig 26 for ⟨p⟩ 0 T eff = 10. Evidently, the evolutions of the quark and gluon spectra in Figs. 25 and 26 are very similar to the quark/antiquark spectra in Figs 20 and 21 obtained in the zero density cases. However, significant differences can be observed for the evolutions of the antiquarks, as for the under-occupied systems of quarks there are no antiquarks present in the initial conditions. Instead, the population of antiquarks observed at later times is produced via gluon splittings g → qq and elastic gg → qq conversions. Hence, the evolutions of the antiquark spectra closely follow the gluon spectra, as can be seen by comparing the upper and lower panels in Figs. 25 and 26.
By comparing the evolutions of the characteristic scales for the zero and finite density systems in Figs. 23 and 27, one finds that the presence of the additional conserved charges does not significantly affect the kinetic equilibration of the system, in accordance with the finding that the evolutions of quark and gluon spectra are essentially unchanged. However, when considering the evolutions of the individual contributions of gluons, quarks and antiquarks to the energy densities in Figs. 24 and 28, one clearly observes that the chemical equilibration associated with the production of antiquarks requires a significant amount of time, with energy densities of gluons and antiquarks only approaching their equilibrium ratios for times ≳ τ R ⟨p⟩ 0 T eq .

V. Equilibration of longitudinally expanding plasmas
So far we have discussed kinetic and chemical equilibrations for homogeneous and isotropic systems. Now in order to address the early time dynamics of high-energy heavy-ion collisions, we will focus on systems which are transversely homogeneous and longitudinally invariant under Lorentz boost, but can feature an anisotropy be- Teq in order to take into account the leading dependence on the initial energy ⟨p⟩0. tween longitudinal and transverse momenta. Denoting the phase-space distribution f (x, p) for a longitudinally boost invariant and transversely homogeneous system can be conveniently expressed in the form f (x, p) = f (τ, p T , p ), where the variable p denotes the longitudinal momentum p = p T sinh(y − η) in the local restframe u µ = (cosh(η), 0, 0, sinh(η)) of the non-equilibrium plasma.
Since the system is homogeneous in transverse coordinates x T and the longitudinal rapidity η, the resulting Boltzmann equation [77]  We note that in comparison to the previous discussion of homogenous and isotropic systems, there are two important physical differences when considering plasmas which are subject to a rapid longitudinal expansion. Due to the expansion, the system will on the one hand become more and more dilute over the course of the evolution, on the other hand the longitudinal expansion tends to reduce the longitudinal momenta in the local rest frame, thereby introducing an anisotropy which can persists on large time scales. We further note that momentum space anisotropic QCD plasmas are generally expected to be unstable [78][79][80][81] due to the non-abelian analogue of the Weibel instability in electrodynamics [82]. While perturbative calculations of the one-loop self- Eq. (48) to represent the different valence quark fractions taking into account a proton to neutron fraction n p (n p + n n ) ≈ 0.4 in heavy nuclei. We present a compact summary of all the simulations performed in Tab. II, where we list the corresponding initial conditions, and coupling strength λ, along with extracted values of the ratio (µ B T ) eq at late times and the shear-viscosity ηT eff (e + p) as discussed below.
A. Early and late time behavior of e and ∆n f Before we present the results of our QCD kinetic theory simulations, it is insightful to analyze the evolution of the energy-momentum tensor and conserved currents at early and late times. Due to the longitudinal expansion, the net-charge ∆n f densities of all flavors are diluted according to indicating that throughout the evolution the net-charge density per unit rapidity τ ∆n f (τ ) = (τ ∆n f ) 0 remains constant. Similarly, the energy density of a homogenous system undergoing a boost-invariant longitudinal expansion decreases according to where in addition to the trivial dilution, the second term on the right hand side characterizes the work performed against the longitudinal expansion [85,86], which is proportional to the longitudinal pressure p L = τ 2 T ηη given by in QCD kinetic theory. Since at early times the system is rapidly expanding in the longitudinal direction, it is unable to maintain a sizeable longitudinal pressure. Early on, one therefore has p L ≪ e , such that initially the energy per unit rapidity τ e(τ ) = (τ e) 0 remains approximately constant. Since initially τ ∆n f (τ ) = (τ ∆n f ) 0 and τ e(τ ) = (τ e) 0 are both constant, this further implies that, for finite density systems, the energy per baryon remains approximately constant at early times. Evidently, this is sharp contrast to the behavior at asymptotically late times, where for an equilibrated QCD plasma the longitudinal pressure becomes p L = e 3 such that τ 4 3 e(τ ) = (τ 4 3 e(τ )) ∞ approaches a constant and the energy per baryon decreases ∝ τ −1 3 . By considering the evolution of e(τ ) along with the ratios of ∆n f (τ ) e(τ ) 3 4 one then finds that the quantity τ 1 3 T eff (τ ) = τ 1 3 T eff ∞ as well as the ratios of the various chemical potentials to the temperature µ f,ldm (τ ) T ldm (τ ) = (µ f T ) ∞ become constant at asymptotically late times.

B. Pressure isotropization and kinetic equilibration
We now turn to the presentation of our QCD kinetic theory results and first analyze the evolution of the bulk anisotropy, characterized by the ratio of the longitudinal pressure p L to the energy density e shown in Fig. 29. Different curves in Fig. 29 show the results for p L e at zero and finite net-baryon density 6 as a function of the scaling variableω = (e + p)τ (4πη), which at zero netbaryon density (µ B T = 0) corresponds to the familiar expressionω = T τ (4πη s) employed in previous works.
Starting from early times, where the system is dominated by the rapid longitudinal expansion and highly anisotropic (p L ≪ e), the longitudinal pressure continuously rises as kinetic interactions become increasingly important. Despite the rapid increase of p L e at early times, the system remains significantly anisotropic throughout the entire evolution shown in Fig. 29 and only approaches an isotropic equilibrium state on much larger time scales.
Nevertheless, starting aroundω ≳ 1 the approach towards equilibrium is described by viscous hydrodynamics, where to leading order in the gradient expansion, the longitudinal pressure is given by Expressing the non-equilibrium correction in terms of the dimensionless ratio ηT eff (e+p) which at zero density 6 Note that instead of characterizing the amount of net baryon density in terms of the initial energy per baryon (eτ ) 0 (∆n B τ ) 0 , the curves at different densities are labeled in terms of the asymptotic ratio of (µ B T )eq extracted from our simulations, and we refer to Tab. II for the corresponding initial state parameters. reduces to the familiar η s, the pressure evolution in hydrodynamics is the determined by By analyzing the late time behavior of the different curves, we can then extract the transport coefficient ηT eff (e + p), who's values are indicated in Tab. II. We note that although in principle ηT eff (e + p) can exhibit a dependence on the chemical potentials µ f T , we find that in the relevant range of µ B T for our simulations, the extracted values only differ by about ten percent, which a posteriori justifies its treatment as a constant when defining the scaling variableω and extracting the values of ηT eff (e + p) based on the late time behavior in Eq. (53). When expressed in terms of the macroscopic scaling variablew, one also observes that the evolution of p L e is rather insensitive to the microscopic coupling strength λ = 5, 10 in Fig. 29, as discussed in detail in [50] for charge neutral plasmas.
By taking into account the (small) µ B T dependence of the transport coefficient ηT eff (e+p), the results for p L e in Fig. 29 are presented such that they all exhibit the same hydrodynamic behavior in Eq. (53) at late times ω ≳ 1 which is indicated by a black dashed line. However, by comparing the results for different net-baryon densities (µ B T ) eq , one clearly observes that at early times ω ≲ 1 the isotropization of the pressure proceeds more slowly for systems with a larger net baryon density. We will show shortly, that this feature can be understood by considering the fact that more baryon rich systems necessarily feature a larger abundance of quarks as compared to the initially gluon dominated zero density systems, which along with the less efficient isotropization of quark and anti-quark distributions leads to a slower build-up of the longitudinal pressure in the system.

C. Kinetic and chemical equilibration of light flavors
Beyond the evolution of the pressure anisotropy, which provides an estimate of the range of applicability of a hydrodynamic description of the QGP, it is also insightful to consider the evolution of the phase-space densities of gluons, quarks and anti-quarks to scrutinize the underlying microscopic dynamics. Our results are compactly summarized in Fig. 30, where we present the evolution of the various distributions for three different values of the net baryon density (µ B T ) eq = 0 (upper panel), (µ B T ) eq = 1.31 (middle panel) and (µ B T ) eq = 2.38 (lower panel). Different rows in each panel correspond to the distributions for different particle species, while different columns show the distributions at four different times, corresponding to the initial conditions in the first column, andω = 0.5, 1.0, 1.5 in the second, third and fourth column. We focus on the evolution of the phase space distributions of gluons (g), up-quarks (u), up-antiquarks (ū) and strange quarks (s), noting that the distributions of strange and anti-strange quarks are identical f s = fs, and that up and down quark distributions exhibit essentially the same features.
Starting from the behavior at zero net-baryon density (µ B T ) eq = 0 depicted in the top panel, where we assume that there are initially no quarks present in the system, one finds that quark/antiquarks of all flavors are democratically produced, and naturally inherit the anisotropy of the gluon distribution. However, the quark/antiquark distributions at intermediate stages of the evolutioñ ω = 0.5, 1.0 exhibit a larger degree of anisotropy as compared to the gluon distribution, indicating the slower isotropization of quarks/antiquarks. By considering the underlying microscopic processes in the bottom up scenario [34], one expects the isotropization of the gluon distribution to be driven by the radiative decay of hard gluons due to collinear g → gg and g → qq processes, followed by gg → gg, gq → gq and gq → gq elastic scatterings which isotropize the momentum distribution of soft gluons, whereas quarks/anti-quarks are pre-dominantly produced via collinear g → qq splittings and to a lesser extent by gg → qq elastic conversions, with the subsequent isotropization of soft quarks/anti-quarks due to qg → qgqg →qg,qq → qq,qq →qq and qq → qq elastic scattering processes. Based on the different color factors for the elastic scattering processes involving quarks and gluons, e.g. M gg gg 2 ∝ C 2 A and M gq gq 2 ∝ C F C A (see Tab. I), it is then natural to expect a faster isotropization of the gluon distribution.
When considering the evolution of the phase-space distributions at finite net-baryon density, shown in the central and bottom panels of Fig. 30 for (µ B T ) eq = 1.31 and 2.38, one finds that the overall behavior of the phasespace distributions at different times is rather similar to the zero density case. However, at finite density, the non-zero values of the conserved u and d charges lead to an overabundance of up and down quarks as compared to anti-quarks of the same flavor. Since at larger netbaryon density, u and d quarks carry a significant fraction of the initial energy, the larger degree of anisotropy of the quark distribution then manifests itself at the level of the bulk anisotropy p L e, seen in Fig. 29.
Besides the dynamics of the up and down flavors, it is also interesting to compare the evolution of the strange quark distribution (f s ) at zero and finite density. While at zero density strange quarks can be efficiently produced via inelastic g → qq and elastic conversions gg → qq conversion, the direct production of ss-pairs from u and d quarks is only possible through quark/anti-quark annihilation qq → qq, which at finite density is suppressed due to the lack of anti-quarks. By comparing the results for f s in the upper and lower panels of Fig. 30, one therefore finds that the strangeness production at finite density is delayed untilω ∼ 1, when strangeness is efficiently produced from inelastic g → qq and elastic conversions gg → qq conversions.
Next in order to further analyze the chemical composition of the QGP, we follow [6] and investigate the fraction of energy e a (τ ) e(τ ) carried by each individual species a during the non-equilibrium evolution. Our results for this quantity, e a (τ ) e total (τ ) are presented in Fig. 31 as a function of the scaling variableω. Different panels in Fig. 31 show the results for different net-baryon densities, with (µ B T ) eq = 0 in the top panel, (µ B T ) eq = 1.31, 1.19, 1.01 in the central panel and (µ B T ) eq = 2.38 in the bottom panel, while different solid, dashed and dotted curves in each panel correspond to the result obtained by varying the chemical composition of the initial state (see Tab. II). Starting with the evolution at zero net-baryon density, we find that for gluon dominated initial conditions (e g,0 e 0 = 1) a large part of the initial energy of gluons is rapidly transferred to quarks and anti-quarks of all flavors. Similarly for quark/anti-quark dominated initial conditions at zero density (e g,0 e 0 = 0), a rapid energy transfer from the quark to the gluon sector occurs, effectively resulting in a memory loss of the initial QGP chemistry on time scalesw ∼ 1. Eventually, forw ≳ 0.5 the zero density plasma becomes gluon dominated, before relaxing towards chemical equilibrium on time scalesω ∼ 1 − 2. Clearly, the situation is different at moderate or large net baryon density shown in the bottom panel of Fig. 31, where u and d quarks carry the dominant fraction of the energy density throughout the evolution. Due to the fact that multiple quark/anti-quark species contribute different amounts, one observes that the evolution of the chemistry of the QGP at moderate and large net-baryon density is significantly more complicated, and the approach towards equilibrium occurs on somewhat larger time scalesω ∼ 1.5 − 2.5, due to the less efficient production of anti-quarks (ū,d) and strangeness (s,s).
We conclude our discussion of equilibration in longitudinally expanding QCD plasmas, by considering once again the evolution of the characteristic scales m 2 D , m 2 Q,u and T * that govern the rates of elastic and inelastic interactions in the plasma. The time evolution of these quantities is presented in Fig. 32, where in order to account for the continuous expansion of the system we have normalized the respective quantities as τ T * ) eq such that forω ≫ 1 all ratios approach unity. By comparing the evolution of the different curves, we find that simulation results at different coupling strength λ = 5, 10 are in good overall agreement when expressing the evolution in terms of the scaling variablew. While the effective temperature T * relaxes towards its equilibrium value on time scales w ∼ 1, the screening masses m 2 D ,m 2 u for gluons and (up-) quarks only approach their equilibrium values at asymptotically late times, indicating residual deviations from local thermal equilibrium on the order of 10%.

VI. Conclusions
We developed a QCD kinetic description of the light flavor QCD degrees of freedom to study near and farfrom equilibrium dynamics of the Quark Gluon Plasma (QGP) at zero and finite density of the conserved baryon number, electric charge and strangeness. Based on numerical solution of the kinetic equations, including all leading order elastic and in-elastic interactions between gluons, quarks and anti-quarks, we exposed the general features of kinetic and chemical equilibration of nonequilibrium QCD plasmas in the perturbative regime at (asymptotically) high energies.
Generally, we find that, albeit the energy transfer between quark and gluon degrees of freedom can take a significant time, kinetic and chemical equilibrations of QCD plasmas occur roughly on the same time scale. By performing detailed investigations of the evolution of the spectra and collision rates, we further established a microscopic understanding of different equilibration processes in QCD plasmas, which generalizes earlier results obtained in pure glue QCD [4,7,46] and QCD at zero density [8,50,84].
Specifically, for over-occupied systems, which initially feature a large number of low energy gluons, we find that the thermalization process proceeds via a self-similar turbulent cascade, before eventually reaching equilibrium on a time scale ∼ 4πη s T eq . Conversely, for underoccupied systems, which initially feature a small number of high energy quarks or gluons, thermalization is achieved via the bottom-up scenario, with a number of interesting features regarding the role of quark and gluon degrees of freedom.
Studies of the equilibration of the QGP in a longitudinally expanding system provide the basis for a realistic matching of the initial state in heavy-ion collisions to initial conditions for the subsequent hydrodynamic evolution. By analyzing the macroscopic evolution of the energy momentum tensor and the microscopic evolution of the phase-space distributions of quarks and gluons, we found that viscous hydrodynamics typically becomes applicable on time scales where (e + p)τ (4πη) ∼ 1; however, isotropization and strangeness production proceed more slowly for finite density systems, and we refer to our companion paper [52] for further discussions of phenomenological consequences.
We finally note that the numerical framework to solve the QCD kinetic equations presented in this paper could be extended in several regards, e.g. by including heavy flavor degrees of freedom or electroweak interactions, to study a variety of aspects regarding the early time dynamics of high-energy heavy-ion collisions and the thermalization of the early universe.

Acknowledgment
We thank Giuliano Giacalone, Aleksi Kurkela, Aleksas Mazeliauskas, Jean-Francois Paquet, Ismail Soudi, Derek Teaney for discussions and collaboration on related projects. We are especially grateful to Aleksas Mazeliauskas for fruitful exchanges on the numerical implementation of QCD kinetic theory. This work is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -project number 315477589 -TRR 211. The authors also gratefully acknowledge computing time provided by the Paderborn Center for Parallel Computing (PC2) and National Energy Research Scientific Computing Center under US Department of Energy.
A. Weight Function Discretization

Weighted Integral: Discretization
We discretize the Boltzmann equation Eq. (6) with the weighted integral of a function F(⃗ p) (particle The weighted integral of distribution function f (⃗ p) reads where w with χ Ωx (x) ∶ Ω x → Z 2 the indicator function on domain Ω x . The completeness relation ensures the summation of weighted integral n(i p , j θ , k φ ) to be the total number of specific particle The weight functions satisfying completeness relation Eq. (A2) can be achieved by decomposing into two parts with left and right sided weights The spectral weight S (p) (p) pi+1 pi needs to be constructed in a form with y(p) an arbitrary function of p so that S (p) (p i ) pi+1 pi = 1, S (p) (p i+1 ) pi+1 pi = 0 and Eq. (A2) satisfied.

Sum Rules
Indeed, the above function y(p ∈ [p i , p i+1 ]) can be reversely expressed with the left and right weights Specifically, we work with a properly choice of the functions for p, cos(θ) and φ y p (p) = p, y θ (θ) = cos(θ), y φ (φ) = 1.
(A9) that provides the way to evaluate energy and longitudinal momentum of the particle in discretized form, following the definition of weighted integral Eq. (A1), completeness relation Eq. (A2) and sum rule Eq. (A8) B. Discretization of Collision Integrals

a. Discretization and Efficient Samplings
The elastic collision integral for particle "a" with momentum ⃗ p 1 in process a, b → c, d (p 1 , p 2 ↔ p 3 , p 4 ) reads: M ab cd (⃗ p 1 , ⃗ p 2 ⃗ p 3 , ⃗ p 4 ) 2 is the matrix element square for process "a, b ↔ c, d" summed over spin and color for all particles, and F ab cd (⃗ p 1 , ⃗ p 2 ⃗ p 3 , ⃗ p 4 ) describes the statistical factor for "a, b ↔ c, d".
As the energy density and longitudinal momentum flux can be directly evaluated from the discretized form in Eq.(A10), energy and longitudinal momentum conservation can be exactly fulfilled by the discretized form of collision integral, as a derivative of distributions.
We take the most complicated process q 1q1 ↔ q 2q2 as an example, other processes follow. According to Eq. (B7), the discretization forms read: (1) For quark q 1 , note that Q q1q1 q2q2 (12 34) = Q q1q1 q2q2 (12 43) We have following conservation laws automatically proved by the discretized collision integral from the completeness relation Eq. (A2) and sum rule Eq. (A8).

Longitudinal Expansion Integrals
In a longitudinally expanding system, there is an additional contribution to the collision integral which can be squares. Once either ω BH , x DQf evolves outside of the grid squares, we recalculate the inelastic rates based on the current scales.

Monte-Carlo Sampling
We perform Monte-Carlo integration of collision integrals for both elastic and inelastic processes.
For the elastic samplings, we first sample q, then −q ≤ ω ≤ q and finally q−ω 2 ≤ p 1 , q+ω 2 ≤ p 2 according to the discussions in Sec. B 1 b. The samplings of angles cos(θ q ), φ q together with φ 1 , φ 2 help us determine the values for p 3 , p 4 . With each of set of samplings for momenta p 1 , p 2 , p 3 , p 4 , we calculate the discretized collision integral according to Eq. (B7) which simultaneously updates the gain and loss terms of all the processes, which by virtue of the sum rules ensures exact energy and particle number conservation as was discussed in Sec. B 1 a. Similarly, the evaluation of inelastic collision integrals are performed by sampling p, z and the angle with respect to the longitudinal direction cos(θ) according to Eq. (B42) which also simultaneously updates the gain and loss terms.
The summation of all relevant processes and all samplings for collision integrals provides the total collision integral in the Boltzmann equation for specific particle.
The sampling numbers are chosen to be N sample,elastic =512 for each specific elastic process and N sample,inelastic =256 for each specific inelastic process.

Adaptive time step
Evolving the particle distributions in the discretized domain, we need an adaptive time step size ∆t to perform a stable increment for each distinct step ∆n a (i p , j θ , k φ , t) = − C 2↔2 a [n](i p , j θ , k φ , t) (C3) +C 1↔2 a [n](i p , j θ , k φ , t) + C z−exp a [n](i p , j θ , k φ , t) ∆t.
In order to do that, we need to make sure essential physics scales are not changing rapidly in each step. Common scales into considerations are total number density n, total energy density e and total longitudinal pressure p L , Debye screening mass square m 2 D quark screening mass square m 2 Qf , effective temperature T * . Some other scales may also be considered. However, more scales will not only increase the stability, but also slow down the evolution with a shorter resulting time step ∆t. According to their expressions listed in Sec. II A, their relative changing rate can be approximated by ∂ t e e = ∫ pd 3 p ν g C A ∂ t f g + ν q C F ∑ f (∂ t f q + ∂ t fq) ∫ pd 3 p ν g C A f g + ν q C F ∑ f (f q + fq) ∂ t p L p L = ∫ pcos 2 (θ)d 3 p ν g C A ∂ t f g + ν q C F ∑ f (∂ t f q + ∂ t fq) ∫ pcos 2 (θ)d 3 p ν g C A f g + ν q C F ∑ f (f q + fq) .