Quark production and thermalization of the quark-gluon plasma

We first assemble a full set of the Boltzmann Equation in Diffusion Approximation (BEDA) for studying thermalization/hydrodynamization as well as the production of massless quarks and antiquarks in out of equilibrium systems. In the BEDA, the time evolution of a generic system is characterized by the following space-time dependent quantities: the jet quenching parameter, the effective temperature, and two more for each quark flavor that describe the conversion between gluons and quarks/antiquarks via the $2\leftrightarrow2$ processes. Out of the latter two quantities, an effective net quark chemical potential is defined, which equals the net quark chemical potential after thermal equilibration. We then study thermalization and the production of three flavors of massless quarks and antiquarks in spatially homogeneous systems initially filled only with gluons. A parametric understanding of thermalization and quark production is obtained for either initially very dense or dilute systems, which are complemented by detailed numerical simulations for intermediate values of initial gluon occupancy $f_0$. For a wide range of $f_0$, the final equilibration time is determined to be about one order of magnitude longer than that in the corresponding pure gluon systems. Moreover, during the final stage of the thermalization process for $f_0\geq 10^{-4}$, gluons are found to thermalize earlier than quarks and antiquarks, undergoing the top-down thermalization.


I. INTRODUCTION
Intricate quantum many-body systems in QCD are observed to be produced in heavy-ion collisions at RHIC and the LHC.So far, a complete description of real-time dynamics of such systems is still beyond the reach of first-principles calculations in QCD.On the other hand, a promising, consistent understanding of this issue in the weak-coupling/high-energy limit is emerging (see refs.[1][2][3] for recent reviews).The validity of such a perturbative description is built upon the theoretical idea of parton saturation, which predicts that the partons initially released in high-energy nuclear collisions are mostly saturated gluons that carry a semi-hard, perturbative scale of the order of Q s , the saturation momentum [4][5][6].
Assuming that Q s lies in the perturbative regime, a unified method in QCD to describe the libration of partons from nuclear wave functions during a heavy-ion collision and the ensuing thermalization/hydrodynamization process is yet to be established.Alternatively, classical field simulations, which sum over tree-level diagrams up to all orders in the strong coupling, have been employed to study the earl-time dynamics in heavy-ion collisions [7,8].However, the classical field approximation is known to break down before the system thermalizes as the classical theory is, in essence, non-renormalizable [9] and the results under this approximation eventually become dependent of the ultra-violet (UV) cutoff [10,11].Based on the quasi-particle approximation, it has been argued that the classical field theory and the Boltzmann equation are equivalent when the gluon occupation number is higher than unity but lower than 1/α s with α s the strong coupling [12].Although such a transition cannot be confirmed by detailed calculations at the lowest order beyond classical diagrams [13,14], switching from classical field simulations to the Boltzmann equation has become a common practice in the studies of pre-equilibrium dynamics [1][2][3].
The Boltzmann equation remains the primary weak-coupling method in QCD for investigating thermalization/hydrodynamization.The leading-order (LO) Boltzmann equation incorporates collision kernels for both the 2 ↔ 2 processes and the 1 ↔ 2 processes.Two variations of the LO Boltzmann equation have been presented in the literature [15,16].In ref. [15], the diffusion approximation [17] has been utilized to the 2 ↔ 2 kernel.This variation is equivalent to that in ref. [16], referred to as Effective Kinetic Theory (EKT), under the leading logarithmic approximation.
The 2 ↔ 2 collision kernel in diffusion approximation [17], after being completed with quantum statistics, has been used to investigate the onset and formation of the Bose-Einstein condensate (BEC) of gluons in initially over-populated systems [18,19].Including only the 2 ↔ 2 kernel, the Boltzmann equation in diffusion approximation (BEDA) admits unique solutions with the gluon distribution composed of a regular distribution for momentum p > 0 and a BEC: a macroscopic accumulation of gluons as a result of non-vanishing gluon flux at p = 0 [19].Such solutions describe the evolution of initially over-populated systems: the excessive number of gluons compared to that can be accommodated by the final thermal distributions eventually form a BEC.Similar solutions were known to exist generically for the Boltzmann equation of bosons when only the 2 ↔ 2 processes are taken into account [20][21][22][23][24].
When the number-changing processes are considered, different theories could behave differently.In the spin-0 scalar theory studied in ref. [25], the BEC could form transiently before it vanishes during the thermalization process.Unexpected from the parametric estimates in ref. [26], this is, however, not the case in QCD.In pure gluon systems, the gluon flux at p = 0 is observed to always vanish in diffusion approximation and no gluons are deposited in a BEC when the 1 ↔ 2 processes are included [27,28].Moreover, more elaborated analyses of infrared modes beyond the scope of the kinetic theory have showed no evidence for a BEC either [29].
Inside a dense QCD medium, the most efficient mechanism for a high-energy parton to lose all its energy is known to be radiative energy loss driven by the Landau-Pomeranchuk-Migdal (LPM) effect [30][31][32].When the 1 → 2 processes due to the LPM effect dominate the evolution of a system, one intriguing phenomenon could arise: comparing to the typical time scale for the 1 → 2 splitting, the thermalization time for soft partons produced from hard ones in the system is so brief that these soft partons establish thermal equilibration among themselves.Then, the thermal bath of soft gluons is gradually heated up by quenching hard partons until the system achieves full thermalization.This phenomenon is referred to as the bottom-up thermalization, which was originally discovered in longitudinally boost-invariant systems of gluons in ref. [15] (and first confirmed by detailed numerical simulations in ref. [33]).In a spatially homogeneous, non-expanding system, the bottomup thermalization also emerges as the last stage of thermalization in initially very dilute cases [28,34].
Allowing the production of quarks and antiquarks delays the thermalization process [19,[35][36][37][38].In spatially homogeneous systems, the equilibration time for the case with three flavors of massless quarks and antiquarks was found to be typically about 5 to 6 times longer than that in the pure gluon case when only the 2 ↔ 2 processes are considered in the BEDA [19].The thermalization process is still delayed when the 1 ↔ 2 processes are additionally taken into account as studied using EKT [36,38].Moreover, gluons are found to achieve kinetic equilibrium among themselves before quarks and antiquarks approach thermal equilibrium [36].In longitudinally boost-invariant cases, the delay due to quark production also occurs, and the system is observed to hydrodynamize before it achieves chemical equilibration [35][36][37][38].
In this paper, we present, for the first time, a full set of BEDA with both the 2 ↔ 2 and 1 ↔ 2 kernels for gluons and massless quarks and antiquarks, and carry out a comprehensive study of thermalization and quark production in spatially homogeneous systems by solving the BEDA.In our studies, we carry out parametrical estimates and obtain numerical solutions of the full theory as well as some analytic results under certain approximations.With this, we obtain a more detailed understanding of both thermalization and quark production in spatially homogeneous systems initially occupied by gluons, as well as the fact that topdown thermalization of gluons always appear due to the slow production rate of quarks and antiquarks if the plasma is not extremely under-occupied initially.Below, we outline the structure of this paper.
In sec.II A, we incorporate the 1 ↔ 2 kernel using the LPM splitting rates [30][31][32]39] into the BEDA while only the 2 ↔ 2 kernel for gluons and N f flavors of massless quarks and antiquarks were included in ref. [19].In this way, we assemble a full set of BEDA as an extension of that for pure gluon systems in ref. [15].Since this paper focuses on investigating spatially homogeneous systems, we streamline the BEDA for this scenario in sec.II B. We also study the low momentum limit of the BEDA in order to examine whether one can consistently impose a boundary condition with a vanishing gluon flux at p = 0 onto the 2 ↔ 2 kernel.This extends previous discussions on the onset of a BEC in pure gluon systems [18,19,[26][27][28][29] to generic systems including additional N f flavors of massless quarks and antiquarks.
By solving the BEDA, the rest part of this paper focuses on the exploration of thermalization and quark production in spatially homogeneous systems.Especially, we investigate whether and to what extent quark production would modify each stage of thermalization in pure gluon systems as previously discussed by the authors in ref. [28].For this purpose, the initial gluon distribution is chosen to be the same as those in [18,19,28] while the distributions for N f = 3 flavors of quarks and antiquarks vanish initially, as detailed in sec.II C. In this case, the time evolution of the system is characterized solely by the jet quenching parameter q and the screening mass squared m 2 D or, equivalently, the effective temperature T * .This subsection also covers the corresponding thermal equilibrium states, and our definition of the equilibration time.
In sec.III A, we carry out parametric estimates for the systems initially composed only of gluons typically carrying a hard momentum Q with an occupation number f 0 ≪ 1.A similar analysis is carried out in sec.IV A for initially over-populated systems of gluons with f 0 ≫ 1.In neither case do we observe any qualitative difference in the thermalization process when compared with the pure gluon scenario [1,15,28,[33][34][35][36][37][38].Complementary quantitative studies are provided in sec.III B for initially under-populated systems with f 0 = 10 −2 , 10 −4 and 10 −6 and in sec.IV B for initially over-populated systems with f 0 = 1, 10 and 100.
Besides, the 1 → 2 splitting rates used in our studies are documented in Appendix A, and the numeric methods for our GPU and CPU simulations are outlined in Appendix B.
The code is available in [40].

II. THE QCD BOLTZMANN EQUATION IN DIFFUSION APPROXIMATION
In this section, we assemble a full set of the Boltzmann equation in diffusion approximation (BEDA) for gluons and N f flavors of massless quarks and antiquarks as an extension of that for pure gluon systems [15,28].

A. The BEDA for gluons, quarks and antiquarks
The leading-order QCD Boltzmann equation includes not only the 2 ↔ 2 processes but also the medium-induced 1 ↔ 2 splitting [15,16] where a = g, q i and qi respectively standing for gluons, quarks and antiquarks of the i th flavor.The distributions are so normalized that the entropy density, number density, and energy density for gluons, quarks and antiquarks are respectively given by s = s g + s q + s q, (2a) where the number of colors ), and the entropy densities of gluons, quarks and antiquarks are respectively defined as In the diffusion approximation the 2 ↔ 2 kernel takes the form [17][18][19]41] where ϵ a = 1 for bosons and ϵ a = −1 for fermions.Here, the jet quenching parameter [42] for species a is defined as qa = C a q, and q ≡ 4πα 2 s L where C a = {C A = N c , C F } respectively for g (adjoint) and q i /q i (fundamental), L ≡ ln and the typical transverse momentum broadening ⟨p 2 t ⟩ = 2 3 ⟨p 2 ⟩ for isotropic systems.The effective temperature is given by [18] T where the screening mass squared is defined as with m 2 D,a standing for the contribution from parton a.Here, m 2 D reduces exactly to the Debye mass squared at O(α s ) for a QGP in thermal equilibrium [43].It is given by two times the effective mass squared m 2 eff,g in ref. [16].
By extending the source terms in ref. [19], one has, for each quark flavor i, where one needs to define two more parameters for each quark flavor: In terms of these two quantities, we define an effective net quark chemical potential for the which reduces to the corresponding net quark potential in thermal equilibrium, denoted by µ i , as a result of the conservation of the net quark number.Besides, the 2 ↔ 2 kernel also conserves the total energy and the total parton number. 1he 1 ↔ 2 processes are described by where ν a is the number of spin times color degrees of freedom for parton a: ν a = 2N c for q or q and ν a = 2(N2 c − 1) for g, and with Here, x is the momentum fraction carried by particle b for the process a → bc, l ≡ xp and k ≡ (1 − x)p.Following ref. [15], the splitting rates in the deep LPM regime [32,39] (see also Appendix A) are used in this paper: where q and q are understood to be of the same flavor, and the quark flavor index has been dropped as the splitting rates are the same for each flavor of massless quarks and antiquarks.
Explicitly, one has It is straightforward to check that the 1 ↔ 2 kernel above conserves the total energy and the net quark number for each flavor.As a result, the BEDA with both the 2 ↔ 2 and 1 ↔ 2 kernels admits the following thermal fixed-point solutions: where T and µ i are the thermal equilibrium temperature and the thermal net quark chemical potential, respectively.In the case where the energy density and the net quark number density for flavor i are conserved, T and µ i can be obtained by the conservation of these two quantities.

B. The BEDA for spatially homogeneous systems
For spatially homogeneous systems, the BEDA reduces to where the overdot and the prime denote the derivatives with respect to t and p, respectively, and we have introduced the current J a as As the 2 ↔ 2 kernel contains second-order derivatives with respect to p, we need to impose two boundary conditions to solve the BEDA in eq. ( 20), as described in detail in Appendix B. One natural choice is to impose For gluons, the boundary condition lim p→0 p 2 J g = 0 may not be imposed consistently: If f g ∝ 1/p at low p, it requires lim p→0 pf g → T * [19].In order to check the consistency of imposing this boundary condition, let us first investigate the low-p behavior of f a below.
Following the pure gluon case [27,28], we expand the right-hand side of eq. ( 20) around p = 0. Here, one needs only to keep the first term on the right-hand side of eq. ( 12) for the 1 ↔ 2 kernel, evaluated with the splitting rates at leading order in x and F at leading order in p.To be more specific, by assuming p ≪ p/x one has where the last two terms on the right-hand side can be neglected for quarks and antiquarks while only the last term is subleading for gluons as f g ∝ 1/p.In this way, one has with µ i * defined in terms of I i c and Īi c in eq. ( 11), and π L).Accordingly, in the limit p → 0, the distribution functions take the form That is, they are the same as thermal distributions with the temperature and the net quark chemical potential respectively given by T * and µ i * .As a result, one can consistently set lim p→0 p 2 J g = 0, and there is no transient formation of a BEC according to [19].
Neglecting the time dependence of the macroscopic quantities in eq. ( 24), the BEDA admits the following approximate solutions C. The initial-state and final-state phase-space distributions In the rest part of this paper, we study thermalization and the production of three flavors of massless quarks (N f = 3) in a spatially homogeneous system initially filled only with gluons.In this case, all the quarks and antiquarks share identical distributions.For brevity, the distribution functions for different species f a are respectively denoted by with the quark flavor index neglected.
1.The distributions at t = 0 As we are interested to study how the whole evolution history in the pure gluon system, as studied in ref. [28], is affected by allowing quark production, the initial gluon distribution is chosen to be the same as that inspired by saturation physics in refs.[18,19,28]: Correspondingly, the initial number density and energy density are respectively given by Since the BEDA is invariant under F ↔ F , the system preserves F = F in the subsequent evolution given F = F at the initial time.Accordingly, one has I c = Īc and That is, in this case the time evolution of the system is governed by the same two quantities: the jet quenching parameter and m2 D as in pure gluon systems.For the initial distributions in eq. ( 29), the approximate solutions in eq. ( 26) can be expressed as As p F * ∼ p A * , we use for the parametric estimates in the following sections,

The thermal equilibrium states and the equilibration time
As µ i = 0, one can straightforwardly use energy conservation: to obtain the thermal equilibrium temperature Then, by using the thermal distributions in eq. ( 19), one has n eq = n eq,g + n eq,q + n eq,q ≡ 2(N 2 c − 1) + s eq = s eq,g + s eq,q + s eq,q ≡ 2(N 2 c − 1) + and And one can check that they satisfy ϵ eq = −P + T s eq + N f i=1 µ i (n eq,q i − n eq,q i ) with P = ϵ eq /3.
Note that for N f = 3 quarks and antiquarks together account for over half of the amount of ϵ, n and s while they only contribute one third of the overall values of qA and m 2 D in thermal equilibrium.
In order to compare quantitatively the thermalization times for N f = 0 and N f = 3, we define the equilibration time t eq in terms of a macroscopic quantity M as with W a constant.Following [19], we choose W = 0.05.

III. THERMALIZATION IN INITIALLY UNDER-POPULATED SYSTEMS
In initially under-populated systems, the total parton number at the initial time, by definition, is less than that in the corresponding thermal equilibrium state.For N f = N c = 3, an initially under-populated system has f 0 < f 0c = 0.308 for the initial distributions in eq. ( 29).
A. Parametric estimates for f 0 ≪ 1 In the pure gluon case, the system establishes thermal equilibrium through three distinct stages for f 0 ≪ 1 [1,28].Below, we investigate whether and how the three stages are modified when quark production is not ignored.
Stage 1. Overheating in the soft sector: During this stage, the properties of the system are dominated by hard gluons with p ∼ Q.

And one has
where the number density of hard gluons is denoted by n h,g and qA is simply denoted by q here and below.Quarks, antiquarks and soft gluons are predominantly generated by hard gluons through the 1 ↔ 2 processes while, as justified below, the 2 ↔ 2 processes play a parametrically negligible role in their production.
Like the pure gluon case [28], soft gluons are mostly produced via the g ↔ gg processes.
From the g → gg splitting, the number density of gluons typically carrying momentum p ≪ Q can be estimated as and, equivalently, The above estimate for f starts to break down when the reverse process: gg → g becomes equally important.According to eq. ( 24), this occurs when pf (p)/T * ∼ 1, corresponding Q.At p ≲ p * , the g → gg splitting and its reverse process balance out, leading to thermalization of soft gluons with an overheated temperature given by T * .This is consistent with the approximate solution in eq. ( 32).Accordingly, the system witnesses a pronounced accumulation of soft gluons typically carrying momentum p s ∼ p * with their number density given by The 2 ↔ 2 processes can be ignored in the above parametric analysis: From multiple elastic scattering, the gluons typically pick up some momentum broadening ∆p ∼ (qt) Quarks and antiquarks can be produced via either the g → q/q conversion in the 2 ↔ 2 scattering or the g → q q process.According to the source terms in eq. ( 9), the g → q/q conversion yields The quarks and antiquarks produced in this process typically carry hard momentum Q as hard gluons dominate m 2 D .From the g → q q splitting, the number density of produced q or q carrying a typical momentum p can be estimated as and, equivalently, The above estimate using single branching starts to break down at p ∼ p * when F (p * ) ∼ 1/(e − µ i * T * + 1) = 1/2.For p ≲ p * , the reverse process q q → g becomes non-negligible and its rivalry with the g → q q process leads to thermalization of this sector: F → 1/2, enforcing the Pauli principle.By taking into account all the typical values of p, one can conclude that the quarks and antiquarks produced by the g → q q splitting are mostly hard with a number density given by That is, in this process a gluon tends to radiate a quark/antiquark as hard as possible while, in contrast, the softest gluon with p ∼ p * is the most likely to be emitted via g → gg.This difference is due to the soft divergence of the latter process, which is absent in the former.
As f 0 ≪ 1, the g → q q splitting gives the dominant contribution, leading to a linear growth of the quark/antiquark number over time: We now verify the dominance of hard gluons during this stage.As the quark number and the soft gluon number are both parametrically smaller than that of hard gluons, q receives dominant contributions from hard gluons only.For m 2 D , soft gluons and soft quarks with typical momentum p s ∼ p * make the following contributions: They are both parametrically smaller than the contributions from hard gluons with Cooling and overcooling in the soft sector: During this stage, most of the partons in the system are still hard gluons, leaving both On the other hand, m 2 D receives dominant contributions from soft gluons whose momenta are pushed up to p s ∼ (qt) 1 2 by multiple elastic scattering.The number density of soft gluons can be estimated from the g → gg splitting: And their contribution to m 2 D takes the form which is parametrically larger than that from hard gluons starting from Qt ∼ α −2 s f 0 .Accordingly, one has The qualitative feature of the gluon distribution at early times.For parametric estimates, the qualitative features of the system in Stages 1 and 2 can be derived from these three distinct portions of f .Similar to Stage 1, f can be roughly divided into three distinct portions illustrated in Fig. 1, which are parametrically the same as pure gluon case [28].In the range of p dominated by the g → gg splitting, one has It holds down to p ∼ p s ∼ (qt) , where p s f (p s )/T * ∼ 1 and the reverse process gg → g becomes equally important.Consequently, detailed balance is established in the soft sector with p ≲ p s .At higher p, the estimate in eq. ( 53) breaks down at p ∼ p with s f where the initial distribution f 0 starts to dominate.
The qualitative features of the quark and antiquark distributions are predominantly determined by the gluon distribution.For Q ≳ p ≳ pF ≡ α Q, quarks and antiquarks are mostly produced by hard gluons via the g → q q splitting with Qt.
For p ≲ p, the intermediate portion of f ∝ p − 7 2 contributes to F as well.However, its contribution is parametrically negligible for p ≳ pF .This can be justified as follows.Since it is the most likely for a gluon to radiate a quark/antiquark carrying about half of its momentum, this portion of f gives a contribution By comparing this estimate with eq. ( 55), one can see that it becomes dominant only for p ≲ pF (note that one always has p s ≲ pF ≪ p throughout this stage).And the above estimate based on single branching breaks down at p ∼ p s as F (p s ) ∼ 1/2.At p ≲ p s , the system maintains a balance between the g → q q process and its reverse process, corresponding to F → 1/2.Accordingly, unlike Stage 1, the system witnesses a pronounced accumulation of soft quarks and antiquarks with p ∼ p s besides hard ones with p ∼ Q: Moreover, starting at 0 Q, soft quarks and antiquarks with p ∼ p s become parametrically more abundant than hard ones.
Let us then look into quark production via the 2 ↔ 2 process.From the conversion term in (9), one can obtain the quark distribution which is negligible compared to the contribution from the g → q q process for p ≳ p s in eq. ( 56).On the other hand, the quark number produced via this process is given by which is parametrically the same as the contribution from the g → q q process.Since the quark number density is parametrically smaller than n s,g for p s ≲ f 0 , the quark and antiquark production does not modify the parametric form of q or m 2 D .In summary, even though soft gluons contribute dominantly to m 2 D right after Qt ∼ α −2 s f 0 , they do not contribute significantly to quark production until Qt ∼ α −2 s .Afterwards, the soft sector of the system becomes overcooled with T * ≲ T , and soft quarks and antiquarks are dominantly produced by soft gluons via both the g → q q and gg → q q processes although there are less soft gluons than hard ones in the system.As a result, one has Stage 3. Reheating and mini-jet quenching: α −2 s f which is consistent with a thermalized QGP of temperature T * ∼ f 0 Q.This is qualitatively different from the previous stages as both q and m 2 D are now determined by soft partons.This can be further justified by the fact that the relaxation time and In the above parametric estimate, one can safely neglect hard quarks and antiquarks as their number density n h,q ∼ α s And the soft sector remains a thermalized QGP with a timedependent temperature T * as the relaxation time when all the relevant quantities have the same parametric forms as those in thermal equilibrium.Therefore, in Stage 3, the system undergoes the bottom-up thermalization, parametrically the same as the final stage in the pure gluon case with [15] or without [28,34] longitudinal expansion.The only distinction is that the soft thermal bath is now a QGP, composed of gluons as well as quarks and antiquarks.

B. Quantitative studies
In this subsection, we carry out some detailed simulations for initially under-populated systems, focusing on features that are not revealed in the above parametric estimates.Fol- D , q and n g with N f = 0 (dotted) and N f = 3 (solid).The right plot shows the result of n q in comparison with the scaling laws obtained by parametric estimates.
lowing [19], L is kept as a constant in our simulations in order to be consistent with the harmonic oscillator approximation used to derive the splitting rates in the deep LPM regime.
And we show the results for L = 1 and α s = 0.1 below. 2   1. Initially dilute systems with f 0 = 0.01 As a direct comparison with the pure gluon case studied in ref. [28], we first investigate the time evolution of the system for f 0 = 0.01.In this case, the four time scales take the following values: Although there is no significant separation between the time scales after Qt ∼ α −2 s , the three stages: overheating, cooling/overcooling, and reheating are still manifest in the time evolution of T * for N f = 0 [28] (see also the left panel of fig.2).Let us examine whether this is still the case for N f = 3.Even though the system is under-populated in terms of the total parton number for N f = 3, there is an initial excess of gluons comparing to that in the corresponding thermal equilibrium state.Therefore, n g has to decrease to its thermal equilibrium value eventually. 2 This choice for α s is arbitrary as long as it remains within the perturbative regime and satisfies α s f 0 ≲ 1.
A change of the value for the coupling constant translates into a re-scaling of the time variable, since both collision kernels are proportional to α 2 s .
However, as shown in the left panel of fig.2, it is observed to increase instead at early times.
Moreover, n g as well as qA , m 2 D and T * are all almost indistinguishable from those in the pure gluon case before Qt ∼ 10.That is, the early-time dynamics of the system is dominated by gluons.This is a consequence of the fact that quarks and antiquarks can only be produced slowly: n q grows linearly over t up to Qt ∼ 1 and its growth accelerates only afterwards, as shown in the right panel of fig. 2. Consequently, the soft sector of the system experiences almost identical overheating and cooling phases before its evolution history starts to deviate from that in the pure gluon system.
For N f = 0, T * undershoots the thermal temperature at the percent level and starts to increase around Qt ∼ 73, similar to that during Stage 3 in the limit f 0 ≪ 1 [28], while the late-time evolution for N f = 3 becomes qualitatively different.At Qt ∼ 73, one has T * ≈ 0.130Q, which is still above its thermal equilibrium value (T = 0.107Q).In the meanwhile, q reduces about 9% compared to that for N f = 0 while m 2 D is less sensitive to N f , only dropping about 3%.Such a trend continues later on, causing T * ∝ q/(α s m 2 D ) to decrease until the system fully thermalizes.Such a qualitative difference from the dilute limit is also manifest in the fact that the growth of n q at later times, although being accelerated, is slower than the estimated scaling behavior of (Qt) Let us delve into the details of thermalization and quark production for N f = 3, as further revealed by the quantities in fig. 3.At Qt ∼ 73, the contributions from quarks and antiquarks to q (the upper left panel) and m 2 D (the upper center panel) are both about one order of magnitude smaller than those from gluons.This means that gluons still dominate the subsequent evolution.Consequently, the system follows the trend in the pure gluon case, generating more gluons despite of the existing excess of gluons, as shown in the lower left panel.This further stokes the growth of q and m 2 D and, hence, increases the likelihood of gluons splitting or converting to quarks and antiquarks in order to curb the unwanted growth of n g .Accordingly, the system witnesses rapid production of quarks and antiquarks via both the 1 ↔ 2 processes (the lower right panel) and the 2 ↔ 2 processes (the lower center panel) around this time.
Around Qt ∼ 87, the gluon number density starts to decrease.Shortly afterwards, the contributions from gluons to q and m 2 D also start to decrease while their contribution to s, as depicted in the upper right panel of fig.3, has already been declining.Consequently, both q and m 2 D develop a peak over time as gluons always give a larger contribution than quarks Nf = 3 : g + q + q Nf = 3 : g Nf = 3 : g + q Nf = 3 : q + q FIG. 3. Contributions of gluons (g), quarks (q) and antiquarks (q) to q, m 2 D , s and n for f 0 = 0.01.Besides, for N f = 3, the time derivatives of the number densities for g, q + q and g + q + q via the 2 ↔ 2 processes and the 1 ↔ 2 processes are presented in the lower center and lower right panels, respectively.The lower right panel also shows ṅg + ṅq = ṅg + ṅq via the 1 ↔ 2 processes.and antiquarks.In contrast, both n and s continuously increases over time toward their thermal equilibrium values to which quarks and antiquarks contribute more significantly than gluons, as shown in eq. ( 37).
The production of quarks and antiquarks is only partially responsible for the rapid decrease in ṅg via the 1 ↔ 2 processes shown in the lower right panel of fig. 3. Another contributing factor to its decrease is the increasing significance of the gg → g process, which counteracts the g → gg process primarily responsible for the earlier increase in n g .Starting around Qt = 200, one has ṅ1↔2 g + ṅ1↔2 q = ṅ1↔2 g + ṅ1↔2 q ≈ 0. This implies that the excess of gluons is dominantly eliminated by the g → q/q conversion in the 2 ↔ 2 kernel and the g → q q splitting.As for the g → gg process, it has nearly established detailed balance with its reverse process (note that all the quantities for N f = 0 shown in fig. 3  reached their thermal equilibrium values).This is consistent with the observation in the left panel of fig.4: the gluon distribution for Qt ≳ 200 can all be well fitted, up to some maximum momentum larger than T * , by the Bose-Einstein distribution with a time-dependent temperature given by T * .In contrast, the quark distribution, presented in the right panel of fig. 4, does not fit the Fermi-Dirac distribution (except that F → 1/2 at low p) until the very end of the thermalization process.That is, gluons undergo the top-down thermalization, meaning that they thermalize first with a decreasing temperature T * before the system fully thermalizes.This observation is qualitatively consistent with the analysis in ref. [35].
Let us now discuss thermalization more quantitatively in terms of the equilibration time defined in eq.(39).For N f = 0, it takes t N f =0 eq ≈ 112Q −1 for the relative differences of T * , n and s from their thermal equilibrium values to all become less than 5% (see the upper left panel of fig.5).As shown in the upper right panel of fig. 5, the gluon distribution matches a Bose-Einstein distribution with a relative error smaller than 5% for all the values of p < 5T * around this time.And it fits perfectly the thermal distribution at still later times though, as exemplified by the result at Qt = 10 3 .
In terms of T * , n and s, one has t eq ≈ 424Q −1 for N f = 3.However, the lower left panel of n g /n eq,g g / eq,g s g /s eq,g n g /n eq,g (T * ) g / eq,g (T * ) s g /s eq,g (T * ) n q /n eq,q q / eq,q s q /s eq,q fig. 5 shows that the energy density, entropy density and number density of gluons and quarks are still quite different from their thermal equilibrium values.The seemingly equilibration of T * , n and s is a result of the contributions from gluons being much higher than their thermal equilibrium values, effectively compensating for the much lower contributions from quarks and antiquarks.On the other hand, at this time, n g , s g and ϵ g are only 3.5%, 5.5% and 7.5% respectively above their thermal equilibrium values with a time-dependent temperature T * .
And it takes about t eq,g = 509Q −1 for all these relative differences to drop below 5%, marking quantitatively the equilibration of gluons according the criteria in eq.(39).Eventually, all the quantities shown in the lower left panel of fig. 5 approach their thermal equilibrium values with a relative difference of less than 5% at a later time t eq = 810Q −1 , which is 7.2 times that for N f = 0.Such a separation between the equilibration times justifies the top-down thermalization of gluons, which is consistent with the observation in EKT that chemical equilibration takes a longer time to achieve [35][36][37][38].
The lower right panel of fig. 5 focuses on the high-p behavior of the gluon distribution at Qt ≥ 200.At Qt = 200, the gluon distribution fits, with a relative deviation being equal to or less than 5%, the Bose-Einstein distribution with temperature T * for all values of p ≤ 3.3T * .At higher momentum, it evidently deviates more from the thermal distribution.
The maximal momentum where f fits the thermal distribution with temperature T * increases over time, as illustrated by the results at Qt = 350 and 10 3 in this plot.In this aspect, the top-down thermalization in the weak-coupling limit of QCD is qualitatively different from that in the strongly-coupled CFT using the AdS/CFT correspondence [47,48], where the UV modes equilibrate first [49][50][51][52][53][54][55][56][57].

Initially very dilute systems
The parametric estimates in sec.III A are valid only if f 0 ≪ 1 such that all the four time scales are well separated.As discussed above, for f 0 = 0.01, gluons, after being coupled to N f = 3 flavors of massless quarks and antiquarks, thermalize in a top-down manner at late times.This is qualitatively different from Stage 3 in the parametric estimates.In this subsection, we investigate quantitatively whether late-time reheating could occur in initially more dilute systems.

Let us start with
0 ) = (0.01, 100, 2154, 3161).The detailed time-evolution of T * , m 2 D , q and n g for N f = 0 and N f = 3 (the left panel) as well as n q for N f = 3 (the right panel) are presented in fig.6.Similar to the case with f 0 = 0.01, these quantities are insensitive to N f at early times, showing that n g /n eq,g g / eq,g s g /s eq,g n g /n eq,g (T * ) g / eq,g (T * ) s g /s eq,g (T * ) n q /n eq,q q / eq,q s q /s eq,q At late times, the top-down thermalization persists for f 0 = 10 −4 .According to their energy density, number density and entropy density shown in the right panel of fig. 7, the gluons first establish thermal equilibration among themselves with a time-dependent temperature T * around t eq,g = 2446Q −1 .Then, the energy density, number density and entropy density of gluons and quarks all converge towards their thermal equilibrium values and start to exhibit a relative deviation of 5% or less around t eq = 3480Q −1 .This qualitatively agrees with the case for f 0 = 0.01 except that the final equilibration time is only about 42% longer than t eq,g .As another indicator of the shortened period of the top-down thermalization, we determine the equilibration time for N f = 0 in terms of T * , n and s and obtain t N f =0 eq = 1489Q −1 .This tells us that the final equilibration time for N f = 3 is only about 2.3 times the value of t N f =0 eq .That is, thermalization is less delayed due to quark production compared to the case with f 0 = 0.01.
The above complex thermalization pattern is also manifest in the evolution of the distribution functions, as presented in fig.8.The left panel shows the gluon distribution at different times.Except for the every early times when the typical momentum of soft gluons p s < p min with p min the infrared momentum cutoff, our numerical results confirm its low-p behavior with pf → T * .Moreover, we find that around Qt = 1700 < t eq,g , the gluon distribution can already be well fitted by the Bose-Einstein distribution with temperature T * for all values of p ≤ 5T * (with a relative error of 5% or less).On the other hand, the quark distribution F , as shown in the right panel, keeps increasing over time throughout the shown range of p.It approaches the thermal distribution from below only at the very late stage of the evolution.
Aiming to illustrate how late-time reheating could emerge for N f = 3, let us study the system with f 0 = 10 −6 , which corresponds to (α 0 ) = (10 −4 , 10 2 , 10 4 , 17783).The late-time behavior of T * is shown in the left panel of fig. 9. 3 It shows that the system overcools irrespective of the values of N f with T * reaching its minimum at Qt ≈ 809 for N f = 0 and at Qt ≈ 1087 for N f = 3.For N f = 0, the system, as expected, reheats towards full thermalization afterwrads.For N f = 3, reheating finally occurs.However, as a vestigial top-down thermalization, T * later on overshoots the thermal n g /n eq,g g / eq,g s g /s eq,g n g /n eq,g (T * ) g / eq,g (T * ) s g /s eq,g (T * ) n q /n eq,q q / eq,q s q /s eq,q N f = 3 equilibrium temperature slightly at the percent level before the system fully thermalizes.
The delay of thermalization due to quark production is much milder for f 0 = 10 −6 , as shown in the right panel of fig. 9.In terms of T * , n and s, one has t for the pure gluon system.For N f = 3, it is no longer evident that gluons thermalize first and undergo the top-down thermalization.And the final equilibration time in terms of the energy density, number density and entropy density of gluons and quarks is determined to be t eq = 1.79 × 10 4 Q −1 , which is only about 50% longer than t N f =0 eq .

IV. THERMALIZATION IN INITIALLY OVER-POPULATED SYSTEMS
In this section, we study initially over-populated systems in which there are a larger number of partons at t = 0 than that can be accommodated by the final thermal equilibrium state.Below we carry out both parametric estimates and numerical simulations for such systems.
A. Parametric estimates for f 0 ≫ 1 In the limit f 0 ≫ 1, one initially has Like the pure gluon case [1,28,34], the system goes through the following two stages towards achieving thermal equilibrium.
Stage 1. Soft gluon radiation and overheating: During this stage, q and m 2 D are dominated by hard gluons, and they are parametrically the same as those at the initial time.Accordingly, the soft sector of the system is the same as that in an overheated system with 0 Q.This can be justified by the following parametric estimates.
The soft sector of the gluon distribution is rapidly populated through radiation off hard gluons.For the range of p dominated by the g → gg splitting, the number density of gluons typically carrying momentum p is given by and, accordingly, the gluon distribution takes the following parametric form The above estimate is valid for Q with the typical soft momentum p s given by the balancing condition: For p ≲ p s ∼ p * , f approaches the thermal distribution with temperature T * .That is, f can be expressed in the same parametric form as that illustrated in fig. 1 (with p s and p given above).Accordingly, besides hard gluons, there is a pronounced accumulation of soft gluons carrying a typical momentum of order p s : Evidently, it is parametrically smaller than the hard gluon number density at Qt ≲ (α s f 0 ) −2 .
Quarks and antiquarks are produced much more slowly than soft gluons.Via the g → q/q conversion in the 2 ↔ 2 processes, one has It is parametrically the same as that produced via g → q q: That is, the 1 ↔ 2 processes are not more efficient than the 2 ↔ 2 processes in producing quarks and antiquarks.On the other hand, for p ≳ p s , the quark distribution is mostly determined by the g → q q splitting: with F (p s ) ∼ 1/2.And at lower p, quarks and antiquarks fill a thermal distribution with During this stage, the contributions from soft gluons, quarks and antiquarks to m 2 D and q are all parametrically smaller than those from hard gluons.At the end of this stage at with qa denoting the contribution from parton a.That is, the properties of the system are still mostly determined by gluons.
Stage 2. Momentum broadening and cooling: As justified below, the system is still dominated by gluons during this stage.In this case, the scaling laws for all the relevant quantities can be derived in the same way as the pure gluon case.That is, solving consistently the following two equations respectively given by momentum broadening due to multiple elastic scattering and energy conservation [28] yields , p ∼ (α s f 0 ) , q ∼ (α s f 0 ) , m 2 D ∼ (α s f 0 ) , T * ∼ (α s f 0 ) Irrespective of the detailed mechanism how the number of excessive gluons are eliminated, the above derivation is valid as long as it does not provide an alternative, faster way than multiple elastic scattering to increase the typical momentum of gluons [26,28,29].
Allowing the production of quarks and antiquarks does not modify the scaling laws derived above.As in Stage 1, the 2 ↔ 2 and 1 ↔ 2 processes are equally efficient to produce quarks and antiquarks, yielding and That is, parametrically, the majority of partons in the system consistently fill thermal dis- At Qt ≲ α , the typical momentum of partons p is always parametrically lower than T * .As a result, q and m 2 D do not receive significant contributions from quarks and antiquarks, and the system is hence dominated by gluons.
At Qt ∼ α n g /n eq,g g / eq,g s g /s eq,g n g /n eq,g (T * ) g / eq,g (T * ) s g /s eq,g (T * ) n q /n eq,q q / eq,q s q /s eq,q N f = 3 panel confirms the early-time linear growth of the quark number density.These observations align with the predicted behaviors in Stage 1 according to the parametric estimates.This is expected since the early-time evolution is universally described by the approximate solutions in eq. ( 26) for all values of f 0 .
The subsequent evolution of q, m 2 D , T * and n g appears qualitatively similar for both N f = 0 and N f = 3: they all steadily decrease towards their thermal equilibrium values.
Alternatively, the evolution history of T * can be read out from the gluon distribution at low momentum with pf → T * , as shown in the left panel of fig.11 for N f = 3.And f approaches the Bose-Einstein distribution with temperature T from above while the quark distribution, shown in the right panel of fig.11, steadily approaches the Fermi-Dirac distribution from below.This is all seemingly consistent with a cooling process in Stage 2. However, q, m 2 D , T * and n g do not fit the scaling laws derived in the limit f 0 ≫ 1.Similarly, the growth of n q is shown to slow down, but it does not exhibit the scaling behavior of (Qt) 3 7 .
The system for N f = 3 also undergoes the top-down thermalization before it achieves fully thermalization.In terms of T * , n and s, one has t eq = 106Q −1 in comparison with D , q and n g for N f = 0 (dotted) and N f = 3 (solid).The right panel shows n q over time in comparison with the scaling in the limit f 0 ≫ 1.

t
N f =0 eq = 23.2Q−1 for N f = 0.However, the left panel of fig. 12 shows that the energy density, number density and entropy density of gluons and quarks are still quite different from their final thermal equilibrium values around this time.That is, the system has not achieved chemical equilibration yet.In contrast, the relative deviations of the energy density, number density and entropy density of gluons from their thermal equilibrium values at temperature T * are much smaller.Moreover, the right panel of fig. 12 shows that the gluon distribution can already fit, with a relative error of 5% or less, the Bose-Einstein distribution with temperature T * up to p = 3.4T * as early as Qt = 100.Then, it takes only a little bit longer for their relative differences to all drop below 5% around t eq,g = 138Q −1 .At the end, in terms of the energy density, number density and entropy density of gluons and quarks, the final equilibration time is given by t eq = 229Q −1 , which is about 9.9 times the equilibration time for N f = 0.

Initially very dense systems
The time evolution of q, m 2 D , T * , n g and n q in the system with f 0 = 10 is shown in the two panels of fig.13.Similar to the cases studied previously, quarks and antiquarks play a n g /n eq,g g / eq,g s g /s eq,g n g /n eq,g (T * ) g / eq,g (T * ) s g /s eq,g (T * ) n q /n eq,q q / eq,q s q /s eq,q negligible role at early times.In terms of q, T * and n g shown in the left panel of fig.13, the system, regardless of the value of N f , universally experiences a rapid cooling stage following a slowly varying overheating one.At Qt ≈ 1.4, q, T * and n g for N f = 3 all start to drop more than 5% below their corresponding values in the pure gluon case while it takes about 10 times as much time for m 2 D to reach a similar level of deviation.In the meanwhile, the production rate of n q , as shown in the right panel of Fig. 13, slows down after a linear growth over time.However, these quantities except m 2 D do not demonstrate convincingly the scaling behaviors in the cooling stage for f 0 ≫ 1.
In the final stage of thermalization, gluons equilibrate first with a time-dependent temperature T * for N f = 3, as qualified by the quantities in the left panel of fig.14.This is qualitatively the same as the cases for f 0 ≥ 10 −4 .Quantitatively, the equilibration time for gluons is determined to be t eq,g = 76.9Q−1 by comparing the energy density, number density and entropy density of gluons to their thermal equilibrium values at temperature T * .
Even before this time, the gluon distribution can be very well fitted by the Bose-Einstein distribution with T * up to some maximum momentum higher than T * , as exemplified in the And solving exactly the BEDA only fills in quantitative details that are not parametrically more important than the above parametric form.Consistent with eq. ( 77), keeping only the dominant terms in the collision kernels under the assumption f ≫ 1 yields the following self-similar solutions with p ≡ (Qt) − 1 7 p/Q.Here, f s is the same as the pure gluon case, which can be obtained based on the same argument as that in refs.[1,26,29].To justify the self-similar form of F , let us assume the ansatz F = (Qt) α F s ((Qt) β p/Q) and take the 1 ↔ 2 kernel as a start.n g /n eq,g g / eq,g s g /s eq,g n g /n eq,g (T * ) g / eq,g (T * ) s g /s eq,g (T * ) n q /n eq,q q / eq,q s q /s eq,q N f = 3 f 0 = 100 F 0 = F0 = 0 Keeping only the dominant term from the g → q q splitting, one has C q 1↔2 = (Qt) −1 Cq (p) with Cq only dependent of p.Then, matching it to the time derivative of F yields α = 0 and β = −1/7.One can easily see that including the 2 ↔ 2 kernel does not change these two exponents.As illustrated in fig.16, our numerical solutions are indeed consistent with the self-similar solutions for some range of t where all the quantities shown in fig.(15) approach the expected scaling behaviors.
After exiting the scaling phase, the pure gluon system equilibrates around t for N f = 0.For N f = 3, gluons undergo the top-down thermalization before achieving full thermalization.This is qualitatively similar to all the cases studied previously for f 0 ≥ 10 −4 .As shown in the left panel of fig.17, gluons first equilibrate around t eq,g = 4479Q −1 .
Even before this time, the gluon distribution can be very well fitted by the Bose-Einstein distribution with T * up to some maximum momentum higher than T * , as exemplified in the right panel of fig.17.Eventually, all the partons equilibrate around t eq = 6491Q −1 , which is about 10.2 times that in the pure gluon case.

V. CONCLUSIONS
We have introduced the Boltzmann Equation in Diffusion Approximation (BEDA) to study the thermalization of a system of massless quarks and gluons with both 2 ↔ 2 and 1 ↔ 2 kernels.First, we have incorporated the 1 ↔ 2 kernel using the LPM splitting rates [30][31][32]39] into the BEDA while only the 2 ↔ 2 kernel for gluons and N f flavors of massless quarks and antiquarks were included in ref. [19].In this way, we have assembled a full set of BEDA as an extension of that for pure gluon systems in ref. [15].Its collision kernels contain the following space-time dependent quantities, calculated from the phase-space distribution functions of partons: the jet quenching parameter q, the effective temperature ) with m 2 D the screening mass squared, and two more for each quark flavor i, denoted by I i c and Īi c in eq. ( 10).In terms of the latter two quantities, a definition of the effective net quark chemical potential µ i * has been given for flavor i in eq. ( 11), which equals the net quark chemical potential after thermal equilibration.
Then, we have studied spatially homogeneous systems of gluons by solving the BEDA.To do so, we have first studied the low momentum limit of the BEDA exploring the possibility of the formation of Bose-Einstein Condensate.Similar to the observations in refs.[27,28], we find that both the gluon and quark distributions in the limit p → 0 are identical to those in thermal equilibrium with the temperature and the net quark chemical potential given by T * and µ i * , respectively.As a result, the gluon flux at p = 0, which is proportional to (lim p→0 pf g − T * ) with f g the gluon distribution [19], always vanishes and, hence, no gluons are deposited into a BEC.
After establishing the framework, we have applied it to investigate the thermalization process and quark production.First, we have focused on the initially underpopulated scenario with the initial occupation number f 0 ≪ 1. Exactly like the pure gluon case [1,28], the system thermalizes through three stages for N f = 3.And we find that quark production does not modify the parametric forms of q and m 2 D during each of these stages.In Stage 1, both q and m 2 D are dominated by hard gluons.And the soft sector is filled with gluons, quarks and antiquarks according to the thermal distributions with an overheated temperature T * ∼ Q and µ i * = 0 up to p ∼ p s .Here, p s is the typical soft momentum much smaller than Q.During this stage, most of quarks and antiquarks are hard and their number density grows linearly over time.And the contributions from soft gluons, quarks and antiquarks to q and m 2 D , albeit increasing over time, are all parametrically smaller than those from hard gluons.In Stage 2, m 2 D receives dominant contributions from soft gluons and keeps increasing over time while q is still dominated by hard gluons, leading to a declining T * .Around the moment when T * cools down to the thermal equilibrium temperature T ∼ f 1/4 0 Q, soft quarks and antiquarks become more abundant than hard ones.Later on, T * keeps cooling while the number density of quarks and antiquarks ∼ (qt) 3/2 grows more rapidly than that of soft gluons ∼ T * qt.This stage ends when there are parametrically equal numbers of soft gluons, hard gluons, quarks and antiquarks.Moreover, the soft partons have enough time to thermalize into a quark-gluon plasma (QGP) with an overcooled temperature T * ∼ f 1/3 0 Q.In Stage 3, both q and m 2 D are predominantly determined by soft partons in the QGP with its temperature T * increasing over time.The system achieves final equilibration after hard partons are fully quenched, causing T * to be reheated up to T .That is, the system in this stage undergoes the bottom-up thermalization, qualitatively similar to the cases studied in refs.[15,28,[33][34][35][36][37][38].
We have also carried out parametric estimates for initially over-populated systems with f 0 ≫ 1.We have found that the two stages preceding thermalization in the pure gluon case [1,28] are not affected by the production of quarks and antiquarks.And their contributions to q and m 2 D are always parametrically smaller than those from gluons before the system achieves thermal equilibration.To sum up, the system during Stage 1 shares some similarities with that in the limit f 0 ≪ 1: q and m 2 D are predominantly determined by hard gluons.Accordingly, the soft sector of the system looks the same as that in thermal equilibrium with an overheated temperature T * ∼ f 0 Q.And the number density of quarks and antiquarks grows linearly over time.This stage ends when the typical momentum of soft partons is pushed up to ∼ Q by multiple elastic scatterings.In Stage 2, the typical momentum of partons continuously grows over time driven by multiple elastic scatterings.
Since there are still parametrically much less quarks and antiquarks compared to gluons, the number density of gluons decreases as mandated by energy conservation.And all the quantities in the system including the number density of quarks and antiquarks evolve according to a set of scaling laws, manifest as a universal scaling solution previously discovered in both classical statistical field simulations [29,58,59] and kinetic theory [26,28,34,38,60].Full thermalization is established when all the quantities exit their scaling regions and approach their corresponding thermal equilibrium values.
In order to complement our parametric descriptions, we have also presented quantitative results for several scenarios with different initial occupancies: f 0 = 10 −6 , 10 −4 , 10 −2 , 1, 10 and 100.Here, we have shown several examples in which the thermalization time becomes longer when quarks and antiquarks are allowed to participate in the system's evolution.
Nevertheless, this is true when one compares directly thermalization times between systems with identical initial conditions.In phenomenological studies, one might be more interested in looking at systems with the same final configuration, that is, the same temperature, viscosity, etc.In this case, one can naively re-scale the time variable with the equilibrium temperature of the system T , such that t th → t th T .The value of T can be computed using eq.( 36), yielding T (N f = 0) ≈ 1.31T (N f = 3).Therefore, this is irrelevant for overoccupied systems, where the thermalization time varies around one order of magnitude.For the under-occupied scenario, this is not relevant either, since this is the case where the thermalization times for N f = 0 and N f = 3 are really close and the scaling factor of T barely changes the result.
Another more complex option has been used in previous works [33,[35][36][37][38] with the main purpose of studying the isotropization of the system.In this case, the re-scaling variable is t → tT 4πη/s .Here, the entropy density s can be extracted from eq. ( 37) and the values of the viscosity η can be calculated for the scenarios that we have studied according to ref. [61].In this way, we observe no significant reduction of the thermalization time for over-populated systems.However, this contribution became more important when the initial system has a small population, where the thermalization time of the scenarios approach faster with the re-scaled time than with our previous set-up.
One universal feature for the cases with f 0 ≥ 10 −4 that could not be revealed by the above parametric estimates is as follows.Before achieving full thermalization, gluons undergo topdown thermalization, meaning that they first equilibrate at temperature T * higher than T , and then the thermal bath of gluons cools down to T via generating quarks and antiquarks.
The whole thermalization process is hence delayed compared to the pure gluon cases.Among the aforementioned values of f 0 , the reheating stage is only observed for f 0 = 10 −6 while the self-similarity in both gluon and quark distributions and the scaling laws are only discernible for f 0 = 100.These observations qualitatively agree with those using EKT in refs.[35][36][37][38], Thus, the BEDA provides us a framework which retains the essential physics of EKT while offering two advantages owing to its simplicity.Firstly, it is more cost-effective to implement numerically, enabling exploration of more realistic geometries.Secondly, the BEDA is easier to handle for analytical estimations.A quantitative comparison between the BEDA and EKT will be presented in an forthcoming publication.
for a uniform grid.Then, the current of the Fokker-Planck term, p 2 J, is easily computed on the p s grid.In terms of J on p s and S on p, one has C a 2↔2,i = p 2 s,i J s,i − p 2 s,i−1 J s,i−1 ∆V i + S i (B4) with ∆V i = 1 3 (p 3 s,i − p 3 s,i−1 ).One also needs to specify two boundary conditions: s,−1 J s,−1 = 0, p 2 s,N J s,N = 0 (B5) with p s,−1 = 0.In initially over-populated cases, the first term can be set without conflict only if the inelastic kernel is included.In this case, the distributions of both quarks and gluons approach an equilibrium profile with temperature T * at small p, as shown in eq. ( 26).
The second boundary condition just approximates the restriction that there can not be particle flux at p = ∞.

Discretization of the 1 ↔ 2 kernel
Evaluation of this kernel consumes most of the computation time since it involves an integral of the form of eq. ( 12) for each momentum grid element p i .That is, we need to compute where K(p i , x) includes the collinear splitting functions and the statistical factors given in eq. ( 14).In our code, the x-integration is carried out by using the trapezoidal rule.
Computing the above integral involves evaluating the distribution functions at momentum p different from the momentum grid points.To handle this, we employ the third order Lagrange interpolation (for pf and F/ F ) by using four adjacent grid points: two before and two after p if it does not lie between the first or the last two grid points.Otherwise, the linear interpolation is used.
The evaluation of C 1↔2 (p) at each time step is amenable to parallel computation.We implement this calculation using both the GPU and CPU programming.Given our current computational resources, the GPU code is found to be about one order of magnitude faster

FIG. 2 .
FIG. 2. The time evolution of different macroscopic quantities for f 0 = 0.01.The left panel shows the results of T * , m 2D , q and n g with N f = 0 (dotted) and N f = 3 (solid).The right plot shows the result of n q in comparison with the scaling laws obtained by parametric estimates.

3 2 ,
as shown in the right panel of fig. 2.

FIG. 4 .
FIG. 4. Gluon (left) and quark (right) distributions at different times for f 0 = 0.01.Starting around Qt = 200, the gluon distribution can be well fitted by a Bose-Einstein distribution (thermal) with a time-dependent temperature T * up to some maximum momentum larger than T * while the quark distribution only approaches a Fermi-Dirac distribution (thermal) towards the end of the process.

FIG. 5 .
FIG.5.Thermalization for f 0 = 0.01.The upper left panel shows the ratios of T * , s and n to their final thermal equilibrium values for N f = 0 and 3.The lower left panel shows the ratios of the energy density, number density and entropy density of g and q to their thermal equilibrium values with the temperature given by T or T * (only their dependence on T * is kept explicit).The two right panels show p 4 f at different times for N f = 0 (upper) and N f = 3 (lower).

FIG. 7 .
FIG. 7. Evolution at late times for f 0 = 10 −4 .The left panel shows an amplified representation of T * for both N f = 0 and N f = 3.The right panel shows thermalization of the same set of quantities as those in the lower left panels of fig. 5.

FIG. 8 .
FIG.8.The gluon (left) and quark (right) distributions at different times for f 0 = 10 −4 .At low momentum, pf follows the complex evolution pattern of T * : overheating, cooling, reheating and cooling except for every early times when the thermalized soft sector is out of the shown range of p.The quark distribution F , on the other hand, keeps increasing over t across the shown range of p towards the thermal distribution (dashed).

F 0 = F0 = 0 FIG. 9 .
FIG. 9. Thermalization for f 0 = 10 −6 .The left panel shows T * as a function of Qt for both N f = 0 and N f = 3.The right panel shows the time evolution of the same set of quantities as those in the lower left panels of fig. 5.

10 FIG. 12 .
FIG. 12.The top-down for f 0 = 1.The left panel shows how the same set of quantities as those in the lower left panels of fig. 5 approach unity.The right panel shows the gluon distribution at different times, comparing to the Bose-Einstein distribution (thermal) with a time-dependent temperature T * .

3 n g /Q 3 f0 = 10 F0FIG. 13 .
FIG. 13.The time evolution of different macroscopic quantities for f 0 = 10.The left panel shows the results of T * , m 2D , q and n g for N f = 0 (dotted) and N f = 3 (solid).The right panel shows n q over time in comparison with the scaling in the limit f 0 ≫ 1.

25 FIG. 14 .
FIG. 14.The top-down thermalization of gluons for f 0 = 10.The left panel shows how the same set of quantities as those in the lower left panel of fig. 5 approach unity.The right panel shows the gluon distribution at different times, comparing to the Bose-Einstein distribution (thermal) with a time-dependent temperature T * .

FIG. 16 .
FIG.16.Self-similar scaling behaviors of the parton distributions for f 0 = 100 and α s = 0.01.The left panels show the gluon (upper) and quark (lower) distributions at three different times in the range of t where all the macroscopic quantities in fig.15follow the scaling laws.The right panels compare the distributions at Qt = 0.3 and the others rescaled according to eq. (78).

FIG. 17 .
FIG. 17.The top-down thermalization of gluons for f 0 = 100 and α s = 0.01.The left panel shows how the same set of quantities as those in the lower left panel of fig. 5 approach unity.The right panel shows the gluon distribution at different times, comparing to the Bose-Einstein distribution (thermal) with a time-dependent temperature T * .
[46]becomes comparable with t, meaning that soft partons have enough time to establish thermal equilibrium among themselves.Later on, each hard parton loses energy by typically radiating one gluon with p br ∼ α 2 s qt 2 [15, 44], which rapidly degrades into soft partons with p s ∼ T * [45] via democratic branching[46]within a time of order t.As the energy density of the QGP, sourced by hard partons, is given by T 4 * ∼ ϵ s ∼ p br n h,g , one has